c
c*****************************************************************************
c*                                                                           *
c*  Copyright Compaq Computer Corporation 1999.         All rights reserved. *
c*                                                                           *
c*****************************************************************************
c
c
c
      PROGRAM EXAMPLE_ITSOL_1
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, USING A
C           USER-DEFINED STORAGE SCHEME.
C       (2) USE OF A USER-DEFINED ROUTINE MATVEC
C       (3) USE OF A USER DEFINED ROUTINE MSTOP
C       (4) USE OF USER DEFINED ROUTINE PCONDL
C       (5) USE OF THE ROUTINE DITSOL_DEFAULTS
C       (6) INFORMATION PRINTED OUT FOR IOLEVEL = 3
C       (7) PRINTING MESSAGES TO FILE. 
C
C 
C ***** SAMPLE OUTPUT FOLLOWS
C
C SOLVING EXAMPLE PROBLEM WITH SPD SPLIT PRECONDITIONED CG
C DIAGONAL PRECONDITIONING USED 
C  method used : cg with spd split preconditioning
C  order of system =    100
C  stopping criterion used :     0
C  maximum iterations allowed:   100
C  tolerance for convergence :   0.10000000E-05
C using user defined stopping criterion
C     iteration =          0    stopping test = 0.20000000E+01
C     iteration =          1    stopping test = 0.92307692E+00
C     iteration =          2    stopping test = 0.58793970E+00
C     iteration =          3    stopping test = 0.63509006E+00
C     iteration =          4    stopping test = 0.42973063E+00
C     iteration =          5    stopping test = 0.48724587E+00
C     iteration =          6    stopping test = 0.41781094E+00
C     iteration =          7    stopping test = 0.50288290E+00
C     iteration =          8    stopping test = 0.17070529E+00
C     iteration =          9    stopping test = 0.38502953E-01
C     iteration =         10    stopping test = 0.16076790E-01
C     iteration =         11    stopping test = 0.75548469E-02
C     iteration =         12    stopping test = 0.15431168E-02
C     iteration =         13    stopping test = 0.14025362E-03
C     iteration =         14    stopping test = 0.42745516E-05
C     iteration =         15    stopping test = 0.63333204E-15
C  solution obtained after          15  iterations
C  normal exit from solver
C  final value of stopping test = 0.63333204E-15
C    TRUE SOLUTION     SOLUTION FROM SOLVER     ABS. DIFFERENCE
C     0.10000000E+01      0.10000000E+01      0.33306691E-15
C     0.10000000E+01      0.10000000E+01      0.44408921E-15
C     0.10000000E+01      0.10000000E+01      0.44408921E-15
C     0.10000000E+01      0.10000000E+01      0.66613381E-15
C                               .
C                               .
C                               .
C                       (edited for brevity)
C                               .
C                               .
C                               .
C     0.10000000E+01      0.10000000E+01      0.88817842E-15
C     0.10000000E+01      0.10000000E+01      0.11102230E-15
C     0.10000000E+01      0.10000000E+01      0.33306691E-15
C     0.10000000E+01      0.10000000E+01      0.33306691E-15
C MAX ERROR IN SOLUTION=  0.23314684E-14
C
C
C

      IMPLICIT NONE 

      INTEGER NMAX
      PARAMETER (NMAX = 100)
C
      REAL*8 X(NMAX), XO(NMAX), RHS(NMAX), QL(NMAX)
      REAL*8  RPARAM(50), RWORK(4*NMAX), DUM, TEMP
      REAL*8 A(NMAX)
C
      INTEGER IPARAM(50), IA(2), IDUM
      INTEGER I,NX, NY, NXNY, IOUNIT, IERROR
C
      REAL*8 A1(NMAX), A2(NMAX), A3(NMAX), A4(NMAX), A5(NMAX)

      COMMON /MATRIX/ A1, A2, A3, A4, A5
C
      EXTERNAL MATVEC, PCONDL, USER_MSTOP, USER_PRINT_ROUTINE
C
C ***** SET UP PROBLEM SIZE
C
      NX = 10
      NY = 10
      NXNY = NX*NY
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 ROUTINE DEFAULTS
C           CHANGE IPCOND TO 4 (FOR SPD SPLIT PRECONDITIONING)
C           CHANGE IOLEVEL TO 3
C           CHANGE ISTOP TO 0 (FOR USER DEFINED MSTOP)
C
      IPARAM(3) = 0
      IPARAM(4) = 4*NXNY
C
      IPARAM(5) = 7
      IPARAM(6) = 3
      IPARAM(7) = 4
      IPARAM(8) = 0
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 SPD SPLIT
     $ PRECONDITIONED CG',/,2X,'DIAGONAL PRECONDITIONING USED ',/)
C 
C ***** GENERATE THE MATRIX (USE IA TO PASS NX AND NY TO SOLVER)
C
      CALL GENMAT(NX, NY, NXNY)
      IA(1) = NX
      IA(2) = NY
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, IA, 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 DIAGONAL PRECONDITIONER
C
      CALL GEN_PCOND (NXNY, QL)
C
C ***** CALL THE SOLVER
C
       CALL DITSOL_PCG ( MATVEC, PCONDL, DUM, USER_MSTOP,
     $                   A, IA, X, RHS,
     $                   NXNY, QL, 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 THE VECTOR X BY THE MATRIX TO OBTAIN VECTOR Y
C       ONLY JOB = 0 NEEDED FOR THIS EXAMPLE
C
      IMPLICIT NONE 
C
    
      REAL*8 X(*), Y(*), A(*), RPARAM(*), W
      INTEGER IA(*),IPARAM(*), N, JOB
       
C
      CALL MULA (IA(1), IA(2), N, X, Y)
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       ONLY JOB = 0 NEEDED FOR THIS EXAMPLE
C
      IMPLICIT NONE 

      REAL*8 X(*), Y(*), A(*), QL(*),RPARAM(*), W
      INTEGER IA(*), IPARAM(*), N, IQL, JOB
C
C ***** DIAGONAL PRECONDITIONING
C
      CALL APPLY_PCOND_DIA (N, X, Y, QL)
C
      RETURN
      END
C
C
C
      SUBROUTINE GEN_PCOND (N, QL)
C
C ***** GENERATE THE LEFT PRECONDITIONER
C
      IMPLICIT NONE 

      INTEGER NMAX
      PARAMETER (NMAX = 100)
      INTEGER N,I

      REAL*8 QL(*), A1(NMAX), A2(NMAX), A3(NMAX), A4(NMAX), A5(NMAX)
      COMMON /MATRIX/ A1, A2, A3, A4, A5

C
      DO I = 1, N
          QL(I) = 1.0D0 / A3(I)
      END DO
C
      RETURN
      END
C
C
C
      SUBROUTINE APPLY_PCOND_DIA (N, X, Y, QL)
C
C ***** APPLY THE DIAGONAL  PRECONDITIONER
C
      IMPLICIT NONE 

      REAL*8 X(*), Y(*), QL(*)
      INTEGER N,I
C
      DO I = 1, N
          Y(I) = X(I) * QL(I)
      END DO
C
      RETURN
      END
C
C
C
      SUBROUTINE GENMAT (NX, NY, NXNY)
C
C ***** GENERATE THE MATRIX FOR THE EXAMPLE
C
      IMPLICIT NONE 

      INTEGER NMAX
      PARAMETER (NMAX = 100)
      REAL*8 A1(NMAX), A2(NMAX), A3(NMAX), A4(NMAX), A5(NMAX)
      COMMON /MATRIX/ A1, A2, A3, A4, A5
      INTEGER NXNY, I, J, K, NX, NY
C
      DO I = 1, NXNY
          A3(I) = 4.0D0
          A1(I) = 0.0D0
          A2(I) = 0.0D0
          A4(I) = 0.0D0
          A5(I) = 0.0D0
      END DO
C
      DO J = 1, NY
          DO I = 2, NX
              K = (J-1)*NX+I
              A2(K) = -1.0D0
          END DO
          DO I = 1, NX-1
              K = (J-1)*NX+I
              A4(K) = -1.0D0
          END DO
      END DO
C
      DO J = 2, NY
          DO I = 1, NX
              K = (J-1)*NX+I
              A1(K) = -1.0D0
          END DO
      END DO
C
      DO J = 1, NY-1
          DO I = 1, NX
              K = (J-1)*NX+I
              A5(K) = -1.0D0
          END DO
      END DO
C
      RETURN
      END

C
C
C

      SUBROUTINE MULA (NX, NY, NXNY, TMP1, TMP2)
C
C ***** TO OBTAIN THE MATRIX VECTOR MULTIPLY
C
      IMPLICIT NONE 

      INTEGER NMAX
      PARAMETER (NMAX = 100)
      REAL*8 TMP1(*), TMP2(*)
      INTEGER NX,NY,NXNY, I

      REAL*8 A1(NMAX), A2(NMAX), A3(NMAX), A4(NMAX), A5(NMAX)
      COMMON /MATRIX/ A1, A2, A3, A4, A5
C
      DO I = 1, NXNY
          TMP2(I) = A3(I)*TMP1(I) 
      END DO
C
      DO I = 1, NXNY-1
          TMP2(I) = TMP2(I) +
     $              A4(I)*TMP1(I+1)
      END DO
C
      DO I = 2, NXNY
          TMP2(I) = TMP2(I) +
     $              A2(I)*TMP1(I-1) 
      END DO
C
      DO I = 1, NXNY-NX
          TMP2(I) = TMP2(I) +
     $              A5(I)*TMP1(I+NX) 
      END DO
C
      DO I = NX+1, NXNY
          TMP2(I) = TMP2(I) +
     $              A1(I)*TMP1(I-NX) 
      END DO
C
      RETURN
      END
C
C
C
      SUBROUTINE USER_MSTOP (IPARAM, RPARAM, X, R, Z, B, N)
C
C ***** ROUTINE FOR THE USER PROVIDED STOPPING CRITERION 
C
      IMPLICIT NONE 

      REAL*8 X(*), B(*), R(*), Z(*), RPARAM(*), STPTST
      INTEGER ITERS, I, IOUNIT, IPARAM(*), N
C
      IOUNIT = IPARAM(5)
      ITERS = IPARAM(10)
      IF (ITERS.EQ.0) THEN
          WRITE(IOUNIT, 100) 
  100     FORMAT(/,2X,'USING USER DEFINED STOPPING CRITERION',/)
      END IF
C
C ***** USER DEFINED STOPPING CRITERION USES MAX NORM OF RESIDUAL FOR
C       STPTST
C
      STPTST = ABS (R(1))
      DO I = 2, N
         STPTST = MAX ( STPTST,ABS(R(I)) )
      END DO
C
      RPARAM(2) = STPTST
C
      RETURN
      END

