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
