Index: tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_stdLinearBlock.tol
===================================================================
--- tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_stdLinearBlock.tol	(revision c6d5b9eb509b88828dc44c2e3e33c1982d66ea2b)
+++ tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_stdLinearBlock.tol	(revision 50d9df26c6ee2bd9649bcf379d62e159e26a6fe1)
@@ -38,4 +38,6 @@
   //Matrix of border of constraining inequations
   VMatrix _.a  = Rand(0,0,0,0);
+  //Transformed restrictions
+  VMatrix _.D = Rand(0,0,0,0);
   //Cholesky decomposition of current information matrix that is not
   //explicitely stored
@@ -48,4 +50,7 @@
   //factorizations
   Real design.is.constant = False;
+  //Matrix of linear system in last iteration
+  VMatrix _old.X  = Rand(0,0,0,0); 
+  Real _changed.X = True;
 
   //Maximum accepted rounding error for Choleski operations
@@ -115,24 +120,34 @@
     If(And(design.is.constant,VRows(_.L)), 0,
     {
-      Real sparse = VNonNullCells(X)/(VRows(X)*VColumns(X));
-      Real numErr = Copy(NError);
-      VMatrix _.L := If(sparse<.20, 
-      {
-        CholeskiFactor(X,"XtX")
-      },
-      {
-        VMatrix X_   = Convert(X,"Blas.R.Dense");
-        VMatrix XtX_ = MtMSqr(X_);
-        VMatrix XtX  = Convert(XtX_,"Cholmod.R.Sparse");
-        CholeskiFactor(XtX,"X")
+      Real _changed.X := If(Or(VColumns(X)!=VColumns(_old.X),
+                               VRows   (X)!=VRows   (_old.X)), -1,
+      {
+        1+VMatMax(Abs(X-_old.X))-1
       });
-      Real If(Copy(NError)>numErr,
-      {
-        Text path = TolAppDataPath+"tmp/"+ToName(_MID)+".oza";
-        WriteLn("Failed factorization matrix will be stored into "+path);
-        Ois.Store([[X]],path)
-      });
-      VMatrix If(_.R, _.Ls := Convert(_.L, "Cholmod.R.Sparse"));
-      1
+      If(!_changed.X ,0,
+      {
+        WriteLn("TRACE"+_MID+" factorize changed="<<_changed.X );
+        Real sparse = VNonNullCells(X)/(VRows(X)*VColumns(X));
+        Real numErr = Copy(NError);
+        VMatrix _.L := If(sparse<.20, 
+        {
+          CholeskiFactor(X,"XtX")
+        },
+        {
+          VMatrix X_   = Convert(X,"Blas.R.Dense");
+          VMatrix XtX_ = MtMSqr(X_);
+          VMatrix XtX  = Convert(XtX_,"Cholmod.R.Sparse");
+          CholeskiFactor(XtX,"X")
+        });
+        Real If(Copy(NError)>numErr,
+        {
+          Text path = TolAppDataPath+"tmp/"+ToName(_MID)+".oza";
+          WriteLn("Failed factorization matrix will be stored into "+path);
+          Ois.Store([[X]],path)
+        });
+        VMatrix If(_.R, _.Ls := Convert(_.L, "Cholmod.R.Sparse"));
+        VMatrix _old.X := X;
+        1
+      })
     }) 
   };
@@ -199,11 +214,11 @@
   {
   //WriteLn("TRACE "+_MID+" draw.constrained 1.1 Max(A*b0-a)="<<VMatMax(_.A*_.blockDef::_.store - a));
-    VMatrix   d   = _.a - _.A*sr::nu;
-    VMatrix   D   = Tra(CholeskiSolve(_.L,_.At,"PtL"));
-    VMatrix   z0  = Tra(Tra(_.blockDef::_.store-sr::nu)*_.Ls);
-  //WriteLn("TRACE "+_MID+" draw.constrained 1.2 Max(D*z0-d)="<<VMatMax(D*z0 - d));
-    VMatrix   z   = TruncStdGaussian(D,d,z0,1,
-                     _.blockDef::_.config::bsr.truncMNormal.gibbsNumIter);
-    Real distBorder = VMatMax(D*z-d);
+    VMatrix d   = _.a - _.A*sr::nu;
+    VMatrix If(_changed.X, _.D  := Tra(CholeskiSolve(_.L,_.At,"PtL")));
+    VMatrix z0  = Tra(Tra(_.blockDef::_.store-sr::nu)*_.Ls);
+  //WriteLn("TRACE "+_MID+" draw.constrained 1.2 Max(D*z0-d)="<<VMatMax(_.D*z0 - d));
+    VMatrix z = TruncStdGaussian(_.D,d,z0,1,
+                 _.blockDef::_.config::bsr.truncMNormal.gibbsNumIter);
+    Real distBorder = VMatMax(_.D*z-d);
     Real If(distBorder<=0,
     {
@@ -213,9 +228,9 @@
     {
       WriteLn(_MID+" Found a point that doesn't match constrains Max(D*z-d)="<<
-        VMatMax(D*z-d)+"\n Trying again");
+        VMatMax(_.D*z-d)+"\n Trying again");
       If(!VRows(_.beta0),0,
       {
         VMatrix z0 := Tra(Tra(_.beta0-sr::nu)*_.Ls);
-        VMatrix z  := TruncStdGaussian(D,d,z0,1,
+        VMatrix z  := TruncStdGaussian(_.D,d,z0,1,
                      _.blockDef::_.config::bsr.truncMNormal.gibbsNumIter);
         1