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
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             {
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
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             {
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 1.1.0 written by Dimitri van Heesch, © 1997-2000