PROGRAM LITESE C C Test easy to use driver for iterative linear solver Good Broyden. C C File: maingbe.f C Further program source files used with this program: C gbit1e.f Easy to use driver for linear solver Good Broyden C gbit1 .f Iterative linear solver Good Broyden C easypack.f SLAP-preconditioner and matrix*vector GIANT C interface subroutines and preconditioners and C matrix times vector subroutine C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(MNPDE=1, MN=961) C PARAMETER(MNZMAX=,KMAX=9) PARAMETER(MNZMAX=4805, KMAX=9) C PARAMETER (LRWK=(KMAX+5)*MN+2*KMAX+50) PARAMETER (LRWK=13522) C PARAMETER(LIWKU=10+2*MNZMAX+4*MN+3+11, LRWKU=2*MNZMAX+MNX+MNY+2) PARAMETER(LIWKU=13478, LRWKU=9674) DOUBLE PRECISION USCAL(MN),U(MN),RHS(MN) DOUBLE PRECISION RWK(LRWK), RWKU(LRWKU) INTEGER IWKU(LIWKU), IOPT(50) DOUBLE PRECISION RTOL INTEGER IERR CHARACTER*100 FILNAM C SUN: C PARAMETER (IUMON=10) C MAC: PARAMETER (IUMON=6) C IF (IUMON.NE.6) THEN WRITE(6,*) 'Monitor output filename ?' READ(5,*) FILNAM OPEN(IUMON,FILE=FILNAM) ENDIF C C Linear solver input parameters and options: C RTOL=1.D-5 IERR=0 C DO 30 I=1,50 IOPT(I)=0 RWK(I)=0.0D0 30 CONTINUE C Print one monitor line for each iterate IOPT(13)=-1 C Monitor output unit IOPT(14)=IUMON C Alternative monitor format (0=std,1=alt) C IOPT(17)=1 C Time monitor on/off (1/0) IOPT(19)=1 C Time monitor output unit IOPT(20)=IUMON C Number of saved vectors IOPT(31)=KMAX C Maximum number of iterations IOPT(32)=1000 C C IOPT(35)=type of preconditioner to be used C 0 = Incomplete LU Decomposition Preconditioner SLAP C 1 = Lower Triangle Preconditioner SLAP C 2 = Diagonal Scaling Preconditioner SLAP C 3 = None C 4 = Block diagonal Scaling Preconditioner ZIB C IOPT(35)=0 IUDATA=21 FILNAM='main_gbit_easy.data' OPEN(IUDATA,FILE=FILNAM,STATUS='OLD') READ(IUDATA,*) N, NELT REWIND (IUDATA) C C Iteration starting values and scaling values C DO 40 I=1,N U(I) = 0.0D0 USCAL(I) = 1.0D0 40 CONTINUE C C Read linear systems matrix, right hand side and iteration C startvector C NFILL = NELT JOB = 3 CALL DTIN( N, NFILL, IWKU(12), IWKU(12+NELT), RWKU, ISYM, U, $ RHS, IUDATA, JOB ) IF (NFILL.NE.NELT) STOP 'Error while reading input data' IWKU(1) = 0 C C Call iterative linear solver C CALL GBIT1E(N,NELT,LRWKU,RWKU,LIWKU,IWKU,RHS,U,USCAL,RTOL,IOPT, $ IERR,LRWK,RWK) C ---------------------------------------------------------------- C Output summary results WRITE (IUMON,20001) RWK(32), RTOL, IERR, IOPT(41), IOPT(43) 20001 FORMAT(' Norm(Solution) = ',D18.10,4X,'Achieved Precision = ', $ D18.10,/,' IERR = ',I6,/ $ ' No. iterations=',I5,4X,'No. MULJAC/PRECON calls=',I5) IF (IUMON.NE.6) CLOSE(IUMON) STOP END SUBROUTINE DTIN( N, NELT, IA, JA, A, ISYM, SOLN, RHS, $ IUNIT, JOB ) C***BEGIN PROLOGUE DTIN C***DATE WRITTEN 871119 (YYMMDD) C***REVISION DATE 881213 (YYMMDD) C***CATEGORY NO. D2A4, D2B4 C***KEYWORDS LIBRARY=SLATEC(SLAP), C TYPE=DOUBLE PRECISION(DTIN-D), C Linear system, SLAP Sparse, Diagnostics C***AUTHOR Seager, Mark K., (LLNL) C Lawrence Livermore National Laboratory C PO BOX 808, L-300 C Livermore, CA 94550 (415) 423-3141 C seager@lll-crg.llnl.gov C***PURPOSE Read in SLAP Triad Format Linear System. C Routine to read in a SLAP Triad format matrix and C right hand side and solution to the system, if known. C***DESCRIPTION C *Usage: C INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, IUNIT, JOB C DOUBLE PRECISION A(NELT), SOLN(N), RHS(N) C C CALL DTIN( N, NELT, IA, JA, A, ISYM, SOLN, RHS, IUNIT, JOB ) C C *Arguments: C N :OUT Integer C Order of the Matrix. C NELT :INOUT Integer. C On input NELT is the maximum number of non-zeros that C can be stored in the IA, JA, A arrays. C On output NELT is the number of non-zeros stored in A. C IA :OUT Integer IA(NELT). C JA :OUT Integer JA(NELT). C A :OUT Double Precision A(NELT). C On output these arrays hold the matrix A in the SLAP C Triad format. See "LONG DESCRIPTION", below. C ISYM :OUT Integer. C Flag to indicate symmetric storage format. C If ISYM=0, all nonzero entries of the matrix are stored. C If ISYM=1, the matrix is symmetric, and only the lower C triangle of the matrix is stored. C SOLN :OUT Double Precision SOLN(N). C The solution to the linear system, if present. This array C is accessed if and only if JOB to read it in, see below. C If the user requests that SOLN be read in, but it is not in C the file, then it is simply zeroed out. C RHS :OUT Double Precision RHS(N). C The right hand side vector. This array is accessed if and C only if JOB is set to read it in, see below. C If the user requests that RHS be read in, but it is not in C the file, then it is simply zeroed out. C IUNIT :IN Integer. C Fortran logical I/O device unit number to write the matrix C to. This unit must be connected in a system dependent fashion C to a file or the console or you will get a nasty message C from the Fortran I/O libraries. C JOB :INOUT Integer. C Flag indicating what I/O operations to perform. C On input JOB indicates what Input operations to try to C perform. C JOB = 0 => Read only the matrix. C = 1 => Read matrix and RHS (if present). C = 2 => Read matrix and SOLN (if present). C = 3 => Read matrix, RHS and SOLN (if present). C On output JOB indicates what operations were actually C performed. C JOB = 0 => Read in only the matrix. C = 1 => Read in the matrix and RHS. C = 2 => Read in the matrix and SOLN. C = 3 => Read in the matrix, RHS and SOLN. C C *Precision: Double Precision C *Portability: C You must make sure that IUNIT is a valid Fortran logical C I/O device unit number and that the unit number has been C associated with a file or the console. This is a system C dependent function. C C***LONG DESCRIPTION C The format for the output is as follows. On the first line C are counters and flags: N, NELT, ISYM, IRHS, ISOLN. N, NELT C and ISYM are described above. IRHS is a flag indicating if C the RHS was written out (1 is yes, 0 is no). ISOLN is a C flag indicating if the SOLN was written out (1 is yes, 0 is C no). The format for the fist line is: 5i10. Then comes the C NELT Triad's IA(I), JA(I) and A(I), I = 1, NELT. The format C for these lines is : 1X,I5,1X,I5,1X,E16.7. Then comes C RHS(I), I = 1, N, if IRHS = 1. Then comes SOLN(I), I = 1, C N, if ISOLN = 1. The format for these lines is: 1X,E16.7. C C =================== S L A P Triad format =================== C This routine requires that the matrix A be stored in the C SLAP Triad format. In this format only the non-zeros are C stored. They may appear in *ANY* order. The user supplies C three arrays of length NELT, where NELT is the number of C non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)). For C each non-zero the user puts the row and column index of that C matrix element in the IA and JA arrays. The value of the C non-zero matrix element is placed in the corresponding C location of the A array. This is an extremely easy data C structure to generate. On the other hand it is not too C efficient on vector computers for the iterative solution of C linear systems. Hence, SLAP changes this input data C structure to the SLAP Column format for the iteration (but C does not change it back). C C Here is an example of the SLAP Triad storage format for a C 5x5 Matrix. Recall that the entries may appear in any order. C C 5x5 Matrix SLAP Triad format for 5x5 matrix on left. C 1 2 3 4 5 6 7 8 9 10 11 C |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21 C |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2 C | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1 C | 0 0 0 44 0| C |51 0 53 0 55| C***REFERENCES (NONE) C***ROUTINES CALLED (NONE) C***END PROLOGUE DTIN IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, JOB DOUBLE PRECISION A(NELT), RHS(N), SOLN(N) C C Local variables. C INTEGER IRHS, ISOLN, I, NELTMAX C C Read in the information heading. C***FIRST EXECUTABLE STATEMENT DTIN NELTMAX = NELT READ(IUNIT,1000) N, NELT, ISYM, IRHS, ISOLN NELT = MIN( NELT, NELTMAX ) C C Read in the matrix non-zeros in Triad format. DO 10 I = 1, NELT READ(IUNIT,1010) IA(I), JA(I), A(I) 10 CONTINUE C C If requested, read in the rhs. JOBRET = 0 IF( JOB.EQ.1 .OR. JOB.EQ.3 ) THEN C C Check to see if rhs is in the file. IF( IRHS.EQ.1 ) THEN JOBRET = 1 READ(IUNIT,1020) (RHS(I),I=1,N) ELSE DO 20 I = 1, N RHS(I) = 0.0D0 20 CONTINUE ENDIF ENDIF C C If requested, read in the soln. IF( JOB.GT.1 ) THEN C C Check to see if soln is in the file. IF( ISOLN.EQ.1 ) THEN JOBRET = JOBRET + 2 READ(IUNIT,1020) (SOLN(I),I=1,N) ELSE DO 30 I = 1, N SOLN(I) = 0.0D0 30 CONTINUE ENDIF ENDIF C JOB = JOBRET RETURN 1000 FORMAT(5I10) 1010 FORMAT(1X,I5,1X,I5,1X,E16.7) 1020 FORMAT(1X,E16.7) C------------- LAST LINE OF DTIN FOLLOWS ---------------------------- END