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

*** empty log message ***

Browse Source
This commit is contained in:
jiri 2008-03-02 16:40:22 +00:00
parent 58d7d5bd50
commit 5f84d14ce2
2 changed files with 11 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
nonclass.cc
View File

@ -178,7 +178,7 @@ void linear_solve(NRMat<double> &A, NRMat<double> *B, double *det, int n)
if(n<=0) n=A.nrows(); //default - whole matrix 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 (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(); 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<double> &A, NRVec<double> &B, double *det, int n) void linear_solve(NRMat<double> &A, NRVec<double> &B, double *det, int n)
@ -186,7 +186,7 @@ void linear_solve(NRMat<double> &A, NRVec<double> &B, double *det, int n)
if(n<=0) n=A.nrows(); //default - whole matrix 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()"); if(n==A.nrows() && A.nrows() != B.size() || n>B.size()||n>A.nrows() ) laerror("incompatible matrices in linear_solve()");
B.copyonwrite(); 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<double> &a, NRVec<double> &wr, NRVec<double> &wi, void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
NRMat<double> *vl, NRMat<double> *vr, const bool corder, int n, NRMat<double> *vl, NRMat<double> *vr, const bool corder, int n,
const int sorttype, const bool biorthonormalize, const int sorttype, const int biorthonormalize,
NRMat<double> *b, NRVec<double> *beta) NRMat<double> *b, NRVec<double> *beta)
{ {
if(n<=0) n = a.nrows(); if(n<=0) n = a.nrows();
@ -470,7 +470,7 @@ void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
else FORNAME(dgeev)(&jobvr, &jobvl, &n, a, &lda, wr, wi, vr?vr[0]:(double *)0, 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); &ldvr, vl?vl[0]:(double *)0, &ldvl, work, &lwork, &r);
delete[] work; delete[] work;
//@@@
//cout <<"TEST dgeev\n"<<wr<<wi<<*vr<<*vl<<endl; //cout <<"TEST dgeev\n"<<wr<<wi<<*vr<<*vl<<endl;
if (r < 0) laerror("illegal argument in ggev/geev in gdiagonalize()"); if (r < 0) laerror("illegal argument in ggev/geev in gdiagonalize()");
@ -487,9 +487,8 @@ void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
//calculate scaling paramter //calculate scaling paramter
double tmp; double tmp;
tmp=cblas_ddot(n,(*vl)[i],1,(*vr)[i], 1); tmp=cblas_ddot(n,(*vl)[i],1,(*vr)[i], 1);
tmp=1./sqrt(abs(tmp)); if(biorthonormalize==1) cblas_dscal(n,1./tmp,(*vl)[i],1);
cblas_dscal(n,tmp,(*vl)[i],1); if(biorthonormalize==2) cblas_dscal(n,1./tmp,(*vr)[i],1);
cblas_dscal(n,tmp,(*vr)[i],1);
i++; i++;
} }
else //complex pair else //complex pair
@ -567,7 +566,7 @@ void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
void gdiagonalize(NRMat<double> &a, NRVec< complex<double> > &w, void gdiagonalize(NRMat<double> &a, NRVec< complex<double> > &w,
NRMat< complex<double> >*vl, NRMat< complex<double> > *vr, NRMat< complex<double> >*vl, NRMat< complex<double> > *vr,
const bool corder, int n, const int sorttype, const bool biorthonormalize, const bool corder, int n, const int sorttype, const int biorthonormalize,
NRMat<double> *b, NRVec<double> *beta) NRMat<double> *b, NRVec<double> *beta)
{ {
if(!corder) laerror("gdiagonalize() corder 0 not implemented"); if(!corder) laerror("gdiagonalize() corder 0 not implemented");

6
nonclass.h
View File

@ -105,11 +105,11 @@ declare_la(complex<double>)
// Separate declarations // Separate declarations
//general nonsymmetric matrix and generalized diagonalization //general nonsymmetric matrix and generalized diagonalization
extern void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi, extern void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
NRMat<double> *vl, NRMat<double> *vr, const bool corder=1, int n=0, const int sorttype=0, const bool biorthonormalize=0, NRMat<double> *vl, NRMat<double> *vr, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
NRMat<double> *b=NULL, NRVec<double> *beta=NULL); NRMat<double> *b=NULL, NRVec<double> *beta=NULL);
extern void gdiagonalize(NRMat<double> &a, NRVec< complex<double> > &w, extern void gdiagonalize(NRMat<double> &a, NRVec< complex<double> > &w,
NRMat< complex<double> >*vl, NRMat< complex<double> > *vr, NRMat< complex<double> >*vl, NRMat< complex<double> > *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<double> *b=NULL, NRVec<double> *beta=NULL); NRMat<double> *b=NULL, NRVec<double> *beta=NULL);
extern NRMat<double> matrixfunction(NRSMat<double> a, double (*f) (double)); extern NRMat<double> matrixfunction(NRSMat<double> a, double (*f) (double));
extern NRMat<double> realmatrixfunction(NRMat<double> a, double (*f) (double)); //a has to by in fact symmetric extern NRMat<double> realmatrixfunction(NRMat<double> a, double (*f) (double)); //a has to by in fact symmetric
@ -143,7 +143,9 @@ const NRMat<T> inverse(NRMat<T> a, T *det=0)
#endif #endif
NRMat<T> result(a.nrows(),a.nrows()); NRMat<T> result(a.nrows(),a.nrows());
result = (T)1.; result = (T)1.;
a.copyonwrite();
linear_solve(a, &result, det); linear_solve(a, &result, det);
result.transposeme(); //tested with noncblas
return result; return result;
} }