bugfix of complex calcinverse and calcsvdinverse in nonclass

nonclass.cc

@@ -266,60 +266,77 @@ linear_solve_do(a,&B[0],1,a.nrows(),det,n);
 
 
 // Roman, complex version of linear_solve()
-extern "C" void FORNAME(zgesv)(const int *N, const int *NRHS, double *A, const int *LDA,
+extern "C" void FORNAME(zgesv)(const FINT *N, const FINT *NRHS, std::complex<double> *A, const FINT *LDA,
-           int *IPIV, double *B, const int *LDB, int *INFO);
+           FINT *IPIV, std::complex<double> *B, const FINT *LDB, FINT *INFO);
 
-void linear_solve(NRMat< std::complex<double> > &A, NRMat< std::complex<double> > *B, std::complex<double> *det, int n)
+// 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
+static void linear_solve_do(NRMat<std::complex<double> > &A, std::complex<double> *B, const int nrhs, const int ldb, std::complex<double> *det, int n)
 {
-        int r, *ipiv;
+	FINT r, *ipiv;
 	int iswap=0;
 
-        if (A.nrows() != A.ncols()) laerror("linear_solve() call for non-square matrix");
+	if(nrhs==0) std::cout<<"Warning: some zgesv implementations might skip LU decomposition when nrhs==0\n";
-        if (B && A.nrows() != B->ncols()) laerror("incompatible matrices in linear_solve()");
+	if (n==A.nrows() && A.nrows() != A.ncols()) laerror("linear_solve() call for non-square matrix");
-        A.copyonwrite();
+	A.copyonwrite();
-        if (B) B->copyonwrite();
+	ipiv = new FINT[n];
-        ipiv = new int[A.nrows()];
+	const int lda=A.ncols();
-  n = A.nrows();
+	std::complex<double> *a= &A(0,0);
-  int nrhs = B ? B->nrows() : 0;
+	for(int i=1; i<n; ++i) for(int j=0; j<i; ++j)  {std::complex<double> t=a[j*lda+i]; a[j*lda+i]=a[i*lda+j]; a[i*lda+j]=t;}
-  int lda = A.ncols();
+	const FINT ntmp=n;
-  int ldb = B ? B->ncols() : A.nrows();
+	const FINT nrhstmp=nrhs;
-        FORNAME(zgesv)(&n, &nrhs, (double *)A[0], &lda, ipiv,
+	const FINT ldatmp=lda;
-                 B ? (double *)(*B)[0] : (double *)0, &ldb, &r);
+	const FINT ldbtmp=ldb;
-        if (r < 0) {
+	FORNAME(zgesv)(&ntmp, &nrhstmp, a, &ldatmp, ipiv, B, &ldbtmp, &r);
-                delete[] ipiv;
+	for(int i=1; i<n; ++i) for(int j=0; j<i; ++j)  {std::complex<double> t=a[j*lda+i]; a[j*lda+i]=a[i*lda+j]; a[i*lda+j]=t;}
-                laerror("illegal argument in lapack_gesv");
+	if (r < 0) {
-        }
+		delete[] ipiv;
-        if (det && r==0) {
+		laerror("illegal argument in lapack_gesv");
-                *det = A[0][0];
+	}
-                for (int i=1; i<A.nrows(); ++i) *det *= A[i][i];
+	if (det && r==0) {
+		*det = 1.;
+		//take into account some numerical instabilities in zgesv for singular matrices
+		for (int i=0; i<n; ++i)  {std::complex<double> 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; i<n; ++i)
+	        for (int i=0; i<n; ++i)
-                        {
+			{
-                        if(ipiv[i]==0) shift=0;
+			if(ipiv[i]==0) shift=0;
-                        if(ipiv[i]<0 || ipiv[i]>n) 
+			if(ipiv[i]<0 || ipiv[i]>n) 
-                                {
+				{
-                                std::cout <<"IPIV["<<i<<"] = "<<ipiv[i]<<std::endl;
+				std::cout <<"IPIV["<<i<<"] = "<<ipiv[i]<<std::endl;
-                                laerror("problem with ipiv in zgesv");
+				laerror("problem with ipiv in zgesv");
-                                }
+				}
-                        }
+			}
 #ifdef IPIV_DEBUG
        std::cout <<"shift = "<<shift<<std::endl;
 #endif
-                //
+		//change sign of det by parity of ipiv permutation
-                //change sign of det by parity of ipiv permutation
+		if(*det!=0.) for (int i=0; i<n; ++i) if(i+shift != ipiv[i]) {*det = -(*det); ++iswap;}
-		if(det) for (int i=0; i<A.nrows(); ++i) if(i+shift != ipiv[i]) {*det = -(*det); ++iswap;}
+	}
-        }
 	if(det && r>0) *det = 0;
 #ifdef IPIV_DEBUG
        std::cout <<"iswap = "<<iswap<<std::endl;
-        
-        std::cout <<"ipiv = ";
-        for (int i=0; i<n; ++i) std::cout <<ipiv[i]<<" ";
-        std::cout <<std::endl;
 #endif
 
-        delete [] ipiv;
+	delete [] ipiv;
-        if (r>0 && B) laerror("singular matrix in zgesv");
+	if (r>0 && B) laerror("singular matrix in zgesv");
+}
+
+void linear_solve(NRMat< std::complex<double> > &A, NRMat< std::complex<double> > *B, std::complex<double> *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<double> > &A, NRVec< std::complex<double> > &B, std::complex<double> *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);
 }
 
 

nonclass.h

@@ -154,8 +154,8 @@ extern T trace2(const NRSMat<T> &a, const NRMat<T> &b, const bool diagscaled=0);
 extern T trace2(const NRMat<T> &a, const NRSMat<T> &b, const bool diagscaled=0);\
 extern void linear_solve(NRMat<T> &a, NRMat<T> *b, T *det=0,int n=0); /*solve Ax^T=b^T (b is nrhs x n) */ \
 extern void linear_solve(NRSMat<T> &a, NRMat<T> *b, T *det=0, int n=0); /*solve Ax^T=b^T (b is nrhs x n) */\
-extern void linear_solve(NRMat<T> &a, NRVec<T> &b, double *det=0, int n=0); \
+extern void linear_solve(NRMat<T> &a, NRVec<T> &b, T *det=0, int n=0); \
-extern void linear_solve(NRSMat<T> &a, NRVec<T> &b, double *det=0, int n=0); \
+extern void linear_solve(NRSMat<T> &a, NRVec<T> &b, T *det=0, int n=0); \
 extern void diagonalize(NRMat<T> &a, NRVec<LA_traits<T>::normtype> &w, const bool eivec=1, const bool corder=1, int n=0, NRMat<T> *b=NULL, const int itype=1); \
 extern void diagonalize(NRSMat<T> &a, NRVec<LA_traits<T>::normtype> &w, NRMat<T> *v, const bool corder=1, int n=0, NRSMat<T> *b=NULL, const int itype=1);\
 extern void singular_decomposition(NRMat<T> &a, NRMat<T> *u, NRVec<LA_traits<T>::normtype> &s, NRMat<T> *v, const bool vnotdagger=0, int m=0, int n=0); \
@@ -285,6 +285,7 @@ const NRMat<T> calcsvdinverse(NRMat<T> a, double thr=0, bool do_absadjoint=false
 		for(int i=0; i<n; ++i) s[i] = (w[i]<thr)? 0. : 1./w[i];
 	v.diagmultr(s);
 	u.transposeme();
+	u.conjugateme(); //A = U W V^H, hence A^-1 = V W^-1 U^H
         return v*u; //could use gemm instead of separate transpose too
 }