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 2005-09-10 16:57:41 +00:00
parent d992dbcae9
commit 301163d965
1 changed files with 12 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
nonclass.cc
View File

@ -29,6 +29,8 @@ INSTANTIZE(int)
INSTANTIZE(short) INSTANTIZE(short)
INSTANTIZE(char) INSTANTIZE(char)
#define EPSDET 1e-300
template <typename T> template <typename T>
void lawritemat(FILE *file,const T *a,int r,int c,const char *form0, void lawritemat(FILE *file,const T *a,int r,int c,const char *form0,
int nodim,int modulo, int issym) int nodim,int modulo, int issym)
@ -127,12 +129,14 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
delete[] ipiv; delete[] ipiv;
laerror("illegal argument in lapack_gesv"); laerror("illegal argument in lapack_gesv");
} }
if (det && r>=0) { if (det && r==0) {
*det = A[0][0]; *det = 1.;
for (int i=1; i<n; ++i) *det *= A[i][i]; //take into account some numerical instabilities in dgesv for singular matrices
for (int i=0; i<n; ++i) {double t=A[i][i]; if(!finite(t) || abs(t) < EPSDET ) {*det=0.; break;} else *det *=t;}
//change sign of det by parity of ipiv permutation //change sign of det by parity of ipiv permutation
for (int i=0; i<n; ++i) *det = -(*det); if(*det) for (int i=0; i<n; ++i) *det = -(*det);
} }
if(det && r>0) *det = 0;
delete [] ipiv; delete [] ipiv;
if (r>0 && B) laerror("singular matrix in lapack_gesv"); if (r>0 && B) laerror("singular matrix in lapack_gesv");
} }
@ -171,11 +175,12 @@ static void linear_solve_do(NRSMat<double> &a, double *b, const int nrhs, const
delete[] ipiv; delete[] ipiv;
laerror("illegal argument in spsv() call of linear_solve()"); laerror("illegal argument in spsv() call of linear_solve()");
} }
if (det && r >= 0) { if (det && r == 0) {
*det = a(0,0); *det = 1.;
for (int i=1; i<n; i++) *det *= a(i,i); for (int i=1; i<n; i++) {double t=a(i,i); if(!finite(t) || abs(t) < EPSDET ) {*det=0.; break;} else *det *= t;}
//do not use ipiv, since the permutation matrix occurs twice in the decomposition and signs thus cancel (man dspsv) //do not use ipiv, since the permutation matrix occurs twice in the decomposition and signs thus cancel (man dspsv)
} }
if (det && r>0) *det = 0;
delete[] ipiv; delete[] ipiv;
if (r > 0 && b) laerror("singular matrix in linear_solve(SMat&, Mat*, double*"); if (r > 0 && b) laerror("singular matrix in linear_solve(SMat&, Mat*, double*");
} }