diff --git a/nonclass.cc b/nonclass.cc index d9f10df..a51837f 100644 --- a/nonclass.cc +++ b/nonclass.cc @@ -178,7 +178,7 @@ void linear_solve(NRMat &A, NRMat *B, 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?(double *)B:NULL,B?B->nrows() : 0, B?B->ncols():A.nrows(), det,n); +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) @@ -186,7 +186,7 @@ void linear_solve(NRMat &A, NRVec &B, 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,(double *)B,1,A.nrows(),det,n); +linear_solve_do(A,&B[0],1,A.nrows(),det,n); } @@ -425,7 +425,7 @@ if(gdbeta) {tmp=gdbeta[i]; gdbeta[i]=gdbeta[j]; gdbeta[j]=tmp;} void gdiagonalize(NRMat &a, NRVec &wr, NRVec &wi, NRMat *vl, NRMat *vr, const bool corder, int n, - const int sorttype, const bool biorthonormalize, + const int sorttype, const int biorthonormalize, NRMat *b, NRVec *beta) { if(n<=0) n = a.nrows(); @@ -470,7 +470,7 @@ void gdiagonalize(NRMat &a, NRVec &wr, NRVec &wi, 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); delete[] work; -//@@@ + //cout <<"TEST dgeev\n"< &a, NRVec &wr, NRVec &wi, //calculate scaling paramter double tmp; tmp=cblas_ddot(n,(*vl)[i],1,(*vr)[i], 1); - tmp=1./sqrt(abs(tmp)); - cblas_dscal(n,tmp,(*vl)[i],1); - cblas_dscal(n,tmp,(*vr)[i],1); + if(biorthonormalize==1) cblas_dscal(n,1./tmp,(*vl)[i],1); + if(biorthonormalize==2) cblas_dscal(n,1./tmp,(*vr)[i],1); i++; } else //complex pair @@ -567,7 +566,7 @@ void gdiagonalize(NRMat &a, NRVec &wr, NRVec &wi, void gdiagonalize(NRMat &a, NRVec< complex > &w, NRMat< complex >*vl, NRMat< complex > *vr, - const bool corder, int n, const int sorttype, const bool biorthonormalize, + const bool corder, int n, const int sorttype, const int biorthonormalize, NRMat *b, NRVec *beta) { if(!corder) laerror("gdiagonalize() corder 0 not implemented"); diff --git a/nonclass.h b/nonclass.h index 38fc336..1cde342 100644 --- a/nonclass.h +++ b/nonclass.h @@ -105,11 +105,11 @@ declare_la(complex) // Separate declarations //general nonsymmetric matrix and generalized diagonalization extern void gdiagonalize(NRMat &a, NRVec &wr, NRVec &wi, - NRMat *vl, NRMat *vr, const bool corder=1, int n=0, const int sorttype=0, const bool biorthonormalize=0, + NRMat *vl, NRMat *vr, 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< complex > &w, NRMat< complex >*vl, NRMat< complex > *vr, - const bool corder=1, int n=0, const int sorttype=0, const bool biorthonormalize=0, + const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0, NRMat *b=NULL, NRVec *beta=NULL); extern NRMat matrixfunction(NRSMat a, double (*f) (double)); extern NRMat realmatrixfunction(NRMat a, double (*f) (double)); //a has to by in fact symmetric @@ -143,7 +143,9 @@ const NRMat inverse(NRMat a, T *det=0) #endif NRMat result(a.nrows(),a.nrows()); result = (T)1.; + a.copyonwrite(); linear_solve(a, &result, det); + result.transposeme(); //tested with noncblas return result; }