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

ipiv diagnostics also for zgesv

Browse Source
This commit is contained in:
Jiri Pittner 2023-05-15 16:53:25 +02:00
parent abdf6cdd73
commit ce2400e703
1 changed files with 22 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
nonclass.cc
View File

@ -258,6 +258,7 @@ extern "C" void FORNAME(zgesv)(const int *N, const int *NRHS, double *A, const i
void linear_solve(NRMat< std::complex<double> > &A, NRMat< std::complex<double> > *B, std::complex<double> *det, int n) void linear_solve(NRMat< std::complex<double> > &A, NRMat< std::complex<double> > *B, std::complex<double> *det, int n)
{ {
int r, *ipiv; int r, *ipiv;
int iswap=0;
if (A.nrows() != A.ncols()) laerror("linear_solve() call for non-square matrix"); if (A.nrows() != A.ncols()) laerror("linear_solve() call for non-square matrix");
if (B && A.nrows() != B->ncols()) laerror("incompatible matrices in linear_solve()"); if (B && A.nrows() != B->ncols()) laerror("incompatible matrices in linear_solve()");
@ -274,12 +275,31 @@ void linear_solve(NRMat< std::complex<double> > &A, NRMat< std::complex<double>
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 = A[0][0];
for (int i=1; i<A.nrows(); ++i) *det *= A[i][i]; for (int i=1; i<A.nrows(); ++i) *det *= A[i][i];
//change sign of det by parity of ipiv permutation int shift=1;
for (int i=0; i<A.nrows(); ++i) *det = -(*det); for (int i=0; i<n; ++i)
{
if(ipiv[i]==0) shift=0;
if(ipiv[i]<0 || ipiv[i]>n) 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
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");
} }