Example: Conjugate Gradient Algorithm :: SAS/IML(R) 12.1 User's Guide

Example: Conjugate Gradient Algorithm

Consider the following small example: where

$A = \left( \begin{array}{cccc} 3 & 1 & 0 & 0 \\ 1 & 4 & 1 & 3 \\ 0 & 1 & 10 & 0 \\ 0 & 3 & 0 & 3 \end{array}\right)$

and the vector of right-hand sides $b = (1\, 1\, 1\, 1)^ T.$ Since the matrix is positive definite and symmetric, you can apply the conjugate gradient algorithm to solve the system. Remember that you must specify only the lower-triangular part of the matrix (so row indices must be greater than or equal to the corresponding column indices.)

The code for this example is as follows:

   /* value   row   col */
   A = { 3     1     1,
         1     2     1,
         4     2     2,
         1     3     2,
         3     4     2,
        10     3     3,
         3     4     4 };

   /* right-hand sides */
   b = {1, 1, 1, 1};

   /* desired solution tolerance (optional) */
   tol = 1e-7;

   /* maximum number of iterations (optional) */
   maxit = 200;

   /* allocate iteration progress (optional) */
   hist = j(50, 1);

   /* provide an initial guess (optional) */
   start = {2, 3, 4, 5};

   /* invoke conjugate gradient method */
   call itsolver (
      x, st, it,          /* output parameters */
      'cg', A, b, 'ic',   /* input parameters */
      tol,                /* optional control parameters */
      maxit,
      start,
      hist
    );

   print  x;  /* print solution */
   print st;  /* print solution tolerance */
   print it;  /* print resultant number of iterations */

Notice that the example used an incomplete Cholesky preconditioner (which is recommended). Here is the program output:

The conjugate gradient method converged successfully within three iterations. You can also print out the hist (iteration progress) array. Different starting points result in different iterative histories.