determinant(): added workaround for openblas dgesv incompatibiility

noncblas.cc

@@ -648,6 +648,7 @@ int clapack_dgesv(const CBLAS_ORDER Order, const int N, const int NRHS,
 		double *A, const int lda, int *ipiv,
 		double *B, const int ldb)
 {
+std::cout <<"In MY clapack_dgesv, N and NRHS = "<< N<<" "<< NRHS<<"\n";
 	FINT INFO=0;
 	if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted");
 	//B should be in the same physical order, just transpose A in place and the LU result on output

nonclass.cc

@@ -140,15 +140,23 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
 	int r, *ipiv;
 	int iswap=0;
 
-	
+	if(nrhs==0) std::cout<<"Warning: some dgesv implementations might skip LU decomposition when nrhs==0\n";	
 	if (n==A.nrows() && A.nrows() != A.ncols()) laerror("linear_solve() call for non-square matrix");
 	A.copyonwrite();
 	ipiv = new int[A.nrows()];
 #ifdef IPIV_DEBUG
 	for(int i=0; i<A.nrows(); ++i) ipiv[i]=123456789;
+	std::cout <<"canary ipiv initialized\n";
+	std::cout <<"A before clapack_dgesv = "<<A<<std::endl;
+	std::cout <<"Active dimension n= "<<n<<std::endl;
 #endif
 	r = clapack_dgesv(CblasRowMajor, n, nrhs, &A(0,0), A.ncols(), ipiv, B , ldb);
-//	std::cout <<"A after clapack_dgesv = "<<A<<std::endl;
+#ifdef IPIV_DEBUG
+	std::cout <<"A after clapack_dgesv = "<<A<<std::endl;
+	std::cout <<"ipiv = ";
+	for (int i=0; i<n; ++i) std::cout <<ipiv[i]<<" ";
+	std::cout <<std::endl;
+#endif
 	if (r < 0) {
 		delete[] ipiv;
 		laerror("illegal argument in lapack_gesv");
@@ -177,10 +185,6 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
 	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;
@@ -212,6 +216,7 @@ extern "C" void FORNAME(dspsv)(const char *UPLO, const FINT *N, const FINT *NRHS
 
 static void linear_solve_do(NRSMat<double> &a, double *b, const int nrhs, const int ldb, double *det, int n)
 {
+ if(nrhs==0) std::cout<<"Warning: some dspsv implementations might skip LU decomposition when nrhs==0\n";
 	FINT r, *ipiv;
 	a.copyonwrite();
 	ipiv = new FINT[n];

nonclass.h

@@ -272,16 +272,9 @@ typename LA_traits<MAT>::normtype MatrixNorm(const MAT  &A, const char norm);
 template<class MAT>
 typename LA_traits<MAT>::normtype CondNumber(const MAT  &A, const char norm);
 
-
+#ifdef HAS_MKL
-//general determinant
+#define NO_OPENBLAS_WORKAROUND
-template<class MAT>
+#endif
-const typename LA_traits<MAT>::elementtype determinant(MAT a)//passed by value
-{
-typename LA_traits<MAT>::elementtype det;
-if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix");
-linear_solve(a,NULL,&det);
-return det;
-}
 
 //general determinant destructive on input 
 template<class MAT>
@@ -289,10 +282,27 @@ const typename LA_traits<MAT>::elementtype determinant_destroy(MAT &a) //passed
 {
 typename LA_traits<MAT>::elementtype det;
 if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix");
+
+//for openblas 0.3.31 we have to fake some RHS otherwise LU decomp. is not performed
+#ifdef NO_OPENBLAS_WORKAROUND
 linear_solve(a,NULL,&det);
+#else
+//fake rhs
+NRVec<typename LA_traits<MAT>::elementtype> b(a.ncols());
+for(int i=0; i<b.size(); ++i) b[i] = (typename LA_traits<MAT>::elementtype)i;
+linear_solve(a,b,&det);
+#endif
 return det;
 }
 
+//general determinant
+template<class MAT>
+const typename LA_traits<MAT>::elementtype determinant(MAT a)//passed by value
+{
+a.copyonwrite();
+return determinant_destroy(a);
+}
+
 
 //------------------------------------------------------------------------------
 // solves set of linear equations using gesvx