1*DECK DSDOMN 2 SUBROUTINE DSDOMN (N, B, X, NELT, IA, JA, A, ISYM, NSAVE, ITOL, 3 + TOL, ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW) 4C***BEGIN PROLOGUE DSDOMN 5C***PURPOSE Diagonally Scaled Orthomin Sparse Iterative Ax=b Solver. 6C Routine to solve a general linear system Ax = b using 7C the Orthomin method with diagonal scaling. 8C***LIBRARY SLATEC (SLAP) 9C***CATEGORY D2A4, D2B4 10C***TYPE DOUBLE PRECISION (SSDOMN-S, DSDOMN-D) 11C***KEYWORDS ITERATIVE PRECONDITION, NON-SYMMETRIC LINEAR SYSTEM SOLVE, 12C SLAP, SPARSE 13C***AUTHOR Greenbaum, Anne, (Courant Institute) 14C Seager, Mark K., (LLNL) 15C Lawrence Livermore National Laboratory 16C PO BOX 808, L-60 17C Livermore, CA 94550 (510) 423-3141 18C seager@llnl.gov 19C***DESCRIPTION 20C 21C *Usage: 22C INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, NSAVE, ITOL, ITMAX 23C INTEGER ITER, IERR, IUNIT, LENW, IWORK(10), LENIW 24C DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR 25C DOUBLE PRECISION RWORK(7*N+3*N*NSAVE+NSAVE) 26C 27C CALL DSDOMN(N, B, X, NELT, IA, JA, A, ISYM, NSAVE, ITOL, TOL, 28C $ ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW ) 29C 30C *Arguments: 31C N :IN Integer. 32C Order of the Matrix. 33C B :IN Double Precision B(N). 34C Right-hand side vector. 35C X :INOUT Double Precision X(N). 36C On input X is your initial guess for solution vector. 37C On output X is the final approximate solution. 38C NELT :IN Integer. 39C Number of Non-Zeros stored in A. 40C IA :IN Integer IA(NELT). 41C JA :IN Integer JA(NELT). 42C A :IN Double Precision A(NELT). 43C These arrays should hold the matrix A in either the SLAP 44C Triad format or the SLAP Column format. See "Description", 45C below. If the SLAP Triad format is chosen, it is changed 46C internally to the SLAP Column format. 47C ISYM :IN Integer. 48C Flag to indicate symmetric storage format. 49C If ISYM=0, all non-zero entries of the matrix are stored. 50C If ISYM=1, the matrix is symmetric, and only the upper 51C or lower triangle of the matrix is stored. 52C NSAVE :IN Integer. 53C Number of direction vectors to save and orthogonalize against. 54C ITOL :IN Integer. 55C Flag to indicate type of convergence criterion. 56C If ITOL=1, iteration stops when the 2-norm of the residual 57C divided by the 2-norm of the right-hand side is less than TOL. 58C If ITOL=2, iteration stops when the 2-norm of M-inv times the 59C residual divided by the 2-norm of M-inv times the right hand 60C side is less than TOL, where M-inv is the inverse of the 61C diagonal of A. 62C ITOL=11 is often useful for checking and comparing different 63C routines. For this case, the user must supply the "exact" 64C solution or a very accurate approximation (one with an error 65C much less than TOL) through a common block, 66C COMMON /DSLBLK/ SOLN( ) 67C If ITOL=11, iteration stops when the 2-norm of the difference 68C between the iterative approximation and the user-supplied 69C solution divided by the 2-norm of the user-supplied solution 70C is less than TOL. 71C TOL :INOUT Double Precision. 72C Convergence criterion, as described above. (Reset if IERR=4.) 73C ITMAX :IN Integer. 74C Maximum number of iterations. 75C ITER :OUT Integer. 76C Number of iterations required to reach convergence, or 77C ITMAX+1 if convergence criterion could not be achieved in 78C ITMAX iterations. 79C ERR :OUT Double Precision. 80C Error estimate of error in final approximate solution, as 81C defined by ITOL. 82C IERR :OUT Integer. 83C Return error flag. 84C IERR = 0 => All went well. 85C IERR = 1 => Insufficient space allocated for WORK or IWORK. 86C IERR = 2 => Method failed to converge in ITMAX steps. 87C IERR = 3 => Error in user input. 88C Check input values of N, ITOL. 89C IERR = 4 => User error tolerance set too tight. 90C Reset to 500*D1MACH(3). Iteration proceeded. 91C IERR = 5 => Preconditioning matrix, M, is not positive 92C definite. (r,z) < 0. 93C IERR = 6 => Breakdown of method detected. 94C (p,Ap) < epsilon**2. 95C IUNIT :IN Integer. 96C Unit number on which to write the error at each iteration, 97C if this is desired for monitoring convergence. If unit 98C number is 0, no writing will occur. 99C RWORK :WORK Double Precision RWORK(LENW). 100C Double Precision array used for workspace. 101C LENW :IN Integer. 102C Length of the double precision workspace, RWORK. 103C LENW >= 7*N+NSAVE*(3*N+1). 104C IWORK :WORK Integer IWORK(LENIW). 105C Used to hold pointers into the RWORK array. 106C LENIW :IN Integer. 107C Length of the integer workspace, IWORK. LENIW >= 10. 108C 109C *Description: 110C This routine is simply a driver for the DOMN routine. It 111C calls the DSDS routine to set up the preconditioning and 112C then calls DOMN with the appropriate MATVEC and MSOLVE 113C routines. 114C 115C The Sparse Linear Algebra Package (SLAP) utilizes two matrix 116C data structures: 1) the SLAP Triad format or 2) the SLAP 117C Column format. The user can hand this routine either of the 118C of these data structures and SLAP will figure out which on 119C is being used and act accordingly. 120C 121C =================== S L A P Triad format =================== 122C 123C In this format only the non-zeros are stored. They may 124C appear in *ANY* order. The user supplies three arrays of 125C length NELT, where NELT is the number of non-zeros in the 126C matrix: (IA(NELT), JA(NELT), A(NELT)). For each non-zero 127C the user puts the row and column index of that matrix 128C element in the IA and JA arrays. The value of the non-zero 129C matrix element is placed in the corresponding location of 130C the A array. This is an extremely easy data structure to 131C generate. On the other hand it is not too efficient on 132C vector computers for the iterative solution of linear 133C systems. Hence, SLAP changes this input data structure to 134C the SLAP Column format for the iteration (but does not 135C change it back). 136C 137C Here is an example of the SLAP Triad storage format for a 138C 5x5 Matrix. Recall that the entries may appear in any order. 139C 140C 5x5 Matrix SLAP Triad format for 5x5 matrix on left. 141C 1 2 3 4 5 6 7 8 9 10 11 142C |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21 143C |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2 144C | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1 145C | 0 0 0 44 0| 146C |51 0 53 0 55| 147C 148C =================== S L A P Column format ================== 149C 150C In this format the non-zeros are stored counting down 151C columns (except for the diagonal entry, which must appear 152C first in each "column") and are stored in the double pre- 153C cision array A. In other words, for each column in the 154C matrix first put the diagonal entry in A. Then put in the 155C other non-zero elements going down the column (except the 156C diagonal) in order. The IA array holds the row index for 157C each non-zero. The JA array holds the offsets into the IA, 158C A arrays for the beginning of each column. That is, 159C IA(JA(ICOL)),A(JA(ICOL)) are the first elements of the ICOL- 160C th column in IA and A, and IA(JA(ICOL+1)-1), A(JA(ICOL+1)-1) 161C are the last elements of the ICOL-th column. Note that we 162C always have JA(N+1)=NELT+1, where N is the number of columns 163C in the matrix and NELT is the number of non-zeros in the 164C matrix. 165C 166C Here is an example of the SLAP Column storage format for a 167C 5x5 Matrix (in the A and IA arrays '|' denotes the end of a 168C column): 169C 170C 5x5 Matrix SLAP Column format for 5x5 matrix on left. 171C 1 2 3 4 5 6 7 8 9 10 11 172C |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35 173C |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3 174C | 0 0 33 0 35| JA: 1 4 6 8 9 12 175C | 0 0 0 44 0| 176C |51 0 53 0 55| 177C 178C *Side Effects: 179C The SLAP Triad format (IA, JA, A) is modified internally to 180C be the SLAP Column format. See above. 181C 182C *Cautions: 183C This routine will attempt to write to the Fortran logical output 184C unit IUNIT, if IUNIT .ne. 0. Thus, the user must make sure that 185C this logical unit is attached to a file or terminal before calling 186C this routine with a non-zero value for IUNIT. This routine does 187C not check for the validity of a non-zero IUNIT unit number. 188C 189C***SEE ALSO DOMN, DSLUOM 190C***REFERENCES (NONE) 191C***ROUTINES CALLED DCHKW, DOMN, DS2Y, DSDI, DSDS, DSMV 192C***REVISION HISTORY (YYMMDD) 193C 890404 DATE WRITTEN 194C 890404 Previous REVISION DATE 195C 890915 Made changes requested at July 1989 CML Meeting. (MKS) 196C 890921 Removed TeX from comments. (FNF) 197C 890922 Numerous changes to prologue to make closer to SLATEC 198C standard. (FNF) 199C 890929 Numerous changes to reduce SP/DP differences. (FNF) 200C 910411 Prologue converted to Version 4.0 format. (BAB) 201C 920407 COMMON BLOCK renamed DSLBLK. (WRB) 202C 920511 Added complete declaration section. (WRB) 203C 921113 Corrected C***CATEGORY line. (FNF) 204C***END PROLOGUE DSDOMN 205C .. Parameters .. 206 INTEGER LOCRB, LOCIB 207 PARAMETER (LOCRB=1, LOCIB=11) 208C .. Scalar Arguments .. 209 DOUBLE PRECISION ERR, TOL 210 INTEGER IERR, ISYM, ITER, ITMAX, ITOL, IUNIT, LENIW, LENW, N, 211 + NELT, NSAVE 212C .. Array Arguments .. 213 DOUBLE PRECISION A(N), B(N), RWORK(LENW), X(N) 214 INTEGER IA(NELT), IWORK(LENIW), JA(NELT) 215C .. Local Scalars .. 216 INTEGER LOCAP, LOCCSA, LOCDIN, LOCDZ, LOCEMA, LOCIW, LOCP, LOCR, 217 + LOCW, LOCZ 218C .. External Subroutines .. 219 EXTERNAL DCHKW, DOMN, DS2Y, DSDI, DSDS, DSMV 220C***FIRST EXECUTABLE STATEMENT DSDOMN 221C 222 IERR = 0 223 IF( N.LT.1 .OR. NELT.LT.1 ) THEN 224 IERR = 3 225 RETURN 226 ENDIF 227C 228C Change the SLAP input matrix IA, JA, A to SLAP-Column format. 229 CALL DS2Y( N, NELT, IA, JA, A, ISYM ) 230C 231C Set up the workspace. 232 LOCIW = LOCIB 233C 234 LOCDIN = LOCRB 235 LOCR = LOCDIN + N 236 LOCZ = LOCR + N 237 LOCP = LOCZ + N 238 LOCAP = LOCP + N*(NSAVE+1) 239 LOCEMA = LOCAP + N*(NSAVE+1) 240 LOCDZ = LOCEMA + N*(NSAVE+1) 241 LOCCSA = LOCDZ + N 242 LOCW = LOCCSA + NSAVE 243C 244C Check the workspace allocations. 245 CALL DCHKW( 'DSDOMN', LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR ) 246 IF( IERR.NE.0 ) RETURN 247C 248 IWORK(4) = LOCDIN 249 IWORK(9) = LOCIW 250 IWORK(10) = LOCW 251C 252C Compute the inverse of the diagonal of the matrix. 253 CALL DSDS(N, NELT, IA, JA, A, ISYM, RWORK(LOCDIN)) 254C 255C Perform the Diagonally Scaled Orthomin iteration algorithm. 256 CALL DOMN(N, B, X, NELT, IA, JA, A, ISYM, DSMV, 257 $ DSDI, NSAVE, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, 258 $ RWORK(LOCR), RWORK(LOCZ), RWORK(LOCP), RWORK(LOCAP), 259 $ RWORK(LOCEMA), RWORK(LOCDZ), RWORK(LOCCSA), 260 $ RWORK, IWORK ) 261 RETURN 262C------------- LAST LINE OF DSDOMN FOLLOWS ---------------------------- 263 END 264