Index: tol/bmath/blinalg/linalg.cpp
===================================================================
--- tol/bmath/blinalg/linalg.cpp	(revision ef245212dde1586694a4984931d62a8888c42f9d)
+++ tol/bmath/blinalg/linalg.cpp	(revision be138d0f4131aeb470d886e60b46f5c1ecd11c88)
@@ -292,4 +292,38 @@
 }
 
+//--------------------------------------------------------------------
+BMatrix<double> CholeskiMinimumResidualsSolve(
+  const BMatrix<double>& A, 
+  const BMatrix<double>& b)
+
+/*! Solves Ax=b using Cholesky decomposition
+ */
+//--------------------------------------------------------------------
+{
+  BMatrix<double> x;
+  BVMat A_, b_, x_;
+  A_.DMat2VMat(A, false,0.50,DEpsilon());
+  b_.DMat2dense(b);
+  if(!BVMat::CholeskiMinRes(A_,b_,x_))
+  {
+    x_.GetDMat(x);
+  }
+  return(x);
+}
+
+//--------------------------------------------------------------------
+BMat CholeskiMinimumResidualsSolve(const BMat& A, const BMat& b)
+
+/*! Solves Ax=b using Cholesky decomposition
+ */
+//--------------------------------------------------------------------
+{
+  BMat x;
+  BMatrix<double>& dA = (BMatrix<double>&)A;
+  BMatrix<double>& db = (BMatrix<double>&)b;
+  BMatrix<double>& dx = (BMatrix<double>&)x;
+  dx = CholeskiMinimumResidualsSolve(dA,db);
+  return(x);
+}
 
 //--------------------------------------------------------------------
@@ -379,5 +413,5 @@
     BMat A = B << Urc.T();
     BMat z(n-r,1); z.SetAllValuesTo(0);
-    BMat v = MinimumResidualsSolve(A, y) ;
+    BMat v = gsl_MinimumResidualsSolve(A, y) ;
     
     for(j=0; j<c; j++)
Index: tol/bmath/blinalg/tol_blinalg.h
===================================================================
--- tol/bmath/blinalg/tol_blinalg.h	(revision ef245212dde1586694a4984931d62a8888c42f9d)
+++ tol/bmath/blinalg/tol_blinalg.h	(revision be138d0f4131aeb470d886e60b46f5c1ecd11c88)
@@ -61,4 +61,12 @@
   int maxIter = -1
 );
+
+BMatrix<double> CholeskiMinimumResidualsSolve(
+  const BMatrix<double>& A, 
+  const BMatrix<double>& b);
+BMat CholeskiMinimumResidualsSolve(
+  const BMat& A, 
+  const BMat& b);
+
 
 BDat LinearLeastSquareSolve