Changeset 3673741 in tol

Timestamp:

Jan 19, 2010, 2:27:45 PM (17 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

1941-docker-image-to-test, devel, master

Children:

4f21bb0

Parents:

d6488fb

Message:

Fixes #846
Fixes #754

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@1983 be4846ae-f051-483a-831f-b64e7d08527f

File:

• tol/btol/vmatrix_type/vmat_iterative.cpp (modified) (11 diffs)

tol/btol/vmatrix_type/vmat_iterative.cpp

@@ -79 +79 @@
   double advance;
   double relAdv;
-  double oldR;
   double R=r.FrobeniusNorm();
+  double oldR=R*2.0;
   BDat xNorm=x.FrobeniusNorm();
   BDat oldxNorm;
@@ -107 +107 @@
     pa = p*a;
     oldMaxpa = maxpa;
-    oldR = R;
-    R = r.FrobeniusNorm();
     oldxNorm = xNorm;
     x = x+pa;
     xNorm = x.FrobeniusNorm();
     xNormAdvance = xNorm-oldxNorm;
-
     maxpa = 0;
     for(i=0; i<cols; i++)
@@ -123 +120 @@
       }
     }
-
+    r = r - q*a;
+    oldR = R;
+    R = r.FrobeniusNorm();
     advance = R-oldR;
     relAdv = advance/oldR;
-    if(k)
-    {
-      if(S.IsUnknown()                    ) { break; }
-      if(maxpa.IsUnknown()                ) { break; }
-      if(advance             >=  0        ) { break; }
-      if(S                   <= tolerance ) { break; }
-      if(1+maxpa-1           == 0         ) { break; }
-/*
-      if(maxpa               <= chop ) { break; }
-      if(fabs(advance)       <= chop ) { break; }
-      if(fabs(relAdv)        <= chop ) { break; }
-      if(fabs(R)             <= chop ) { break; }
-      if(oldMaxpa            <= tol  ) { break; }
-      if(maxpa               <= tol  ) { break; }
-      if(fabs(advance)       <= tol  ) { break; }
-      if(xNormAdvance.Abs()  <= tol  ) { break; }
-*/
-    }
-    r = r - q*a;
   //s  = At*r;
     s = (r.T()*A).T();
     oldS = S;
     S  = s.SquaredSum();
+    if(S                   <= tolerance ) { break; }
+    if(S.IsUnknown()                    ) { break; }
+    if(maxpa.IsUnknown()                ) { break; }
+    if(advance             >=  0        ) { break; }
+    if(1+maxpa-1           == 0         ) { break; }
+/*
+    if(maxpa               <= chop ) { break; }
+    if(fabs(advance)       <= chop ) { break; }
+    if(fabs(relAdv)        <= chop ) { break; }
+    if(fabs(R)             <= chop ) { break; }
+    if(oldMaxpa            <= tol  ) { break; }
+    if(maxpa               <= tol  ) { break; }
+    if(fabs(advance)       <= tol  ) { break; }
+    if(xNormAdvance.Abs()  <= tol  ) { break; }
+*/
     b  = S.Value()/oldS.Value();
     p = p*b+s;
@@ -238 +233 @@
   int i, k;
   BVMat x, r;
-  if(!CheckOprProd(A,funName)) { return(x); }
   r.Convert(b0,ESC_blasRdense);
   if(!x0.Rows())
@@ -246 +240 @@
   else
   {
-    x .Convert(x0,ESC_blasRdense);
+    x.Convert(x0,ESC_blasRdense);
   }
   BVMat q;
@@ -259 +253 @@
   int rows = q.Rows();
   int cols = p.Rows();
-  //PRECONDITIONS: ENSURING DIMENSIONS 
+  /* PRECONDITIONS: ENSURING DIMENSIONS */
   if(!rows||!cols||!r.Columns()||!x.Columns()||
      (rows<cols)||(rows!=r.Rows())||(cols!=x.Rows())) 
@@ -265 +259 @@
     return(x); 
   }
-  
+
   if(chop<0) { chop = sqrt(DEpsilon()/cols); }
   double tolerance = chop*chop;
   if((rows<cols)||(rows!=r.Rows())) { return(x); }
-
+//BVMat At = A.T();
+//BVMat s = At*r;
   BVMat pa;
   BVMat rpa;
@@ -277 +272 @@
   double advance;
   double relAdv;
-  double oldR;
   double R=r.FrobeniusNorm();
+  double oldR=R*2.0;
   BDat xNorm=x.FrobeniusNorm();
   BDat oldxNorm;
@@ -290 +285 @@
   for(k=0; k<maxIter;  k++)
   {
+/*
     if(k && !(k%cols))
     {
@@ -298 +294 @@
         " advance="+advance);
     }
-    OprProd(A,p,false,q);
+*/
+    OprProd(A,x,false, q);
     q2 = q.SquaredSum();
     a = S.Value()/q2;
     pa = p*a;
     oldMaxpa = maxpa;
-    oldR = R;
-    R = r.FrobeniusNorm();
     oldxNorm = xNorm;
     x = x+pa;
     xNorm = x.FrobeniusNorm();
     xNormAdvance = xNorm-oldxNorm;
-
     maxpa = 0;
     for(i=0; i<cols; i++)
@@ -319 +313 @@
       }
     }
-
+    r = r - q*a;
+    oldR = R;
+    R = r.FrobeniusNorm();
     advance = R-oldR;
     relAdv = advance/oldR;
-    if(k)
-    {
-      if(S.IsUnknown()    ) 
-      { break; }
-      if(maxpa.IsUnknown()) 
-      { break; }
-      if(S <= tolerance   ) 
-      { break; }
-/* */
-      if(1+maxpa-1 == 0   ) 
-      { break; }
-      if(advance >   0    ) 
-      { break; }
-/* */
-    }
-    r = r - q*a;
   //s  = At*r;
     OprProd(A,r,true,s);
     oldS = S;
     S  = s.SquaredSum();
+    if(S                   <= tolerance ) { break; }
+    if(S.IsUnknown()                    ) { break; }
+    if(maxpa.IsUnknown()                ) { break; }
+    if(advance             >=  0        ) { break; }
+    if(1+maxpa-1           == 0         ) { break; }
+/*
+    if(maxpa               <= chop ) { break; }
+    if(fabs(advance)       <= chop ) { break; }
+    if(fabs(relAdv)        <= chop ) { break; }
+    if(fabs(R)             <= chop ) { break; }
+    if(oldMaxpa            <= tol  ) { break; }
+    if(maxpa               <= tol  ) { break; }
+    if(fabs(advance)       <= tol  ) { break; }
+    if(xNormAdvance.Abs()  <= tol  ) { break; }
+*/
     b  = S.Value()/oldS.Value();
     p = p*b+s;