jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

*** empty log message ***

Browse Source
This commit is contained in:
jiri 2009-05-28 14:14:12 +00:00
parent 46170fcca9
commit 8a060e9abf
1 changed files with 15 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
nonclass.cc
View File

@ -286,10 +286,10 @@ void diagonalize(NRMat<double> &a, NRVec<double> &w, const bool eivec,
extern "C" void FORNAME(zheev)(const char *JOBZ, const char *UPLO, const int *N, extern "C" void FORNAME(zheev)(const char *JOBZ, const char *UPLO, const int *N,
complex<double> *A, const int *LDA, double *W, double *WORK, const int *LWORK, int *INFO); complex<double> *A, const int *LDA, double *W, complex<double> *WORK, const int *LWORK, double *RWORK, int *INFO);
extern "C" void FORNAME(zhegv)(const int *ITYPE, const char *JOBZ, const char *UPLO, const int *N, extern "C" void FORNAME(zhegv)(const int *ITYPE, const char *JOBZ, const char *UPLO, const int *N,
complex<double> *A, const int *LDA, complex<double> *B, const int *LDB, double *W, double *WORK, const int *LWORK, int *INFO); complex<double> *A, const int *LDA, complex<double> *B, const int *LDB, double *W, complex<double> *WORK, const int *LWORK, double *RWORK, int *INFO);
// a will contain eigenvectors (columns if corder==1), w eigenvalues // a will contain eigenvectors (columns if corder==1), w eigenvalues
@ -313,21 +313,23 @@ void diagonalize(NRMat<complex<double> > &a, NRVec<double> &w, const bool eivec,
char vectors = 'V'; char vectors = 'V';
if (!eivec) vectors = 'N'; if (!eivec) vectors = 'N';
int LWORK = -1; int LWORK = -1;
double WORKX; complex<double> WORKX;
int ldb=0; if(b) ldb=b->ncols(); int ldb=0; if(b) ldb=b->ncols();
// First call is to determine size of workspace // First call is to determine size of workspace
if(b) FORNAME(zhegv)(&itype,&vectors, &U, &n, a, &m, *b, &ldb, w, &WORKX, &LWORK, &r ); double *RWORK = new double(3*n);
else FORNAME(zheev)(&vectors, &U, &n, a, &m, w, &WORKX, &LWORK, &r ); if(b) FORNAME(zhegv)(&itype,&vectors, &U, &n, a, &m, *b, &ldb, w, &WORKX, &LWORK, RWORK, &r );
LWORK = (int)WORKX; else FORNAME(zheev)(&vectors, &U, &n, a, &m, w, &WORKX, &LWORK, RWORK, &r );
double *WORK = new double[LWORK]; LWORK = (int)WORKX.real();
if(b) FORNAME(zhegv)(&itype,&vectors, &U, &n, a, &m, *b,&ldb, w, WORK, &LWORK, &r ); complex<double> *WORK = new complex<double>[LWORK];
else FORNAME(zheev)(&vectors, &U, &n, a, &m, w, WORK, &LWORK, &r ); if(b) FORNAME(zhegv)(&itype,&vectors, &U, &n, a, &m, *b,&ldb, w, WORK, &LWORK, RWORK, &r );
else FORNAME(zheev)(&vectors, &U, &n, a, &m, w, WORK, &LWORK, RWORK, &r );
delete[] WORK; delete[] WORK;
delete[] RWORK;
if (vectors == 'V' && corder) a.transposeme(n); if (vectors == 'V' && corder) a.transposeme(n);
if (r < 0) laerror("illegal argument in sygv/syev in diagonalize()"); if (r < 0) laerror("illegal argument in hegv/heev in diagonalize()");
if (r > 0) laerror("convergence problem in sygv/syev in diagonalize()"); if (r > 0) laerror("convergence problem in hegv/heev in diagonalize()");
} }
@ -407,8 +409,8 @@ void diagonalize(NRSMat<complex<double> > &a, NRVec<double> &w, NRMat<complex<do
delete[] RWORK; delete[] RWORK;
if (v && corder) v->transposeme(n); if (v && corder) v->transposeme(n);
if (r < 0) laerror("illegal argument in spgv/spev in diagonalize()"); if (r < 0) laerror("illegal argument in hpgv/hpev in diagonalize()");
if (r > 0) laerror("convergence problem in spgv/spev in diagonalize()"); if (r > 0) laerror("convergence problem in hpgv/hpev in diagonalize()");
} }