% test_conjgrad.m  simple case 4 by 4  matrix
%                  x computed  using  b = 1,2,3,4
function test_conjgrad()
  diary('test_conjgrad_m.out');
  disp('test_conjgrad_m running with simeq and conjgrad')
  disp('solve Ax=b for x given A and b')

  A = [ 1 2 5 9 ; 2 6 9 10 ; 5 9 11 13 ; 9 10 13 19 ]
  b = [1 ; 2 ; 3 ; 4 ] % must be column
  x = A \ b;  % solve A x = b

  disp('simeq solution x=')
  disp(x)

  disp('initial guess for conjgrad')
  x0 = [0 ; 0 ; 0 ; 0 ]

  xc = conjgrad(A, b, x0, 1.0e-10, 9);

  disp('conjgrad solution xc=')
  disp(xc)

  disp('check A det, eig')
  Ad = det(A)
  [v e] = eig(A);
  disp(e)
  disp('l u lu(A)')
  [l u] = lu(A)
  disp('ichol(A)')
  disp('C = ichol(A) % fails because A is not sparse')
  disp('C = ichol(sparse(A))')
  disp('% still fails, non positive pivot')
  disp('test_conjgrad.m finished')
end % test_conjgrad
