DGEES sample using the select function

[srao12]

Hi, I am using C, to call the LAPACK function DGEES to obtain the Schur decomposition of a matrix. However, it always put the eigenvalue that has a real part 1 and imaginary part 0, on the (1,1) position, which I don't want to allow. I tried using the sort option, but that causes a crash ( the second time I call dgeees). I couldn't find any sample example for doing this, and used the documentation in http://www.netlib.org/lapack/explore-ht ... 36b284b438 for reference. Below is my code. Could someone please guide me as to where I am going wrong? It worked fine when I had the sort option set to "N".

Code: Select all

#include "lapacke.h"
#include <stdio.h>
#include <stdlib.h>
int schur(double *val, int nrows, double *U, double *T)
{
   int lwork = -1;
   int   info;
   int sdim = nrows-1;
   double wkopt;
   double *wr = NULL, *wi = NULL, *Z = NULL;
   double   *work = NULL;
   extern  LAPACK_D_SELECT2          SELECT;
   lapack_logical    *bwork;
   
   wr = (double *)calloc(nrows,sizeof(double));
   wi = (double *)calloc(nrows,sizeof(double));
   Z  = (double *)calloc(nrows,sizeof(double));
   bwork  = (lapack_logical *)calloc(nrows,sizeof(lapack_logical));
   
   LAPACK_dgees ( "V", "S", SELECT, &nrows, val, &nrows, &sdim, wr, wi, Z, &nrows, &wkopt, &lwork, bwork, &info ); /*SELECT is never references because SORT option is set to "N", can ignore warning about unitialized variable*/ 
    lwork = (int)wkopt;
    work = (double*)malloc( lwork*sizeof(double) );
   LAPACK_dgees ( "V", "S", SELECT, &nrows, val, &nrows, &sdim, wr, wi, Z, &nrows, work, &lwork, bwork, &info ); /*SELECT is never references because SORT option is set to "N", can ignore warning about unitialized variable*/   

    if( info > 0 )
   {
      printf( "The algorithm computing Schur Decompsition failed to converge.\n" );
        exit( 1 );
    }
   free(wr);
   free(wi);
   free(Z);
   free(work);

   return(0);
}


lapack_logical SELECT ( double* ER, double* EI)
{
   if ((*ER==1) && (*EI==0))
      return(0);
   else
      return(1);
}