c
c
c*****************************************************************************
c*                                                                           *
c*  Copyright Compaq Computer Corporation 1999.         All rights reserved. *
c*                                                                           *
c*****************************************************************************
c
c
c
      PROGRAM EXAMPLE_ITSOL_3
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 LEAST SQUARES CONJUGATE GRADIENT METHOD, 
C           WITH LEFT PRECONDITIONING USING ILU FACTORIZATION. THE 
C           MATRIX IS STORED USING GENERAL STORAGE BY ROWS.
C       (2) USE OF ROUTINE MATVEC (DMATVEC_GENR)
C       (3) USE OF ROUTINE PCONDL TO CALL ROUTINES FOR APPLYING  
C           THE PRECONDITIONER (DAPPLY_ILU_GENR_L AND DAPPLY_ILU_GENR_U)
C       (4) USE OF THE ROUTINE DITSOL_DEFAULTS
C
C 
C ***** SAMPLE OUTPUT FOLLOWS 
C
C SOLVING EXAMPLE PROBLEM WITH LEFT PRECONDITIONED
C  LEAST SQUARES CG. ILU PRECONDITIONING USED 
C MATRIX STORED USING GENERAL SPARSE STORAGE BY ROWS
C  method used : least squares cg with left 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          17  iterations
C  normal exit from solver
C  final value of stopping test = 0.25683815E-06
C    TRUE SOLUTION     SOLUTION FROM SOLVER     ABS. DIFFERENCE
C     0.10000000E+01      0.99999999E+00      0.65643619E-08
C     0.10000000E+01      0.10000000E+01      0.14181074E-08
C     0.10000000E+01      0.10000000E+01      0.79716975E-08
C     0.10000000E+01      0.10000000E+01      0.48284889E-08
C                               .
C                               .
C                               .
C                       (edited for brevity)
C                               .
C                               .
C                               .
C     0.10000000E+01      0.10000000E+01      0.44230040E-08
C     0.10000000E+01      0.10000000E+01      0.45262150E-08
C     0.10000000E+01      0.10000000E+01      0.42126541E-08
C     0.10000000E+01      0.99999999E+00      0.78054114E-08
C MAX ERROR IN SOLUTION =  0.21752978E-07
C
C
C
      IMPLICIT NONE
C
      INTEGER NMAX
      PARAMETER (NMAX = 100)
      INTEGER NDIM
      PARAMETER (NDIM = 100)
      INTEGER NZ
      PARAMETER (NZ = 5*NMAX)
      INTEGER MAX_NRWK 
      PARAMETER (MAX_NRWK = 5*NMAX)
C
      REAL*8 X(NMAX), XO(NMAX), RHS(NMAX)
      REAL*8 RPARAM(50), RWORK(MAX_NRWK), DUM
      REAL*8 A_GENR(NZ), P_ILU(NZ), TEMP
C
      INTEGER INDEX_GENR(NMAX+1+NZ)
      INTEGER IPARAM(50), IDUM
      INTEGER I, NX, NY, NXNY, IOUNIT, NZEROS, IERROR
C
      EXTERNAL MATVEC, PCONDL
      EXTERNAL 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 DEFAULTS
C              CHANGE IOUNIT TO 7
C              CHANGE IPCOND TO 1 (FOR LEFT PRECONDITIONING)
C              CHANGE IOLEVEL TO 2
C
      IPARAM(3) = 0
      IPARAM(4) = MAX_NRWK
C
      IPARAM(5) = 7
      IPARAM(6) = 2
      IPARAM(7) = 1
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 LEFT
     $ PRECONDITIONED',/,2X,' LEAST SQUARES CG.',
     $ ' ILU PRECONDITIONING USED ',/
     $ 2X,'MATRIX STORED USING GENERAL SPARSE STORAGE BY ROWS',/)
C 
C ***** GENERATE THE MATRIX
C
      CALL GENMAT(NX, NY, NXNY, A_GENR, INDEX_GENR, NZEROS)
      IPARAM(31) = NZEROS
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_GENR, INDEX_GENR, 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_GENR (A_GENR, INDEX_GENR(1), INDEX_GENR(NXNY+2),
     $                      NZEROS, P_ILU, NXNY)
C
C ***** CALL THE SOLVER
C
       CALL DITSOL_PLSCG ( MATVEC, PCONDL, DUM, DUM, 
     $                     A_GENR, INDEX_GENR, X, RHS, NXNY, 
     $                     P_ILU, INDEX_GENR, 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(*), NZEROS, JOB, IPARAM(*), N
C
      NZEROS = IPARAM(31)
C
      CALL DMATVEC_GENR (JOB, A, IA(1), IA(N+2), 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 LEFT PRECONDITIONER 
C
      IMPLICIT NONE

      INTEGER NMAX
      PARAMETER (NMAX = 100)
      REAL*8 X(*), Y(*), RPARAM(*), QL(*), A(*), TMP(NMAX), W
      INTEGER JOB, IPARAM(*), IQL(*), IA(*), IPCOND, NZEROS, N
      
C
C ***** ILU PRECONDITIONING
C
      IPCOND = IPARAM(7)
      NZEROS = IPARAM(31)
C
      IF (JOB.EQ.0)  THEN
C
          CALL DAPPLY_ILU_GENR_L (JOB, QL, IQL(1), IQL(N+2), NZEROS,
     $                           X, TMP, N)
          CALL DAPPLY_ILU_GENR_U (JOB, QL, IQL(1), IQL(N+2), NZEROS,
     $                           TMP, Y, N)
C
      ELSE IF (JOB.EQ.1) THEN
C
          CALL DAPPLY_ILU_GENR_U (JOB, QL, IQL(1), IQL(N+2), NZEROS,
     $                           X, TMP, N)
          CALL DAPPLY_ILU_GENR_L (JOB, QL, IQL(1), IQL(N+2), NZEROS,
     $                           TMP, Y, N)
C
      END IF
C
      RETURN
      END
C
C
C
      SUBROUTINE GENMAT (NX, NY, NXNY, A, IA, NZEROS)
C
C ***** GENERATE THE MATRIX FOR THE EXAMPLE USING THE GENERAL STORAGE 
C       BY ROWS FORM. JA IS STORED IN IA STARTING FROM N+2 ONWARDS.
C
      IMPLICIT NONE

      REAL*8 A(*)
      INTEGER IA(*), JBEG, NXNY, NZEROS, NX, NY, J, I, K
C
      JBEG = NXNY+1
      IA(1) = 1
      NZEROS = 0
C          
      DO J = 1, NY
          DO I = 1, NX
              K = (J-1)*NX+I
              IF (J.GT.1) THEN
                  NZEROS = NZEROS+1
                  A(NZEROS) = -1.0D0
                  IA(JBEG+NZEROS) = K-NX        
              END IF
              IF (I.GT.1) THEN
                  NZEROS = NZEROS+1
                  A(NZEROS) = -1.0D0
                  IA(JBEG+NZEROS) = K-1              
              END IF
              IF (J.LT.NY) THEN
                  NZEROS = NZEROS+1
                  A(NZEROS) = -1.0D0
                  IA(JBEG+NZEROS) = K+NX              
              END IF
              IF (I.LT.NX) THEN
                  NZEROS = NZEROS+1
                  A(NZEROS) = -1.0D0
                  IA(JBEG+NZEROS) = K+1              
              END IF
              NZEROS = NZEROS+1
              A(NZEROS) = 4.0D0
              IA(JBEG+NZEROS) = K
C
              IA(K+1) = NZEROS+1
C
          END DO
C
      END DO
C
      RETURN
      END
