For  Math. 128 B                         19 April 2004

Assignment:  Program the solution of  A4*z = b  by
Conjugate Gradient Iteration  without ever storing
the matrix  A4 .  Inherit its dimensions from  b .

Then reprogram the iteration after preconditioning
the system  A4*z = b  by the matrix implicit in 
invd2 ,  which is presumed very roughly proportional
to the inverse of  A4 .

Incidentally,  why does  invd2  invert  d2 ?
How much good does preconditioning by  invd2  do?



function  y = a4(x)
%  y = a4(x)  is the matrix product  A4*x  of 
%  a symmetric positive definite matrix  A4 
%  and a column vector  x(:)  whose  length(x) 
%  determines the dimensions of ...
%          [ 4  -1  -1                          ]
%          [-1   4  -1  -1                      ]
%   A4  =  [-1  -1   4  -1  -1                  ]
%          [    -1  -1   4  -1  -1              ]
%          [        -1  -1   4  -1  -1          ]
%          ...          ... ... ... ... ... ... ]
n = length(x) ;
y = conv([-1 -1 4 -1 -1], x(:)) ;
y = y(3:n+2) ;




function Y = d2(X)
%  y = d2(x) = -diff(diff([0;x;0]))  for any column  x .
%  Y = d2(X)  applies  d2  to every column of matrix  X .
%                [ 2  -1                     ]
%                [-1   2  -1                 ]
%  Matrix  d2 =  [    -1   2  -1             ]
%                [        -1  ... ...        ]
%                [            ...   2  -1    ]
%                [                 -1   2  -1]
%                [                     -1   2]
%  is symmetric positive definite tridiagonal.
[m,n] = size(X) ;  on = zeros(1,n) ;
Y = -diff(diff( [on;X;on] )) ;




function  X = invd2(B)
%  X = invd2(d2(X)) = d2(invd2(X))
[N,m] = size(B) ;  uN = ones(N,1) ; 
S = cumsum(B) ;  s = sum(S)/(N+1) ;
X = flipud(cumsum(flipud(S - uN*s))) ;


                             Prof. W. Kahan