c
c
c*****************************************************************************
c*                                                                           *
c*  Copyright Compaq Computer Corporation 1999.         All rights reserved. *
c*                                                                           *
c*****************************************************************************
c
c
c
      PROGRAM EXAMPLE_ITSOL_4
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 GMRES METHOD METHOD, WITH SPLIT INCOMPLETE
C           CHOLESKY PRECONDITIONING. THE MATRIX IS STORED USING
C           THE UNSYMMETRIC DIAGONAL FORMAT
C       (2) USE OF ROUTINE MATVEC (DMATVEC_UDIA)
C       (3) USE OF ROUTINES PCONDL AND PCONDR TO CALL ROUTINES FOR APPLYING 
C           THE PRECONDITIONER (DAPPLY_ILU_UDIA_L AND DAPPLY_ILU_UDIA_U)
C       (5) USE OF THE ROUTINE DITSOL_DEFAULTS
C
C 
C ***** SAMPLE OUTPUT FOLLOWS
C
C SOLVING EXAMPLE PROBLEM WITH SPLIT PRECONDITIONED GMRES
C ILU PRECONDITIONING USED 
C MATRIX STORED IN UNSYMMETRIC DIAGONAL FORMAT
C  method used : gmres with split preconditioning
C           3  previous residual vectors are stored
C  order of system =    100
C  stopping criterion used :     3
C  maximum iterations allowed:   100
C  tolerance for convergence :   0.10000000E-05
C     iteration =          0    stopping test = 0.27403910E+01
C     iteration =          1    stopping test = 0.83454209E+00
C     iteration =          2    stopping test = 0.44870733E+00
C     iteration =          3    stopping test = 0.17433646E+00
C     iteration =          4    stopping test = 0.72856695E-01
C     iteration =          5    stopping test = 0.38408003E-01
C     iteration =          6    stopping test = 0.14164437E-01
C     iteration =          7    stopping test = 0.56991076E-02
C     iteration =          8    stopping test = 0.28790502E-02
C     iteration =          9    stopping test = 0.13411353E-02
C     iteration =         10    stopping test = 0.63725219E-03
C     iteration =         11    stopping test = 0.32213502E-03
C     iteration =         12    stopping test = 0.13084728E-03
C     iteration =         13    stopping test = 0.54223092E-04
C     iteration =         14    stopping test = 0.27502204E-04
C     iteration =         15    stopping test = 0.13133042E-04
C     iteration =         16    stopping test = 0.63404914E-05
C     iteration =         17    stopping test = 0.32280646E-05
C     iteration =         18    stopping test = 0.13311315E-05
C     iteration =         19    stopping test = 0.55695871E-06
C  solution obtained after          19  iterations
C  normal exit from solver
C  final value of stopping test = 0.55695871E-06
C    TRUE SOLUTION     SOLUTION FROM SOLVER     ABS. DIFFERENCE
C     0.10000000E+01      0.99999997E+00      0.32325871E-07
C     0.10000000E+01      0.99999987E+00      0.13396025E-06
C     0.10000000E+01      0.99999978E+00      0.21953694E-06
C     0.10000000E+01      0.99999975E+00      0.24533450E-06
C                               .
C                               .
C                               .
C                       (edited for brevity)
C                               .
C                               . 
C                               .
C     0.10000000E+01      0.99999972E+00      0.28341823E-06
C     0.10000000E+01      0.99999974E+00      0.25671993E-06
C     0.10000000E+01      0.99999983E+00      0.16690401E-06
C     0.10000000E+01      0.99999995E+00      0.48797753E-07
C MAX ERROR IN SOLUTION =  0.78315035E-06
C
C
C
C
      IMPLICIT NONE

      INTEGER NMAX
      PARAMETER (NMAX = 100)
      INTEGER NDIM
      PARAMETER (NDIM = 100)
      INTEGER KPREV_MAX 
      PARAMETER (KPREV_MAX = 5)
      INTEGER NWK 
      PARAMETER (NWK = NMAX*(KPREV_MAX+2) + KPREV_MAX*(KPREV_MAX+5)+1 )
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, NX, NY, NXNY, IOUNIT, IERROR, NZEROS 
C
      EXTERNAL MATVEC, PCONDL, PCONDR, USER_PRINT_ROUTINE
C
C ***** SET UP PROBLEM SIZE
C
      NX = 10
      NY = 10
      NXNY = NX*NY
      NZEROS = 5
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 3 (FOR SPLIT PRECONDITIONING)
C              CHANGE IOLEVEL TO 3
C              CHANGE ISTOP TO 3 (ONLY ISTOP=3 AND 4 ALLOWED)
C
      IPARAM(3) = 0
      IPARAM(4) = NWK
C
      IPARAM(5) = 7
      IPARAM(6) = 3
      IPARAM(7) = 3
      IPARAM(8) = 3
C
C ***** ASSIGN VALUE TO KPREV (NUMBER OF PREVIOUS RESIDUALS STORED)
C       NOTE THAT THE SIZE OF RWORK ALLOWS A MAXIMUM VALUE OF KPREV = 5
C
      IPARAM(34) = 3
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 WITH SPLIT
     $ PRECONDITIONED GMRES',/,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
C
      CALL MATVEC (0, IPARAM, RPARAM, A_UDIA, INDEX_UDIA, 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_PGMRES ( MATVEC, PCONDL, PCONDR, 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 MATVEC(JOB, IPARAM, RPARAM, A, IA, DUM, 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
      INTEGER IA(*), JOB, NDIM, N, NZEROS, IPARAM(*)
C
      NZEROS = IPARAM(31)
      NDIM = IPARAM(32)
C
      CALL DMATVEC_UDIA (JOB, A, IA, NDIM, NZEROS, DUM, X, Y, N)
C
      RETURN
      END
C
C
C
      SUBROUTINE PCONDL(JOB, IPARAM, RPARAM, QL, IQL, A, IA, 
     $                  W, X, Y, N)
C
C ***** CALL THE LEFT PRECONDITIONER 
C
      IMPLICIT NONE

      INTEGER NMAX
      PARAMETER (NMAX = 100)

      REAL*8 X(*), Y(*), RPARAM(*), QL(*), A(*), W, TMP(NMAX)
      INTEGER JOB, IPARAM(*), IQL(*), IA(*), IPCOND, NZEROS, N, NDIM
C
C ***** ILU PRECONDITIONING
C
      IPCOND = IPARAM(7)
      NZEROS = IPARAM(31)
      NDIM = IPARAM(32)
C
      CALL DAPPLY_ILU_UDIA_L (JOB, QL, IQL, NDIM, NZEROS,
     $                           X, Y, N)
C
      RETURN
      END
C
C
C
      SUBROUTINE PCONDR(JOB, IPARAM, RPARAM, QR, IQR, A, IA, 
     $                  W, X, Y, N)
C
C ***** CALL THE RIGHT PRECONDITIONER 
C
      IMPLICIT NONE

      REAL*8 X(*), Y(*), RPARAM(*), QR(*), A(*), W
      INTEGER JOB, IPARAM(*), IQR(*), IA(*)
      INTEGER IPCOND, NZEROS, NDIM, N
C
C ***** ILU PRECONDITIONING
C
      IPCOND = IPARAM(7)
      NZEROS = IPARAM(31)
      NDIM = IPARAM(32)
C
      CALL DAPPLY_ILU_UDIA_U (JOB, QR, IQR, NDIM, NZEROS,
     $                           X, Y, N)
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 UNSYMMETRIC DIAGONAL
C       FORM
C
      IMPLICIT NONE
       
      INTEGER NDIM, NXNY, NX, NY, J, 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
