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

Compare commits

base: jiri:ce2400e70341bf022f333afded5a27eede571996
Branches Tags
jiri:master
..
compare: jiri:5c494684ab116abc5c9c9ca9bb7dcff11f305e46
Branches Tags
jiri:master

No commits in common. "ce2400e70341bf022f333afded5a27eede571996" and "5c494684ab116abc5c9c9ca9bb7dcff11f305e46" have entirely different histories.
ce2400e703 ... 5c494684ab

2 changed files with 4 additions and 30 deletions
Show Stats Expand all files Collapse all files

2
laerror.cc
View File

@ -174,7 +174,6 @@ extern "C" {
#include "cblas.h"
}
#include <stdarg.h>
#ifndef AVOID_DUPLICATE_CBLAS_XERBLA
extern "C" void cblas_xerbla(int p, const char *rout, const char *form, ...)
{
va_list argptr;
@ -190,7 +189,6 @@ extern "C" void cblas_xerbla(int p, const char *rout, const char *form, ...)
va_end(argptr);
laerror("terminating in cblas_xerbla");
}
#endif
extern "C" int cblas_errprn(int ierr, int info, char *form, ...) {
char msg0[1024], *msg;

32
nonclass.cc
View File

@ -26,7 +26,7 @@
#include "qsort.h"
#include "fortran.h"
#undef IPIV_DEBUG
#define IPIV_DEBUG
namespace LA {
@ -155,11 +155,7 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
for (int i=0; i<n; ++i) {double t=A[i][i]; if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *=t;}
//find out whether ipiv are numbered from 0 or from 1
int shift=1;
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 clapack_dgesv");
}
for (int i=0; i<n; ++i) if(ipiv[i]==0) shift=0;
#ifdef IPIV_DEBUG
std::cout <<"shift = "<<shift<<std::endl;
#endif
@ -258,7 +254,6 @@ 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)
{
int r, *ipiv;
int iswap=0;
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()");
@ -275,31 +270,12 @@ void linear_solve(NRMat< std::complex<double> > &A, NRMat< std::complex<double>
delete[] ipiv;
laerror("illegal argument in lapack_gesv");
}
if (det && r==0) {
if (det && r>=0) {
*det = A[0][0];
for (int i=1; i<A.nrows(); ++i) *det *= A[i][i];
int shift=1;
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;}
for (int i=0; i<A.nrows(); ++i) *det = -(*det);
}
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;
if (r>0 && B) laerror("singular matrix in zgesv");
}