diff --git a/la_traits.h b/la_traits.h index 788db4e..226b01a 100644 --- a/la_traits.h +++ b/la_traits.h @@ -456,6 +456,11 @@ generate_traits_smat(NRSMat_from1) template struct LA_traits : LA_traits_aux::scalar_type> {}; + +//precision-dependent default thresholds (the double values are those used so far) +template struct LA_precision {static T svd_threshold() {return (T)1e-12;}}; +template<> struct LA_precision {static float svd_threshold() {return 1e-5f;}}; + }//namespace #endif diff --git a/mat.cc b/mat.cc index fcaedd9..c022fda 100644 --- a/mat.cc +++ b/mat.cc @@ -5348,24 +5348,48 @@ NRMat NRMat::inverse() const return calcinverse(*this); } +template<> +NRMat NRMat::inverse() const +{ +return calcinverse(*this); +} + template<> NRMat > NRMat >::inverse() const { return calcinverse(*this); } +template<> +NRMat > NRMat >::inverse() const +{ +return calcinverse(*this); +} + template<> NRMat NRMat::svdinverse(const double thr) const { return calcsvdinverse(*this,thr,false); } +template<> +NRMat NRMat::svdinverse(const float thr) const +{ +return calcsvdinverse(*this,thr,false); +} + template<> NRMat > NRMat >::svdinverse(const double thr) const { return calcsvdinverse(*this,thr,false); } +template<> +NRMat > NRMat >::svdinverse(const float thr) const +{ +return calcsvdinverse(*this,thr,false); +} + @@ -5375,12 +5399,24 @@ NRMat NRMat::absadjoint(const double thr) const return calcsvdinverse(*this,thr,true); } +template<> +NRMat NRMat::absadjoint(const float thr) const +{ +return calcsvdinverse(*this,thr,true); +} + template<> NRMat > NRMat >::absadjoint(const double thr) const { return calcsvdinverse(*this,thr,true); } +template<> +NRMat > NRMat >::absadjoint(const float thr) const +{ +return calcsvdinverse(*this,thr,true); +} + diff --git a/noncblas.h b/noncblas.h index 30ab416..137f77f 100644 --- a/noncblas.h +++ b/noncblas.h @@ -689,5 +689,8 @@ int cblas_errprn(int ierr, int info, char *form, ...); int clapack_dgesv(const CBLAS_ORDER Order, const int N, const int NRHS, double *A, const int lda, int *ipiv, double *B, const int ldb); +int clapack_sgesv(const CBLAS_ORDER Order, const int N, const int NRHS, + float *A, const int lda, int *ipiv, + float *B, const int ldb); #endif #endif diff --git a/nonclass.cc b/nonclass.cc index d6a3dae..ade946a 100644 --- a/nonclass.cc +++ b/nonclass.cc @@ -195,6 +195,62 @@ static void linear_solve_do(NRMat &A, double *B, const int nrhs, const i if (r>0 && B) laerror("singular matrix in lapack_gesv"); } +static void linear_solve_do(NRMat &A, float *B, const int nrhs, const int ldb, float *det, int n) +{ + int r, *ipiv; + int iswap=0; + + if(nrhs==0) std::cout<<"Warning: some dgesv implementations might skip LU decomposition when nrhs==0\n"; + if (n==A.nrows() && A.nrows() != A.ncols()) laerror("linear_solve() call for non-square matrix"); + A.copyonwrite(); + ipiv = new int[A.nrows()]; +#ifdef IPIV_DEBUG + for(int i=0; in) + { + std::cout <<"IPIV["< &A, NRMat *B, double *det, int n) { if(n<=0) n=A.nrows(); //default - whole matrix @@ -203,6 +259,14 @@ 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(NRMat &A, NRMat *B, float *det, int n) +{ +if(n<=0) n=A.nrows(); //default - whole matrix +if (n==A.nrows() && B && A.nrows() != B->ncols() || B && n>B->ncols() ||n>A.nrows()) laerror("incompatible matrices in 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(NRMat &A, NRVec &B, double *det, int n) { if(n<=0) n=A.nrows(); //default - whole matrix @@ -211,12 +275,22 @@ B.copyonwrite(); linear_solve_do(A,&B[0],1,A.nrows(),det,n); } +void linear_solve(NRMat &A, NRVec &B, float *det, int n) +{ +if(n<=0) n=A.nrows(); //default - whole matrix +if(n==A.nrows() && A.nrows() != B.size() || n>B.size()||n>A.nrows() ) laerror("incompatible matrices in linear_solve()"); +B.copyonwrite(); +linear_solve_do(A,&B[0],1,A.nrows(),det,n); +} + // Next routines are not available in clapack, fotran ones will be used with an // additional swap/transpose of outputs when needed extern "C" void FORNAME(dspsv)(const char *UPLO, const FINT *N, const FINT *NRHS, double *AP, FINT *IPIV, double *B, const FINT *LDB, FINT *INFO); +extern "C" void FORNAME(sspsv)(const char *UPLO, const FINT *N, const FINT *NRHS, + float *AP, FINT *IPIV, float *B, const FINT *LDB, FINT *INFO); static void linear_solve_do(NRSMat &a, double *b, const int nrhs, const int ldb, double *det, int n) { @@ -254,6 +328,42 @@ static void linear_solve_do(NRSMat &a, double *b, const int nrhs, const if (r > 0 && b) laerror("singular matrix in linear_solve(SMat&, Mat*, double*"); } +static void linear_solve_do(NRSMat &a, float *b, const int nrhs, const int ldb, float *det, int n) +{ + if(nrhs==0) std::cout<<"Warning: some dspsv implementations might skip LU decomposition when nrhs==0\n"; + FINT r, *ipiv; + a.copyonwrite(); + ipiv = new FINT[n]; + char U = LAPACK_FORTRANCASE('u'); +#ifdef FORINT + const FINT ntmp=n; + const FINT nrhstmp=nrhs; + const FINT ldbtmp=ldb; + FORNAME(sspsv)(&U, &ntmp, &nrhstmp, a, ipiv, b, &ldbtmp,&r); +#else + FORNAME(sspsv)(&U, &n, &nrhs, a, ipiv, b, &ldb,&r); +#endif +// std::cout <<"A after dspsv = "<0) *det = 0; + delete[] ipiv; + if (r > 0 && b) laerror("singular matrix in linear_solve(SMat&, Mat*, float*"); +} + void linear_solve(NRSMat &a, NRMat *B, double *det, int n) { @@ -264,6 +374,15 @@ if(n<=0) n=a.nrows(); linear_solve_do(a,B?(*B)[0]:NULL,B?B->nrows() : 0, B?B->ncols():a.nrows(),det,n); } +void linear_solve(NRSMat &a, NRMat *B, float *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 symmetric 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, double *det, int n) { @@ -274,10 +393,21 @@ if(n<=0) n=a.nrows(); linear_solve_do(a,&B[0],1,a.nrows(),det,n); } +void linear_solve(NRSMat &a, NRVec &B, float *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 symmetric linear_solve()"); + B.copyonwrite(); +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); +extern "C" void FORNAME(chpsv)(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) { @@ -314,6 +444,41 @@ static void linear_solve_do(NRSMat > &a, std::complex 0 && b) laerror("singular matrix in linear_solve(SMat&, Mat*, complex*"); } +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) { @@ -324,6 +489,15 @@ if(n<=0) n=a.nrows(); linear_solve_do(a,B?(*B)[0]:NULL,B?B->nrows() : 0, B?B->ncols():a.nrows(),det,n); } +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) { @@ -334,10 +508,21 @@ if(n<=0) n=a.nrows(); linear_solve_do(a,&B[0],1,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); +extern "C" void FORNAME(cgesv)(const FINT *N, const FINT *NRHS, std::complex *A, const FINT *LDA, + FINT *IPIV, std::complex *B, const FINT *LDB, FINT *INFO); // A will be overwritten, B will contain the solutions, A is nxn, B is rhs x n // Fortran zgesv sees the row-major A transposed, so A is transposed in place before and after, like in clapack_dgesv in noncblas.cc @@ -393,6 +578,58 @@ static void linear_solve_do(NRMat > &A, std::complex0 && B) laerror("singular matrix in zgesv"); } +static void linear_solve_do(NRMat > &A, std::complex *B, const int nrhs, const int ldb, std::complex *det, int n) +{ + FINT r, *ipiv; + int iswap=0; + + if(nrhs==0) std::cout<<"Warning: some zgesv implementations might skip LU decomposition when nrhs==0\n"; + if (n==A.nrows() && A.nrows() != A.ncols()) laerror("linear_solve() call for non-square matrix"); + A.copyonwrite(); + ipiv = new FINT[n]; + const int lda=A.ncols(); + std::complex *a= &A(0,0); + for(int i=1; i t=a[j*lda+i]; a[j*lda+i]=a[i*lda+j]; a[i*lda+j]=t;} + const FINT ntmp=n; + const FINT nrhstmp=nrhs; + const FINT ldatmp=lda; + const FINT ldbtmp=ldb; + FORNAME(cgesv)(&ntmp, &nrhstmp, a, &ldatmp, ipiv, B, &ldbtmp, &r); + for(int i=1; i t=a[j*lda+i]; a[j*lda+i]=a[i*lda+j]; a[i*lda+j]=t;} + if (r < 0) { + delete[] ipiv; + laerror("illegal argument in lapack_gesv"); + } + if (det && r==0) { + *det = 1.; + //take into account some numerical instabilities in zgesv for singular matrices + for (int i=0; i t=A[i][i]; if(!std::isfinite(t.real()) || !std::isfinite(t.imag()) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *=t;} + //find out whether ipiv are numbered from 0 or from 1 + int shift=1; + for (int i=0; in) + { + std::cout <<"IPIV["< > &A, NRMat< std::complex > *B, std::complex *det, int n) { if(n<=0) n=A.nrows(); //default - whole matrix @@ -401,6 +638,14 @@ 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(NRMat< std::complex > &A, NRMat< std::complex > *B, std::complex *det, int n) +{ +if(n<=0) n=A.nrows(); //default - whole matrix +if (n==A.nrows() && B && A.nrows() != B->ncols() || B && n>B->ncols() ||n>A.nrows()) laerror("incompatible matrices in 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(NRMat< std::complex > &A, NRVec< std::complex > &B, std::complex *det, int n) { if(n<=0) n=A.nrows(); //default - whole matrix @@ -409,13 +654,23 @@ B.copyonwrite(); linear_solve_do(A,&B[0],1,A.nrows(),det,n); } +void linear_solve(NRMat< std::complex > &A, NRVec< std::complex > &B, std::complex *det, int n) +{ +if(n<=0) n=A.nrows(); //default - whole matrix +if(n==A.nrows() && A.nrows() != B.size() || n>B.size()||n>A.nrows() ) laerror("incompatible matrices in linear_solve()"); +B.copyonwrite(); +linear_solve_do(A,&B[0],1,A.nrows(),det,n); +} + //other version of linear solver based on gesvx //------------------------------------------------------------------------------ extern "C" void FORNAME(zgesvx)(const char *fact, const char *trans, const FINT *n, const FINT *nrhs, std::complex *A, const FINT *lda, std::complex *AF, const FINT *ldaf, const FINT *ipiv, char *equed, double *R,double *C, std::complex *B, const FINT *ldb, std::complex *X, const FINT *ldx, double *rcond, double *ferr, double *berr, std::complex *work, double *rwork, FINT *info); +extern "C" void FORNAME(cgesvx)(const char *fact, const char *trans, const FINT *n, const FINT *nrhs, std::complex *A, const FINT *lda, std::complex *AF, const FINT *ldaf, const FINT *ipiv, char *equed, float *R,float *C, std::complex *B, const FINT *ldb, std::complex *X, const FINT *ldx, float *rcond, float *ferr, float *berr, std::complex *work, float *rwork, FINT *info); extern "C" void FORNAME(dgesvx)(const char *fact, const char *trans, const FINT *n, const FINT *nrhs, double *A, const FINT *lda, double *AF, const FINT *ldaf, const FINT *ipiv, char *equed, double *R,double *C, double *B, const FINT *ldb, double *X, const FINT *ldx, double *rcond, double *ferr, double *berr, double *work, FINT *iwork, FINT *info); +extern "C" void FORNAME(sgesvx)(const char *fact, const char *trans, const FINT *n, const FINT *nrhs, float *A, const FINT *lda, float *AF, const FINT *ldaf, const FINT *ipiv, char *equed, float *R,float *C, float *B, const FINT *ldb, float *X, const FINT *ldx, float *rcond, float *ferr, float *berr, float *work, FINT *iwork, FINT *info); //------------------------------------------------------------------------------ // solves set of linear equations using dgesvx // input: @@ -489,6 +744,66 @@ int linear_solve_x(NRMat &_A, double *_B, const int _rhsCount, c } return (int)info; } + +template<> +int linear_solve_x(NRMat &_A, float *_B, const int _rhsCount, const int _eqCount, const bool _eq, const bool _saveA, float *_rcond){ + const int A_rows = _A.nrows(); + const int A_cols = _A.ncols(); + + const char fact = LAPACK_FORTRANCASE(_eq?'E':'N'); + const char trans = LAPACK_FORTRANCASE('T');//because of c-order + char equed = LAPACK_FORTRANCASE('B');//if fact=='N' then equed is an output argument, therefore not declared as const + + if(_eqCount < 0 || _eqCount > A_rows || _eqCount > A_cols || _rhsCount < 0){ + laerror("linear_solve_x: invalid input matrices"); + } + + float *A; + if(_eq && !_saveA) _A.copyonwrite(); //A will be overwritten by the equilibrated matrix, detach before taking the pointer + float * const _A_data = (float*)_A; + + FINT info; + const FINT nrhs = _rhsCount; + const FINT n = _eqCount; + FINT lda = A_cols; + const FINT ldaf = lda; + + float rcond; + float ferr[nrhs], berr[nrhs], work[4*n]; + float R[n], C[n]; + + FINT *const iwork = new FINT[n]; + FINT *const ipiv = new FINT[n]; + + float *X = new float[n*nrhs]; + float *AF = new float[ldaf*n]; + + A = _A_data; + if(_eq){ + if(_saveA){//store the corresponding submatrix of _A (not needed provided fact=='N') + A = new float[n*n]; + int offset1 = 0;int offset2 = 0; + for(register int i=0;i >(NRMat > &_A, std: } return (int)info; } + +template<> +int linear_solve_x >(NRMat > &_A, std::complex *_B, const int _rhsCount, const int _eqCount, const bool _eq, const bool _saveA, float *_rcond){ + const int A_rows = _A.nrows(); + const int A_cols = _A.ncols(); + + const char fact = LAPACK_FORTRANCASE(_eq?'E':'N'); + const char trans = LAPACK_FORTRANCASE('T');//because of c-order + char equed = LAPACK_FORTRANCASE('B');//if fact=='N' then equed is an output argument, therefore not declared as const + + if(_eqCount < 0 || _eqCount > A_rows || _eqCount > A_cols || _rhsCount < 0){ + laerror("linear_solve_x: invalid input matrices"); + } + + 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; + const FINT nrhs = _rhsCount; + const FINT n = _eqCount; + FINT lda = A_cols; + const FINT ldaf = lda; + + float rcond; + float ferr[nrhs], berr[nrhs]; + float R[n], C[n], rwork[2*n]; + std::complex work[2*n]; + + FINT *const ipiv = new FINT[n]; + + std::complex *X = new std::complex[n*nrhs]; + std::complex *AF = new std::complex[ldaf*n]; + + A = _A_data; + if(_eq){ + if(_saveA){//store the corresponding submatrix of _A (not needed provided fact=='N') + A = new std::complex[n*n]; + int offset1 = 0;int offset2 = 0; + for(register int i=0;i B^TX^T = A^T @@ -618,6 +994,50 @@ int multiply_by_inverse(NRMat &_A, NRMat &_B, bool _useE return (int)info; } + +template<> +int multiply_by_inverse(NRMat &_A, NRMat &_B, bool _useEq, float *_rcond){ + + const FINT n = _A.nrows(); + const FINT m = _A.ncols(); + if(n != m || n != _B.nrows() || n != _B.ncols()){ + laerror("multiply_by_inverse: incompatible matrices"); + } + + const char fact = _useEq?'E':'N'; + 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 + const int n2 = n*n; + + _B.copyonwrite();//even if fact='N', call copyonwrite because the solution is going to be stored in _B; before taking the pointer + float * const B = (float*)_B; + //_A is the right hand side for gesvx, which scales it if equilibrated: work on a copy then, _A is input only + float * const A = _useEq ? new float[n2] : (float*)_A; + if(_useEq) memcpy(A,(float*)_A,n2*sizeof(float)); + + FINT info; + float rcond; + float ferr[n], berr[n], work[4*n]; + float R[n], C[n]; + + FINT *const iwork = new FINT[n]; + FINT *const ipiv = new FINT[n]; + + float *X = new float[n2]; + float *AF = new float[n2]; + + FORNAME(sgesvx)(&fact, &trans, &n, &n, B, &n, AF, &n, &ipiv[0], &equed, &R[0], &C[0], A, &n, X, &n, &rcond, ferr, berr, work, iwork, &info); + + + if(_rcond)*_rcond = rcond; + cblas_scopy(n2, X, 1, B, 1);//store the solution + + delete[] iwork;delete[] ipiv; + delete[] AF;delete[] X; + if(_useEq) delete[] A; + + return (int)info; +} //------------------------------------------------------------------------------ // for given square matrices A, B computes X = AB^{-1} as follows // XB = A => B^TX^T = A^T @@ -673,15 +1093,63 @@ int multiply_by_inverse >(NRMat > &_A, return (int)info; } + +template<> +int multiply_by_inverse >(NRMat > &_A, NRMat > &_B, bool _useEq, float *_rcond){ + + const FINT n = _A.nrows(); + const FINT m = _A.ncols(); + if(n != m || n != _B.nrows() || n != _B.ncols()){ + laerror("multiply_by_inverse: incompatible matrices"); + } + const int n2 = n*n; + + const char fact = _useEq?'E':'N'; + 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 + + _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; + //_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; + float rcond; + float ferr[n], berr[n]; + float R[n], C[n], rwork[2*n]; + std::complex work[2*n]; + + FINT *const ipiv = new FINT[n]; + + std::complex *X = new std::complex[n2]; + std::complex *AF = new std::complex[n2]; + + FORNAME(cgesvx)(&fact, &trans, &n, &n, B, &n, AF, &n, &ipiv[0], &equed, &R[0], &C[0], A, &n, X, &n, &rcond, ferr, berr, work, rwork, &info); + + + if(_rcond)*_rcond = rcond; + cblas_ccopy(n2, X, 1, B, 1);//store the solution + + delete[] ipiv; + delete[] AF;delete[] X; + if(_useEq) delete[] A; + + return (int)info; +} //------------------------------------------------------------------------------ extern "C" void FORNAME(dsyev)(const char *JOBZ, const char *UPLO, const FINT *N, double *A, const FINT *LDA, double *W, double *WORK, const FINT *LWORK, FINT *INFO); +extern "C" void FORNAME(ssyev)(const char *JOBZ, const char *UPLO, const FINT *N, + float *A, const FINT *LDA, float *W, float *WORK, const FINT *LWORK, FINT *INFO); extern "C" void FORNAME(dsygv)(const FINT *ITYPE, const char *JOBZ, const char *UPLO, const FINT *N, double *A, const FINT *LDA, double *B, const FINT *LDB, double *W, double *WORK, const FINT *LWORK, FINT *INFO); +extern "C" void FORNAME(ssygv)(const FINT *ITYPE, const char *JOBZ, const char *UPLO, const FINT *N, + float *A, const FINT *LDA, float *B, const FINT *LDB, float *W, float *WORK, const FINT *LWORK, FINT *INFO); // a will contain eigenvectors (columns if corder==1), w eigenvalues @@ -738,13 +1206,70 @@ void diagonalize(NRMat &a, NRVec &w, const bool eivec, if (r > 0) laerror("convergence problem in sygv/syev in diagonalize()"); } +void diagonalize(NRMat &a, NRVec &w, const bool eivec, + const bool corder, int n, NRMat *b, const int itype) +{ + FINT m = a.nrows(); + if (m != a.ncols()) laerror("diagonalize() call with non-square matrix"); + if (a.nrows() != w.size()) + laerror("inconsistent dimension of eigenvalue vector in diagonalize()"); + if(n==0) n=m; + if(n<0||n>m) laerror("actual dimension out of range in diagonalize"); + if(b) if(n>b->nrows() || n> b->ncols()) laerror("wrong B matrix dimension in diagonalize"); + + a.copyonwrite(); + w.copyonwrite(); + if(b) b->copyonwrite(); + + FINT r = 0; + char U =LAPACK_FORTRANCASE('u'); + char vectors = LAPACK_FORTRANCASE('v'); + if (!eivec) vectors = LAPACK_FORTRANCASE('n'); + FINT LWORK = -1; + float WORKX; + FINT ldb=0; if(b) ldb=b->ncols(); + +#ifdef FORINT + const FINT itypetmp = itype; + FINT ntmp = n; + // First call is to determine size of workspace + if(b) FORNAME(ssygv)(&itypetmp,&vectors, &U, &ntmp, a, &m, *b, &ldb, w, &WORKX, &LWORK, &r ); + else FORNAME(ssyev)(&vectors, &U, &ntmp, a, &m, w, &WORKX, &LWORK, &r ); +#else + // First call is to determine size of workspace + if(b) FORNAME(ssygv)(&itype,&vectors, &U, &n, a, &m, *b, &ldb, w, &WORKX, &LWORK, &r ); + else FORNAME(ssyev)(&vectors, &U, &n, a, &m, w, &WORKX, &LWORK, &r ); +#endif + + LWORK = (FINT)WORKX; + float *WORK = new float[LWORK]; + +#ifdef FORINT + if(b) FORNAME(ssygv)(&itypetmp,&vectors, &U, &ntmp, a, &m, *b, &ldb, w, WORK, &LWORK, &r ); + else FORNAME(ssyev)(&vectors, &U, &ntmp, a, &m, w, WORK, &LWORK, &r ); +#else + if(b) FORNAME(ssygv)(&itype,&vectors, &U, &n, a, &m, *b,&ldb, w, WORK, &LWORK, &r ); + else FORNAME(ssyev)(&vectors, &U, &n, a, &m, w, WORK, &LWORK, &r ); +#endif + + delete[] WORK; + if (LAPACK_FORTRANCASE(vectors) == LAPACK_FORTRANCASE('v') && corder) a.transposeme(n); + + if (r < 0) laerror("illegal argument in sygv/syev in diagonalize()"); + if (r > 0) laerror("convergence problem in sygv/syev in diagonalize()"); +} + extern "C" void FORNAME(zheev)(const char *JOBZ, const char *UPLO, const FINT *N, std::complex *A, const FINT *LDA, double *W, std::complex *WORK, const FINT *LWORK, double *RWORK, FINT *INFO); +extern "C" void FORNAME(cheev)(const char *JOBZ, const char *UPLO, const FINT *N, + std::complex *A, const FINT *LDA, float *W, std::complex *WORK, const FINT *LWORK, float *RWORK, FINT *INFO); extern "C" void FORNAME(zhegv)(const FINT *ITYPE, const char *JOBZ, const char *UPLO, const FINT *N, std::complex *A, const FINT *LDA, std::complex *B, const FINT *LDB, double *W, std::complex *WORK, const FINT *LWORK, double *RWORK, FINT *INFO); +extern "C" void FORNAME(chegv)(const FINT *ITYPE, const char *JOBZ, const char *UPLO, const FINT *N, + std::complex *A, const FINT *LDA, std::complex *B, const FINT *LDB, float *W, std::complex *WORK, const FINT *LWORK, float *RWORK, FINT *INFO); // a will contain eigenvectors (columns if corder==1), w eigenvalues @@ -802,13 +1327,71 @@ void diagonalize(NRMat > &a, NRVec &w, const bool e if (r > 0) laerror("convergence problem in hegv/heev in diagonalize()"); } +void diagonalize(NRMat > &a, NRVec &w, const bool eivec, + const bool corder, int n, NRMat > *b, const int itype) +{ + FINT m = a.nrows(); + if (m != a.ncols()) laerror("diagonalize() call with non-square matrix"); + if (a.nrows() != w.size()) + laerror("inconsistent dimension of eigenvalue vector in diagonalize()"); + if(n==0) n=m; + if(n<0||n>m) laerror("actual dimension out of range in diagonalize"); + if(b) if(n>b->nrows() || n> b->ncols()) laerror("wrong B matrix dimension in diagonalize"); + + a.copyonwrite(); + w.copyonwrite(); + if(b) b->copyonwrite(); + + FINT r = 0; + char U =LAPACK_FORTRANCASE('U'); + char vectors = LAPACK_FORTRANCASE('V'); + if (!eivec) vectors = LAPACK_FORTRANCASE('n'); + FINT LWORK = -1; + std::complex WORKX; + FINT ldb=0; if(b) ldb=b->ncols(); + + // First call is to determine size of workspace + float *RWORK = new float[3*n+2]; +#ifdef FORINT + const FINT itypetmp = itype; + FINT ntmp = n; + if(b) FORNAME(chegv)(&itypetmp,&vectors, &U, &ntmp, a, &m, *b, &ldb, w, &WORKX, &LWORK, RWORK, &r ); + else FORNAME(cheev)(&vectors, &U, &ntmp, a, &m, w, &WORKX, &LWORK, RWORK, &r ); +#else + if(b) FORNAME(chegv)(&itype,&vectors, &U, &n, a, &m, *b, &ldb, w, &WORKX, &LWORK, RWORK, &r ); + else FORNAME(cheev)(&vectors, &U, &n, a, &m, w, &WORKX, &LWORK, RWORK, &r ); +#endif + + LWORK = (FINT)WORKX.real(); + std::complex *WORK = new std::complex[LWORK]; + +#ifdef FORINT + if(b) FORNAME(chegv)(&itypetmp,&vectors, &U, &ntmp, a, &m, *b, &ldb, w, WORK, &LWORK, RWORK, &r ); + else FORNAME(cheev)(&vectors, &U, &ntmp, a, &m, w, WORK, &LWORK, RWORK, &r ); +#else + if(b) FORNAME(chegv)(&itype,&vectors, &U, &n, a, &m, *b,&ldb, w, WORK, &LWORK, RWORK, &r ); + else FORNAME(cheev)(&vectors, &U, &n, a, &m, w, WORK, &LWORK, RWORK, &r ); +#endif + + delete[] WORK; + delete[] RWORK; + if (LAPACK_FORTRANCASE(vectors) == LAPACK_FORTRANCASE('v') && corder) {a.transposeme(n); a.conjugateme();} + + if (r < 0) laerror("illegal argument in hegv/heev in diagonalize()"); + if (r > 0) laerror("convergence problem in hegv/heev in diagonalize()"); +} + extern "C" void FORNAME(dspev)(const char *JOBZ, const char *UPLO, const FINT *N, double *AP, double *W, double *Z, const FINT *LDZ, double *WORK, FINT *INFO); +extern "C" void FORNAME(sspev)(const char *JOBZ, const char *UPLO, const FINT *N, + float *AP, float *W, float *Z, const FINT *LDZ, float *WORK, FINT *INFO); extern "C" void FORNAME(dspgv)(const FINT *ITYPE, const char *JOBZ, const char *UPLO, const FINT *N, double *AP, double *BP, double *W, double *Z, const FINT *LDZ, double *WORK, FINT *INFO); +extern "C" void FORNAME(sspgv)(const FINT *ITYPE, const char *JOBZ, const char *UPLO, const FINT *N, + float *AP, float *BP, float *W, float *Z, const FINT *LDZ, float *WORK, FINT *INFO); // v will contain eigenvectors, w eigenvalues @@ -849,12 +1432,53 @@ void diagonalize(NRSMat &a, NRVec &w, NRMat *v, if (r > 0) laerror("convergence problem in spgv/spev in diagonalize()"); } +void diagonalize(NRSMat &a, NRVec &w, NRMat *v, + const bool corder, int n, NRSMat *b, const int itype) +{ + if(n<=0) n = a.nrows(); + if (v) if (v->nrows() != v ->ncols() || n > v->nrows() || n > a.nrows()) + laerror("diagonalize() call with inconsistent dimensions"); + if (n==a.nrows() && n != w.size() || n>w.size()) laerror("inconsistent dimension of eigenvalue vector"); + + if(b) if(n>b->nrows() || n> b->ncols()) laerror("wrong B matrix dimension in diagonalize"); + + a.copyonwrite(); + w.copyonwrite(); + if(v) v->copyonwrite(); + if(b) b->copyonwrite(); + + FINT r = 0; + char U = LAPACK_FORTRANCASE('u'); + char job = LAPACK_FORTRANCASE(v ? 'v' : 'n'); + + float *WORK = new float[3*n]; + FINT ldv=v?v->ncols():n; +#ifdef FORINT + const FINT itypetmp = itype; + FINT ntmp = n; + if(b) FORNAME(sspgv)(&itypetmp,&job, &U, &ntmp, a, *b, w, v?(*v)[0]:(float *)0, &ldv, WORK, &r ); + else FORNAME(sspev)(&job, &U, &ntmp, a, w, v?(*v)[0]:(float *)0, &ldv, WORK, &r ); +#else + if(b) FORNAME(sspgv)(&itype,&job, &U, &n, a, *b, w, v?(*v)[0]:(float *)0, &ldv, WORK, &r ); + else FORNAME(sspev)(&job, &U, &n, a, w, v?(*v)[0]:(float *)0, &ldv, WORK, &r ); +#endif + delete[] WORK; + if (v && corder) v->transposeme(n); + + if (r < 0) laerror("illegal argument in spgv/spev in diagonalize()"); + if (r > 0) laerror("convergence problem in spgv/spev in diagonalize()"); +} + extern "C" void FORNAME(zhpev)(const char *JOBZ, const char *UPLO, const FINT *N, std::complex *AP, double *W, std::complex *Z, const FINT *LDZ, std::complex *WORK, double *RWORK, FINT *INFO); +extern "C" void FORNAME(chpev)(const char *JOBZ, const char *UPLO, const FINT *N, + std::complex *AP, float *W, std::complex *Z, const FINT *LDZ, std::complex *WORK, float *RWORK, FINT *INFO); extern "C" void FORNAME(zhpgv)(const FINT *ITYPE, const char *JOBZ, const char *UPLO, const FINT *N, std::complex *AP, std::complex *BP, double *W, std::complex *Z, const FINT *LDZ, std::complex *WORK, double *RWORK, FINT *INFO); +extern "C" void FORNAME(chpgv)(const FINT *ITYPE, const char *JOBZ, const char *UPLO, const FINT *N, + std::complex *AP, std::complex *BP, float *W, std::complex *Z, const FINT *LDZ, std::complex *WORK, float *RWORK, FINT *INFO); // v will contain eigenvectors, w eigenvalues @@ -897,11 +1521,53 @@ void diagonalize(NRSMat > &a, NRVec &w, NRMat 0) laerror("convergence problem in hpgv/hpev in diagonalize()"); } +void diagonalize(NRSMat > &a, NRVec &w, NRMat > *v, + const bool corder, int n, NRSMat > *b, const int itype) +{ + if(n<=0) n = a.nrows(); + if (v) if (v->nrows() != v ->ncols() || n > v->nrows() || n > a.nrows()) + laerror("diagonalize() call with inconsistent dimensions"); + if (n==a.nrows() && n != w.size() || n>w.size()) laerror("inconsistent dimension of eigenvalue vector"); + + if(b) if(n>b->nrows() || n> b->ncols()) laerror("wrong B matrix dimension in diagonalize"); + + a.copyonwrite(); + w.copyonwrite(); + if(v) v->copyonwrite(); + if(b) b->copyonwrite(); + + FINT r = 0; + char U = LAPACK_FORTRANCASE('u'); + char job = LAPACK_FORTRANCASE(v ? 'v' : 'n'); + + std::complex *WORK = new std::complex[2*n]; + float *RWORK = new float[3*n]; + FINT ldv=v?v->ncols():n; +#ifdef FORINT + const FINT itypetmp = itype; + FINT ntmp = n; + if(b) FORNAME(chpgv)(&itypetmp,&job, &U, &ntmp, a, *b, w, v?(*v)[0]:(std::complex *)0, &ldv, WORK, RWORK, &r ); + else FORNAME(chpev)(&job, &U, &ntmp, a, w, v?(*v)[0]:(std::complex *)0, &ldv, WORK, RWORK, &r ); +#else + if(b) FORNAME(chpgv)(&itype,&job, &U, &n, a, *b, w, v?(*v)[0]:(std::complex *)0, &ldv, WORK, RWORK, &r ); + else FORNAME(chpev)(&job, &U, &n, a, w, v?(*v)[0]:(std::complex *)0, &ldv, WORK, RWORK, &r ); +#endif + delete[] WORK; + delete[] RWORK; + 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()"); +} + extern "C" void FORNAME(dgesvd)(const char *JOBU, const char *JOBVT, const FINT *M, const FINT *N, double *A, const FINT *LDA, double *S, double *U, const FINT *LDU, double *VT, const FINT *LDVT, double *WORK, const FINT *LWORK, FINT *INFO ); +extern "C" void FORNAME(sgesvd)(const char *JOBU, const char *JOBVT, const FINT *M, + const FINT *N, float *A, const FINT *LDA, float *S, float *U, const FINT *LDU, + float *VT, const FINT *LDVT, float *WORK, const FINT *LWORK, FINT *INFO ); //normally in v returns vtransposed for default vnotdagger=0 void singular_decomposition(NRMat &a, NRMat *u, NRVec &s, @@ -960,6 +1626,62 @@ void singular_decomposition(NRMat &a, NRMat *u, NRVec &s if (r > 0) laerror("convergence problem in gesvd() of singular_decomposition()"); } +void singular_decomposition(NRMat &a, NRMat *u, NRVec &s, + NRMat *v, const bool vnotdagger, int m, int n) +{ + FINT m0 = a.nrows(); + FINT n0 = a.ncols(); + if(m<=0) m=(int)m0; + if(n<=0) n=(int)n0; + if(n>n0 || m>m0) laerror("bad dimension in singular_decomposition"); + if (u) if (m > u->nrows() || s.size()> u->ncols()) + laerror("inconsistent dimension of U Mat in singular_decomposition()"); + if (s.size() < m && s.size() < n) + laerror("inconsistent dimension of S Vec in singular_decomposition()"); + if (v) if (s.size() > v->nrows() || n > v->ncols()) + laerror("inconsistent dimension of V Mat in singular_decomposition()"); + + a.copyonwrite(); + s.copyonwrite(); + if (u) u->copyonwrite(); + if (v) v->copyonwrite(); + + // C-order (transposed) input and swap u,v matrices, + // v should be transposed at the end + char jobu = u ? (u->nrows()==u->ncols() ? 'A' : 'S') : 'N'; + char jobv = v ? (v->nrows()==v->ncols() ? 'A' : 'S') : 'N'; + + float work0; + FINT lwork = -1; + FINT r; + + FINT lda=a.ncols(); + FINT ldu= u ? u->ncols():0; + FINT ldv= v ? v->ncols():0; +#ifdef FORINT + FINT ntmp = n; + FINT mtmp = m; + FORNAME(sgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, &work0, &lwork, &r); +#else + FORNAME(sgesvd)(&jobv, &jobu, &n, &m, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, &work0, &lwork, &r); +#endif + + lwork = (FINT) work0; + float *work = new float[lwork]; + +#ifdef FORINT + FORNAME(sgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, work, &lwork, &r); +#else + FORNAME(sgesvd)(&jobv, &jobu, &n, &m, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, work, &lwork, &r); +#endif + + delete[] work; + if (v && vnotdagger) v->transposeme(n); + + if (r < 0) laerror("illegal argument in gesvd() of singular_decomposition()"); + if (r > 0) laerror("convergence problem in gesvd() of singular_decomposition()"); +} + @@ -968,6 +1690,9 @@ void singular_decomposition(NRMat &a, NRMat *u, NRVec &s extern "C" void FORNAME(zgesvd)(const char *JOBU, const char *JOBVT, const FINT *M, const FINT *N, std::complex *A, const FINT *LDA, double *S, std::complex *U, const FINT *LDU, std::complex *VT, const FINT *LDVT, std::complex *WORK, const FINT *LWORK, double *RWORK, FINT *INFO ); +extern "C" void FORNAME(cgesvd)(const char *JOBU, const char *JOBVT, const FINT *M, + const FINT *N, std::complex *A, const FINT *LDA, float *S, std::complex *U, const FINT *LDU, + std::complex *VT, const FINT *LDVT, std::complex *WORK, const FINT *LWORK, float *RWORK, FINT *INFO ); void singular_decomposition(NRMat > &a, NRMat > *u, NRVec &s, NRMat > *v, const bool vnotdagger, int m, int n) @@ -1028,6 +1753,65 @@ void singular_decomposition(NRMat > &a, NRMat 0) laerror("convergence problem in gesvd() of singular_decomposition()"); } +void singular_decomposition(NRMat > &a, NRMat > *u, NRVec &s, + NRMat > *v, const bool vnotdagger, int m, int n) +{ + FINT m0 = a.nrows(); + FINT n0 = a.ncols(); + if(m<=0) m=(int)m0; + if(n<=0) n=(int)n0; + if(n>n0 || m>m0) laerror("bad dimension in singular_decomposition"); + if (u) if (m > u->nrows() || s.size()> u->ncols()) + laerror("inconsistent dimension of U Mat in singular_decomposition()"); + if (s.size() < m && s.size() < n) + laerror("inconsistent dimension of S Vec in singular_decomposition()"); + if (v) if (s.size() > v->nrows() || n > v->ncols()) + laerror("inconsistent dimension of V Mat in singular_decomposition()"); + + int nmin = ncopyonwrite(); + if (v) v->copyonwrite(); + + // C-order (transposed) input and swap u,v matrices, + // v should be transposed at the end + char jobu = u ? (u->nrows()==u->ncols() ? 'A' : 'S') : 'N'; + char jobv = v ? (v->nrows()==v->ncols() ? 'A' : 'S') : 'N'; + + std::complex work0; + FINT lwork = -1; + FINT r; + float *rwork = new float[5*nmin]; + + FINT lda=a.ncols(); + FINT ldu= u ? u->ncols():0; + FINT ldv= v ? v->ncols():0; +#ifdef FORINT + FINT ntmp = n; + FINT mtmp = m; + FORNAME(cgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, &work0, &lwork, rwork, &r); +#else + FORNAME(cgesvd)(&jobv, &jobu, &n, &m, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, &work0, &lwork, rwork, &r); +#endif + + lwork = (FINT) work0.real(); + std::complex *work = new std::complex[lwork]; + +#ifdef FORINT + FORNAME(cgesvd)(&jobv, &jobu, &ntmp, &mtmp, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, work, &lwork, rwork, &r); +#else + FORNAME(cgesvd)(&jobv, &jobu, &n, &m, a, &lda, s, v?(*v)[0]:0, &ldv, u?(*u)[0]:0, &ldu, work, &lwork, rwork, &r); +#endif + + delete[] work; + delete[] rwork; + if (v && vnotdagger) {v->transposeme(n); v->conjugateme();} + + if (r < 0) laerror("illegal argument in gesvd() of singular_decomposition()"); + if (r > 0) laerror("convergence problem in gesvd() of singular_decomposition()"); +} + //QR decomposition //extern "C" void FORNAME(dgeqrf)(const int *M, const int *N, double *A, const int *LDA, double *TAU, double *WORK, int *LWORK, int *INFO); @@ -1036,24 +1820,40 @@ void singular_decomposition(NRMat > &a, NRMat *A, const FINT *LDA, std::complex *W, std::complex *VL, const FINT *LDVL, std::complex *VR, const FINT *LDVR, std::complex *WORK, const FINT *LWORK, double *RWORK, FINT *INFO ); +extern "C" void FORNAME(cgeev)(const char *JOBVL, const char *JOBVR, const FINT *N, + std::complex *A, const FINT *LDA, std::complex *W, std::complex *VL, const FINT *LDVL, + std::complex *VR, const FINT *LDVR, std::complex *WORK, const FINT *LWORK, + float *RWORK, FINT *INFO ); extern "C" void FORNAME(zggev)(const char *JOBVL, const char *JOBVR, const FINT *N, std::complex *A, const FINT *LDA, std::complex *B, const FINT *LDB, std::complex *W, std::complex *WBETA, std::complex *VL, const FINT *LDVL, std::complex *VR, const FINT *LDVR, std::complex *WORK, const FINT *LWORK, double *RWORK, FINT *INFO ); +extern "C" void FORNAME(cggev)(const char *JOBVL, const char *JOBVR, const FINT *N, + std::complex *A, const FINT *LDA, std::complex *B, const FINT *LDB, std::complex *W, std::complex *WBETA, + std::complex *VL, const FINT *LDVL, std::complex *VR, const FINT *LDVR, + std::complex *WORK, const FINT *LWORK, float *RWORK, FINT *INFO ); //tridiagonal extern "C" void FORNAME(dsteqr)(const char *compz, const FINT *n, double *d, double *d1, double *z, const FINT *ldz, double *work, FINT *info); +extern "C" void FORNAME(ssteqr)(const char *compz, const FINT *n, float *d, float *d1, float *z, const FINT *ldz, float *work, FINT *info); void diagonalize3(NRVec &d, NRVec &d1, NRMat *v, const bool corder, int n0) { @@ -1078,47 +1878,76 @@ void diagonalize3(NRVec &d, NRVec &d1, NRMat *v, const b if (info > 0) laerror("convergence problem in dsteqr() fo diagonalize3"); } +void diagonalize3(NRVec &d, NRVec &d1, NRMat *v, const bool corder, int n0) +{ + FINT n = n0? n0: d.size(); + if(d1.size()nrows()||n!=v->ncols()) laerror("inconsistent dimensions in diagonalize3");} + d.copyonwrite(); + d1.copyonwrite(); + if(v) v->copyonwrite(); + + FINT ldz=n; + FINT info; + int worksize=2*n-2; + if(worksize<1) worksize=1; + float *work = new float[worksize]; + char job=v?'i':'n'; + FORNAME(ssteqr)(&job,&n,&d[0],&d1[0],(v?&(*v)(0,0):NULL),&ldz,work,&info); + delete[] work; + if (v && corder) v->transposeme(); -//statics for sorting + if (info < 0) laerror("illegal argument in dsteqr() fo diagonalize3"); + if (info > 0) laerror("convergence problem in dsteqr() fo diagonalize3"); +} + + +//statics for sorting, one set per floating type (variable templates) static int *gdperm; -static double *gdwr, *gdwi, *gdbeta; +template static T *gdwr; +template static T *gdwi; +template static T *gdbeta; //compare methods -static double realonly(const int i, const int j) +template +static T realonly(const int i, const int j) { -if(gdbeta) +if(gdbeta) { - if(gdbeta[i]==0. && gdbeta[j]!=0) return 1.; - if(gdbeta[j]==0. && gdbeta[i]!=0) return -1.; - if(gdbeta[i]==0. && gdbeta[j]==0) return 0.; - double tmp = gdwr[i]/gdbeta[i]-gdwr[j]/gdbeta[j]; + if(gdbeta[i]==0. && gdbeta[j]!=0) return 1.; + if(gdbeta[j]==0. && gdbeta[i]!=0) return -1.; + if(gdbeta[i]==0. && gdbeta[j]==0) return 0.; + T tmp = gdwr[i]/gdbeta[i]-gdwr[j]/gdbeta[j]; if(tmp) return tmp; - return gdwi[j]/gdbeta[j]-gdwi[i]/gdbeta[i]; + return gdwi[j]/gdbeta[j]-gdwi[i]/gdbeta[i]; } //else -double tmp = gdwr[i]-gdwr[j]; +T tmp = gdwr[i]-gdwr[j]; if(tmp) return tmp; -return gdwi[j]-gdwi[i]; +return gdwi[j]-gdwi[i]; } -static double realfirst(const int i, const int j) +template +static T realfirst(const int i, const int j) { -if(gdwi[i] && ! gdwi[j]) return 1.; -if(!gdwi[i] && gdwi[j]) return -1.; -return realonly(i,j); +if(gdwi[i] && ! gdwi[j]) return 1.; +if(!gdwi[i] && gdwi[j]) return -1.; +return realonly(i,j); } -static double (* gdcompar[2])(const int, const int) = {&realonly, &realfirst}; +template +static T (* gdcompar[2])(const int, const int) = {&realonly, &realfirst}; //swap method +template static void gdswap(const int i, const int j) { -double tmp; +T tmp; int itmp; itmp=gdperm[i]; gdperm[i]=gdperm[j]; gdperm[j]=itmp; -tmp=gdwr[i]; gdwr[i]=gdwr[j]; gdwr[j]=tmp; -tmp=gdwi[i]; gdwi[i]=gdwi[j]; gdwi[j]=tmp; -if(gdbeta) {tmp=gdbeta[i]; gdbeta[i]=gdbeta[j]; gdbeta[j]=tmp;} +tmp=gdwr[i]; gdwr[i]=gdwr[j]; gdwr[j]=tmp; +tmp=gdwi[i]; gdwi[i]=gdwi[j]; gdwi[j]=tmp; +if(gdbeta) {tmp=gdbeta[i]; gdbeta[i]=gdbeta[j]; gdbeta[j]=tmp;} } @@ -1273,9 +2102,9 @@ void gdiagonalize(NRMat &a, NRVec &wr, NRVec &wi, NRVec perm(n); for(int i=0; i = *beta; else gdbeta = NULL; + gdwr = wr, gdwi = wi; + genqsort(0,n-1,gdcompar[sorttype-1],gdswap); if(vl) { for(int i=0; i &a, NRVec &wr, NRVec &wi, } +void gdiagonalize(NRMat &a, NRVec &wr, NRVec &wi, + NRMat *vl, NRMat *vr, const bool corder, int n, + const int sorttype, const int biorthonormalize, + NRMat *b, NRVec *beta) +{ + if(n<=0) {n = a.nrows(); if(a.ncols()!=a.nrows() ) laerror("gdiagonalize() call for a non-square matrix");} + if (n > a.ncols() || n>a.nrows() ) laerror("gdiagonalize() of too big submatrix"); + if (n > wr.size()) + laerror("inconsistent dimension of eigen vector in gdiagonalize()"); + if (vl) if (n > vl->nrows() || n > vl->ncols()) + laerror("inconsistent dimension of vl in gdiagonalize()"); + if (vr) if (n > vr->nrows() || n > vr->ncols()) + laerror("inconsistent dimension of vr in gdiagonalize()"); + if (beta) if(n > beta ->size()) laerror("inconsistent dimension of beta in gdiagonalize()"); + if(b) if(n > b->nrows() || n > b->ncols()) + laerror("inconsistent dimension of b in gdiagonalize()"); + if(b && !beta || beta && !b) laerror("missing array for generalized diagonalization"); + + a.copyonwrite(); + wr.copyonwrite(); + wi.copyonwrite(); + if (vl) vl->copyonwrite(); + if (vr) vr->copyonwrite(); + if (beta) beta->copyonwrite(); + if (b) b->copyonwrite(); + + char jobvl = LAPACK_FORTRANCASE(vl ? 'v' : 'n'); + char jobvr = LAPACK_FORTRANCASE(vr ? 'v' : 'n'); + float work0; + FINT lwork = -1; + FINT r; + FINT lda=a.ncols(); + FINT ldb=0; + if(b) ldb=b->ncols(); + FINT ldvl= vl?vl->ncols():lda; + FINT ldvr= vr?vr->ncols():lda; + +#ifdef FORINT + FINT ntmp = n; + if(b) FORNAME(sggev)(&jobvr, &jobvl, &ntmp, a, &lda, *b, &ldb, wr, wi, *beta, vr?(*vr)[0]:(float *)0, + &ldvr, vl?(*vl)[0]:(float *)0, &ldvl, &work0, &lwork, &r); + else FORNAME(sgeev)(&jobvr, &jobvl, &ntmp, a, &lda, wr, wi, vr?(*vr)[0]:(float *)0, + &ldvr, vl?(*vl)[0]:(float *)0, &ldvl, &work0, &lwork, &r); +#else + if(b) FORNAME(sggev)(&jobvr, &jobvl, &n, a, &lda, *b, &ldb, wr, wi, *beta, vr?(*vr)[0]:(float *)0, + &ldvr, vl?(*vl)[0]:(float *)0, &ldvl, &work0, &lwork, &r); + else FORNAME(sgeev)(&jobvr, &jobvl, &n, a, &lda, wr, wi, vr?(*vr)[0]:(float *)0, + &ldvr, vl?(*vl)[0]:(float *)0, &ldvl, &work0, &lwork, &r); +#endif + + lwork = (FINT) work0; + float *work = new float[lwork]; + +#ifdef FORINT + if(b) FORNAME(sggev)(&jobvr, &jobvl, &ntmp, a, &lda, *b, &ldb, wr, wi, *beta, vr?(*vr)[0]:(float *)0, + &ldvr, vl?(*vl)[0]:(float *)0, &ldvl, work, &lwork, &r); + else FORNAME(sgeev)(&jobvr, &jobvl, &ntmp, a, &lda, wr, wi, vr?(*vr)[0]:(float *)0, + &ldvr, vl?(*vl)[0]:(float *)0, &ldvl, work, &lwork, &r); +#else + if(b) FORNAME(sggev)(&jobvr, &jobvl, &n, a, &lda, *b, &ldb, wr, wi, *beta, vr?(*vr)[0]:(float *)0, + &ldvr, vl?(*vl)[0]:(float *)0, &ldvl, work, &lwork, &r); + else FORNAME(sgeev)(&jobvr, &jobvl, &n, a, &lda, wr, wi, vr?(*vr)[0]:(float *)0, + &ldvr, vl?(*vl)[0]:(float *)0, &ldvl, work, &lwork, &r); +#endif + delete[] work; + + + if (r < 0) laerror("illegal argument in ggev/geev in gdiagonalize()"); + if (r > 0) laerror("convergence problem in ggev/geev in gdiagonalize()"); + +//std::cout <<"TEST dgeev\n"<0) + { + NRVec perm(n); + for(int i=0; i = *beta; else gdbeta = NULL; + gdwr = wr, gdwi = wi; + genqsort(0,n-1,gdcompar[sorttype-1],gdswap); + if(vl) + { + for(int i=0; itransposeme(n); + if (vr) vr->transposeme(n); + } + +} + //most general complex routine @@ -1402,6 +2404,109 @@ void gdiagonalize(NRMat > &a, NRVec< std::complex > } +template<> +void gdiagonalize(NRMat > &a, NRVec< std::complex > &w, + NRMat< std::complex >*vl, NRMat< std::complex > *vr, + const bool corder, int n, const int sorttype, const int biorthonormalize, + NRMat > *b, NRVec > *beta) +{ + + if(n<=0) {n = a.nrows(); if(a.ncols()!=a.nrows() ) laerror("gdiagonalize() call for a non-square matrix");} + if (n > a.ncols() || n>a.nrows() ) laerror("gdiagonalize() of too big submatrix"); + if (n > w.size()) + laerror("inconsistent dimension of eigen vector in gdiagonalize()"); + if (vl) if (n > vl->nrows() || n > vl->ncols()) + laerror("inconsistent dimension of vl in gdiagonalize()"); + if (vr) if (n > vr->nrows() || n > vr->ncols()) + laerror("inconsistent dimension of vr in gdiagonalize()"); + if (beta) if(n > beta ->size()) laerror("inconsistent dimension of beta in gdiagonalize()"); + if(b) if(n > b->nrows() || n > b->ncols()) + laerror("inconsistent dimension of b in gdiagonalize()"); + if(b && !beta || beta && !b) laerror("missing array for generalized diagonalization"); + + a.copyonwrite(); + w.copyonwrite(); + if (vl) vl->copyonwrite(); + if (vr) vr->copyonwrite(); + if (beta) beta->copyonwrite(); + if (b) b->copyonwrite(); + + char jobvl = LAPACK_FORTRANCASE(vl ? 'v' : 'n'); + char jobvr = LAPACK_FORTRANCASE(vr ? 'v' : 'n'); + std::complex work0; + FINT lwork = -1; + FINT r; + FINT lda=a.ncols(); + FINT ldb=0; + if(b) ldb=b->ncols(); + FINT ldvl= vl?vl->ncols():lda; + FINT ldvr= vr?vr->ncols():lda; + + float *rwork = new float[n*(b?8:2)]; + +#ifdef FORINT + FINT ntmp = n; + if(b) FORNAME(cggev)(&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(cgeev)(&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(cggev)(&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(cgeev)(&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(cggev)(&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(cgeev)(&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(cggev)(&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(cgeev)(&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; + delete[] rwork; + +//std::cout <<"TEST zg(g|e)ev\n"< 0) laerror("convergence problem in ggev/geev in gdiagonalize()"); + + if(biorthonormalize && vl && vr) + { + if(b || beta) laerror("@@@ biorthonormalize not implemented yet for generalized non-hermitian eigenproblem");//metric b would be needed + for(int i=0; i tmp; + cblas_cdotc_sub(n,(*vr)[i],1,(*vl)[i], 1, &tmp); + tmp = 1.f/tmp; + if(biorthonormalize==1) cblas_cscal(n,&tmp,(*vl)[i],1); + if(biorthonormalize==2) cblas_cscal(n,&tmp,(*vr)[i],1); + } + } + + if(sorttype>0) + { + laerror("sorting not implemented in complex gdiagonalize"); + } + + + if (corder) { + if (vl) {vl->transposeme(n); vl->conjugateme();} + if (vr) {vr->transposeme(n); vr->conjugateme();} + } + +} + template<> void gdiagonalize(NRMat &a, NRVec< std::complex > &w, @@ -1472,6 +2577,75 @@ void gdiagonalize(NRMat &a, NRVec< std::complex > &w, if (rvr) delete rvr; } +template<> +void gdiagonalize(NRMat &a, NRVec< std::complex > &w, + NRMat< std::complex >*vl, NRMat< std::complex > *vr, + const bool corder, int n, const int sorttype, const int biorthonormalize, + NRMat *b, NRVec *beta) +{ + if(n<=0) {n = a.nrows(); if(a.ncols()!=a.nrows() ) laerror("gdiagonalize() call for a non-square matrix");} + if(n> a.nrows() || n == a.nrows() && n != a.ncols()) laerror("gdiagonalize() of too big submatrix"); + + NRVec wr(n), wi(n); + NRMat *rvl = 0; + NRMat *rvr = 0; + if (vl) rvl = new NRMat(n, n); + if (vr) rvr = new NRMat(n, n); + gdiagonalize(a, wr, wi, rvl, rvr, 0, n, sorttype, biorthonormalize, b, beta); + + //process the results into complex matrices + int i; + for (i=0; i(wr[i], wi[i]); + if (rvl || rvr) { + i = 0; + while (i < n) { + if (wi[i] == 0) { + if(corder) + { + if (vl) for (int j=0; j((*rvl)[i][j], (*rvl)[i+1][j]); + (*vl)[j][i+1] = std::complex((*rvl)[i][j], -(*rvl)[i+1][j]); + } + else + { + (*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) + for (int j=0; j((*rvr)[i][j], (*rvr)[i+1][j]); + (*vr)[j][i+1] = std::complex((*rvr)[i][j], -(*rvr)[i+1][j]); + } + else + { + (*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; + } + } + } + if (rvl) delete rvl; + if (rvr) delete rvr; +} + //for compatibility in davidson void gdiagonalize(NRMat > &a, NRVec &wr, NRVec &wi, @@ -1490,6 +2664,22 @@ for(int i=0; i > &a, NRVec &wr, NRVec &wi, + NRMat > *vl, NRMat > *vr, const bool corder, int n, const int sorttype, const int biorthonormalize, + NRMat > *b, NRVec > *beta) +{ +if(wr.size()!=wi.size()) laerror("length mismatch in gdiagonalize"); +NRVec > w(wr.size()); +gdiagonalize(a,w,vl,vr,corder,n,sorttype,biorthonormalize,b,beta); +wr.copyonwrite(); +wi.copyonwrite(); +for(int i=0; i const SparseMat realpart > >(const SparseMat< std::complex > &a) { @@ -1532,6 +2722,24 @@ const NRMat realpart > >(const NRMat< std::c return result; } +template<> +const NRMat realpart > >(const NRMat< std::complex > &a) +{ + NRMat result(a.nrows(), a.ncols()); + +#ifdef CUDALA + if(a.location == cpu){ +#endif +// NRMat result(a.nrows(), a.ncols()); + cblas_scopy(a.nrows()*a.ncols(), (const float *)a[0], 2, result, 1); +#ifdef CUDALA + }else{ + laerror("not implemented for cuda yet"); + } +#endif + return result; +} + template<> const NRMat imagpart > >(const NRMat< std::complex > &a) { @@ -1551,6 +2759,25 @@ const NRMat imagpart > >(const NRMat< std::c return result; } +template<> +const NRMat imagpart > >(const NRMat< std::complex > &a) +{ + NRMat result(a.nrows(), a.ncols()); + +#ifdef CUDALA + if(a.location == cpu){ +#endif + +// NRMat result(a.nrows(), a.ncols()); + cblas_scopy(a.nrows()*a.ncols(), (const float *)a[0]+1, 2, result, 1); +#ifdef CUDALA + }else{ + laerror("not implemented for cuda yet"); + } +#endif + return result; +} + template<> const NRSMat realpart > >(const NRSMat< std::complex > &a) { @@ -1568,6 +2795,23 @@ const NRSMat realpart > >(const NRSMat< std return result; } +template<> +const NRSMat realpart > >(const NRSMat< std::complex > &a) +{ + NRSMat result(a.nrows()); + +#ifdef CUDALA + if(a.location == cpu){ +#endif + cblas_scopy(a.nrows()*(a.ncols()+1)/2, (const float *)(&a[0]), 2, result, 1); +#ifdef CUDALA + }else{ + laerror("not implemented for cuda yet"); + } +#endif + return result; +} + template<> const NRSMat imagpart > >(const NRSMat< std::complex > &a) { @@ -1586,6 +2830,24 @@ const NRSMat imagpart > >(const NRSMat< std return result; } +template<> +const NRSMat imagpart > >(const NRSMat< std::complex > &a) +{ + NRSMat result(a.nrows()); + +#ifdef CUDALA + if(a.location == cpu){ +#endif + + cblas_scopy(a.nrows()*(a.ncols()+1)/2, (const float *)(&a[0])+1, 2, result, 1); +#ifdef CUDALA + }else{ + laerror("not implemented for cuda yet"); + } +#endif + return result; +} + template<> const NRVec realpart > >(const NRVec< std::complex > &a) { @@ -1603,6 +2865,23 @@ const NRVec realpart > >(const NRVec< std::c return result; } +template<> +const NRVec realpart > >(const NRVec< std::complex > &a) +{ + NRVec result(a.size()); + +#ifdef CUDALA + if(a.location == cpu){ +#endif + cblas_scopy(a.size(), (const float *)&a[0], 2, result, 1); +#ifdef CUDALA + }else{ + laerror("not implemented for cuda yet"); + } +#endif + return result; +} + template<> const NRVec imagpart > >(const NRVec< std::complex > &a) { @@ -1621,6 +2900,24 @@ const NRVec imagpart > >(const NRVec< std::c return result; } +template<> +const NRVec imagpart > >(const NRVec< std::complex > &a) +{ + NRVec result(a.size()); + +#ifdef CUDALA + if(a.location == cpu){ +#endif + + cblas_scopy(a.size(), (const float *)&a[0]+1, 2, result, 1); +#ifdef CUDALA + }else{ + laerror("not implemented for cuda yet"); + } +#endif + return result; +} + template<> @@ -1645,6 +2942,28 @@ const NRMat< std::complex > realmatrix > (const NRMat +const NRMat< std::complex > realmatrix > (const NRMat &a) +{ + + NRMat > result(a.nrows(), a.ncols()); + +#ifdef CUDALA + if(a.location == cpu){ +#endif + + +// NRMat > result(a.nrows(), a.ncols()); + cblas_scopy(a.nrows()*a.ncols(), a, 1, (float *)result[0], 2); +#ifdef CUDALA + }else{ + laerror("not implemented for cuda yet"); + } +#endif + + return result; +} + template<> const NRMat< std::complex > imagmatrix > (const NRMat &a) { @@ -1665,6 +2984,26 @@ const NRMat< std::complex > imagmatrix > (const NRMat +const NRMat< std::complex > imagmatrix > (const NRMat &a) +{ + NRMat< std::complex > result(a.nrows(), a.ncols()); + +#ifdef CUDALA + if(a.location == cpu){ +#endif + +// NRMat< std::complex > result(a.nrows(), a.ncols()); + cblas_scopy(a.nrows()*a.ncols(), a, 1, (float *)result[0]+1, 2); +#ifdef CUDALA + }else{ + laerror("not implemented for cuda yet"); + } +#endif + + return result; +} + template<> const NRMat< std::complex > complexmatrix > (const NRMat &re, const NRMat &im) { @@ -1675,6 +3014,16 @@ const NRMat< std::complex > complexmatrix > (const NRMat +const NRMat< std::complex > complexmatrix > (const NRMat &re, const NRMat &im) +{ + if(re.nrows()!=im.nrows() || re.ncols() != im.ncols()) laerror("incompatible sizes of real and imaginary parts"); + NRMat< std::complex > result(re.nrows(), re.ncols()); + cblas_scopy(re.nrows()*re.ncols(), re, 1, (float *)result[0], 2); + cblas_scopy(re.nrows()*re.ncols(), im, 1, (float *)result[0]+1, 2); + return result; +} + template<> const SparseSMat< std::complex > complexmatrix >(const SparseSMat &re, const SparseSMat &im) { if(re.nrows()!=im.nrows() || re.ncols() != im.ncols()) laerror("incompatible sizes of real and imaginary parts"); @@ -1743,6 +3092,20 @@ NRMat realmatrixfunction(NRMat a, double (*f) (const double)) return r; } +NRMat realmatrixfunction(NRMat a, float (*f) (const float)) +{ + int n = a.nrows(); + NRVec w(n); + diagonalize(a, w, true, false); + + for (int i=0; i u = a; + a.diagmultl(w); + NRMat r(n, n); + r.gemm(0.0, u, 't', a, 'n', 1.0); + return r; +} + NRMat > complexmatrixfunction(NRMat a, double (*fre) (const double), double (*fim) (const double)) { @@ -1763,6 +3126,25 @@ NRMat > complexmatrixfunction(NRMat a, double (*fre return r; } +NRMat > complexmatrixfunction(NRMat a, float (*fre) (const float), float (*fim) (const float)) +{ + int n = a.nrows(); + NRVec wre(n),wim(n); + diagonalize(a, wre, true, false); + for (int i=0; i u = a; + NRMat b = a; + a.diagmultl(wre); + b.diagmultl(wim); + NRMat t(n,n),tt(n,n); + t.gemm(0.0, u, 't', a, 'n', 1.0); + tt.gemm(0.0, u, 't', b, 'n', 1.0); + NRMat > r(n, n); + for (int i=0; i(t(i,j),tt(i,j)); + return r; +} + @@ -1773,64 +3155,124 @@ std::complex myccopy (const std::complex &x) return x; } +std::complex myccopy (const std::complex &x) +{ + return x; +} + double mycopy (const double x) { return x; } +float mycopy (const float x) +{ + return x; +} + std::complex myclog (const std::complex &x) { return log(x); } +std::complex myclog (const std::complex &x) +{ + return log(x); +} + std::complex mycexp (const std::complex &x) { return std::exp(x); } +std::complex mycexp (const std::complex &x) +{ + return std::exp(x); +} + std::complex sqrtinv (const std::complex &x) { return 1./std::sqrt(x); } +std::complex sqrtinv (const std::complex &x) +{ + return 1.f/std::sqrt(x); +} + double sqrtinv (const double x) { return 1./std::sqrt(x); } +float sqrtinv (const float x) +{ + return 1./std::sqrt(x); +} + NRMat log(const NRMat &a) { return matrixfunction(a, &myclog); } +NRMat log(const NRMat &a) +{ + return matrixfunction(a, &myclog); +} + NRMat > log(const NRMat > &a) { return matrixfunction(a, &myclog); } +NRMat > log(const NRMat > &a) +{ + return matrixfunction(a, &myclog); +} + NRMat exp0(const NRMat &a) { return matrixfunction(a, &mycexp); } +NRMat exp0(const NRMat &a) +{ + return matrixfunction(a, &mycexp); +} + NRMat > exp0(const NRMat > &a) { return matrixfunction(a, &mycexp); } +NRMat > exp0(const NRMat > &a) +{ + return matrixfunction(a, &mycexp); +} + NRMat > copytest(const NRMat > &a) { return matrixfunction(a, &myccopy); } +NRMat > copytest(const NRMat > &a) +{ + return matrixfunction(a, &myccopy); +} + NRMat copytest(const NRMat &a) { return matrixfunction(a, &myccopy); } +NRMat copytest(const NRMat &a) +{ + return matrixfunction(a, &myccopy); +} + @@ -1852,6 +3294,23 @@ const NRVec diagofproduct(const NRMat &a, const NRMat &b return result; } +const NRVec diagofproduct(const NRMat &a, const NRMat &b, + bool trb, bool conjb) +{ + if (trb && (a.nrows() != b.nrows() || a.ncols() != b.ncols()) || + !trb && (a.nrows() != b.ncols() || a.ncols() != b.nrows())) + laerror("incompatible Mats in diagofproduct()"); + NRVec result(a.nrows()); + if (trb) + for(int i=0; i > diagofproduct(const NRMat< std::complex > &a, const NRMat< std::complex > &b, bool trb, bool conjb) @@ -1880,6 +3339,33 @@ const NRVec< std::complex > diagofproduct(const NRMat< std::complex > diagofproduct(const NRMat< std::complex > &a, + const NRMat< std::complex > &b, bool trb, bool conjb) +{ + if (trb && (a.nrows() != b.nrows() || a.ncols() != b.ncols()) || + !trb && (a.nrows() != b.ncols() || a.ncols() != b.nrows())) + laerror("incompatible Mats in diagofproduct()"); + NRVec< std::complex > result(a.nrows()); + if (trb) { + if (conjb) { + for(int i=0; i &a, const NRMat &b, bool trb) { @@ -1895,6 +3381,20 @@ double trace2(const NRMat &a, const NRMat &b, bool trb) return sum; } +float trace2(const NRMat &a, const NRMat &b, bool trb) +{ + if (trb && (a.nrows() != b.nrows() || a.ncols() != b.ncols()) || + !trb && (a.nrows() != b.ncols() || a.ncols() != b.nrows())) + laerror("incompatible Mats in trace2()"); + if (trb) return cblas_sdot(a.nrows()*a.ncols(), a, 1, b, 1); + + float sum = 0.0; + for (int i=0; i trace2(const NRMat > &a, const NRMat > &b, bool adjb) { @@ -1914,6 +3414,24 @@ std::complex trace2(const NRMat > &a, const NRMat trace2(const NRMat > &a, const NRMat > &b, bool adjb) +{ + if (adjb && (a.nrows() != b.nrows() || a.ncols() != b.ncols()) || + !adjb && (a.nrows() != b.ncols() || a.ncols() != b.nrows())) + laerror("incompatible Mats in trace2()"); + + std::complex dot; + if (adjb) { cblas_cdotc_sub(a.nrows()*a.ncols(), b, 1, a, 1, &dot); return dot; } + + std::complex sum = std::complex(0.,0.); + for (int i=0; i &a, const NRSMat &b, const bool diagscaled) { @@ -1927,6 +3445,19 @@ double trace2(const NRSMat &a, const NRSMat &b, return r; } +float trace2(const NRSMat &a, const NRSMat &b, + const bool diagscaled) +{ + if (a.nrows() != b.nrows()) laerror("incompatible SMats in trace2()"); + + //float r = 0; for (int i=0; i &a, const NRMat &b, const bool diagscaled) { if (a.nrows() != b.nrows()||b.nrows()!=b.ncols()) laerror("incompatible SMats in trace2()"); @@ -1936,6 +3467,18 @@ double r=0; for (j=0; j<=i;j++) r += a[k++] * (b[i][j] + (i!=j||diagscaled ? b[j][i] : 0 )); +return r; +} + +float trace2(const NRSMat &a, const NRMat &b, const bool diagscaled) +{ + if (a.nrows() != b.nrows()||b.nrows()!=b.ncols()) laerror("incompatible SMats in trace2()"); +float r=0; + int i, j, k=0; + for (i=0; i &a, const NRSMat &b, const bool diagsc return trace2(b,a,diagscaled); } +float 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, @@ -1959,6 +3507,19 @@ std::complex trace2(const NRSMat > &a, const NRSMat return r; } +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_cdotc_sub(a.nrows()*(a.nrows()+1)/2, b, 1, a, 1, &dot); //sum over the lower triangle of a_ij conj(b_ij) + float 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()"); @@ -1969,15 +3530,32 @@ std::complex r=0; return r; } +std::complex 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); } +std::complex trace2(const NRMat > &a, const NRSMat > &b, const bool diagscaled) +{ +return trace2(b,a,diagscaled); +} + //Cholesky interface extern "C" void FORNAME(dpotrf)(const char *UPLO, const FINT *N, double *A, const FINT *LDA, FINT *INFO); +extern "C" void FORNAME(spotrf)(const char *UPLO, const FINT *N, float *A, const FINT *LDA, FINT *INFO); extern "C" void FORNAME(zpotrf)(const char *UPLO, const FINT *N, std::complex *A, const FINT *LDA, FINT *INFO); +extern "C" void FORNAME(cpotrf)(const char *UPLO, const FINT *N, std::complex *A, const FINT *LDA, FINT *INFO); void cholesky(NRMat &a, bool upper) { @@ -1997,6 +3575,24 @@ else } +void cholesky(NRMat &a, bool upper) +{ +if(a.nrows()!=a.ncols()) laerror("matrix must be square in Cholesky"); +FINT lda=a.ncols(); +FINT n=a.nrows(); +char uplo= LAPACK_FORTRANCASE(upper?'u':'l'); +FINT info; +a.copyonwrite(); +FORNAME(spotrf)(&uplo, &n, a, &lda, &info); +if(info) {std::cerr << "Lapack error "< > &a, bool upper) { @@ -2016,10 +3612,30 @@ else for(int i=0; i > &a, bool upper) +{ +if(a.nrows()!=a.ncols()) laerror("matrix must be square in Cholesky"); +FINT lda=a.ncols(); +FINT n=a.nrows(); +char uplo= LAPACK_FORTRANCASE(upper?'u':'l'); +FINT info; +a.copyonwrite(); +a.transposeme();//switch to Fortran order +FORNAME(cpotrf)(&uplo, &n, a, &lda, &info); +if(info) {std::cerr << "Lapack error "< *A, const FINT *LDA, double *WORK); +extern "C" float FORNAME(clange)( const char *NORM, const FINT *M, const FINT *N, const std::complex *A, const FINT *LDA, float *WORK); extern "C" double FORNAME(dlange)( const char *NORM, const FINT *M, const FINT *N, const double *A, const FINT *LDA, double *WORK); +extern "C" float FORNAME(slange)( const char *NORM, const FINT *M, const FINT *N, const float *A, const FINT *LDA, float *WORK); //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) @@ -2046,6 +3662,17 @@ LA_traits > >::normtype MatrixNorm(const NRMat +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(); + float work[M]; + const float ret = FORNAME(clange)(&TypNorm, &M, &N, A[0], &M, &work[0]); + return ret; +} + template<> LA_traits >::normtype MatrixNorm(const NRMat &A, const char norm) { @@ -2057,13 +3684,28 @@ LA_traits >::normtype MatrixNorm(const NRMat &A, const cha return ret; } +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(); + float work[M]; + const float ret = FORNAME(slange)(&TypNorm, &M, &N, A[0], &M, &work[0]); + return ret; +} + //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(cgetrf)(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(sgetrf)(const FINT *M, const FINT *N, float *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(cgecon)( const char *norm, const FINT *n, const std::complex *A, const FINT *LDA, const float *anorm, float *rcond, std::complex *work, float *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); +extern "C" void FORNAME(sgecon)( const char *norm, const FINT *n, const float *A, const FINT *LDA, const float *anorm, float *rcond, float *work, FINT *iwork, FINT *info); template<> LA_traits > >::normtype CondNumber(const NRMat > &A, const char norm) @@ -2093,6 +3735,34 @@ LA_traits > >::normtype CondNumber(const NRMat::infinity() : 1./rcond; } +template<> +LA_traits > >::normtype CondNumber(const NRMat > &A, const char norm) +{ + 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(); + if(N != A.ncols()) laerror("nonsquare matrix in CondNumber"); + float Norma(0.0), rcond(0.0); + FINT info; + + NRMat > LU(A); + LU.copyonwrite(); + FINT *ipiv = new FINT[N]; + FORNAME(cgetrf)(&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]; + float *rwork = new float[2*N]; + Norma = MatrixNorm(A, norm); + FORNAME(cgecon)(&TypNorm, &N, LU[0], &N, &Norma, &rcond, work, rwork, &info); + delete[] work; + delete[] rwork; + if(info<0) laerror("illegal argument in zgecon in CondNumber"); + return rcond==0. ? std::numeric_limits::infinity() : 1./rcond; +} + template<> LA_traits >::normtype CondNumber(const NRMat &A, const char norm) { @@ -2121,6 +3791,34 @@ LA_traits >::normtype CondNumber(const NRMat &A, const cha return rcond==0. ? std::numeric_limits::infinity() : 1./rcond; } +template<> +LA_traits >::normtype CondNumber(const NRMat &A, const char norm) +{ + 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(); + if(N != A.ncols()) laerror("nonsquare matrix in CondNumber"); + float Norma(0.0), rcond(0.0); + FINT info; + + NRMat LU(A); + LU.copyonwrite(); + FINT *ipiv = new FINT[N]; + FORNAME(sgetrf)(&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 + + float *work = new float[4*N]; + FINT *iwork = new FINT[N]; + Norma = MatrixNorm(A, norm); + FORNAME(sgecon)(&TypNorm, &N, LU[0], &N, &Norma, &rcond, work, iwork, &info); + delete[] work; + delete[] iwork; + if(info<0) laerror("illegal argument in dgecon in CondNumber"); + return rcond==0. ? std::numeric_limits::infinity() : 1./rcond; +} + #ifdef obsolete void gendiagonalize(NRMat &a, NRVec &w, NRMat b, int n) @@ -2222,5 +3920,29 @@ if(nchange&1)//still adjust to get determinant=1 } } +void adjustphases(NRMat &v) +{ +int n=v.nrows(); +float det=determinant(v); +int nchange=0; +for(int i=0; i &a, NRMat *u, NRVec: declare_la(double) declare_la(std::complex) +declare_la(float) +declare_la(std::complex) //svd_solve without contamination from nullspace, b is n x nrhs, result returns in b, A is destroyed template -NRVec svd_solve(NRMat &a, const NRVec &b, double thres=1e-12, int *nullity = NULL) +NRVec svd_solve(NRMat &a, const NRVec &b, typename LA_traits::normtype thres=LA_precision::normtype >::svd_threshold(), int *nullity = NULL) { if(b.size()!=a.nrows()) laerror("size mismatch in svd_solve"); if(nullity) *nullity=0; -NRVec w(a.ncols()); +NRVec::normtype > w(a.ncols()); NRMat u(a.nrows(),a.ncols()); NRMat vt(a.ncols(),a.ncols()); singular_decomposition(a,&u,w,&vt,false); @@ -191,11 +193,11 @@ return res; template -NRMat svd_solve(NRMat &a, const NRMat &b, double thres=1e-12, int *nullity = NULL) +NRMat svd_solve(NRMat &a, const NRMat &b, typename LA_traits::normtype thres=LA_precision::normtype >::svd_threshold(), int *nullity = NULL) { if(b.nrows()!=a.nrows()) laerror("size mismatch in svd_solve"); if(nullity) *nullity=0; -NRVec w(a.ncols()); +NRVec::normtype > w(a.ncols()); NRMat u(a.nrows(),a.ncols()); NRMat vt(a.ncols(),a.ncols()); singular_decomposition(a,&u,w,&vt,false); @@ -221,10 +223,13 @@ return res; 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: for a complex pair (i,i+1) the eigenvector of w_i is Re + i Im = vector i + i*vector i+1 +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); template -extern void gdiagonalize(NRMat &a, NRVec< std::complex > &w, - NRMat< std::complex >*vl=NULL, NRMat< std::complex > *vr=NULL, +extern void gdiagonalize(NRMat &a, NRVec::complextype > &w, + NRMat::complextype > *vl=NULL, NRMat::complextype > *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0, NRMat *b=NULL, NRVec *beta=NULL); //eigenvectors are stored in complex matrices for T both double and complex @@ -232,9 +237,13 @@ extern void gdiagonalize(NRMat &a, NRVec< std::complex > &w, 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); +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); //diagonalization of symmetric real tridiagonal matrix extern void diagonalize3(NRVec &d, NRVec &d1, NRMat *v, const bool corder=1, int n0=0); +extern void diagonalize3(NRVec &d, NRVec &d1, NRMat *v, const bool corder=1, int n0=0); //complex,real,imaginary parts of various entities template @@ -251,6 +260,8 @@ extern const typename LA_traits::complextype complexmatrix (const T&, const T //Cholesky decomposition extern void cholesky(NRMat &a, bool upper=1); extern void cholesky(NRMat > &a, bool upper=1); +extern void cholesky(NRMat &a, bool upper=1); +extern void cholesky(NRMat > &a, bool upper=1); //inverse by means of linear solve, pass by value to preserve argument intact template @@ -268,13 +279,13 @@ const NRMat calcinverse(NRMat a, T *det=NULL) //inverse by means of SVD , pass by value to preserve argument intact template -const NRMat calcsvdinverse(NRMat a, double thr=0, bool do_absadjoint=false) +const NRMat calcsvdinverse(NRMat a, typename LA_traits::normtype thr=0, bool do_absadjoint=false) { if(a.nrows()!=a.ncols()) laerror("svdinverse() for non-square matrix"); int n=a.nrows(); a.copyonwrite(); NRMat u(n,n),v(n,n); - NRVec w(n),s(n); + NRVec::normtype > w(n),s(n); singular_decomposition(a,&u,w,&v,true); if(do_absadjoint) //absadjoint matrix - missing det(U)*det(V) to be full adjoint - a sign or complex phase { @@ -348,7 +359,7 @@ return determinant_destroy(a); // the info parameter of gesvx is returned (see man dgesvx) //------------------------------------------------------------------------------ template -int linear_solve_x(NRMat &A, T *B, const int rhsCount, const int eqCount, const bool eq, const bool saveA, double *rcond); +int linear_solve_x(NRMat &A, T *B, const int rhsCount, const int eqCount, const bool eq, const bool saveA, typename LA_traits::normtype *rcond); //------------------------------------------------------------------------------ @@ -364,7 +375,7 @@ int linear_solve_x(NRMat &A, T *B, const int rhsCount, const int eqCount, con // the info parameter of dgesvx is returned (see man dgesvx) //------------------------------------------------------------------------------ template -int multiply_by_inverse(NRMat &A, NRMat &B, bool useEq=false, double *rcond=NULL); +int multiply_by_inverse(NRMat &A, NRMat &B, bool useEq=false, typename LA_traits::normtype *rcond=NULL); //general submatrix, INDEX will typically be NRVec or even int* @@ -414,6 +425,7 @@ return r; //auxiliary routine to adjust eigenvectors to guarantee real logarithm extern void adjustphases(NRMat &v); +extern void adjustphases(NRMat &v); //declaration of template interface to cblas routines with full options available @@ -443,6 +455,24 @@ inline double xdot(int n, const double *x, int incx, const double *y, in return cblas_ddot(n,x,incx,y,incy); } +template<> +inline void xcopy (int n, const float *x, int incx, float *y, int incy) +{ +cblas_scopy(n, x, incx, y, incy); +} + +template<> +inline void xaxpy(int n, const float &a, const float *x, int incx, float *y, int incy) +{ +cblas_saxpy(n, a, x, incx, y, incy); +} + +template<> +inline float xdot(int n, const float *x, int incx, const float *y, int incy) +{ +return cblas_sdot(n,x,incx,y,incy); +} + //debugging aid: reconstruct an explicit matrix from the implicit version @@ -468,12 +498,14 @@ return r; extern NRMat realmatrixfunction(NRMat a, double (*f) (double)); //a has to by in fact symmetric extern NRMat > complexmatrixfunction(NRMat a, double (*fre) (double), double (*fim) (double)); //a has to by in fact symmetric +extern NRMat realmatrixfunction(NRMat a, float (*f) (float)); //a has to by in fact symmetric +extern NRMat > complexmatrixfunction(NRMat a, float (*fre) (float), float (*fim) (float)); //a has to by in fact symmetric template -NRMat matrixfunction(NRSMat a, double (*f) (double)) //of symmetric/hermitian matrix +NRMat matrixfunction(NRSMat a, typename LA_traits::normtype (*f) (typename LA_traits::normtype)) //of symmetric/hermitian matrix { int n = a.nrows(); - NRVec w(n); + NRVec::normtype > w(n); NRMat v(n, n); diagonalize(a, w, &v, 0); @@ -488,18 +520,18 @@ NRMat matrixfunction(NRSMat a, double (*f) (double)) //of symmetric/hermit template -extern NRMat matrixfunction(NRMat a, std::complex (*f)(const std::complex &)) //of a general real/complex matrix +extern NRMat matrixfunction(NRMat a, typename LA_traits::complextype (*f)(const typename LA_traits::complextype &)) //of a general real/complex matrix { int n = a.nrows(); - NRVec > w(n); - NRMat > u(n,n),v(n,n); + NRVec::complextype > w(n); + NRMat::complextype > u(n,n),v(n,n); #ifdef debugmf -NRMat > a0=a; +NRMat::complextype > a0=a; #endif gdiagonalize(a, w, &u, &v, false,n,0,false,NULL,NULL);//a gets destroyed, eigenvectors are rows - NRVec< std::complex > z = diagofproduct(u, v, 1, 1); + NRVec< typename LA_traits::complextype > z = diagofproduct(u, v, 1, 1); #ifdef debugmf std::cout <<"TEST matrixfunction\n"< > wz(n); +NRVec< typename LA_traits::complextype > wz(n); for (int i=0; i > r(n, n); + NRMat< typename LA_traits::complextype > r(n, n); r.gemm(0.0, v, 'c', u, 'n', 1.0); return (NRMat) r; //convert back to real if applicable by the explicit decomplexifying constructor; it is NOT checked to which accuracy the imaginary part is actually zero } @@ -528,6 +560,8 @@ std::cout <<"TEST matrixfunction4 "<< a0<< v.transpose(true)*diagonalmatrix(wz)* extern std::complex sqrtinv(const std::complex &); extern double sqrtinv(const double); +extern std::complex sqrtinv(const std::complex &); +extern float sqrtinv(const float); //functions on matrices inline NRMat sqrt(const NRSMat &a) { return matrixfunction(a,&std::sqrt); } @@ -541,6 +575,17 @@ extern NRMat > exp0(const NRMat > &a); extern NRMat > copytest(const NRMat > &a); extern NRMat copytest(const NRMat &a); extern NRMat exp0(const NRMat &a); +inline NRMat sqrt(const NRSMat &a) { return matrixfunction(a,&std::sqrt); } +inline NRMat sqrtinv(const NRSMat &a) { return matrixfunction(a,&sqrtinv); } +inline NRMat realsqrt(const NRMat &a) { return realmatrixfunction(a,&std::sqrt); } +inline NRMat realsqrtinv(const NRMat &a) { return realmatrixfunction(a,&sqrtinv); } +inline NRMat log(const NRSMat &a) { return matrixfunction(a,&std::log); } +extern NRMat log(const NRMat &a); +extern NRMat > log(const NRMat > &a); +extern NRMat > exp0(const NRMat > &a); +extern NRMat > copytest(const NRMat > &a); +extern NRMat copytest(const NRMat &a); +extern NRMat exp0(const NRMat &a); }//namespace diff --git a/smat.cc b/smat.cc index af17b97..0755ee7 100644 --- a/smat.cc +++ b/smat.cc @@ -1324,9 +1324,15 @@ void NRSMat::storesubmatrix(const NRVec &selection, const NRSMat &rhs) template<> NRSMat NRSMat::inverse() {return NRSMat(NRMat(*this).inverse());} +template<> +NRSMat NRSMat::inverse() {return NRSMat(NRMat(*this).inverse());} + template<> NRSMat > NRSMat >::inverse() {return NRSMat >(NRMat >(*this).inverse());} +template<> +NRSMat > NRSMat >::inverse() {return NRSMat >(NRMat >(*this).inverse());} + /***************************************************************************//**