c
c
c*****************************************************************************
c*                                                                           *
c*  Copyright Compaq Computer Corporation 1999.         All rights reserved. *
c*                                                                           *
c*****************************************************************************
c
c
c
      PROGRAM EXAMPLE_ITSOL_5
C
C ***** TO DEMONSTRATE THE USE OF THE SPARSE ITERATIVE SOLVER.
C
C ***** THIS PROGRAM ILLUSTRATES THE FOLLOWING:
C
C       (1) USE OF THE DRIVER ROUTINE TO SOLVE THE TEST PROBLEM VIA
C           PRECONDITIONED CONJUGATE GRADIENT METHOD, WITH SPD SPLIT 
C           INCOMPLETE CHOLESKY PRECONDITIONING. THE MATRIX IS STORED USING
C           THE UNSYMMETRIC DIAGONAL FORMAT
C       (2) USE OF ROUTINE DMATVEC_UDIA
C       (3) USE OF ROUTINE DCREATE_ILU_UDIA TO CREATE THE PRECONDITIONER 
C       (4) USE OF THE ROUTINE DITSOL_DEFAULTS
C
C 
C ***** SAMPLE OUTPUT FOLLOWS
C
C  SOLVING EXAMPLE PROBLEM USING DRIVER ROUTINE WITH  SPLIT PRECONDITIONED CG
C  ILU PRECONDITIONING USED 
C  MATRIX STORED IN UNSYMMETRIC DIAGONAL FORMAT
C
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      iteration =          0    stopping test = 0.69282032E+01
C      iteration =          1    stopping test = 0.19194193E+01
C      iteration =          2    stopping test = 0.12100937E+01
C      iteration =          3    stopping test = 0.52439623E+00
C      iteration =          4    stopping test = 0.84029860E-01
C      iteration =          5    stopping test = 0.20539881E-01
C      iteration =          6    stopping test = 0.34309306E-02
C      iteration =          7    stopping test = 0.47063334E-03
C      iteration =          8    stopping test = 0.16605002E-03
C      iteration =          9    stopping test = 0.45072557E-04
C      iteration =         10    stopping test = 0.72087304E-05
C      iteration =         11    stopping test = 0.88047573E-06
C
C   solution obtained after          11  iterations
C   normal exit from solver
C   final value of stopping test = 0.88047573E-06
C
C     TRUE SOLUTION     SOLUTION FROM SOLVER     ABS. DIFFERENCE
C
C      0.10000000E+01      0.10000000E+01      0.23690334E-07
C      0.10000000E+01      0.99999998E+00      0.24240992E-07
C      0.10000000E+01      0.10000000E+01      0.18777453E-07
C      0.10000000E+01      0.99999999E+00      0.12495111E-07
C                                .
C                                .
C                        (edited for brevity)
C                                .
C                                .
C      0.10000000E+01      0.99999999E+00      0.14116978E-07
C      0.10000000E+01      0.10000000E+01      0.29831236E-07
C      0.10000000E+01      0.99999998E+00      0.20860417E-07
C      0.10000000E+01      0.10000000E+01      0.25380622E-07
C
C  MAX ERROR IN SOLUTION =  0.56260082E-07
C
C
C
C
      IMPLICIT NONE

      INTEGER NMAX
      PARAMETER (NMAX = 100)
      INTEGER NDIM
      PARAMETER (NDIM = 100)
      INTEGER NWK
      PARAMETER (NWK = 5*NMAX)
C
      REAL*8 X(NMAX), XO(NMAX), RHS(NMAX)
      REAL*8 RPARAM(50), RWORK(NWK), DUM
      REAL*8 A_UDIA(NDIM,5), P_ILU(NDIM,5), TEMP
      INTEGER IP_ILU(5), INDEX_UDIA(5)
C
      INTEGER IPARAM(50), IDUM
      INTEGER I, J, NX, NY, NXNY, IERROR, IOUNIT, NZEROS
C
      EXTERNAL DMATVEC_DRIVER, DPCONDL_DRIVER, DPCONDR_DRIVER
      EXTERNAL USER_PRINT_ROUTINE
C
C ***** SET UP PROBLEM SIZE
C
      NX = 10
      NY = 10
      NXNY = NX*NY
      NZEROS = 5
C
C ***** Initialize P_ILU
C 
      DO I = 1, NDIM 
         DO J = 1, 5
          P_ILU(I,J) = 0.0D0
         ENDDO
      ENDDO
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 IOUNIT TO 7
C              CHANGE IPCOND TO 4 (FOR SPD SPLIT PRECONDITIONING)
C              CHANGE IOLEVEL TO 3
C
      IPARAM(3) = 0
      IPARAM(4) = NWK
      IPARAM(35) = 3
      IPARAM(36) = 3
      IPARAM(37) = 1
C
      IPARAM(5) = 7
      IPARAM(6) = 3
      IPARAM(7) = 4
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)

      WRITE (IOUNIT,101)
101   FORMAT (/,2X,'SOLVING EXAMPLE PROBLEM USING DRIVER ROUTINE WITH
     $  SPLIT PRECONDITIONED CG',/,2X,'ILU PRECONDITIONING USED ',/
     $ 2X,'MATRIX STORED IN UNSYMMETRIC DIAGONAL FORMAT',/)
C 
C ***** GENERATE THE MATRIX
C
      CALL GENMAT(NX, NY, NXNY, A_UDIA, INDEX_UDIA, 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 USING DMATVEC_UDIA
C
      CALL DMATVEC_UDIA (0, A_UDIA, INDEX_UDIA, NDIM, NZEROS, 
     $                   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 ILU PRECONDITIONER
C
      CALL DCREATE_ILU_UDIA (A_UDIA, INDEX_UDIA, NDIM, NZEROS,
     $                      P_ILU, IP_ILU, NXNY)
C
C ***** CALL THE SOLVER
C
       CALL DITSOL_DRIVER ( DMATVEC_DRIVER, DPCONDL_DRIVER,
     $                      DPCONDR_DRIVER, DUM, 
     $                      A_UDIA, INDEX_UDIA, X, RHS, NXNY, 
     $                      P_ILU, IP_ILU, P_ILU, IP_ILU, 
     $                      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 GENMAT (NX, NY, NXNY, A, IA, NDIM, NZEROS)
C
C ***** GENERATE THE MATRIX FOR THE EXAMPLE IN THE UNSYMMETRIC DIAGONAL
C       FORM
C
      IMPLICIT NONE

      INTEGER NDIM, J, NXNY, NX, NY, I, K, NZEROS
      REAL*8 A(NDIM,*)
      INTEGER IA(*)
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-1
          DO I = 1, NX
              K = (J-1)*NX+I
              A(K,2) = -1.0D0
          END DO
      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 = 1, NY
          DO I = 1, NX-1
              K = (J-1)*NX+I
              A(K,4) = -1.0D0
          END DO
      END DO
C
      DO J = 2, NY
          DO I = 1, NX
              K = (J-1)*NX+I
              A(K,5) = -1.0D0
          END DO
      END DO
C
      IA(1) = 0
      IA(2) = NX
      IA(3) = -1
      IA(4) = 1
      IA(5) = -NX
C
      RETURN
      END


