*
      SUBROUTINE SSPDCHK( X, LDX, B, X0, WORK, LDW, PFORM, MATVEC, 
     $                    MATVECTRANS, PSOLVE, PSOLVETRANS, PSOLVEQ,
     $                    PSOLVETRANSQ, BACKSOLVE, TOL, SCALEDTOL,
     $                    LTEST, SPDRES, NUMTESTS, NUMSUSP, CRITERR )
*
*     .. Scalar Arguments ..
      INTEGER            LDW, LDX, NUMTESTS, NUMSUSP, CRITERR
      REAL               TOL, SCALEDTOL
      LOGICAL            SPDRES
*     ..
*     .. Array Arguments ..
      REAL               X( LDX,* ), B( * ), X0( * ), WORK( * )
      LOGICAL            LTEST( * )
      CHARACTER*5        PFORM( * )
*     ..
*     .. External Functions ..
      EXTERNAL           MATVEC, MATVECTRANS, PSOLVE, PSOLVETRANS,
     $                   PSOLVEQ, PSOLVETRANSQ, BACKSOLVE
*     ..
*     .. Common Blocks ..
*     MAXDIM2 = MAXDIM*MAXDIM.
*
      INTEGER             MAXDIM, MAXDIM2
      PARAMETER         ( MAXDIM = 200, MAXDIM2 = 40000 )
*
      INTEGER             N, LDA
      REAL                A, M
      CHARACTER           CURPFORM*5
*
      COMMON            / SYSTEM / A( MAXDIM2 ), M( MAXDIM )
      COMMON            / MATDIM / N, LDA
      COMMON            / FORMS  / CURPFORM
*
*  Purpose
*  =======
*
*  Subroutine to test the performance of the template kernels
*  on symmetric positivie definite matrices.
*
*  Generates, solves, and checks accuracy of linear systems.
*
*  Algorithms tested:
*
*     1. CG
*     2. Cheby
*     3. SOR
*     4. BiCG
*     5. CGS
*     6. BiCGSTAB
*     7. GMRES
*     8. QMR
*     9. Jacobi ( for diagonally dominant Poisson matrices only )
*
*  Various systems are generated. Each method attempts to solve the 
*  system to the input TOL in MAXIT iterations. Each method iterates 
*  using various preconditioners. 
*
*  The result is compared with the normalized scaled residual 
*
*     || b - A*x || / ( ||A||||x||*N*TOL ). 
*
*  In order to do this, the solution vectors are stored in  matrix 
*  X( LDX,* ). Column j contains the solution vector for algorithm j,
*  j as defined above.
*
*  =================================================================
*
*     .. Parameters ..
*
      REAL               ONE
      PARAMETER        ( ONE = 1.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J, K, NX, NY, NZ, MATFORM, NSYS, NPFORMS,
     $                   MAXIT, RESTRT, IVAL, IINDX, IPOINTR
      REAL               ANORM, MATNORM,
     $                   ONEFUN, ZEROFUN, NEGONEFUN, HENKFUN, HENKDFUN,
     $                   THOUSFUN, THOUSXFUN, TEN5X2FUN, NEGTHOUSXFUN
      LOGICAL            TSTCG, TSTCHEBY, TSTSOR, TSTBICG, TSTCGS,
     $                   TSTBICGS, TSTGMRES, TSTQMR, TSTJACOBI,
     $                   FLAG, LSAMEN
      CHARACTER          AFORM*4, RHS*4, INITIALGUESS*4
*     ..
*     .. Local Arrays ..
      INTEGER            INFO( 9 ), ITER( 9 )
      REAL               RESID( 9 )
*     ..
*     .. External Functions ..
*     PDE Coefficient functions.
      EXTERNAL           ONEFUN, ZEROFUN, NEGONEFUN, HENKFUN, HENKDFUN,
     $                   THOUSFUN, THOUSXFUN, TEN5X2FUN, NEGTHOUSXFUN
*     ..
*     .. Executable Statements ..
*
      NPFORMS  = 2
      NUMTESTS = 0
      NUMSUSP  = 0
      CRITERR  = 0
*
      READ(9,*) NSYS
*
*     Check for quick return.
*
      IF ( NSYS.LT.0 ) THEN
         WRITE(*,*) 'ERROR IN NONSYMMETRIC TESTER: NUMBER OF SYSTEMS TO
     $BE GENERATED IS LESS THAN 0'
         RETURN
      ENDIF
*
      FLAG = .FALSE.
      DO 5 I = 1, 9
         IF ( LTEST( I ) ) FLAG = .TRUE.
    5 CONTINUE
      IF ( .NOT.FLAG ) RETURN
*
      TSTCG     = LTEST( 1 )
      TSTCHEBY  = LTEST( 2 )
      TSTSOR    = LTEST( 3 )
      TSTBICG   = LTEST( 4 )
      TSTCGS    = LTEST( 5 )
      TSTBICGS  = LTEST( 6 )
      TSTGMRES  = LTEST( 7 )
      TSTQMR    = LTEST( 8 )
      TSTJACOBI = LTEST( 9 )
*
   10 CONTINUE
*
      DO 60 MATFORM = 1, NSYS
*
         READ(9,*) AFORM, NX, NY, NZ, RHS, INITIALGUESS
*
*        The following two matrices are generated using a 5- or 7-point 
*        stencil using centered differences on a 1d, 2d, or 3d grid, 
*        with Dirichlet boundary conditions.
*
*        The last 7 arguments to this routine are the coefficient 
*        functions for the PDE:
*
*           delx( a delx u ) + dely ( b dely u) + delz ( c delz u ) +
*           delx ( d u ) + dely (e u) + delz( f u ) + g u
*
         IF ( LSAMEN( 4, AFORM,'F2SH' ) ) THEN
*
*           -u_xx - u_yy
*
            N = NX*NY*NZ
            IPOINTR = 1
            IINDX   = IPOINTR + N+1
            IVAL    = IINDX   + 5*N
            CALL GEN57PT( NX, NY, NZ, WORK(IVAL), WORK(IINDX),
     $                    WORK(IPOINTR), NEGONEFUN, NEGONEFUN,
     $                    ZEROFUN, ZEROFUN, ZEROFUN, ZEROFUN, ZEROFUN )
            CALL COMP2DENSE( WORK(IVAL), WORK(IPOINTR), 
     $                       WORK(IINDX), N, A, LDA,'ROW', INFO(1) )
            IF ( INFO(1).NE.0 ) THEN
               WRITE(*,81) N, INFO(1)
               WRITE(10,81) N, INFO(1)
               GO TO 60
            ENDIF
*
         ELSE IF ( LSAMEN( 4, AFORM,'F3SH' ) ) THEN
*
*           -u_xx - u_yy - u_zz
*
            N = NX*NY*NZ
            IPOINTR = 1
            IINDX   = IPOINTR + N+1
            IVAL    = IINDX   + 10*N
            CALL GEN57PT( NX, NY, NZ, WORK(IVAL), WORK(IINDX),
     $                    WORK(IPOINTR), NEGONEFUN, NEGONEFUN,
     $                    NEGONEFUN, ZEROFUN, ZEROFUN, ZEROFUN,
     $                    ZEROFUN )
            CALL COMP2DENSE( WORK(IVAL), WORK(IPOINTR),
     $                       WORK(IINDX), N, A, LDA,'ROW', INFO(1) )
            IF ( INFO(1).NE.0 ) THEN
               WRITE(*,82) N, INFO(1)
               WRITE(10,82) N, INFO(1)
               GO TO 60
            ENDIF
         ELSE IF ( LSAMEN( 4, AFORM,'WATHEN' ) ) THEN
*
*           Form Wathen matrix. 
*           ( Matrix order will be 3*NX*NY + 2*NX + 2*NY + 1 )
*
            K = 0
            CALL WATHEN( NX, NY, K, N, A, LDA, WORK, LDW, INFO(1) )
            IF ( INFO(1).NE.0 ) THEN
               WRITE(*,83)  3*NX*NY + 2*NX + 2*NY + 1, INFO(1)
               WRITE(10,83) 3*NX*NY + 2*NX + 2*NY + 1, INFO(1)
               GO TO 60
            ENDIF
         ELSE
            WRITE(*,*) AFORM, 'IS AN UNKNOWM MATRIX TYPE'
            GO TO 60
         ENDIF
*
         MAXIT = 4*N
         ANORM = MATNORM( N, A, LDA )
*
*        Form RHS and initial guess vectors.
*
         CALL VECGEN( RHS, N, A, LDA, B, INFO(2) )
         CALL VECGEN( INITIALGUESS, N, A, LDA, X0, INFO(3) )
*
         IF ( INFO(2).NE.0 ) THEN
            WRITE(*,92)  RHS
            WRITE(10,92) RHS
            GO TO 60
         ELSE IF ( INFO(3).NE.0 ) THEN
            WRITE(*,93)  INITIALGUESS
            WRITE(10,93) INITIALGUESS
            GO TO 60
         ENDIF
*
   20    CONTINUE
*
*        Solve system using the various algorithms, using no 
*        preconditioning, then diagonal preconditioning.
*
         DO 50 I = 1, NPFORMS
*
            DO 30 J = 1, 9
               INFO( J ) = 0
               ITER( J )  = MAXIT
               RESID( J ) = TOL
               CALL SCOPY( N, X0, 1, X(1,J), 1 )
   30       CONTINUE
*
            CALL PRECON( N, A, LDA, PFORM( I ), M, INFO(1) )
            IF ( INFO(1).NE.0 ) THEN
               WRITE(*,94) PFORM(I)
               WRITE(10,94) PFORM(I)
               GO TO 50
            ENDIF
            CURPFORM = PFORM( I )
*
            IF ( TSTCG )
     $         CALL CG( N, B, X(1,1), WORK, LDW, ITER(1), RESID(1),
     $                  MATVEC, PSOLVE, INFO(1) )
*
            IF ( TSTCHEBY )
     $         CALL CHEBY( N, B, X(1,2), WORK, LDW, ITER(2), RESID(2),
     $                     MATVEC, PSOLVE, INFO(2) )
*
            IF  ( TSTSOR ) THEN
*
*              Set OMEGA
*
               IF ( MATFORM.EQ.1 ) THEN
*
*                 Guass-Seidel
*
                  WORK(1) = ONE
               ELSE
                  WORK(1) = 1.2E+0
               ENDIF
               CALL SOR( N, B, X(1,3), WORK, LDW, ITER(3),
     $                   RESID(3), MATVEC, BACKSOLVE, INFO(3) )
            ENDIF
*
            IF ( TSTBICG )
     $         CALL BICG( N, B, X(1,4), WORK, LDW, ITER(4), RESID(4),
     $                    MATVEC, MATVECTRANS, PSOLVE, PSOLVETRANS,
     $                    INFO(4) )
*
            IF ( TSTCGS )
     $         CALL CGS( N, B, X(1,5), WORK, LDW, ITER(5), RESID(5),
     $                   MATVEC, PSOLVE, INFO(5) )
*
            IF ( TSTBICGS )
     $         CALL BICGSTAB( N, B, X(1,6), WORK, LDW, ITER(6),
     $                        RESID(6), MATVEC, PSOLVE, INFO(6) )
*
            IF ( TSTGMRES ) THEN
*
*              For the symmetric case, restarts = N.
*
               RESTRT = N
               CALL GMRES( N, B, X(1,7), RESTRT, WORK, LDW,
     $                     WORK((6+RESTRT)*LDW+1), LDW, ITER(7),
     $                     RESID(7), MATVEC, PSOLVE, INFO(7) )
            ENDIF
*
            IF ( TSTQMR )
     $          CALL QMR( N, B, X(1,8), WORK, LDW, ITER(8), RESID(8),
     $                    MATVEC, MATVECTRANS, PSOLVEQ, PSOLVETRANSQ, 
     $                    INFO(8))
*
            IF ( TSTJACOBI ) THEN
               IF ( .NOT.LSAMEN( 3, AFORM,'WATH' ).AND.( I.EQ.1 ) ) THEN
*
*                 Since preconditioning does not apply to Jacobi, it is
*                 only called once per test matrix. (Wathen matrix is 
*                 not diagonally dominant.)
*
                  CALL JACOBI( N, B, X(1,9), WORK, LDW, ITER(9), 
     $                         RESID(9), MATVEC, INFO(9) )
                  IF ( INFO( 9 ).NE.0 ) NUMSUSP = NUMSUSP + 1
               ELSE
*
*                 Flag not to check accuracy.
*
                  INFO( 9 ) = 100
               ENDIF
            ENDIF
*
*           Check for convergence.
*
            DO 40 J = 1, 8
               IF ( INFO( J ).NE.0 ) NUMSUSP = NUMSUSP + 1 
   40       CONTINUE
*
*           Write results to file.
*
            CALL RESULT( N, A, LDA, X, LDX, B, WORK(1), 'SPD', PFORM(I),
     $                   ITER, RESID, TOL, INFO, AFORM, ANORM,
     $                   LTEST, SCALEDTOL, SPDRES, CRITERR )
            NUMTESTS = NUMTESTS + 1
*
   50    CONTINUE
*
   60 CONTINUE
*
   81 FORMAT('WARNING: COULD NOT FORM 2-D POISSON (N=', I4,'),INFO=',I2)
   82 FORMAT('WARNING: COULD NOT FORM 3-D POISSON (N=', I4,'),INFO=',I2)
   83 FORMAT('WARNING: COULD NOT FORM ORDER', I4,' WATHEN MATRIX')
*
   92 FORMAT('ERROR: RHS FOR TEST MATRIX', A4, ' NOT FORMED')
   93 FORMAT('ERROR: INITIAL GUESS FOR TEST MATRIX', A4, ' NOT FORMED')
   94 FORMAT('ERROR: PRECONDITIONER', A5, ' NOT FORMED')
*
  200 CONTINUE
*
      RETURN
*
*     -- End of SSPDCHK
*
      END
