Index: tol/stdlib/math/linalg/_solve_precond_unsym.tol
===================================================================
--- tol/stdlib/math/linalg/_solve_precond_unsym.tol	(revision a9abf540831dafe29f879dba90d09bf6ff972464)
+++ tol/stdlib/math/linalg/_solve_precond_unsym.tol	(revision a9abf540831dafe29f879dba90d09bf6ff972464)
@@ -0,0 +1,179 @@
+//////////////////////////////////////////////////////////////////////////////
+// FILE    : _solve_precond_unsym.tol
+// PURPOSE : Soves an unsymmetric overdetermined linear systems by diagonal 
+//           expansion preconditioning of previously normalized system
+//////////////////////////////////////////////////////////////////////////////
+
+Text _.autodoc.member.SolvePrecondUnsym = 
+"Soves an unsymmetric overdetermined linear system, that could be large, "
+"sparse and ill-conditioned, in a robust and fast way by applying diagonal "
+"additive preconditioner of previously normalized system foilowing these steps:\n"
+"1. Original unsymmetric system M*z = y is previously normalized by non null "
+"elements of root squares of diagonal elements of M'M that will called D:\n"
+"  H*u = y where H=M*Di; Di=D^-1; u=D*x;\n"
+"2. A digonal expansion preconditioner Hc=(H<<c*I)^-1 will be built to be \n"
+"  Hc'*Hc=(S+c*I)^-1*S; where S=H'*H \n"
+"close to identity and number c is found to be the lesser than H<<c*I "
+"admits Choleski decomposition S+c*I = Lc*Lc'. \n"
+"3. Preconditioned system is solved by Lanczos iterative method with "
+"partial reorthogonalization MinimumResidualsSolve. Preconditioner is "
+"stored as a Code operator using just sparse matrix Lc to avoid store "
+"((S+c*I)^-1)*S that is not neccessary so sparse as S nor L.\n"
+"4. Solution is denormalized to solve original problem.\n"
+"Using instructions to solve M*x=y :\n"
+"  //Builds preconditioning handler\n"
+"  NameBlock SPDP = SymPosDefPrecond(H); \n"
+"  //Finds optimal preconditioner\n"
+"  Real SPDP::ScaleOptimize(1.E-7); \n"
+"  //Solves the system\n"
+"  NameBlock solution = SPDP::Solve(y,1.E-10); \n" 
+"  //Takes the solution matrix\n"
+"  VMatrix x = solution::_.x;\n"
+"\n";
+
+NameBlock SolvePrecondUnsym(VMatrix M, Real do.normalize)
+{[[
+  //Original s  //Original symmetric system to be solved M*x = y
+  VMatrix _.M = M;
+  Real _.do.normalize = do.normalize;
+  Real _.m = VRows(_.M);
+  Real _.n = VColumns(_.M);
+
+  //Returns a diagonal normalizer VMatrix for a symmetric VMatrix
+  VMatrix Normalizer(VMatrix M)
+  {
+    Real m = VRows(M);
+    Real n = VColumns(M);
+    If(!_.do.normalize, Eye(n), {
+      VMatrix di = Mat2VMat(SetCol(For(1,n,Real(Real j) 
+      { 
+        VMatrix Mj = SubCol(M,[[j]]);
+        Real d = Sqrt(VMatDat(MtMSqr(Mj),1,1));
+        If(d>0,1/d,1)
+      })));
+      Eye(n,n,0,di)
+    })
+  };
+
+  VMatrix _.Di = Normalizer(_.M);
+
+  //Normalized matrix
+  VMatrix _.H = If(!_.do.normalize, _.M, _.M*_.Di );
+
+  //Scale parameter of diagonal preconditioner
+  Real    _.c   = 0;
+  //Tolerance of c parameter
+  Real    _.tol.scale = 1.E-7;
+  //Preparing to store Sc = S+c*I
+  VMatrix _.Hc  = Rand(0,0,0,0);
+  //Preparing to store Lc matching Sc = Lc*Tra(Lc)
+  VMatrix _.Lc  = Rand(0,0,0,0);
+  //Intenta factorizar S+c*I y si lo consigue devuelve 1 y si no 0
+  Real ScaleEval(Real c)
+  {
+    Text _MID = "[SolvePrecondUnsym()::ScaleEval()] ";
+    VMatrix Dc = Eye(_.n,_.n,0,c);
+    VMatrix Hc = _.H << Dc;
+    Real old = Show(0,"ALL");
+    VMatrix Lc = CholeskiFactor(Hc,"XtX");
+    Real Show(old,"ALL");
+    Real ok = VRows(Lc)==_.n;
+    Real If(ok,
+    {
+      VMatrix _.Hc := Hc;
+      VMatrix _.Lc := Lc;
+      1
+    });
+    WriteLn(_MID+"("<<c+") = "<<ok); 
+    ok
+  };
+  //Aplica el método recursivo de la bisectriz en un intervalo buscando el
+  //menor valor de c para el que S+c*I es definida positiva
+  Real ScaleBisect(Real a, Real b)
+  {
+    If(b-a<=_.tol.scale, b,
+    {
+      Real c = (a+b)/2;
+      Real cOk = ScaleEval(c);
+      If(cOk, ScaleBisect(a,c), ScaleBisect(c,b))
+    })
+  };
+  //Busca un intervalo inicial [a,b] tal que S+a*I no sea definida positiva y
+  //S+b*I sí lo sea y luego aplica el método de la bisectriz para encontrar la
+  //frontera 
+  Real ScaleOptimize(Real tolerance) 
+  {
+    Real _.tol.scale := tolerance;
+    Real a = 0;
+    Real aOk = ScaleEval(a);
+    Real _.c := If(aOk, a,
+    {
+      Real b = _.tol.scale;
+      Real bOk = ScaleEval(b);
+      While(!bOk,
+      {
+        a := b;
+        b :=2*b;
+        bOk := ScaleEval(b)
+      });
+      ScaleBisect(a,b)
+    });
+    _.c
+  };
+
+  //Linear operator of preconditioned overdetermined system
+  //VBR: MRS doesnŽt works fine with this unsymmetric system
+  VMatrix PreCndLinOpr(VMatrix u, Real traspose)
+  {
+    If(!traspose,
+    {
+      _.H*CholeskiSolve(_.Lc,u,"LtP")
+    },
+    {
+      VMatrix Htu = Tra(Tra(u)*_.H);
+      CholeskiSolve(_.Lc,Htu,"LtP")
+    })
+  };
+
+  //Linear operator of preconditioned symmetric system HŽH*x=HŽy
+  VMatrix PreCndSymLinOpr(VMatrix u, Real traspose)
+  {
+     VMatrix v = Tra(Tra(_.H*u)*_.H);
+     CholeskiSolve(_.Lc,v,"PtLLtP") 
+  };
+
+  NameBlock Solve(VMatrix y, Real tolerance, Real forze.lanczos)
+  {
+    [[
+      Text _MID = "[SolvePrecondUnsym()::Solve()] ";
+      VMatrix _.y = y;
+      Real    _.tol.solver = tolerance;
+      Real    _.tol.mrs    = (tolerance^2)/(_.m*_.n);
+      VMatrix _.v = Tra(Tra(y)*_.H);
+      VMatrix _.u0 = CholeskiSolve(_.Lc,_.v,"PtLLtP");
+      VMatrix _.u  = If(And(!forze.lanczos, !_.c), _.u0,
+      {
+        MinimumResidualsSolve(PreCndSymLinOpr,_.u0,_.tol.mrs,_.u0)
+      });
+    //VMatrix _.u  = MinimumResidualsSolve(PreCndLinOpr,_.y,_.tol.mrs,_.u1);
+      VMatrix _.x = If(!_.do.normalize,_.u,_.Di*_.u);
+      VMatrix _.residuals = _.M*_.x-_.y;
+      Real    _.stdErr = Sqrt(VMatDat(MtMSqr(_.residuals),1,1)/_.m);
+      Real    _.maxErr = VMatMax(Abs(_.residuals));
+      Real    _.isOk = 
+      {
+        Real ok.1 = VRows(_.x)==_.n;
+        Real ok.2 = _.maxErr <= _.tol.solver;
+        If(!ok.1,WriteLn(_MID+"Wrong dimensions" ,"W"));
+        If(!ok.2,WriteLn(_MID+"Solution has maximum absolute error "<<_.maxErr+
+          " great than specified tolerance "<<_.tol.solver,"W"));
+        And(ok.1,ok.2)
+      }
+    ]]
+  }
+]]};
+
+/////////////////////////////////////////////////////////////////////////////
+PutDescription(_.autodoc.member.SolvePrecondUnsym,SolvePrecondUnsym);
+/////////////////////////////////////////////////////////////////////////////
+