*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*     File symmlqtest.f
*
*     These routines are for testing  SYMMLQ.
*
*     04 Sep 1991: logical checkA added to parameter list.
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine Aprod ( n, x, y )

      implicit           double precision (a-h,o-z)
      double precision   x(n), y(n)

*     ------------------------------------------------------------------
*     Aprod  computes  y = A*x  for some matrix  A.
*     This is a simple example for testing  SYMMLQ.
*     ------------------------------------------------------------------

      do 10 i = 1, n
         d    = i * 1.1
         d    = d / n
         y(i) = d * x(i)
   10 continue

*     end of Aprod
      end

      subroutine Msolve( n, x, y )

      implicit           double precision (a-h,o-z)
      double precision   x(n), y(n)

*     ------------------------------------------------------------------
*     Msolve  solves  M*y = x  for some symmetric positive-definite
*     matrix  M.
*     This is a simple example for testing  SYMMLQ.
*     shiftm will be the same as shift in SYMMLQ.
*
*     If pertbn = 0, the preconditioner will be exact, so
*     SYMMLQ should require either one or two iterations,
*     depending on whether (A - shift*I) is positive definite or not.
*
*     If pertbn is nonzero, somewhat more iterations will be required.
*     ------------------------------------------------------------------

      intrinsic          abs, mod
      common    /mshift/ shiftm, pertm

      do 10 i = 1, n
         d    = i * 1.1
         d    = d / n
         d    = abs( d - shiftm )
         if (mod(i,10) .eq. 0) d = d + pertm
         y(i) = x(i) / d
   10 continue

*     end of Msolve
      end

      subroutine test  ( n, goodb, precon, shift, pertbn )

      implicit           double precision (a-h,o-z)
      logical            goodb, precon

*     ------------------------------------------------------------------
*     test   solves sets up and solves a system (A - shift * I)x = b,
*     using Aprod to define A and Msolve to define a preconditioner.
*     ------------------------------------------------------------------

      common    /mshift/ shiftm, pertm
                                                            
      intrinsic          abs
      external           aprod,  msolve
      logical            checkA
      double precision   b(100), r1(100), r2(100), v(100), w(100)
      double precision   x(100), y(100),  xtrue(100)

      parameter        ( one = 1.0d+0,  two = 2.0d+0 )

      shiftm = shift
      pertm  = abs( pertbn )
      nout   = 6
      write(nout, 1000) shift, pertbn

*     Set the true solution and the rhs
*     so that  (A - shift*I) * xtrue = b.

      do 100 i = 1, n
         xtrue(i) = n + 1 - i
  100 continue

      call aprod ( n, xtrue, b )
      call daxpy ( n, (- shift), xtrue, 1, b, 1 )

*     Set other parameters and solve.

      checkA = .true.
      itnlim = n * 2
      rtol   = 1.0E-12

      call symmlq( n, b, r1, r2, v, w, x, y,
     $             aprod, msolve, checkA, goodb, precon, shift,
     $             nout , itnlim, rtol,
     $             istop, itn, anorm, acond, rnorm, ynorm )

*     Compute the final residual,  r1 = b - (A - shift*I)*x.

      call Aprod ( n, x, y )
      call daxpy ( n, (- shift), x, 1, y, 1 )

      do 400 i = 1, n
         r1(i) = b(i) - y(i)
  400 continue

      r1norm = dnrm2 ( n, r1, 1 )
      write(nout, 2000) r1norm

*     Print the solution and some clue about whether it is OK.

      write(nout, 2100) (j, x(j), j = 1, n)

      do 500 j = 1, n
         w(j)  = x(j) - xtrue(j)
  500 continue

      enorm    = dnrm2 ( n, w, 1 ) / dnrm2 ( n, xtrue, 1 )
      etol     = 1.0e-5
      if (enorm .le. etol) then
          write(nout, 3000) enorm
      else
          write(nout, 3100) enorm
      end if

 1000 format(//' ------------------------------------------------------'
     $       / ' Test of  SYMMLQ.'
     $       / ' ------------------------------------------------------'
     $       / ' shift =', f12.4, 6x, 'pertbn =', f12.4)
 2000 format(/ ' Final residual =', 1p, e8.1)
 2100 format(/ ' Solution  x' / 1p, 4(i6, e14.6))
 3000 format(/ ' SYMMLQ  appears to be successful.',
     $         '  Relative error in x =', 1p, e8.1)
 3100 format(/ ' SYMMLQ  appears to have failed.  ',
     $         '  Relative error in x =', 1p, e8.1)
                  
*     end of test
      end

*     ------------------------------------------------------------------
*     Main Program.
*     ------------------------------------------------------------------
      implicit           double precision (a-h,o-z)
      logical            goodb, normal, precon

      parameter        ( zero = 0.0d+0,  one = 1.0d+0 )
                                                    
      normal = .false.
      precon = .true.
      shift  = 0.25
      pertbn = 0.1

      do 100 ntest = 1, 2
         if (ntest .eq. 1) goodb = .false.
         if (ntest .eq. 2) goodb = .true.

*        Test the unlikely tiny cases that often trip us up.
      
         call test( 1, goodb, normal, zero , zero   )
         call test( 2, goodb, normal, zero , zero   )
         call test( 1, goodb, precon, zero , zero   )
         call test( 2, goodb, precon, zero , zero   )
      
*        Test small positive-definite and indefinite systems with
*        exact preconditioners.  SYMMLQ should take 1 and 2 iterations
*        respectively.
      
         n      = 5
         call test( n, goodb, precon, zero , zero   )
         call test( n, goodb, precon, shift, zero   )
      
*        Test other combinations on larger systems.
*        With no preconditioning, SYMMLQ will take about n iterations.
      
         n      = 50
         call test( n, goodb, normal, zero , zero   )
         call test( n, goodb, normal, shift, zero   )
      
*        PERTBN makes the preconditioners incorrect in n/10 entries.
*        SYMMLQ should take about n/10 iterations.
      
         call test( n, goodb, precon, zero , pertbn )
         call test( n, goodb, precon, shift, pertbn )
  100 continue

*     end of Main Program for testing SYMMLQ
      end
