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

added debug testing ipiv in linear_solve_do

Browse Source
This commit is contained in:
Jiri Pittner 2023-05-15 15:47:33 +02:00
parent 475de4869f
commit 5c494684ab
2 changed files with 14 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
nonclass.cc
View File

@ -26,6 +26,8 @@
#include "qsort.h" #include "qsort.h"
#include "fortran.h" #include "fortran.h"
#define IPIV_DEBUG
namespace LA { namespace LA {
@ -154,17 +156,20 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
//find out whether ipiv are numbered from 0 or from 1 //find out whether ipiv are numbered from 0 or from 1
int shift=1; int shift=1;
for (int i=0; i<n; ++i) if(ipiv[i]==0) shift=0; for (int i=0; i<n; ++i) if(ipiv[i]==0) shift=0;
#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) for (int i=0; i<n; ++i) if(i+shift != ipiv[i]) {*det = -(*det); ++iswap;} if(*det) for (int i=0; i<n; ++i) if(i+shift != ipiv[i]) {*det = -(*det); ++iswap;}
} }
/*
std::cout <<"iswap = "<<iswap<<std::endl;
if(det && r>0) *det = 0; if(det && r>0) *det = 0;
#ifdef IPIV_DEBUG
std::cout <<"iswap = "<<iswap<<std::endl;
std::cout <<"ipiv = "; std::cout <<"ipiv = ";
for (int i=0; i<n; ++i) std::cout <<ipiv[i]<<" "; for (int i=0; i<n; ++i) std::cout <<ipiv[i]<<" ";
std::cout <<std::endl; std::cout <<std::endl;
*/ #endif
delete [] ipiv; delete [] ipiv;
if (r>0 && B) laerror("singular matrix in lapack_gesv"); if (r>0 && B) laerror("singular matrix in lapack_gesv");

10
t.cc
View File

@ -1237,7 +1237,7 @@ cout <<r;
} }
if(0) if(1)
{ {
int n,m; int n,m;
cin>>n >>m; cin>>n >>m;
@ -1254,10 +1254,12 @@ for(int i=0;i<n;++i) for(int j=0;j<m;++j)
for(int i=0;i<n;++i) b[i] = i; for(int i=0;i<n;++i) b[i] = i;
NRVec<double> bb=b; NRVec<double> bb=b;
//cout <<a; cout <<a;
//cout <<b; //cout <<b;
NRMat<double>aa=a; NRMat<double>aa=a;
linear_solve(aa,bb); double d;
linear_solve(aa,bb,&d);
cout <<"det = "<<d<<endl;
//cout <<bb; //cout <<bb;
gmres(a,b,x,1,1e-10,100,1,0,1,0); gmres(a,b,x,1,1e-10,100,1,0,1,0);
//conjgrad(a,b,x,1,1e-10,200,1,0,1); //conjgrad(a,b,x,1,1e-10,200,1,0,1);
@ -2575,7 +2577,7 @@ cout <<test;
cout <<"Error = "<<(expitszsz-test).norm()<<endl; cout <<"Error = "<<(expitszsz-test).norm()<<endl;
} }
if(1) if(0)
{ {
NRVec<double> x({1,2,3}); NRVec<double> x({1,2,3});
NRVec<double> y({4,5,6}); NRVec<double> y({4,5,6});