Index: tol/btol/vmatrix_type/vmat_chol.cpp
===================================================================
--- tol/btol/vmatrix_type/vmat_chol.cpp	(revision 1d83266c3400e5129f0be97d280b45e2089f5774)
+++ tol/btol/vmatrix_type/vmat_chol.cpp	(revision 23b4ae2a8d3203c09aa5f13c08f49e7323cc8518)
@@ -109,18 +109,21 @@
 ////////////////////////////////////////////////////////////////////////////////
   int BVMat::CholeskiFactor(const BVMat& X_, BVMat& L_, 
-                            ECholFacOri ori,  bool checkSymmetric)
+                            ECholFacOri ori,  
+                            bool checkSymmetric,
+                            bool forceNaturalOrder)
 ////////////////////////////////////////////////////////////////////////////////
 {
   if(!X_.CheckDefined("CholeskiFactor")) { return(-1); }
-
+  bool old_force = force_natural_order(forceNaturalOrder);
   BVMat* X__, *L__;
-  convertIfNeeded_cRt2cRs(X_,L_,X__,L__,"CholeskiSolve");
+  convertIfNeeded_cRt2cRs(X_,L_,X__,L__,"CholeskiFactor");
   BVMat &X = *X__,  &L = *L__;
   int r = X.Rows();
   int c = X.Columns();
+  int result = 0;
   if((ori==ECFO_X)&&(r!=c))
   {
     err_cannot_apply("CholeskiFactor",I2("non square","no cuadrada"),X);
-    return(-1);
+    result = -1;
   }
   else if((ori==ECFO_XtX)&&(r<c))
@@ -129,5 +132,5 @@
       I2("row deficient to build",
          "no tiene bastantes filas para construir")+ " S = X' X",X);
-    return(-2);
+    result =-2;
   }
   else if((ori==ECFO_XXt)&&(r>c))
@@ -136,39 +139,49 @@
       I2("column deficient to build",
          "no tiene bastantes columnas para construir")+" S = X X'",X);
-    return(-3);
-  }
-  if(ori!=ECFO_X) { checkSymmetric=false; }
-  bool isNotSymm = (checkSymmetric)?!X.IsSymmetric():false;
-  bool isOk = true;
-  bool isNotPosDef=false;
-  if(isNotSymm)
-  {
-    err_cannot_apply("CholeskiFactor",
-      I2("non symmetric","no simetrica"),X);
-    return(-4);
-  }
-  if(checkSymmetric && !isNotSymm) { ((BVMat&)X).CompactSymmetric(false); }
-  const StrCholFac* cholFac = FindCholFac(X.code_, ori);
-  if(cholFac)
-  {
-    (*cholFac->fun)(X, L, isOk, isNotPosDef);
+    result =-3;
   }
   else
   {
-    err_invalid_subtype("CholeskiFactor",X); 
-    isOk = false;
-  }
-  if(!isOk)
-  {
-    if(isNotPosDef)
+    if(ori!=ECFO_X) { checkSymmetric=false; }
+    bool isNotSymm = (checkSymmetric)?!X.IsSymmetric():false;
+    bool isOk = true;
+    bool isNotPosDef=false;
+    if(isNotSymm)
     {
       err_cannot_apply("CholeskiFactor",
-                       I2("non positive definite",
-                          "no definida positiva"),X);
-      L.Delete();
-      return(-5);
+        I2("non symmetric","no simetrica"),X);
+      result = -4;
     }
-    return(-6);
-  }    
+    else
+    {
+      if(checkSymmetric && !isNotSymm) { ((BVMat&)X).CompactSymmetric(false); }
+      const StrCholFac* cholFac = FindCholFac(X.code_, ori);
+      if(cholFac)
+      {
+        (*cholFac->fun)(X, L, isOk, isNotPosDef);
+      }
+      else
+      {
+        err_invalid_subtype("CholeskiFactor",X); 
+        isOk = false;
+      }
+      if(!isOk)
+      {
+        if(isNotPosDef)
+        {
+          err_cannot_apply("CholeskiFactor",
+                           I2("non positive definite",
+                              "no definida positiva"),X);
+          L.Delete();
+          result = -5;
+        }
+        else
+        {
+          result = -6;
+        }
+      } 
+    }
+  }
+  force_natural_order(old_force);
   return(0);
 }
@@ -176,5 +189,7 @@
 ////////////////////////////////////////////////////////////////////////////////
   int BVMat::CholeskiFactor(const BVMat& X, BVMat& L,
-                            const BText& oriName, bool checkSymmetric)
+                            const BText& oriName, 
+                            bool checkSymmetric,
+                            bool forceNaturalOrder)
 ////////////////////////////////////////////////////////////////////////////////
 {
@@ -190,24 +205,28 @@
     return(-3);
   }
-  return(CholeskiFactor(X,L,ori,checkSymmetric));
-}
-  
-////////////////////////////////////////////////////////////////////////////////
-  BVMat BVMat::CholeskiFactor(ECholFacOri ori, bool checkSymmetric)
+  return(CholeskiFactor(X,L,ori,checkSymmetric,forceNaturalOrder));
+}
+  
+////////////////////////////////////////////////////////////////////////////////
+  BVMat BVMat::CholeskiFactor(ECholFacOri ori, 
+                              bool checkSymmetric,
+                              bool forceNaturalOrder)
 //Matrix algebra operator
 ////////////////////////////////////////////////////////////////////////////////
 {
   BVMat aux;
-  CholeskiFactor(*this, aux, ori, checkSymmetric);
+  CholeskiFactor(*this, aux, ori, checkSymmetric, forceNaturalOrder);
   return(aux);
 };
 
 ////////////////////////////////////////////////////////////////////////////////
-  BVMat BVMat::CholeskiFactor(const BText& oriName, bool checkSymmetric)
+  BVMat BVMat::CholeskiFactor(const BText& oriName, 
+                              bool checkSymmetric,
+                              bool forceNaturalOrder)
 //Matrix algebra operator
 ////////////////////////////////////////////////////////////////////////////////
 {
   BVMat aux;
-  CholeskiFactor(*this, aux, oriName, checkSymmetric);
+  CholeskiFactor(*this, aux, oriName, checkSymmetric,forceNaturalOrder);
   return(aux);
 }
@@ -357,4 +376,19 @@
   x.s_.blasRdense_ = cholmod_solve
   (
+    CHOLMOD_Pt,
+    L.s_.chlmRfactor_,
+    b.s_.blasRdense_,
+    common_
+  );
+  return(0);    
+};
+
+////////////////////////////////////////////////////////////////////////////////
+  int BVMat::cRs_bRd_cholSolPt(const BVMat& L, const BVMat& b, BVMat& x)
+////////////////////////////////////////////////////////////////////////////////
+{
+  x.code_ = ESC_blasRdense;
+  x.s_.blasRdense_ = cholmod_solve
+  (
     CHOLMOD_P,
     L.s_.chlmRfactor_,
@@ -366,19 +400,4 @@
 
 ////////////////////////////////////////////////////////////////////////////////
-  int BVMat::cRs_bRd_cholSolPt(const BVMat& L, const BVMat& b, BVMat& x)
-////////////////////////////////////////////////////////////////////////////////
-{
-  x.code_ = ESC_blasRdense;
-  x.s_.blasRdense_ = cholmod_solve
-  (
-    CHOLMOD_Pt,
-    L.s_.chlmRfactor_,
-    b.s_.blasRdense_,
-    common_
-  );
-  return(0);    
-};
-
-////////////////////////////////////////////////////////////////////////////////
   int BVMat::cRs_cRs_cholSolPtLLtP(const BVMat& factor, const BVMat& b, BVMat& x)
 ////////////////////////////////////////////////////////////////////////////////
@@ -440,9 +459,13 @@
   cholmod_sparse* aux = cholmod_spsolve
   (
-    CHOLMOD_Pt,
+    CHOLMOD_P,
     factor.s_.chlmRfactor_,
     b_,
     common_
   );
+  if(b_ != b.s_.chlmRsparse_)
+  {
+    cholmod_free_sparse(&b_, common_); 
+  }
   x.s_.chlmRsparse_= cholmod_spsolve
   (
@@ -453,8 +476,4 @@
   );
   cholmod_free_sparse(&aux, common_); 
-  if(b_ != b.s_.chlmRsparse_)
-  {
-    cholmod_free_sparse(&b_, common_); 
-  }
   return(0);    
 };
@@ -477,7 +496,11 @@
     common_
   );
+  if(b_ != b.s_.chlmRsparse_)
+  {
+    cholmod_free_sparse(&b_, common_); 
+  }
   x.s_.chlmRsparse_= cholmod_spsolve
   (
-    CHOLMOD_P,
+    CHOLMOD_Pt,
     factor.s_.chlmRfactor_,
     aux,
@@ -485,8 +508,4 @@
   );
   cholmod_free_sparse(&aux, common_); 
-  if(b_ != b.s_.chlmRsparse_)
-  {
-    cholmod_free_sparse(&b_, common_); 
-  }
   return(0);    
 };
@@ -552,5 +571,5 @@
   x.s_.chlmRsparse_ = cholmod_spsolve
   (
-    CHOLMOD_P,
+    CHOLMOD_Pt,
     L.s_.chlmRfactor_,
     b_,
@@ -576,5 +595,5 @@
   x.s_.chlmRsparse_ = cholmod_spsolve
   (
-    CHOLMOD_Pt,
+    CHOLMOD_P,
     L.s_.chlmRfactor_,
     b_,