% Testing program for gmres_x, a version of gmres organized as suggested by
% R.J. Hanson and D.R. Kincaid,"Notes on GMRES Algorithm Organization",
% Technical Report TR-05-05, March 2005,
% Computer Sciences Department, University of Texas at Austin.
% Date = 18 April 2005.

disp ' Solve for <Gabriel''s pillow>: Del u + 1 = 0, u = 0 on boundary of unit square'
disp ' To start try 32 grid points and 16 GMRES iterations.'

while 1
  disp ' '
  while 1
    n = input(' Enter number ( > 1) of inner grid points (in x or y direction) = ');
    if (n > 1)
      break
      end
  end
  h = 1/(n+1);

  K = n^2;
  % This allocates the sparse array that holds the difference equation
  % matrix.  To first order, there are 5 non-zero bands.
  Z = spalloc(K,K,5*K);
  % The loops for writing the sparse difference equations are broken into
  % three groups.  The main one writes the equations whose entries do not
  % touch the boundary.  The alternate sets traverse the edges and
  % then the corners.

  % South edge: (i==1)
  for j=2:n-1
    p = j;
    Z(p,p-1) = 1;
    Z(p,p+1) = 1;
    Z(p+n,p) = 1;
  end
  % East edge: (j==n)
  for i=2:n-1
    p = (i-1)*n+n;
    Z(p,p-1) = 1;
    Z(p-n,p) = 1;
    Z(p+n,p) = 1;
  end
  % North edge : (i==n)
  for j=2:n-1
    p = (n-1)*n+j;
    Z(p,p+1) = 1;
    Z(p,p-1) = 1;
    Z(p-n,p) = 1;
  end
  % West edge: (j==1)
  for i=2:n-1
    p = (i-1)*n+1;
    Z(p,p+1) = 1;
    Z(p-n,p) = 1;
    Z(p+n,p) = 1;
  end

  % SW corner: (i==1 and j==1)
  p = 1;
  Z(p,p+1) = 1; Z(p+n,p) = 1;
  % SE corner: (i==1 and j==n)
  p = n;
  Z(p,p-1) = 1; Z(p+n,p) = 1;
  % NE corner: (i==n, j==n)
  p = n^2;
  Z(p-n,p) = 1; Z(p,p-1) = 1;
  % NW corner: (i==n, j==1)
  p = n*(n-1)+1;
  Z(p,p+1) = 1; Z(p-n,p) = 1;
  
% All nodes:
  for i=1:n
    for j=1:n
      p = (i-1)*n+j;
      Z(p,p) = -4;
    end
  end

  % Inside:
  for j=2:n-1
    for i=2:n-1
      p = (i-1)*n+j;
      Z(p,p-1) = 1;
      Z(p,p+1) = 1;
      Z(p-n,p) = 1;
      Z(p+n,p) = 1;
    end
  end
  % Use tridiagonal principal minor of matrix as factor for pre-conditioner.
  W = spalloc(K,K,3*K);

  for j=2:K
    W(j-1,j) = Z(j-1,j);
    W(j,j) = Z(j,j);
    W(j,j-1) = Z(j,j-1);
  end
  W(1,1) = Z(1,1);
  [L,U,P] = lu(W);

  % Right hand side with sign change and scaling applied:
  b = ones(K,1);
  b = -h^2*b;
  tol = 1.E-4;
  Test_num = 2;
    
  [z0,sigma_sequence,iterations] = gmres_x(Z,b,L,U,P,tol);
  out=my_summary(n,z0,sigma_sequence,iterations,'test_y');

  n = input('Enter space for more computing or 0 to stop... ');
  if(n==0)
    break
  end
  clear sigma_sequence;
end