SOLVERS

The ``${\tt damped\_solver}$'' routine presented above calls the routine ${\tt <type\gt arraysolver}$. This routine could solve the problem by any of a wide number of methods. The author of the operator class and the author who sets up the particular problem do not have to modify any of the solver code, and the author of a solver routine does not need to assume anything about the operator that will be used other than they must be derived from one of the base operator classes.

The design of the operator and space classes separates the jobs of implementing a particular operator and writing a routine that solves the linear problem. The scientists can concentrate on writing code specific to their domain of expertise and then ask the local scientific computing specialist to write a spiffy solver routine for them. In contrast, when writing in Fortran we typically need to write a new solver routine for each new operator that we implement.

Routines such as solvers and dot-product tests can be written once, and then just called with operators and spaces as arguments. Below is a Hestenes conjugate gradient solver to solve the equation ${\bf y} \approx {\bf A}\,{\bf x}$ for x, where y is a known ${\tt <type\gt space}$ and A is of class ${\tt <type\gt op}$.Note that the only ${\tt <type\gt space}$ member functions used here are the functions that scale by a float, add another ${\tt <type\gt space}$ and take the norm-squared.

 // The hestenes routine is passed an initial output space "x"
 // the rhs space "y", the operator "A" and the maximum number of
 // iterations to use.
 // It returns an estimate of the solution.
 // N.B. this routine does not use the initial solution other than
 // using it to define the model space.

 <type>space <type>hestenes(<type>space x,   
                            <type>space & y,
                            <type>op & A,
                            int maxiter ) {

   <type>space step,cstep,grad,resid;
   float gamma,gamma1,alpha,beta,ynorm;

   x = 0;
   ynorm = y.norm2();
   resid = y;
   grad = A.Conjugate(resid);
   step = grad;
   gamma1 = step.norm2();

   int i=0;

   while ((i++<maxiter)  && (resid.norm() / ynorm > .00001)  ) {
      cstep = A.Forward(step);
      alpha = gamma1 / cstep.norm2();
      x += step * alpha;
      resid = resid - cstep * alpha;
      grad = A.Conjugate(resid);
      gamma = grad.norm2();
      beta = gamma / gamma1;
      gamma1 = gamma;
      step = step * beta + grad;

   }
   return x;
}

One major task that lies ahead is to implement more types of solver routines. We wish to compare different algorithms for non-Hermitian problems, the use of preconditioners and for nonlinear optimizations. Since these routines can be written in terms of abstract operator and space classes we will only have to write these routines once, they will then be available for use on any problem. In an ideal world we would hope to collaborate with members of the computer science department in writing these routines. We would be able to provide many examples of real world problems and a framework in which to use them and they would be able to provide new algorithms for solving the problems.