Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members  

Solve.cc

Go to the documentation of this file.
00001 /*
00002     File:           Solve.cc
00003 
00004     Function:       Implements Solve.h
00005 
00006     Author(s):      Andrew Willmott
00007 
00008     Copyright:      (c) 1995-2000, Andrew Willmott
00009 
00010     Notes:          
00011 
00012 */
00013 
00014 
00015 #include "vl/Solve.h"
00016 
00035 TMReal SolveOverRelax(
00036             const TMat      &A,
00037             TVec            &x,
00038             const TVec      &b,
00039             TMReal          epsilon,
00040             TMReal          omega,
00041             Int             *steps
00042         )
00043 {
00044     Int     i, j, jMax;
00045     TMReal  sum;
00046     TMReal  diagonal, xOld, error;
00047 
00048     Assert(A.IsSquare(), "(SolveOverRelax) Matrix not square");
00049     j = 0;
00050     if (steps)
00051         jMax = *steps;
00052     else
00053         jMax = VL_MAX_STEPS;
00054         
00055     do
00056     {
00057         error = 0.0;
00058         for (i = 0; i < A.Rows(); i++)
00059         {
00060             sum = b[i] - dot(A[i], x);      
00061             diagonal = A[i][i];
00062             sum += diagonal * x[i];     
00063             // I.e., sum = b[i] - (A[i] * x - A[i,i])
00064             xOld = x[i];
00065             
00066             if (diagonal == 0)
00067                 _Warning("(SolveOverRelax) diagonal element = 0");
00068             else if (omega == 1.0)  // Gauss-Seidel 
00069                 x[i] = sum / diagonal;
00070             else                    // Overrelax
00071                 x[i] = Mix(xOld, sum / diagonal, (Real) omega); 
00072                 
00073             sum -= diagonal * xOld;
00074             error += sqr(sum);
00075         }
00076         j++;
00077     }
00078     while (error > epsilon && j < jMax);
00079 
00080     if (steps)
00081         *steps = j;
00082 
00083     return(error);
00084 }
00085 
00091 TMReal SolveOverRelax(
00092             const TSparseMat    &A,
00093             TVec                &x,
00094             const TVec          &b,
00095             TMReal              epsilon,
00096             TMReal              omega,
00097             Int                 *steps
00098         )
00099 {
00100     Int     i, j, jMax;
00101     TMReal  sum;
00102     TMReal  diagonal, xOld, error;
00103 
00104     Assert(A.IsSquare(), "(SolveOverRelax) Matrix not square");
00105     j = 0;
00106     if (steps)
00107         jMax = *steps;
00108     else
00109         jMax = VL_MAX_STEPS;
00110         
00111     do
00112     {
00113         error = 0.0;
00114         for (i = 0; i < A.Rows(); i++)
00115         {
00116             // sum = b[i] - (A[i] dot x - A[i,i])            
00117         
00118             sum = b[i] - dot(A[i], x);      
00119             diagonal = A[i].Get(i);
00120             sum += diagonal * x[i];     
00121             
00122             xOld = x[i];
00123             
00124             if (diagonal == 0)
00125                 _Warning("(SolveOverRelax) diagonal element = 0");
00126             else if (omega == 1)
00127                 x[i] = sum / diagonal;  // Gauss-Seidel 
00128             else
00129                 x[i] = Mix(xOld, sum / diagonal, (Real) omega); 
00130                 
00131             sum -= diagonal * xOld;
00132             error += sqr(sum);
00133         }
00134         j++;
00135     }
00136     while (error > epsilon && j < jMax);
00137 
00138     if (steps)
00139         *steps = j;
00140         
00141     return(error);
00142 }
00143 
00144 
00159 TMReal SolveConjGrad(
00160                 const TMat      &A,             // Solve Ax = b.
00161                 TVec            &x,
00162                 const TVec      &b,
00163                 TMReal          epsilon,        // how low should we go?
00164                 Int             *steps          // iterations to converge.
00165             )
00166 {
00167     Int     i, iMax;
00168     TMReal  alpha, beta, rSqrLen, rSqrLenOld, u;
00169     TVec    r(A.Rows());        // Residual vector, b - Ax 
00170     TVec    d(A.Rows());        // Direction to move x at each step 
00171     TVec    t(A.Rows());        // temp!
00172 
00173     Assert(A.IsSquare(), "(SolveConjGrad) Matrix not square");
00174     
00175     r = b;
00176     r -= A * x;             
00177     rSqrLen = sqrlen(r);
00178     d = r;                      
00179     i = 0;
00180     if (steps)
00181         iMax = *steps;
00182     else
00183         iMax = VL_MAX_STEPS;
00184         
00185     if (rSqrLen < epsilon)
00186         while (i < iMax)    
00187         {   
00188             i++;
00189             t = A * d;      
00190             u = dot(d, t); 
00191             
00192             if (u == 0)
00193             {
00194                 _Warning("(SolveConjGrad) d'Ad = 0");
00195                 break;
00196             }
00197             
00198             alpha = rSqrLen / u;        // How far should we go?            
00199             x += alpha * d;             // Take a step along direction d
00200             if (i & 0x3F)
00201                 r -= alpha * t; 
00202             else
00203             {
00204                 r = b;      // For stability, correct r every 64th iteration
00205                 r -= A * x;
00206             }
00207                 
00208             rSqrLenOld = rSqrLen;
00209             rSqrLen = sqrlen(r); 
00210             
00211             if (rSqrLen <= epsilon)
00212                 break;                  // Converged! Let's get out of here
00213             
00214             beta = rSqrLen/rSqrLenOld;
00215             d *= beta;                  // Change direction: d = r + beta * d
00216             d += r;
00217         }
00218         
00219     if (steps)
00220         *steps = i;
00221     
00222     return(rSqrLen);
00223 }
00224 
00230 TMReal SolveConjGrad(
00231                 const TSparseMat    &A,
00232                 TVec                &x,
00233                 const TVec          &b,
00234                 TMReal              epsilon,
00235                 Int                 *steps
00236             )
00237 {
00238     Int     i, iMax;
00239     TMReal  alpha, beta, rSqrLen, rSqrLenOld, u;
00240     TVec    r(b.Elts());        // Residual vector, b - Ax 
00241     TVec    d(b.Elts());        // Direction to move x at each step 
00242     TVec    t(b.Elts());
00243 
00244     Assert(A.IsSquare(), "(SolveConjGrad) Matrix not square");
00245     r = b;
00246     r -= A * x;             
00247     rSqrLen = sqrlen(r);
00248     d = r;                      
00249     i = 0;
00250     if (steps)
00251         iMax = *steps;      
00252     else
00253         iMax = VL_MAX_STEPS;
00254         
00255     if (rSqrLen > epsilon)              // If we haven't already converged...
00256         while (i < iMax)
00257         {   
00258             i++;
00259             t = A * d;      
00260             u = dot(d, t);
00261             
00262             if (u == 0.0)
00263             {
00264                 _Warning("(SolveConjGrad) d'Ad = 0");
00265                 break;
00266             }
00267             
00268             alpha = rSqrLen / u;        // How far should we go?
00269             x += alpha * d;             // Take a step along direction d
00270             r -= alpha * t; 
00271             rSqrLenOld = rSqrLen;
00272             rSqrLen = sqrlen(r); 
00273             
00274             if (rSqrLen <= epsilon)
00275                 break;                  // Converged! Let's get out of here
00276             
00277             beta = rSqrLen / rSqrLenOld;
00278             d = r + beta * d;           //  Change direction
00279         }
00280 
00281     if (steps)
00282         *steps = i;
00283     
00284     return(rSqrLen);
00285 }
00286 

Generated at Sat Aug 5 00:16:47 2000 for Class Library by doxygen 1.1.0 written by Dimitri van Heesch, © 1997-2000