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