SUBROUTINE GBIT1E(N,NELT,LA,A,LINDAW,INDAWK,RHS,X,XSCAL,RTOL,IOPT, $ IERR,LRWK,RWK) C* Begin Prologue GBIT1E INTEGER N, NELT, LA DOUBLE PRECISION A(NELT) INTEGER LINDAW, INDAWK(LINDAW) DOUBLE PRECISION RHS(N),X(N),XSCAL(N),RTOL INTEGER IOPT(50),IERR INTEGER LRWK DOUBLE PRECISION RWK(LRWK) C ------------------------------------------------------------ C C* Title C C Good Broyden easy to use subroutine - C Iterative solution of a linear system with matrix given in C SLAP Triad or SLAP column format. C C* Written by U. Nowak, L. Weimann C* Purpose Iterative solution of large scale systems of C linear equations C* Method Secant method Good Broyden with adapted C linesearch C* Category D2a. - Large Systems of Linear Equations C* Keywords Linear Equations; Large Systems; C* Iterative Methods C* Version 1.1 C* Revision March 1991 C* Latest Change March 1991 C* Library CodeLib C* Code Fortran 77, Double Precision C* Environment Standard Fortran 77 environment on PC's, C workstations and hosts. C* Copyright (c) Konrad-Zuse-Zentrum fuer C Informationstechnik Berlin (ZIB) C Takustrasse 7, D-14195 Berlin-Dahlem C phone : + 49/30/84185-0 C fax : + 49/30/84185-125 C* Contact Bodo Erdmann C ZIB, Division Scientific Computing, C Department Numerical Analysis and Modelling C phone : + 49/30/84185-185 C fax : + 49/30/84185-107 C e-mail: erdmann@zib.de C C* References C C /1/ P. Deuflhard, R. Freund, A. Walter: C Fast Secant Methods for the Iterative Solution of C Large Nonsymmetric Linear Systems. C ZIB, Preprint SC 90-5 (July 1990) C C /2/ U. Nowak, L. Weimann: C GIANT - A Software Package for the Numerical Solution C of Very Large Systems of Highly Nonlinear Equations. C ZIB, Technical Report TR 90-11 (December 1990) C C --------------------------------------------------------------- C C* Licence C You may use or modify this code for your own non commercial C purposes for an unlimited time. C In any case you should not deliver this code without a special C permission of ZIB. C In case you intend to use the code commercially, we oblige you C to sign an according licence agreement with ZIB. C C* Warranty C This code has been tested up to a certain level. Defects and C weaknesses, which may be included in the code, do not establish C any warranties by ZIB. ZIB does not take over any liabilities C which may follow from aquisition or application of this code. C C* Software status C This code is under care of ZIB and belongs to ZIB software class 1. C C --------------------------------------------------------------- C C C* Parameters list description (* marks inout parameters) C ====================================================== C C* Input parameters (* marks inout parameters) C =========================================== C C N Int Dimension of linear system C NELT Int Number of nonzero elements of the matrix C LA Int Length of the array A. The required minimum C depends on the selected preconditioner C ( see IPRE=IOPT(35) below ). C Required minimum amount is: C for IPRE=0: 2*NELT C for IPRE=1: NELT + INT( (NELT+N)/2) C for IPRE=2: NELT + N C for IPRE=3: NELT C for IPRE=4: NELT + N*NBSIZE C ( for NBSIZE see IOPT(36) below) C * A Dble The first NELT positions hold the real values C of the matrix stored in SLAP Triad or Column C format. The locations up from A(NELT+1) are C used as workspace by the preconditioner sub- C routines. C LINDAW Int Length of the array INDAWK. The required C minimum depends on the selected preconditioner C ( see IPRE=IOPT(35) below ). C Required minimum amount is: C max( 10+2*NELT , LINDPR), where C for IPRE=0: LINDPR = 2*NELT+4*N+14 C for IPRE=1: LINDPR = INT((3*NELT+N)/2)+2*N+13 C for IPRE=2: LINDPR = NELT+N+12 C for IPRE=3: LINDPR = NELT+N+12 C for IPRE=4: LINDPR = NELT+N*(NBSIZE+1)+12 C ( for NBSIZE see IOPT(36) below) C * INDAWK Int The first 11 positions of INDAWK are used as C workspace by EASYPACK. The following elements C holds the index arrays IA and JA of the SLAP C format (see Note 4 for details on the formats). C If INDAWK(1) .EQ. 0 on input, C IA must be stored from INDAWK(12) upto C INDAWK(11+NELT), and JA must be stored from C INDAWK(11+NELT+1) upto INDAWK(11+2*NELT). C The locations up from INDAWK(11+2*NELT+1) are C used as workspace by the preconditioner sub- C routines. C If INDAWK(1)=NZMAX .NE. 0 on input, C IA must be stored from INDAWK(12) upto C INDAWK(11+NELT), and JA must be stored from C INDAWK(11+NZMAX+1) upto INDAWK(11+NZMAX+NELT). C The elements representing JA are internally C restored such that the same storage scheme C results as requested with INDAWK(1) .EQ. 0 . C RHS(N) Dble The right hand side vector C * X(N) Dble The iteration starting vector C * XSCAL(N) Dble Scaling vector C * RTOL Dble Required (relative) precision C * IOPT(50) Int Options and integer workspace array C (see below for usage) C C* Output parameters C ================= C C * X(N) Dble The solution or final iterate vector C * XSCAL(N) Dble Scaling vector. C Modified, if it containes bad input values: C If (XSCAL(I).LT. SMALL) XSCAL(I) = SMALL , C If (XSCAL(I).GT. GREAT) XSCAL(I) = GREAT . C For SMALL and GREAT, see section machine C constants below and regard note 2. C * RTOL Dble Achieved (relative) precision (see note 1.) C * IOPT(50) Int Options and integer workspace array C (see below for usage) C IERR Int Error code - a zero return signals ok. C See below for nonzero error codes. C C* Workspace parameters C ==================== C C LRWK Int Declared dimension of real workspace C Required minimum: (KMAX+5)*N+2*KMAX+50 C RWK(LRWK) Dble Real Workspace C C Note 1. C The iteration stops, if the criterium for convergence is C RHO*sqrt(sig_mid/N) <= RTOL*norm(X_iter) C is satisfied, where X_iter denotes the actual iterate X and C sig_mid the mean sigma value C sig_mid := (sig_iter-2+2*sig_iter-1+sig_iter)/4, C with the raw error estimate sig_j for the iterate j. The C factor RHO may be modified - see RWK(21). C The norm used herein is defined by C norm(Z) := sqrt( ( Sum(i=1,N) (Z(i)/XSCAL(i))**2 )/N ). C C Note 2. C The machine dependent values SMALL, GREAT and EPMACH are C gained from calls of the machine constants function ZIBCONST. C As delivered, this function is adapted to use constants C suitable for all machines with IEEE arithmetic. If you use C another type of machine, you may change ZIBCONST. C C* Options IOPT: C ============= C C Pos. Name Default Meaning C C 1 QSUCC 0 Indicator for the iteration modus: C =0 : A new iteration will be started C =1 : A previously terminated iteration will C be continued - in general without a C restart C (usally with a smaller tolerance RTOL C prescribed as before or a larger C maximum iteration number count). C 2..12 Reserved C 13 MPRMON 0 Output level of print monitor C = 0 : no output will be written C = 1 : only a summary output will be written C > 1 : reserved for future use C =-j : A summary output and additionally each C j-th iterates statistics will be C written C 14 LUMON 6 Logical unit number for print monitor C 15..16 Reserved C 17 MPROPT 0 Print monitor option: C = 0: Standard print monitor C = 1: Test print monitor for special purposes C 18 Reserved C 19 MPRTIM 0 Output level for the time monitor C = 0 : no time measurement and no output C = 1 : time measurement will be done and C summary output will be written - C regard note 3. C 20 LUTIM 6 Logical output unit for time monitor C 21 MPRSPE 0 Output level for special information C =0: No special output to be done C =1: Special information will be written C The special information may be used to C generate a graphic display of error- and C solution norms and of restart behaviour. C 22 LUSPE 6 Logical output unit for special information C 23..30 Reserved C 31 KMAX 9 Maximum number of latest iterates to be saved C and used by the iterative linear solver. C An input <=-2 means KMAX=0 will be used. C A "-1" input means, that there no limit C applies, e.g. the core routine will be C served with the maximum possible value C allowed by the total workspace amount. C A zero input means KMAX=9 (or smaller C downto a reasonable minimum value, C if workspace is not sufficient). C 32 LITMAX 100 Maximum number of iterations allowed C 33..34 Reserved C 35 IPRE 0 Type of preconditioner, valid values are: C 0: The SLAP preconditioner ILU (Incomplete LU- C decomposition) will be used, C 1: The SLAP preconditioner lower triangle C will be used, C 2: The SLAP preconditioner diagonal scaling C will be used, C 3: No preconditioner will be used, C 4: The ZIB supplied preconditioner block diagonal C scaling will be used. C 36 NBSIZE 2 Only, if IPRE=4: the size of the diagonal C blocks. If IPRE=4 is selected and N is not a C multiple of NBSIZE, an error termination C occurs. C 37..40 Reserved C 41 NITER (Int-Ws) Number of iteration steps totally done C 42 K (Int-Ws) Number of iteration steps done since C latest internal restart due to restrictions C implied by KMAX, TAUMIN and TAUMAX. C 43 NMULJ (Int-Ws) Number of (pairwise) calls of MULJAC and C PRECON C 44 NRW (Int-Ws) Amount of real workspace used C 45 NRWKFR (Int-Ws) First element of RWK which is free to be used C as workspace between successive calls of C GBIT1. C 46..50 Reserved C C Note 3. C The integrated time monitor calls the machine dependent C subroutine SECOND to get the actual time stamp in form C of a real number (Single precision). As delivered, this C subroutine always return 0.0 as time stamp value. Refer C to the compiler- or library manual of the FORTRAN compiler C which you actually use to find out how to get the actual C time stamp on your machine. C C Note 4. C Description of the Jacobian SLAP format C (as extracted from the SLAP package description): 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 C =================== S L A P Column format ================== C This routine requires that the matrix A be stored in the C SLAP Column format. In this format the non-zeros are stored C counting down columns (except for the diagonal entry, which C must appear first in each "column") and are stored in the C double precision array A. In other words, for each column C in the matrix put the diagonal entry in A. Then put in the C other non-zero elements going down the column (except the C diagonal) in order. The IA array holds the row index for C each non-zero. The JA array holds the offsets into the IA, C A arrays for the beginning of each column. That is, C IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the C ICOL-th column in IA and A. IA(JA(ICOL+1)-1), C A(JA(ICOL+1)-1) points to the end of the ICOL-th column. C Note that we always have JA(N+1) = NELT+1, where N is the C number of columns in the matrix and NELT is the number of C non-zeros in the matrix. C C Here is an example of the SLAP Column storage format for a C 5x5 Matrix (in the A and IA arrays '|' denotes the end of a C column): C C 5x5 Matrix SLAP Column 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: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35 C |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3 C | 0 0 33 0 35| JA: 1 4 6 8 9 12 C | 0 0 0 44 0| C |51 0 53 0 55| C C* Optional REAL input/output in RWK: C ==================================== C C Pos. Name Meaning C C 1..20 Reserved C 21 RHO IN Security factor in error estimate (RHO.GE.1). C Default: 4.0D0 C 22 TAUMIN IN Minimum accepted stepsize factor tk. C Default 1.0D-8 C 23 TAUMAX IN Maximum accepted stepsize factor tk. C Default: 1.0D2 C 24 TAUEQU IN Parameter for "near equal" cycle check of C rejected tau values. Default: 1.0D-2 C 25..30 Reserved C 31 DIFNRM OUT The norm2 of the latest correction of the C iterate. C 32 SOLNRM OUT The norm2 of the final iterate. C 33..40 Reserved C 41 SIGMAK INTERN Norm2-square of latest raw correction DELTAK. C 42 TAUN1 INTERN First kept tau value for tau-cycle check. C 43 TAUN2 INTERN Second kept tau value for tau-cycle check. C 44 SOLNRA INTERN Saves norm2 of the previous iterate. C 45..47 C SIGH INTERN Storage to hold raw error estimates of latest C three iterates for mean error calculation. C 48..50 Reserved C C Error codes (Stored to IERR): C ============================= C C -1 Iteration number limit (as set by LITMAX) exceeded C -2 Iteration diverges: The error estimate exceeds an C internal maximum allowed value (DIFLIM). C -3 Iteration diverges: Bad tau values and norms of C corrections oscillate in a two-iterates-cycle. C -4 Iteration cannot continue: bad tau value immediate C at start or restart. C -10 Insufficient workspace C -20 Bad input to dimensional parameter N C -21 Nonpositive value for RTOL supplied C -22 Negative scaling value via vector XSCAL supplied C -25 Block diagonal scaling selected, and N is not a multiple C of the specified diagonal blocks size NBSIZE. C -30 One or more fields specified in IOPT are invalid C (for more information, see error-printout) C C* Machine dependent constants used: C ================================= C C GREAT = squareroot of maxreal divided by 10 C SMALL = squareroot of "smallest positive machine number C divided by relative machine precision" C DOUBLE PRECISION GREAT,SMALL (used in GBPCHK) C ------------------------------------------------------------ C* End Prologue EXTERNAL MULJAC, PRECON INTEGER MPRMON, LUMON, NBSIZE C LUMON = IOPT(14) IF (LUMON.LT.1 .OR. LUMON.GT.99) LUMON=6 MPRMON = IOPT(13) C NZMAX = INDAWK(1) IF (NZMAX.NE.0) THEN IF (NZMAX.LT.NELT) THEN IERR = -26 IF (MPRMON.GT.0) WRITE(LUMON, 10001) NZMAX, NELT 10001 FORMAT(1X,'Invalid NZMAX (=INDAWK(1)) specified: ',I7, $ ' - must be greater than NELT = ',I8,/ $ ' or zero if equal to NELT') RETURN ELSE IF (NZMAX.GT.NELT) THEN DO 10 I=1,NELT INDAWK(11+NELT+I) = INDAWK(11+NZMAX+I) 10 CONTINUE ENDIF ENDIF INDAWK(1) = LA INDAWK(2) = LINDAW INDAWK(3) = 0 INDAWK(4) = 1 INDAWK(5) = NELT INDAWK(7) = LUMON INDAWK(8) = MPRMON INDAWK(9) = IOPT(35) IF ( IOPT(35) .EQ. 4 ) THEN NBSIZE = IOPT(36) IF (NBSIZE.EQ.0) NBSIZE = 2 IF ( MOD(N,NBSIZE) .NE. 0 ) THEN IERR = -25 IF (MPRMON.GT.0) WRITE(LUMON, 10002) N, NBSIZE 10002 FORMAT(1X,'Incompatible size N of the linear system and ', $ /,'block size NBSIZE for block diagonal scaling',/, $ 'N must be a multiple of NBSIZE - ', $ 'N = ',I7,' NBSIZE = ',I4) RETURN ENDIF INDAWK(10) = NBSIZE ENDIF INDAWK(11) = NELT C Setup preconditioner by EASYPACK subroutine EJAC CALL EJAC(DUMMY,N,X,XSCAL,RDUMMY,A,INDAWK,IDUMMY,IERR) IF (IERR.NE.0) RETURN C Call Good Broyden IOPT(35) = 0 IOPT(36) = 0 CALL GBIT1(N,MULJAC,PRECON,RHS,X,XSCAL,RTOL,IOPT,IERR, $ LRWK,RWK,LINDAW,INDAWK,LA,A) RETURN END C C SUBROUTINE JAC1(FCN,N,U,UWGT,F,NELT,IDUMMY,A,IA,JA,NFILL, $ RWKU,IWKU,NJAC,IFAIL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) EXTERNAL FCN INTEGER N DOUBLE PRECISION U(N),UWGT(N),F(N) INTEGER NELT,IDUMMY DOUBLE PRECISION A(NELT) INTEGER IA(NELT),JA(NELT) INTEGER NFILL DOUBLE PRECISION RWKU(*) INTEGER IWKU(*) INTEGER NJAC,IFAIL C C Pass number of nonzero matrix elements to EASYPACK NFILL = IWKU(1) RETURN END C C SUBROUTINE FCN1() C C Dummy external to satisfy reference from EASYPACK (not called) C RETURN END