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