function [outstring] = driveGmres_Fom(instring);
%
%  Usage  disp(driveGmres_Fom(instring));
%
%  Demonstrates GMRES convergence on various types of matrices
%
%          for the linear system     Ax = b
%
%  Graphical displays are given of 
%
%          norm(b - Ax) convergence history
%          for Gmres and Fom
%
%          placement of roots of GMRES and FOM polynomials
%
%          After each 10 iterations a surface and contour
%          plot of log(abs(p(z)))  where p(z) is the GMRES
%          polynomial. 
%
%  Cases are specified with 
%
%  Input: A character string
%
%         instring  =  'random'
%
%                        specifies order 50 random system
%                        shows slow convergence
%
%         instring  =  'non-normal'
%
%                        specifies an order 120  non-normal system
%                        with eigenvalues clustered at lambda = 1,2,3
%
%         instring  =  'normal'
%
%                        specifies an order 120  normal system
%                        with eigenvalues clustered at lambda = 1,2,3
%                        having the same distribution as the non-normal
%                        system
%
%         instring  =  'symmetric_3'
%
%                        specifies an order 120  symmetric system
%                        with eigenvalues at lambda = 1,2,3
%
%         instring  =  'symmetric_cluster'
%
%                        specifies an order 120  symmetric system
%                        with clustered eigenvalues at lambda = 1,2,3
%
%         instring  =  'symmetric_random'
%
%                        specifies an order 120  symmetric 
%                        positive definite system
%                        with eigenvalues generated by rand 
%
%
%
%   D.C. Sorensen
%
%   Fall 2002
%   User interface added Oct. 2003
%   Name Change to Gmres_Fom.m  on 27 Oct 2008
%   FOM residual and graphics update added 27 Oct 2008
%

    MaxIter = 4; Itr = 1;

        switch lower(instring)  
          case 'random'
            MaxIter = 5;
            n = 50;
            A = randn(n); b = randn(n,1);
            s = eig(A);
            outstring = 'Random order 50 system';

          case {'normal', 'non-normal'}
            d = [ones(20,1) ; 2*ones(60,1) ; 3*ones(40,1)];
            n = length(d);
            b = randn(n,1);
            [Q,R] = qr(randn(n)); 
%
%            change value of tau here to vary non-normality
%             larger tau gives higher non-normality
%
%           tau =  0;
%           tau = .05;
            tau = .5;

            R = R*tau;
            for j = 1:n, R(j,j) = d(j); end

            A = Q*R*Q';
            s = eig(A);

            switch lower(instring)
               case 'normal'
                  A = Q*diag(s)*Q';
                  outstring = 'Normal system with 3 clusters of eigenvalues';
               otherwise
                 outstring = 'Non-normal system with 3 clusters of eigenvalues';
               end

         case 'symmetric_3'

            MaxIter = 2;
            d = [ones(20,1) ; 2*ones(60,1) ; 3*ones(40,1)];
            n = length(d);
            b = randn(n,1);
            [Q,R] = qr(randn(n)); 
            A = Q*diag(d)*Q'; s = d;
            outstring = 'Symmetric system with 3 eigenvalues';

         case 'symmetric_cluster'

            d = [ones(20,1) ; 2*ones(60,1) ; 3*ones(40,1)];
            n = length(d);
            b = randn(n,1);
            [Q,R] = qr(randn(n)); 
            d = d + .05*randn(n,1);
            A = Q*diag(d)*Q'; s = d;
            outstring = 'Symmetric system with 3 clusters of eigenvalues';

         case 'symmetric_random'

            d = rand(120,1);
            n = length(d);
            b = randn(n,1);
            [Q,R] = qr(randn(n)); 
            A = Q*diag(d)*Q'; s = d;
            outstring = 'Symmetric positive definite system';

        otherwise
            outstring = 'Unknown case'; return;
        end
%
%  Begin computing and displaying results
%
   r = [];
   rf = [];
   Rho = [];
   Rhof = [];
   h1 = 1;
   h2 = 2;
   h3 = 3;
   t1 = [1:12];
   t2 = (1.0e-8)*ones(12,1);
   t2(1) = 10; t2(12) = 10;
   for k = 1:1:n,
       
      [x,rho,xf,rhof,ritz,hritz] = Gmres_Fom(A,b,k);

      r = [r, norm(b - A*x)];
      rf = [rf, norm(b - A*xf)];
      Rho = [Rho,rho];
      Rhof = [Rhof,rhof];
    
      figure(h1)
      clf
      semilogy(t1,t2,'k.')
      hold on
      [hand1] = semilogy(r,'linewidth',1.5)
      semilogy(r,'*')
      [hand2] = semilogy(rf,'r-','linewidth',1.5)
      semilogy(rf,'r*')
      legend([hand1, hand2], 'GMRES Resid','FOM Resid') 
      title('GMRES residual vs FOM residual','fontsize',20)
      xlabel(' Iterate Number k ','fontsize',14)
      ylabel(' Residual  norm( b - A x_k ) at Iterate k ','fontsize',14)
      hold off

      figure(h2)
      clf
      plot(real(s),imag(s),'k+')
      hold
      t = hritz;
      plot(real(t),imag(t),'o','linewidth',1.5)
      t = ritz;
      plot(real(t),imag(t),'r*')
      legend('eigs','h-Ritz','Ritz')
      title('GMRES roots (h-Ritz) vs FOM roots (Ritz)','fontsize',20 )
      hold
      pause(.7)
      if k == 1, 
            disp(' Arrange Graphs and Strike a Key')
            pause 
      end

      if (mod(k,10) == 0 ),
          disp('Strike a Key') 
          pause 
          logP = 0;
          t = hritz;
%
%         to see FOM poly put t = ritz
%
%         t = ritz;
%
          figure(h2)
          x = axis;
          hx = (x(2)-x(1))/101;
          hy = (x(4)-x(3))/101;
          xx = x(1)+hx:hx:x(2)-hx;
          yy = x(3)+hx:hy:x(4)-hy;
          yy = x(3)+hy:hy:x(4)-hy;
          m = length(yy);
          n = length(xx);
          e = ones(n,1);
          X = e*xx + i*yy'*e';
          P = zeros(m,n);
          p0 = sum(log(abs(t)));
          for jj = 1:length(t),
              P = P + log(abs(X - t(jj)));
          end     
          P = P - p0*ones(m,n);
%
%         P should have the values of log(abs(p(z)/p(0)))
%           for z in region of axis where p is GMRES poly
%
          hold
             [cc,hh] = contour(xx,yy,P); 
%            clabel(cc,hh);
          hold
          figure(h3)
          clf
          surfc(P);
          shading interp
          colorbar
          title('Abs. Value GMRES Poly over Complex Plane','fontsize',18)
          disp('Strike a Key')
          Itr = Itr + 1;
      pause 
      if (Itr > MaxIter), return, end;
      end
   
      Norm_Resid_Diff = norm(Rho-r)
      Norm_Resid_Diff_Fom = norm(Rhof-rf) 
   end