From 8b55c3d0c1cdc6553cfcb7102571cb4672cce22c Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Tue, 27 Oct 2026 15:15:10 +0100 Subject: [PATCH] bugfixes in nonclass, mainly complex-related --- nonclass.cc | 336 +++++++++++++++++++++++++++++++++++++--------------- nonclass.h | 13 +- 2 files changed, 247 insertions(+), 102 deletions(-) diff --git a/nonclass.cc b/nonclass.cc index e004de6..21620dc 100644 --- a/nonclass.cc +++ b/nonclass.cc @@ -236,8 +236,14 @@ static void linear_solve_do(NRSMat &a, double *b, const int nrhs, const } if (det && r == 0) { *det = 1.; - for (int i=0; i0) *det = 0; delete[] ipiv; @@ -265,6 +271,66 @@ linear_solve_do(a,&B[0],1,a.nrows(),det,n); } +// complex hermitian version: lower triangle stored, upper triangle is its complex conjugate +extern "C" void FORNAME(zhpsv)(const char *UPLO, const FINT *N, const FINT *NRHS, + std::complex *AP, FINT *IPIV, std::complex *B, const FINT *LDB, FINT *INFO); + +static void linear_solve_do(NRSMat > &a, std::complex *b, const int nrhs, const int ldb, std::complex *det, int n) +{ + if(nrhs==0) std::cout<<"Warning: some zhpsv implementations might skip LU decomposition when nrhs==0\n"; + FINT r, *ipiv; + a.copyonwrite(); + ipiv = new FINT[n]; + char U = LAPACK_FORTRANCASE('u'); + //the packed C-order lower triangle is the fortran upper triangle of the transposed = complex conjugated hermitian matrix, + //solve conj(A) conj(x) = conj(b) + for(int i=0; i0) *det = 0; + delete[] ipiv; + if (r > 0 && b) laerror("singular matrix in linear_solve(SMat&, Mat*, complex*"); +} + + +void linear_solve(NRSMat > &a, NRMat > *B, std::complex *det, int n) +{ +if(n<=0) n=a.nrows(); + if (n==a.nrows() && B && a.nrows() != B->ncols() || B && n>B->ncols() || n>a.nrows()) + laerror("incompatible matrices in hermitian linear_solve()"); + if (B) B->copyonwrite(); +linear_solve_do(a,B?(*B)[0]:NULL,B?B->nrows() : 0, B?B->ncols():a.nrows(),det,n); +} + + +void linear_solve(NRSMat > &a, NRVec > &B, std::complex *det, int n) +{ +if(n<=0) n=a.nrows(); + if (n==a.nrows() && a.nrows()!= B.size() || n>B.size() || n>a.nrows()) + laerror("incompatible matrices in hermitian linear_solve()"); + B.copyonwrite(); +linear_solve_do(a,&B[0],1,a.nrows(),det,n); +} + + // Roman, complex version of linear_solve() extern "C" void FORNAME(zgesv)(const FINT *N, const FINT *NRHS, std::complex *A, const FINT *LDA, FINT *IPIV, std::complex *B, const FINT *LDB, FINT *INFO); @@ -360,7 +426,8 @@ extern "C" void FORNAME(dgesvx)(const char *fact, const char *trans, const FINT // solution is stored in _B // the info parameter of dgesvx is returned (see man dgesvx) //------------------------------------------------------------------------------ -int linear_solve_x(NRMat &_A, double *_B, const int _rhsCount, const int _eqCount, const bool _eq, const bool _saveA, double *_rcond){ +template<> +int linear_solve_x(NRMat &_A, double *_B, const int _rhsCount, const int _eqCount, const bool _eq, const bool _saveA, double *_rcond){ const int A_rows = _A.nrows(); const int A_cols = _A.ncols(); @@ -373,6 +440,7 @@ int linear_solve_x(NRMat &_A, double *_B, const int _rhsCount, const int } double *A; + if(_eq && !_saveA) _A.copyonwrite(); //A will be overwritten by the equilibrated matrix, detach before taking the pointer double * const _A_data = (double*)_A; FINT info; @@ -402,8 +470,6 @@ int linear_solve_x(NRMat &_A, double *_B, const int _rhsCount, const int offset2 += n; } lda = n;//!!! - }else{ - _A.copyonwrite(); } } @@ -414,7 +480,7 @@ int linear_solve_x(NRMat &_A, double *_B, const int _rhsCount, const int delete[] iwork;delete[] ipiv; delete[] AF;delete[] X; - if(_saveA){ + if(_eq && _saveA){//the copy is made only with equilibration delete[] A; } return (int)info; @@ -433,7 +499,8 @@ int linear_solve_x(NRMat &_A, double *_B, const int _rhsCount, const int // solution is stored in _B // the info parameter of dgesvx is returned (see man dgesvx) //------------------------------------------------------------------------------ -int linear_solve_x(NRMat > &_A, std::complex *_B, const int _rhsCount, const int _eqCount, const bool _eq, const bool _saveA, double *_rcond){ +template<> +int linear_solve_x >(NRMat > &_A, std::complex *_B, const int _rhsCount, const int _eqCount, const bool _eq, const bool _saveA, double *_rcond){ const int A_rows = _A.nrows(); const int A_cols = _A.ncols(); @@ -446,6 +513,7 @@ int linear_solve_x(NRMat > &_A, std::complex *_B, c } std::complex *A; + if(_eq && !_saveA) _A.copyonwrite(); //A will be overwritten by the equilibrated matrix, detach before taking the pointer std::complex * const _A_data = (std::complex*)_A; FINT info; @@ -475,8 +543,6 @@ int linear_solve_x(NRMat > &_A, std::complex *_B, c offset2 += n; } lda = n;//!!! - }else{ - _A.copyonwrite(); } } @@ -488,7 +554,7 @@ int linear_solve_x(NRMat > &_A, std::complex *_B, c delete[] ipiv; delete[] AF;delete[] X; - if(_saveA){ + if(_eq && _saveA){//the copy is made only with equilibration delete[] A; } return (int)info; @@ -519,9 +585,11 @@ int multiply_by_inverse(NRMat &_A, NRMat &_B, bool _useE char equed = 'B';//if fact=='N' then equed is an output argument, therefore not declared as const const int n2 = n*n; - double * const A = (double*)_A; + _B.copyonwrite();//even if fact='N', call copyonwrite because the solution is going to be stored in _B; before taking the pointer double * const B = (double*)_B; - _B.copyonwrite();//even if fact='N', call copyonwrite because the solution is going to be stored in _B + //_A is the right hand side for gesvx, which scales it if equilibrated: work on a copy then, _A is input only + double * const A = _useEq ? new double[n2] : (double*)_A; + if(_useEq) memcpy(A,(double*)_A,n2*sizeof(double)); FINT info; double rcond; @@ -542,6 +610,7 @@ int multiply_by_inverse(NRMat &_A, NRMat &_B, bool _useE delete[] iwork;delete[] ipiv; delete[] AF;delete[] X; + if(_useEq) delete[] A; return (int)info; } @@ -571,9 +640,11 @@ int multiply_by_inverse >(NRMat > &_A, const char trans = 'N';//because of c-order char equed = 'B';//if fact=='N' then equed is an output argument, therefore not declared as const - std::complex * const A = (std::complex*)_A; + _B.copyonwrite();//even if fact='N', call copyonwrite because the solution is going to be stored in _B; before taking the pointer std::complex * const B = (std::complex*)_B; - _B.copyonwrite();//even if fact='N', call copyonwrite because the solution is going to be stored in _B + //_A is the right hand side for gesvx, which scales it if equilibrated: work on a copy then, _A is input only + std::complex * const A = _useEq ? new std::complex[n2] : (std::complex*)_A; + if(_useEq) memcpy(A,(std::complex*)_A,n2*sizeof(std::complex)); FINT info; double rcond; @@ -594,6 +665,7 @@ int multiply_by_inverse >(NRMat > &_A, delete[] ipiv; delete[] AF;delete[] X; + if(_useEq) delete[] A; return (int)info; } @@ -648,8 +720,8 @@ void diagonalize(NRMat &a, NRVec &w, const bool eivec, double *WORK = new double[LWORK]; #ifdef FORINT - if(b) FORNAME(dsygv)(&itypetmp,&vectors, &U, &ntmp, a, &m, *b, &ldb, w, &WORKX, &LWORK, &r ); - else FORNAME(dsyev)(&vectors, &U, &ntmp, a, &m, w, &WORKX, &LWORK, &r ); + if(b) FORNAME(dsygv)(&itypetmp,&vectors, &U, &ntmp, a, &m, *b, &ldb, w, WORK, &LWORK, &r ); + else FORNAME(dsyev)(&vectors, &U, &ntmp, a, &m, w, WORK, &LWORK, &r ); #else if(b) FORNAME(dsygv)(&itype,&vectors, &U, &n, a, &m, *b,&ldb, w, WORK, &LWORK, &r ); else FORNAME(dsyev)(&vectors, &U, &n, a, &m, w, WORK, &LWORK, &r ); @@ -711,8 +783,8 @@ void diagonalize(NRMat > &a, NRVec &w, const bool e std::complex *WORK = new std::complex[LWORK]; #ifdef FORINT - if(b) FORNAME(zhegv)(&itypetmp,&vectors, &U, &ntmp, a, &m, *b, &ldb, w, &WORKX, &LWORK, RWORK, &r ); - else FORNAME(zheev)(&vectors, &U, &ntmp, a, &m, w, &WORKX, &LWORK, RWORK, &r ); + if(b) FORNAME(zhegv)(&itypetmp,&vectors, &U, &ntmp, a, &m, *b, &ldb, w, WORK, &LWORK, RWORK, &r ); + else FORNAME(zheev)(&vectors, &U, &ntmp, a, &m, w, WORK, &LWORK, RWORK, &r ); #else 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 ); @@ -815,7 +887,7 @@ void diagonalize(NRSMat > &a, NRVec &w, NRMattransposeme(n); + if (v && corder) {v->transposeme(n); v->conjugateme();} if (r < 0) laerror("illegal argument in hpgv/hpev in diagonalize()"); if (r > 0) laerror("convergence problem in hpgv/hpev in diagonalize()"); @@ -1086,30 +1158,30 @@ void gdiagonalize(NRMat &a, NRVec &wr, NRVec &wi, #ifdef FORINT FINT ntmp = n; - if(b) FORNAME(dggev)(&jobvr, &jobvl, &ntmp, a, &lda, *b, &ldb, wr, wi, *beta, vr?vr[0]:(double *)0, - &ldvr, vl?vl[0]:(double *)0, &ldvl, &work0, &lwork, &r); - else FORNAME(dgeev)(&jobvr, &jobvl, &ntmp, a, &lda, wr, wi, vr?vr[0]:(double *)0, - &ldvr, vl?vl[0]:(double *)0, &ldvl, &work0, &lwork, &r); + if(b) FORNAME(dggev)(&jobvr, &jobvl, &ntmp, a, &lda, *b, &ldb, wr, wi, *beta, vr?(*vr)[0]:(double *)0, + &ldvr, vl?(*vl)[0]:(double *)0, &ldvl, &work0, &lwork, &r); + else FORNAME(dgeev)(&jobvr, &jobvl, &ntmp, a, &lda, wr, wi, vr?(*vr)[0]:(double *)0, + &ldvr, vl?(*vl)[0]:(double *)0, &ldvl, &work0, &lwork, &r); #else - if(b) FORNAME(dggev)(&jobvr, &jobvl, &n, a, &lda, *b, &ldb, wr, wi, *beta, vr?vr[0]:(double *)0, - &ldvr, vl?vl[0]:(double *)0, &ldvl, &work0, &lwork, &r); - else FORNAME(dgeev)(&jobvr, &jobvl, &n, a, &lda, wr, wi, vr?vr[0]:(double *)0, - &ldvr, vl?vl[0]:(double *)0, &ldvl, &work0, &lwork, &r); + if(b) FORNAME(dggev)(&jobvr, &jobvl, &n, a, &lda, *b, &ldb, wr, wi, *beta, vr?(*vr)[0]:(double *)0, + &ldvr, vl?(*vl)[0]:(double *)0, &ldvl, &work0, &lwork, &r); + else FORNAME(dgeev)(&jobvr, &jobvl, &n, a, &lda, wr, wi, vr?(*vr)[0]:(double *)0, + &ldvr, vl?(*vl)[0]:(double *)0, &ldvl, &work0, &lwork, &r); #endif lwork = (FINT) work0; double *work = new double[lwork]; #ifdef FORINT - if(b) FORNAME(dggev)(&jobvr, &jobvl, &ntmp, a, &lda, *b, &ldb, wr, wi, *beta, vr?vr[0]:(double *)0, - &ldvr, vl?vl[0]:(double *)0, &ldvl, work, &lwork, &r); - else FORNAME(dgeev)(&jobvr, &jobvl, &ntmp, a, &lda, wr, wi, vr?vr[0]:(double *)0, - &ldvr, vl?vl[0]:(double *)0, &ldvl, work, &lwork, &r); + if(b) FORNAME(dggev)(&jobvr, &jobvl, &ntmp, a, &lda, *b, &ldb, wr, wi, *beta, vr?(*vr)[0]:(double *)0, + &ldvr, vl?(*vl)[0]:(double *)0, &ldvl, work, &lwork, &r); + else FORNAME(dgeev)(&jobvr, &jobvl, &ntmp, a, &lda, wr, wi, vr?(*vr)[0]:(double *)0, + &ldvr, vl?(*vl)[0]:(double *)0, &ldvl, work, &lwork, &r); #else - if(b) FORNAME(dggev)(&jobvr, &jobvl, &n, a, &lda, *b, &ldb, wr, wi, *beta, vr?vr[0]:(double *)0, - &ldvr, vl?vl[0]:(double *)0, &ldvl, work, &lwork, &r); - else FORNAME(dgeev)(&jobvr, &jobvl, &n, a, &lda, wr, wi, vr?vr[0]:(double *)0, - &ldvr, vl?vl[0]:(double *)0, &ldvl, work, &lwork, &r); + if(b) FORNAME(dggev)(&jobvr, &jobvl, &n, a, &lda, *b, &ldb, wr, wi, *beta, vr?(*vr)[0]:(double *)0, + &ldvr, vl?(*vl)[0]:(double *)0, &ldvl, work, &lwork, &r); + else FORNAME(dgeev)(&jobvr, &jobvl, &n, a, &lda, wr, wi, vr?(*vr)[0]:(double *)0, + &ldvr, vl?(*vl)[0]:(double *)0, &ldvl, work, &lwork, &r); #endif delete[] work; @@ -1178,6 +1250,19 @@ void gdiagonalize(NRMat &a, NRVec &wr, NRVec &wi, } } + //LAPACK saw the transposed matrix and left and right were swapped, so for a complex pair it returns the eigenvector of w_i as Re - i Im; + //negate Im to obtain the dgeev convention Re + i Im (for both left and right eigenvectors; biorthonormality is preserved) + { + int i=0; + while(i0) { @@ -1252,30 +1337,30 @@ void gdiagonalize(NRMat > &a, NRVec< std::complex > #ifdef FORINT FINT ntmp = n; - if(b) FORNAME(zggev)(&jobvr, &jobvl, &ntmp, a, &lda, *b, &ldb, w, *beta, vr?vr[0]:(std::complex *)0, - &ldvr, vl?vl[0]:(std::complex *)0, &ldvl, &work0, &lwork, rwork, &r); - else FORNAME(zgeev)(&jobvr, &jobvl, &ntmp, a, &lda, w, vr?vr[0]:(std::complex *)0, - &ldvr, vl?vl[0]:(std::complex *)0, &ldvl, &work0, &lwork, rwork, &r); + if(b) FORNAME(zggev)(&jobvr, &jobvl, &ntmp, a, &lda, *b, &ldb, w, *beta, vr?(*vr)[0]:(std::complex *)0, + &ldvr, vl?(*vl)[0]:(std::complex *)0, &ldvl, &work0, &lwork, rwork, &r); + else FORNAME(zgeev)(&jobvr, &jobvl, &ntmp, a, &lda, w, vr?(*vr)[0]:(std::complex *)0, + &ldvr, vl?(*vl)[0]:(std::complex *)0, &ldvl, &work0, &lwork, rwork, &r); #else - if(b) FORNAME(zggev)(&jobvr, &jobvl, &n, a, &lda, *b, &ldb, w, *beta, vr?vr[0]:(std::complex *)0, - &ldvr, vl?vl[0]:(std::complex *)0, &ldvl, &work0, &lwork, rwork, &r); - else FORNAME(zgeev)(&jobvr, &jobvl, &n, a, &lda, w, vr?vr[0]:(std::complex *)0, - &ldvr, vl?vl[0]:(std::complex *)0, &ldvl, &work0, &lwork, rwork, &r); + if(b) FORNAME(zggev)(&jobvr, &jobvl, &n, a, &lda, *b, &ldb, w, *beta, vr?(*vr)[0]:(std::complex *)0, + &ldvr, vl?(*vl)[0]:(std::complex *)0, &ldvl, &work0, &lwork, rwork, &r); + else FORNAME(zgeev)(&jobvr, &jobvl, &n, a, &lda, w, vr?(*vr)[0]:(std::complex *)0, + &ldvr, vl?(*vl)[0]:(std::complex *)0, &ldvl, &work0, &lwork, rwork, &r); #endif lwork = (FINT) work0.real(); std::complex *work = new std::complex[lwork]; #ifdef FORINT - if(b) FORNAME(zggev)(&jobvr, &jobvl, &ntmp, a, &lda, *b, &ldb, w, *beta, vr?vr[0]:(std::complex *)0, - &ldvr, vl?vl[0]:(std::complex *)0, &ldvl, work, &lwork, rwork, &r); - else FORNAME(zgeev)(&jobvr, &jobvl, &ntmp, a, &lda, w, vr?vr[0]:(std::complex *)0, - &ldvr, vl?vl[0]:(std::complex *)0, &ldvl, work, &lwork, rwork, &r); + if(b) FORNAME(zggev)(&jobvr, &jobvl, &ntmp, a, &lda, *b, &ldb, w, *beta, vr?(*vr)[0]:(std::complex *)0, + &ldvr, vl?(*vl)[0]:(std::complex *)0, &ldvl, work, &lwork, rwork, &r); + else FORNAME(zgeev)(&jobvr, &jobvl, &ntmp, a, &lda, w, vr?(*vr)[0]:(std::complex *)0, + &ldvr, vl?(*vl)[0]:(std::complex *)0, &ldvl, work, &lwork, rwork, &r); #else - if(b) FORNAME(zggev)(&jobvr, &jobvl, &n, a, &lda, *b, &ldb, w, *beta, vr?vr[0]:(std::complex *)0, - &ldvr, vl?vl[0]:(std::complex *)0, &ldvl, work, &lwork, rwork, &r); - else FORNAME(zgeev)(&jobvr, &jobvl, &n, a, &lda, w, vr?vr[0]:(std::complex *)0, - &ldvr, vl?vl[0]:(std::complex *)0, &ldvl, work, &lwork, rwork, &r); + if(b) FORNAME(zggev)(&jobvr, &jobvl, &n, a, &lda, *b, &ldb, w, *beta, vr?(*vr)[0]:(std::complex *)0, + &ldvr, vl?(*vl)[0]:(std::complex *)0, &ldvl, work, &lwork, rwork, &r); + else FORNAME(zgeev)(&jobvr, &jobvl, &n, a, &lda, w, vr?(*vr)[0]:(std::complex *)0, + &ldvr, vl?(*vl)[0]:(std::complex *)0, &ldvl, work, &lwork, rwork, &r); #endif delete[] work; @@ -1295,7 +1380,6 @@ void gdiagonalize(NRMat > &a, NRVec< std::complex > std::complex tmp; cblas_zdotc_sub(n,(*vr)[i],1,(*vl)[i], 1, &tmp); tmp = 1./tmp; - std::cout <<"scaling by "< &a, NRVec< std::complex > &w, } else { - (*vl)[i][j] = std::complex((*rvl)[i][j], (*rvl)[i+1][j]); - (*vl)[i+1][j] = std::complex((*rvl)[i][j], -(*rvl)[i+1][j]); + (*vl)[i][j] = std::complex((*rvl)[i][j], -(*rvl)[i+1][j]); //rows are the conjugated eigenvectors + (*vl)[i+1][j] = std::complex((*rvl)[i][j], (*rvl)[i+1][j]); } } if (vr) @@ -1372,8 +1456,8 @@ void gdiagonalize(NRMat &a, NRVec< std::complex > &w, } else { - (*vr)[i][j] = std::complex((*rvr)[i][j], (*rvr)[i+1][j]); - (*vr)[i+1][j] = std::complex((*rvr)[i][j], -(*rvr)[i+1][j]); + (*vr)[i][j] = std::complex((*rvr)[i][j], -(*rvr)[i+1][j]); //rows are the conjugated eigenvectors + (*vr)[i+1][j] = std::complex((*rvr)[i][j], (*rvr)[i+1][j]); } } i += 2; @@ -1851,7 +1935,37 @@ double r=0; return r; } -inline double trace2(const NRMat &a, const NRSMat &b, const bool diagscaled) +double trace2(const NRMat &a, const NRSMat &b, const bool diagscaled) +{ +return trace2(b,a,diagscaled); +} + +//complex hermitian versions: lower triangle stored, upper triangle is its complex conjugate +//Tr(A B) = sum_ij A_ij B_ji = sum_ij A_ij conj(B_ij), real for hermitian A, B +std::complex trace2(const NRSMat > &a, const NRSMat > &b, + const bool diagscaled) +{ + if (a.nrows() != b.nrows()) laerror("incompatible SMats in trace2()"); + + std::complex dot; + cblas_zdotc_sub(a.nrows()*(a.nrows()+1)/2, b, 1, a, 1, &dot); //sum over the lower triangle of a_ij conj(b_ij) + double r = 2.0*dot.real(); + if (diagscaled) return r; + for (int i=0; i trace2(const NRSMat > &a, const NRMat > &b, const bool diagscaled) +{ + if (a.nrows() != b.nrows()||b.nrows()!=b.ncols()) laerror("incompatible SMats in trace2()"); +std::complex r=0; + int i, j, k=0; + for (i=0; i trace2(const NRMat > &a, const NRSMat > &b, const bool diagscaled) { return trace2(b,a,diagscaled); } @@ -1900,24 +2014,40 @@ else //various norms -extern "C" double FORNAME(zlange)( const char *NORM, const FINT *M, const FINT *N, std::complex *A, const FINT *LDA, double *WORK); -extern "C" double FORNAME(dlange)( const char *NORM, const FINT *M, const FINT *N, double *A, const FINT *LDA, double *WORK); +extern "C" double FORNAME(zlange)( const char *NORM, const FINT *M, const FINT *N, const std::complex *A, const FINT *LDA, double *WORK); +extern "C" double FORNAME(dlange)( const char *NORM, const FINT *M, const FINT *N, const double *A, const FINT *LDA, double *WORK); -double MatrixNorm(NRMat > &A, const char norm) +//LAPACK sees the C-order matrix transposed: 1-norm and infinity norm are swapped, Frobenius and max-abs norms are invariant +static char lapack_norm_type(const char norm) { - const char TypNorm = (tolower(norm) == 'o')?'I':'O'; //switch c-order/fortran-order - const FINT M = A.nrows(); - const FINT N = A.ncols(); +switch(tolower(norm)) + { + case 'o': case '1': return 'I'; + case 'i': return 'O'; + case 'f': case 'e': return 'F'; + case 'm': return 'M'; + default: laerror("unknown norm type in MatrixNorm/CondNumber"); + } +return 0; +} + +template<> +LA_traits > >::normtype MatrixNorm(const NRMat > &A, const char norm) +{ + const char TypNorm = lapack_norm_type(norm); + const FINT M = A.ncols(); //fortran-order view of the C-order matrix is its transpose + const FINT N = A.nrows(); double work[M]; const double ret = FORNAME(zlange)(&TypNorm, &M, &N, A[0], &M, &work[0]); return ret; } -double MatrixNorm(NRMat &A, const char norm) +template<> +LA_traits >::normtype MatrixNorm(const NRMat &A, const char norm) { - const char TypNorm = (tolower(norm) == 'o')?'I':'O'; //switch c-order/fortran-order - const FINT M = A.nrows(); - const FINT N = A.ncols(); + const char TypNorm = lapack_norm_type(norm); + const FINT M = A.ncols(); //fortran-order view of the C-order matrix is its transpose + const FINT N = A.nrows(); double work[M]; const double ret = FORNAME(dlange)(&TypNorm, &M, &N, A[0], &M, &work[0]); return ret; @@ -1925,54 +2055,66 @@ double MatrixNorm(NRMat &A, const char norm) -//condition number -extern "C" void FORNAME(zgecon)( const char *norm, const FINT *n, std::complex *A, const FINT *LDA, const double *anorm, double *rcond, std::complex *work, double *rwork, FINT *info); -extern "C" void FORNAME(dgecon)( const char *norm, const FINT *n, double *A, const FINT *LDA, const double *anorm, double *rcond, double *work, double *rwork, FINT *info); +//condition number ||A|| ||A^{-1}|| in the 1-norm or infinity norm, from the LU factorization of a copy of A +extern "C" void FORNAME(zgetrf)(const FINT *M, const FINT *N, std::complex *A, const FINT *LDA, FINT *IPIV, FINT *INFO); +extern "C" void FORNAME(dgetrf)(const FINT *M, const FINT *N, double *A, const FINT *LDA, FINT *IPIV, FINT *INFO); +extern "C" void FORNAME(zgecon)( const char *norm, const FINT *n, const std::complex *A, const FINT *LDA, const double *anorm, double *rcond, std::complex *work, double *rwork, FINT *info); +extern "C" void FORNAME(dgecon)( const char *norm, const FINT *n, const double *A, const FINT *LDA, const double *anorm, double *rcond, double *work, FINT *iwork, FINT *info); -double CondNumber(NRMat > &A, const char norm) +template<> +LA_traits > >::normtype CondNumber(const NRMat > &A, const char norm) { - const char TypNorm = (tolower(norm) == 'o')?'I':'O'; //switch c-order/fortran-order + const char TypNorm = lapack_norm_type(norm); + if(TypNorm!='O' && TypNorm!='I') laerror("CondNumber implemented only for the 1-norm and infinity norm"); const FINT N = A.nrows(); - double Norma(0.0), ret(0.0); + if(N != A.ncols()) laerror("nonsquare matrix in CondNumber"); + double Norma(0.0), rcond(0.0); FINT info; - std::complex *work; - double *rwork; - if(N != A.ncols()){ - laerror("nonsquare matrix in zgecon"); - return 0.0; - } - work = new std::complex[2*N]; - rwork = new double[2*N]; + NRMat > LU(A); + LU.copyonwrite(); + FINT *ipiv = new FINT[N]; + FORNAME(zgetrf)(&N, &N, LU[0], &N, ipiv, &info); + delete[] ipiv; + if(info<0) laerror("illegal argument in zgetrf in CondNumber"); + if(info>0) return std::numeric_limits::infinity(); //exactly singular + std::complex *work = new std::complex[2*N]; + double *rwork = new double[2*N]; Norma = MatrixNorm(A, norm); - FORNAME(zgecon)(&TypNorm, &N, A[0], &N, &Norma, &ret, &work[0], &rwork[0], &info); + FORNAME(zgecon)(&TypNorm, &N, LU[0], &N, &Norma, &rcond, work, rwork, &info); delete[] work; delete[] rwork; - return ret; + if(info<0) laerror("illegal argument in zgecon in CondNumber"); + return rcond==0. ? std::numeric_limits::infinity() : 1./rcond; } -double CondNumber(NRMat &A, const char norm) +template<> +LA_traits >::normtype CondNumber(const NRMat &A, const char norm) { - const char TypNorm = (tolower(norm) == 'o')?'I':'O'; //switch c-order/fortran-order + const char TypNorm = lapack_norm_type(norm); + if(TypNorm!='O' && TypNorm!='I') laerror("CondNumber implemented only for the 1-norm and infinity norm"); const FINT N = A.nrows(); - double Norma(0.0), ret(0.0); + if(N != A.ncols()) laerror("nonsquare matrix in CondNumber"); + double Norma(0.0), rcond(0.0); FINT info; - double *work; - double *rwork; - if(N != A.ncols()){ - laerror("nonsquare matrix in zgecon"); - return 0.0; - } - work = new double[2*N]; - rwork = new double[2*N]; + NRMat LU(A); + LU.copyonwrite(); + FINT *ipiv = new FINT[N]; + FORNAME(dgetrf)(&N, &N, LU[0], &N, ipiv, &info); + delete[] ipiv; + if(info<0) laerror("illegal argument in dgetrf in CondNumber"); + if(info>0) return std::numeric_limits::infinity(); //exactly singular + double *work = new double[4*N]; + FINT *iwork = new FINT[N]; Norma = MatrixNorm(A, norm); - FORNAME(dgecon)(&TypNorm, &N, A[0], &N, &Norma, &ret, &work[0], &rwork[0], &info); + FORNAME(dgecon)(&TypNorm, &N, LU[0], &N, &Norma, &rcond, work, iwork, &info); delete[] work; - delete[] rwork; - return ret; + delete[] iwork; + if(info<0) laerror("illegal argument in dgecon in CondNumber"); + return rcond==0. ? std::numeric_limits::infinity() : 1./rcond; } diff --git a/nonclass.h b/nonclass.h index e03d4d1..952606b 100644 --- a/nonclass.h +++ b/nonclass.h @@ -177,13 +177,16 @@ NRVec w(a.ncols()); NRMat u(a.nrows(),a.ncols()); NRMat vt(a.ncols(),a.ncols()); singular_decomposition(a,&u,w,&vt,false); -NRVec utb = b*u; +NRVec utb(a.ncols()); +utb.gemv((T)0,u,'c',(T)1,b); //U^H b for(int i=0; ithres) utb[i] /= w[i]; else {utb[i]=0; if(nullity) ++*nullity;} } -return utb*vt; +NRVec res(a.ncols()); +res.gemv((T)0,vt,'c',(T)1,utb); //V W^-1 U^H b, vt contains V^H +return res; } @@ -217,7 +220,7 @@ return res; //corder =0 ... C rows are eigenvectors, =1 ... C columns are eigenvectors extern void gdiagonalize(NRMat &a, NRVec &wr, NRVec &wi, NRMat *vl=NULL, NRMat *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0, - NRMat *b=NULL, NRVec *beta=NULL); //this used real storage of eigenvectors like dgeev + NRMat *b=NULL, NRVec *beta=NULL); //this used real storage of eigenvectors like dgeev: for a complex pair (i,i+1) the eigenvector of w_i is Re + i Im = vector i + i*vector i+1 template extern void gdiagonalize(NRMat &a, NRVec< std::complex > &w, @@ -290,11 +293,11 @@ const NRMat calcsvdinverse(NRMat a, double thr=0, bool do_absadjoint=false } -//several matrix norms +//several matrix norms: norm = '1' or 'O' (max column sum), 'I' (max row sum), 'F' or 'E' (Frobenius), 'M' (max abs element) template typename LA_traits::normtype MatrixNorm(const MAT &A, const char norm); -//condition number +//condition number ||A|| ||A^-1|| for norm '1'/'O' or 'I' (LAPACK estimate: a lower bound, usually exact; infinity for a singular matrix) template typename LA_traits::normtype CondNumber(const MAT &A, const char norm);