added LINBCG and SPRSTP

src/libinit_ac.cpp

@@ -86,6 +86,9 @@ void LibInit_ac()
   new DLibFunRetNew(lib::sprsax_fun,string("SPRSAX"),2,sprsaxKey);
   new DLibFunRetNew(lib::fulstr_fun,string("FULSTR"),1);
   new DLibFunRetNew(lib::sprstp_fun,string("SPRSTP"),1);
+  const string linbcgKey[]={"ITER", "DOUBLE", KLISTEND};
+  const string linbcgIgnoreKey[]={"ITOL", "TOL", "ITMAX", KLISTEND};
+  new DLibFunRetNew(lib::linbcg_fun,string("LINBCG"),3,linbcgKey,linbcgIgnoreKey);
 
   const string spl1Key[]={"YP0","YPN_1","YP1","DOUBLE","HELP",KLISTEND}; //YP1 is old value for YP0
   new DLibFunRetNew(lib::spl_init_fun,string("SPL_INIT"),2,spl1Key);

src/sparse_matrix.cpp

@@ -188,7 +188,6 @@ namespace lib {
     SizeT nParam = e->NParam(1);
 	SPMATRowMajDbl *Mat=getFromPtr(e, 0);
 	SPMATRowMajDbl *res=new SPMATRowMajDbl((*Mat).transpose());
-	delete Mat;
 	return convertToPtr(res);
   }
   BaseGDL* sprsax_fun(EnvT* e) {
@@ -208,6 +207,36 @@ namespace lib {
 	return convertToGDL(Mat3);
   }
 
+  BaseGDL* linbcg_fun(EnvT* e) {
+	SizeT nParam = e->NParam(3);
+	SPMATRowMajDbl* A = getFromPtr(e, 0);
+	BaseGDL* p1 = e->GetParDefined(1); // Right Hand B 
+	DType varType = p1->Type();
+	if (p1->Dim().Rank() != 1) e->Throw("Argument " + e->GetString(1) + " must be a vector.");
+	if (varType == GDL_STRING) e->Throw("Argument " + e->GetString(1) + " must not be of STRING type.");
+	int m = p1->Dim(0);
+	if (m != A->cols()) e->Throw("Argument " + e->GetString(1) + " does not have correct size.");
+	DDoubleGDL* Bgdl = e->GetParAs<DDoubleGDL>(1);
+	Eigen::Map<Eigen::VectorXd> B(&(*Bgdl)[0], m);
+	BaseGDL* p2 = e->GetParDefined(2); // Initial Guess
+	varType = p2->Type();
+	if (p2->Dim().Rank() != 1) e->Throw("Argument " + e->GetString(2) + " must be a vector.");
+	if (varType == GDL_STRING) e->Throw("Argument " + e->GetString(2) + " must not be of STRING type.");
+	int n = p2->Dim(0);
+	if (n != m) e->Throw("Argument " + e->GetString(2) + " does not have correct size.");
+	DDoubleGDL* Xgdl = e->GetParAs<DDoubleGDL>(2);
+	Eigen::Map<Eigen::VectorXd> X(&(*Xgdl)[0], n); 
+	//solve ax=b
+	Eigen::SparseLU<Eigen::SparseMatrix<double>> sparseLU;
+    sparseLU.compute(*A);
+    if(sparseLU.info()!=Eigen::Success) e->Throw("Matrix decomposition failed.");
+	X=sparseLU.solve(B);
+    if(sparseLU.info()!=Eigen::Success) e->Throw("No solution found.");
+	DDoubleGDL* resD =new DDoubleGDL(n, BaseGDL::NOZERO);
+	Eigen::Map<Eigen::VectorXd>(&(*resD)[0], n) = X;
+	return resD;
+  }
+
   BaseGDL* fulstr_fun(EnvT* e){
 	SizeT nParam = e->NParam(1);
 //    SPMATRowMajDbl Mat=getFromIndexInMatVec(e, 0);

src/sparse_matrix.hpp

@@ -28,6 +28,6 @@ namespace lib {
   BaseGDL* sprsab_fun(EnvT* e);
   BaseGDL* sprstp_fun(EnvT* e);
   BaseGDL* fulstr_fun(EnvT* e);
-
+  BaseGDL* linbcg_fun(EnvT* e);
 }
 #endif

testsuite/test_sparse.pro

@@ -0,0 +1,20 @@
+; no a test, really. Just a suggestion for test, plu sthe fact
+; that it exercises the coverage.
+pro test_sparse
+
+  a = [[ 2.0,  1.0,  1.0],[ 4.0, -6.0,  0.0],[-2.0,  7.0,  2.0]]
+  z=sprsin(a, thresh=0.5)
+  zz=sprstp(z)
+  q=fulstr(z)
+  if abs(total(a-q)) gt 1E-6 then exit,status=1
+  res=sprsab(z,zz)
+  result=fulstr(res)
+  if total(result) ne 29 then exit,status=1
+  if total(sprsax(z,[1,1,1]))  ne 9 then exit,status=1
+  aludc=a
+  LUDC, aludc, index
+  b = [3,-8.0,10]
+  x = LUSOL(aludc, index, b)
+  r= LINBCG(SPRSIN(a), B, X)
+  if abs(total(r-x)) gt 1E-6 then exit,status=1
+end