c
c
c*****************************************************************************
c*                                                                           *
c*  Copyright Compaq Computer Corporation 1999.         All rights reserved. *
c*                                                                           *
c*****************************************************************************
c
c
c
      PROGRAM EXAMPLE_ITSOL_2
C
C ***** TO DEMONSTRATE THE USE OF THE SPARSE ITERATIVE SOLVER.
C
C ***** THIS PROGRAM ILLUSTRATES THE FOLLOWING:
C
C       (1) USE OF THE SOLVER TO SOLVE THE TEST PROBLEM VIA
C           PRECONDITIONED CONJUGATE GRADIENT METHOD, WITH POLYNOMIAL
C           PRECONDITIONING (OF DEGREE=3). THE MATRIX IS STORED USING THE
C           SYMMETRIC DIAGONAL FORMAT
C       (2) USE OF ROUTINE MATVEC (DMATVEC_SDIA)
C       (3) USE OF ROUTINE PCONDL TO CALL ROUTINES FOR APPLYING THE 
C           PRECONDITIONER (DAPPLY_POLY_SDIA)
C       (5) USE OF THE ROUTINE DITSOL_DEFAULTS
C       (6) INFORMATION PRINTED OUT FOR IOLEVEL = 2
C
C
C ***** SAMPLE OUTPUT FOLLOWS
C
C SOLVING EXAMPLE PROBLEM WITH SPD SPLIT PRECONDITIONED CG
C POLYNOMIAL PRECONDITIONING USED 
C MATRIX STORED IN SYMMETRIC DIAGONAL FORMAT
C  method used : cg with spd split preconditioning
C  order of system =    100
C  stopping criterion used :     1
C  maximum iterations allowed:   100
C  tolerance for convergence :   0.10000000E-05
C  solution obtained after          10  iterations
C  normal exit from solver
C  final value of stopping test = 0.20293204E-06
C    TRUE SOLUTION     SOLUTION FROM SOLVER     ABS. DIFFERENCE
C     0.10000000E+01      0.10000000E+01      0.14674600E-08
C     0.10000000E+01      0.10000000E+01      0.28075924E-08
C     0.10000000E+01      0.10000000E+01      0.33202063E-08
C     0.10000000E+01      0.10000000E+01      0.41068805E-08
C                               .
C                               .
C                               .
C                       (edited for brevity)
C                               .
C                               .
C                               .
C     0.10000000E+01      0.10000000E+01      0.41068805E-08
C     0.10000000E+01      0.10000000E+01      0.33202063E-08
C     0.10000000E+01      0.10000000E+01      0.28075924E-08
C     0.10000000E+01      0.10000000E+01      0.14674602E-08
C MAX ERROR IN SOLUTION =  0.11405409E-07
C
C
C
C
      IMPLICIT NONE

      INTEGER NMAX
      PARAMETER (NMAX = 100)
      INTEGER NDIM
      PARAMETER (NDIM = 100)
C
      REAL*8 X(NMAX), XO(NMAX), RHS(NMAX)
      REAL*8 RPARAM(50), RWORK(4*NMAX), DUM
      REAL*8 A_SDIA(NDIM,3), P_POLY(3*NMAX), TEMP
      INTEGER INDEX_SDIA(3)
C
      INTEGER IPARAM(50), IDUM
      INTEGER I, NX, NY, NXNY, NZEROS, IOUNIT, IERROR
C
      EXTERNAL MATVEC, PCONDL
      EXTERNAL USER_PRINT_ROUTINE
C
C ***** SET UP PROBLEM SIZE
C
      NX = 10
      NY = 10
      NXNY = NX*NY
      NZEROS = 3
C
C ***** SET THE PARAMETERS (INTEGER AND REAL)
C
      CALL DITSOL_DEFAULTS (IPARAM, RPARAM)
C
C ***** CHANGE ANY VALUES THAT ARE DIFFERENT FROM THE DEFAULT
C       ASSIGN VALUES TO PARAMETERS NOT SET BY DEFAULTS
C              CHANGE IPCOND TO 4 (FOR SPD SPLIT PRECONDITIONING)
C              CHANGE IOLEVEL TO 2
C
      IPARAM(3) = 0
      IPARAM(4) = 4*NXNY
C
      IPARAM(5) = 7
      IPARAM(6) = 2
      IPARAM(7) = 4
C
C ***** SET THE DEGREE OF THE POLYNOMIAL
C
      IPARAM(33) = 2
C
C ***** SETUP OUTPUT FILE
C
      IOUNIT = IPARAM(5) 
      OPEN (UNIT=IOUNIT,FILE='OUTPUT.DATA',STATUS='UNKNOWN')
      REWIND IOUNIT

C
C *****  SETUP PRINT ROUTINE
C
      CALL CXML_ITSOL_SET_PRINT_ROUTINE(USER_PRINT_ROUTINE,
     $         IOUNIT, IPARAM)
C
      WRITE (IOUNIT,101)
101   FORMAT (/,2X,'SOLVING EXAMPLE PROBLEM WITH SPD SPLIT
     $ PRECONDITIONED CG',/,2X,'POLYNOMIAL PRECONDITIONING USED ',/
     $ 2X,'MATRIX STORED IN SYMMETRIC DIAGONAL FORMAT',/)
C 
C ***** GENERATE THE MATRIX
C
      CALL GENMAT(NX, NY, NXNY, A_SDIA, INDEX_SDIA, NDIM, NZEROS)
      IPARAM(31) = NZEROS
      IPARAM(32) = NDIM
C
C ***** GENERATE XO, THE TRUE SOLUTION
C
      DO I = 1, NXNY
          XO(I) = 1.0D0
      END DO
C
C ***** OBTAIN THE RIGHT HAND SIDE
C
      CALL MATVEC (0, IPARAM, RPARAM, A_SDIA, INDEX_SDIA, DUM, XO, 
     $             RHS, NXNY)
C
C ***** OBTAIN INITIAL GUESS (ALL ZEROS)
C
      DO I = 1, NXNY
         X(I) = 0.0D0
      END DO
C
C ***** GENERATE THE POLYNOMIAL PRECONDITIONER
C
      CALL DCREATE_POLY_SDIA (A_SDIA, INDEX_SDIA, NDIM, NZEROS,
     $                      P_POLY, NXNY)
C
C ***** CALL THE SOLVER
C
       CALL DITSOL_PCG ( MATVEC, PCONDL, DUM, DUM, 
     $                   A_SDIA, INDEX_SDIA, X, RHS, NXNY, 
     $                   P_POLY, IDUM, DUM, IDUM, 
     $                   IPARAM, RPARAM, IDUM, RWORK, IERROR)
C
C ***** PRINT OUT THE SOLUTION
C
      WRITE (IOUNIT,102)
102   FORMAT (/,5X,'TRUE SOLUTION',5X,'SOLUTION FROM SOLVER',
     $5X,'ABS. DIFFERENCE'/)
      WRITE (IOUNIT,103) (XO(I),X(I),ABS(XO(I)-X(I)),I=1,NXNY)
103   FORMAT (/,3(5X,E15.8))
C
C ***** FIND MAX ERROR IN SOLUTION
C
      TEMP = ABS(XO(1) - X(1))
      DO I =2,NXNY
         TEMP = MAX( TEMP, ABS(XO(I) - X(I)) )
      END DO
      WRITE(IOUNIT,104) TEMP
104   FORMAT(/,2X,'MAX ERROR IN SOLUTION = ', E15.8,/)
C
      STOP
      END
C
C
C
      SUBROUTINE USER_PRINT_ROUTINE( CONTEXT, BUF, SIZE)

      INTEGER CONTEXT
      INTEGER SIZE
      INTEGER BUF(*)

      INCLUDE 'cxml_itsol_print.for'

      CHARACTER*(CXML_ITSOL_MAX_MESSAGE_LENGTH) NEW_BUF

      CALL CXML_FORMAT_STRING(BUF, SIZE, NEW_BUF)
      WRITE(CONTEXT, '(A)') NEW_BUF(1:SIZE)

      END
C
C
C
      SUBROUTINE MATVEC(JOB, IPARAM, RPARAM, A, IA, W, X, Y, N)
C
C  MULTIPLY N VECTOR X BY MATRIX TO OBTAIN Y
C

      IMPLICIT NONE
C
      REAL*8 X(*), Y(*), A(*), RPARAM(*), DUM, W
      INTEGER IA(*), JOB, NDIM,NZEROS, IPARAM(*), N 
C
      NZEROS = IPARAM(31)
      NDIM = IPARAM(32)
C
      CALL DMATVEC_SDIA (JOB, A, IA, NDIM, NZEROS, DUM, X, Y, N, IPARAM)
C
      RETURN
      END
C
C
C
      SUBROUTINE PCONDL(JOB, IPARAM, RPARAM, QL, IQL, A, IA, 
     $                  W, X, Y, N)
C
C ***** CALL THE PRECONDITIONER (LEFT)
C
      IMPLICIT NONE

      REAL*8 X(*), Y(*), RPARAM(*), QL(*), A(*), W
      INTEGER JOB, NDIM,NDEG, IPARAM(*), IQL(*), IA(*)
      INTEGER NZEROS, N
C
C ***** POLYNOMIAL PRECONDITIONING
C
      NZEROS = IPARAM(31)
      NDIM = IPARAM(32)
      NDEG = IPARAM(33)
C
      CALL DAPPLY_POLY_SDIA (0, QL, A, IA, NDIM, NZEROS, X, Y, 
     $                      NDEG, N, IPARAM)
C
      RETURN
      END
C
C
C
      SUBROUTINE GENMAT (NX, NY, NXNY, A, IA, NDIM, NZEROS)
C
C ***** GENERATE THE MATRIX FOR THE EXAMPLE IN THE SYMMETRIC DIAGONAL
C       FORM, STORING THE LOWER TRIANGLULAR PART.
C
      IMPLICIT NONE
     
      INTEGER NDIM, NXNY, NZEROS, NY, NX  
      REAL*8 A(NDIM,*)
      INTEGER IA(*), I, J, K  
C
      DO J = 2, NZEROS
          DO I = 1, NXNY
              A(I,J) = 0.0D0
          END DO
      END DO
C
      DO I = 1, NXNY
          A(I,1) = 4.0D0
      END DO
C
      DO J = 1, NY
          DO I = 2, NX
              K = (J-1)*NX+I
              A(K,3) = -1.0D0
          END DO
      END DO
C
      DO J = 2, NY
          DO I = 1, NX
              K = (J-1)*NX+I
              A(K,2) = -1.0D0
          END DO
      END DO
C
      IA(1) = 0
      IA(3) = -1
      IA(2) = -NX
C
      RETURN
      END
