1*DECK DLLTI2 2 SUBROUTINE DLLTI2 (N, B, X, NEL, IEL, JEL, EL, DINV) 3C***BEGIN PROLOGUE DLLTI2 4C***PURPOSE SLAP Backsolve routine for LDL' Factorization. 5C Routine to solve a system of the form L*D*L' X = B, 6C where L is a unit lower triangular matrix and D is a 7C diagonal matrix and ' means transpose. 8C***LIBRARY SLATEC (SLAP) 9C***CATEGORY D2E 10C***TYPE DOUBLE PRECISION (SLLTI2-S, DLLTI2-D) 11C***KEYWORDS INCOMPLETE FACTORIZATION, ITERATIVE PRECONDITION, SLAP, 12C SPARSE, SYMMETRIC LINEAR SYSTEM SOLVE 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, NEL, IEL(NEL), JEL(NEL) 23C DOUBLE PRECISION B(N), X(N), EL(NEL), DINV(N) 24C 25C CALL DLLTI2( N, B, X, NEL, IEL, JEL, EL, DINV ) 26C 27C *Arguments: 28C N :IN Integer 29C Order of the Matrix. 30C B :IN Double Precision B(N). 31C Right hand side vector. 32C X :OUT Double Precision X(N). 33C Solution to L*D*L' x = b. 34C NEL :IN Integer. 35C Number of non-zeros in the EL array. 36C IEL :IN Integer IEL(NEL). 37C JEL :IN Integer JEL(NEL). 38C EL :IN Double Precision EL(NEL). 39C IEL, JEL, EL contain the unit lower triangular factor of 40C the incomplete decomposition of the A matrix stored in 41C SLAP Row format. The diagonal of ones *IS* stored. This 42C structure can be set up by the DS2LT routine. See the 43C "Description", below for more details about the SLAP Row 44C format. 45C DINV :IN Double Precision DINV(N). 46C Inverse of the diagonal matrix D. 47C 48C *Description: 49C This routine is supplied with the SLAP package as a routine 50C to perform the MSOLVE operation in the SCG iteration routine 51C for the driver routine DSICCG. It must be called via the 52C SLAP MSOLVE calling sequence convention interface routine 53C DSLLI. 54C **** THIS ROUTINE ITSELF DOES NOT CONFORM TO THE **** 55C **** SLAP MSOLVE CALLING CONVENTION **** 56C 57C IEL, JEL, EL should contain the unit lower triangular factor 58C of the incomplete decomposition of the A matrix stored in 59C SLAP Row format. This IC factorization can be computed by 60C the DSICS routine. The diagonal (which is all one's) is 61C stored. 62C 63C ==================== S L A P Row format ==================== 64C 65C This routine requires that the matrix A be stored in the 66C SLAP Row format. In this format the non-zeros are stored 67C counting across rows (except for the diagonal entry, which 68C must appear first in each "row") and are stored in the 69C double precision array A. In other words, for each row in 70C the matrix put the diagonal entry in A. Then put in the 71C other non-zero elements going across the row (except the 72C diagonal) in order. The JA array holds the column index for 73C each non-zero. The IA array holds the offsets into the JA, 74C A arrays for the beginning of each row. That is, 75C JA(IA(IROW)),A(IA(IROW)) are the first elements of the IROW- 76C th row in JA and A, and JA(IA(IROW+1)-1), A(IA(IROW+1)-1) 77C are the last elements of the IROW-th row. Note that we 78C always have IA(N+1) = NELT+1, where N is the number of rows 79C in the matrix and NELT is the number of non-zeros in the 80C matrix. 81C 82C Here is an example of the SLAP Row storage format for a 5x5 83C Matrix (in the A and JA arrays '|' denotes the end of a row): 84C 85C 5x5 Matrix SLAP Row format for 5x5 matrix on left. 86C 1 2 3 4 5 6 7 8 9 10 11 87C |11 12 0 0 15| A: 11 12 15 | 22 21 | 33 35 | 44 | 55 51 53 88C |21 22 0 0 0| JA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3 89C | 0 0 33 0 35| IA: 1 4 6 8 9 12 90C | 0 0 0 44 0| 91C |51 0 53 0 55| 92C 93C With the SLAP Row format the "inner loop" of this routine 94C should vectorize on machines with hardware support for 95C vector gather/scatter operations. Your compiler may require 96C a compiler directive to convince it that there are no 97C implicit vector dependencies. Compiler directives for the 98C Alliant FX/Fortran and CRI CFT/CFT77 compilers are supplied 99C with the standard SLAP distribution. 100C 101C***SEE ALSO DSICCG, DSICS 102C***REFERENCES (NONE) 103C***ROUTINES CALLED (NONE) 104C***REVISION HISTORY (YYMMDD) 105C 871119 DATE WRITTEN 106C 881213 Previous REVISION DATE 107C 890915 Made changes requested at July 1989 CML Meeting. (MKS) 108C 890922 Numerous changes to prologue to make closer to SLATEC 109C standard. (FNF) 110C 890929 Numerous changes to reduce SP/DP differences. (FNF) 111C 910411 Prologue converted to Version 4.0 format. (BAB) 112C 920511 Added complete declaration section. (WRB) 113C 921113 Corrected C***CATEGORY line. (FNF) 114C 930701 Updated CATEGORY section. (FNF, WRB) 115C***END PROLOGUE DLLTI2 116C .. Scalar Arguments .. 117 INTEGER N, NEL 118C .. Array Arguments .. 119 DOUBLE PRECISION B(N), DINV(N), EL(NEL), X(N) 120 INTEGER IEL(NEL), JEL(NEL) 121C .. Local Scalars .. 122 INTEGER I, IBGN, IEND, IROW 123C***FIRST EXECUTABLE STATEMENT DLLTI2 124C 125C Solve L*y = b, storing result in x. 126C 127 DO 10 I=1,N 128 X(I) = B(I) 129 10 CONTINUE 130 DO 30 IROW = 1, N 131 IBGN = IEL(IROW) + 1 132 IEND = IEL(IROW+1) - 1 133 IF( IBGN.LE.IEND ) THEN 134CLLL. OPTION ASSERT (NOHAZARD) 135CDIR$ IVDEP 136CVD$ NOCONCUR 137CVD$ NODEPCHK 138 DO 20 I = IBGN, IEND 139 X(IROW) = X(IROW) - EL(I)*X(JEL(I)) 140 20 CONTINUE 141 ENDIF 142 30 CONTINUE 143C 144C Solve D*Z = Y, storing result in X. 145C 146 DO 40 I=1,N 147 X(I) = X(I)*DINV(I) 148 40 CONTINUE 149C 150C Solve L-trans*X = Z. 151C 152 DO 60 IROW = N, 2, -1 153 IBGN = IEL(IROW) + 1 154 IEND = IEL(IROW+1) - 1 155 IF( IBGN.LE.IEND ) THEN 156CLLL. OPTION ASSERT (NOHAZARD) 157CDIR$ IVDEP 158CVD$ NOCONCUR 159CVD$ NODEPCHK 160 DO 50 I = IBGN, IEND 161 X(JEL(I)) = X(JEL(I)) - EL(I)*X(IROW) 162 50 CONTINUE 163 ENDIF 164 60 CONTINUE 165C 166 RETURN 167C------------- LAST LINE OF DLLTI2 FOLLOWS ---------------------------- 168 END 169