The TELEMAC-MASCARET system  trunk
sd_cdrv.f
Go to the documentation of this file.
1 ! ******************
2  SUBROUTINE sd_cdrv
3 ! ******************
4 !
5  &(n, r,c,ic, ia,ja,a, b, z, nsp,isp,rsp,esp, path, flag)
6 !
7 !***********************************************************************
8 ! BIEF V6P3 30/06/2013
9 !***********************************************************************
10 !
11 !brief DRIVER FOR SUBROUTINES TO SOLVE SPARSE NONSYMMETRICAL
12 !+ SYSTEMS OF LINEAR EQUATIONS (UNCOMPRESSED POINTER STORAGE).
13 !code
14 !+ PARAMETERS
15 !+ CLASS ABBREVIATIONS ARE --
16 !+ N - INTEGER VARIABLE
17 !+ F - REAL VARIABLE
18 !+ V - SUPPLIES A VALUE TO THE DRIVER
19 !+ R - RETURNS A RESULT FROM THE DRIVER
20 !+ I - USED INTERNALLY BY THE DRIVER
21 !+ A - ARRAY
22 !+
23 !+ CLASS \ PARAMETER
24 !+ ------+----------
25 !+ \
26 !+ THE NONZERO ENTRIES OF THE COEFFICIENT MATRIX M ARE STORED
27 !+ ROW-BY-ROW IN THE ARRAY A. TO IDENTIFY THE INDIVIDUAL NONZERO
28 !+ ENTRIES IN EACH ROW, WE NEED TO KNOW IN WHICH COLUMN EACH ENTRY
29 !+ LIES. THE COLUMN INDICES WHICH CORRESPOND TO THE NONZERO ENTRIES
30 !+ OF M ARE STORED IN THE ARRAY JA; I.E., IF A(K) = M(I,J), THEN
31 !+ JA(K) = J. IN ADDITION, WE NEED TO KNOW WHERE EACH ROW STARTS AND
32 !+ HOW LONG IT IS. THE INDEX POSITIONS IN JA AND A WHERE THE ROWS OF
33 !+ M BEGIN ARE STORED IN THE ARRAY IA; I.E., IF M(I,J) IS THE FIRST
34 !+ NONZERO ENTRY (STORED) IN THE I-TH ROW AND A(K) = M(I,J), THEN
35 !+ IA(I) = K. MOREOVER, THE INDEX IN JA AND A OF THE FIRST LOCATION
36 !+ FOLLOWING THE LAST ELEMENT IN THE LAST ROW IS STORED IN IA(N+1).
37 !+ THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS GIVEN BY
38 !+ IA(I+1) - IA(I), THE NONZERO ENTRIES OF THE I-TH ROW ARE STORED
39 !+ CONSECUTIVELY IN
40 !+ A(IA(I)), A(IA(I)+1), ..., A(IA(I+1)-1),
41 !+ AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN
42 !+ JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1).
43 !+ FOR EXAMPLE, THE 5 BY 5 MATRIX
44 !+ ( 1. 0. 2. 0. 0.)
45 !+ ( 0. 3. 0. 0. 0.)
46 !+ M = ( 0. 4. 5. 6. 0.)
47 !+ ( 0. 0. 0. 7. 0.)
48 !+ ( 0. 0. 0. 8. 9.)
49 !+ WOULD BE STORED AS
50 !+ \ 1 2 3 4 5 6 7 8 9
51 !+ ---+--------------------------
52 !+ IA \ 1 3 4 7 8 10
53 !+ JA \ 1 3 2 2 3 4 4 4 5
54 !+ A \ 1. 2. 3. 4. 5. 6. 7. 8. 9. .
55 !+
56 !+ NV \ N - NUMBER OF VARIABLES/EQUATIONS.
57 !+ FVA \ A - NONZERO ENTRIES OF THE COEFFICIENT MATRIX M, STORED
58 !+ \ BY ROWS.
59 !+ \ SIZE = NUMBER OF NONZERO ENTRIES IN M.
60 !+ NVA \ IA - POINTERS TO DELIMIT THE ROWS IN A.
61 !+ \ SIZE = N+1.
62 !+ NVA \ JA - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF A.
63 !+ \ SIZE = SIZE OF A.
64 !+ FVA \ B - RIGHT-HAND SIDE B; B AND Z CAN THE SAME ARRAY.
65 !+ \ SIZE = N.
66 !+ FRA \ Z - SOLUTION X; B AND Z CAN BE THE SAME ARRAY.
67 !+ \ SIZE = N.
68 !+
69 !+ THE ROWS AND COLUMNS OF THE ORIGINAL MATRIX M CAN BE
70 !+ REORDERED (E.G., TO REDUCE FILLIN OR ENSURE NUMERICAL STABILITY)
71 !+ BEFORE CALLING THE DRIVER. IF NO REORDERING IS DONE, THEN SET
72 !+ R(I) = C(I) = IC(I) = I FOR I=1,...,N. THE SOLUTION Z IS RETURNED
73 !+ IN THE ORIGINAL ORDER.
74 !+
75 !+ NVA \ R - ORDERING OF THE ROWS OF M.
76 !+ \ SIZE = N.
77 !+ NVA \ C - ORDERING OF THE COLUMNS OF M.
78 !+ \ SIZE = N.
79 !+ NVA \ IC - INVERSE OF THE ORDERING OF THE COLUMNS OF M; I.E.,
80 !+ \ IC(C(I)) = I FOR I=1,...,N.
81 !+ \ SIZE = N.
82 !+
83 !+ THE SOLUTION OF THE SYSTEM OF LINEAR EQUATIONS IS DIVIDED INTO
84 !+ THREE STAGES --
85 !+ NSF -- THE MATRIX M IS PROCESSED SYMBOLICALLY TO DETERMINE WHERE
86 !+ FILLIN WILL OCCUR DURING THE NUMERIC FACTORIZATION.
87 !+ NNF -- THE MATRIX M IS FACTORED NUMERICALLY INTO THE PRODUCT LDU
88 !+ OF A UNIT LOWER TRIANGULAR MATRIX L, A DIAGONAL MATRIX D,
89 !+ AND A UNIT UPPER TRIANGULAR MATRIX U, AND THE SYSTEM
90 !+ MX = B IS SOLVED.
91 !+ NNS -- THE LINEAR SYSTEM MX = B IS SOLVED USING THE LDU
92 !+ OR FACTORIZATION FROM NNF.
93 !+ NNT -- THE TRANSPOSED LINEAR SYSTEM MT X = B IS SOLVED USING
94 !+ THE LDU FACTORIZATION FROM NNF.
95 !+ FOR SEVERAL SYSTEMS WHOSE COEFFICIENT MATRICES HAVE THE SAME
96 !+ NONZERO STRUCTURE, NSF NEED BE DONE ONLY ONCE (FOR THE FIRST
97 !+ SYSTEM); THEN NNF IS DONE ONCE FOR EACH ADDITIONAL SYSTEM. FOR
98 !+ SEVERAL SYSTEMS WITH THE SAME COEFFICIENT MATRIX, NSF AND NNF NEED
99 !+ BE DONE ONLY ONCE (FOR THE FIRST SYSTEM); THEN NNS OR NNT IS DONE
100 !+ ONCE FOR EACH ADDITIONAL RIGHT-HAND SIDE.
101 !+
102 !+ NV \ PATH - PATH SPECIFICATION; VALUES AND THEIR MEANINGS ARE --
103 !+ \ 1 PERFORM NSF AND NNF.
104 !+ \ 2 PERFORM NNF ONLY (NSF IS ASSUMED TO HAVE BEEN
105 !+ \ DONE IN A MANNER COMPATIBLE WITH THE STORAGE
106 !+ \ ALLOCATION USED IN THE DRIVER).
107 !+ \ 3 PERFORM NNS ONLY (NSF AND NNF ARE ASSUMED TO
108 !+ \ HAVE BEEN DONE IN A MANNER COMPATIBLE WITH THE
109 !+ \ STORAGE ALLOCATION USED IN THE DRIVER).
110 !+ \ 4 PERFORM NNT ONLY (NSF AND NNF ARE ASSUMED TO
111 !+ \ HAVE BEEN DONE IN A MANNER COMPATIBLE WITH THE
112 !+ \ STORAGE ALLOCATION USED IN THE DRIVER).
113 !+ \ 5 PERFORM NSF ONLY.
114 !+
115 !+ VARIOUS ERRORS ARE DETECTED BY THE DRIVER AND THE INDIVIDUAL
116 !+ SUBROUTINES.
117 !+
118 !+ NR \ FLAG - ERROR FLAG; VALUES AND THEIR MEANINGS ARE --
119 !+ \ 0 NO ERRORS DETECTED
120 !+ \ N+K NULL ROW IN A -- ROW = K
121 !+ \ 2N+K DUPLICATE ENTRY IN A -- ROW = K
122 !+ \ 3N+K INSUFFICIENT STORAGE IN NSF -- ROW = K
123 !+ \ 4N+1 INSUFFICIENT STORAGE IN NNF
124 !+ \ 5N+K NULL PIVOT -- ROW = K
125 !+ \ 6N+K INSUFFICIENT STORAGE IN NSF -- ROW = K
126 !+ \ 7N+1 INSUFFICIENT STORAGE IN NNF
127 !+ \ 8N+K ZERO PIVOT -- ROW = K
128 !+ \ 10N+1 INSUFFICIENT STORAGE IN NDRV
129 !+ \ 11N+1 ILLEGAL PATH SPECIFICATION
130 !+
131 !+ WORKING STORAGE IS NEEDED FOR THE FACTORED FORM OF THE MATRIX
132 !+ M PLUS VARIOUS TEMPORARY VECTORS. THE ARRAYS ISP AND RSP SHOULD BE
133 !+ EQUIVALENCED; INTEGER STORAGE IS ALLOCATED FROM THE BEGINNING OF
134 !+ ISP AND REAL STORAGE FROM THE END OF RSP.
135 !+
136 !+ NV \ NSP - DECLARED DIMENSION OF RSP; NSP GENERALLY MUST
137 !+ \ BE LARGER THAN 5N+3 + 2K (WHERE K = (NUMBER OF
138 !+ \ NONZERO ENTRIES IN M)).
139 !+ NVIRA \ ISP - INTEGER WORKING STORAGE DIVIDED UP INTO VARIOUS ARRAYS
140 !+ \ NEEDED BY THE SUBROUTINES; ISP AND RSP SHOULD BE
141 !+ \ EQUIVALENCED.
142 !+ \ SIZE = LRATIO*NSP
143 !+ FVIRA \ RSP - REAL WORKING STORAGE DIVIDED UP INTO VARIOUS ARRAYS
144 !+ \ NEEDED BY THE SUBROUTINES; ISP AND RSP SHOULD BE
145 !+ \ EQUIVALENCED.
146 !+ \ SIZE = NSP.
147 !+ NR \ ESP - IF SUFFICIENT STORAGE WAS AVAILABLE TO PERFORM THE
148 !+ \ SYMBOLIC FACTORIZATION (NSF), THEN ESP IS SET TO THE
149 !+ \ AMOUNT OF EXCESS STORAGE PROVIDED (NEGATIVE IF
150 !+ \ INSUFFICIENT STORAGE WAS AVAILABLE TO PERFORM THE
151 !+ \ NUMERIC FACTORIZATION (NNF)).
152 !
153 !note IMPORTANT : INSPIRED FROM PACKAGE CMLIB3 - YALE UNIVERSITE-YSMP
154 !note
155 !+
156 !+ TO CONVERT THESE ROUTINES FOR DOUBLE PRECISION ARRAYS, SIMPLY USE
157 !+
158 !+ THE DOUBLE PRECISION DECLARATIONS IN PLACE OF THE REAL DECLARATIONS
159 !+
160 !+ IN EACH SUBPROGRAM; IN ADDITION, THE DATA VALUE OF THE INTEGER
161 !+
162 !+ VARIABLE LRATIO MUST BE SET AS INDICATED IN SUBROUTINE NDRV.
163 !
164 !history E. RAZAFINDRAKOTO (LNH)
165 !+ 18/02/08
166 !+ V5P9
167 !+
168 !
169 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
170 !+ 13/07/2010
171 !+ V6P0
172 !+ Translation of French comments within the FORTRAN sources into
173 !+ English comments
174 !
175 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
176 !+ 21/08/2010
177 !+ V6P0
178 !+ Creation of DOXYGEN tags for automated documentation and
179 !+ cross-referencing of the FORTRAN sources
180 !
181 !history J-M HERVOUET (LNHE)
182 !+ 10/08/2012
183 !+ V6P2
184 !+ LRATIO set to 1 instead of 2 and suppressed.
185 !
186 !history C.PEYRARD (LNHE)
187 !+ 30/06/2013
188 !+ V6P3
189 !+ SD_NDRV => SD_CDRV
190 !
191 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
192 !| A |-->|NONZERO ENTRIES OF THE COEFFICIENT MATRIX M,
193 !| | |STORED BY ROWS
194 !| B |-->|RIGHT-HAND SIDE B ;
195 !| C |-->|ORDERING OF THE COLUMNS OF MATRIX
196 !| ESP |---|INTEGER DIMENSION STORAGE : IF SUFFICIENT STORAGE
197 !| | |WAS AVAILABLE TO PERFORM THE SYMBOLIC
198 !| | |FACTORIZATION (NSF), THEN ESP IS SET TO THE AMOUNT
199 !| | | OF EXCESS STORAGE PROVIDED (NEGATIVE IF
200 !| | |INSUFFICIENT STORAGE WAS AVAILABLE TO PERFORM
201 !| | |THE NUMERIC FACTORIZATION (NNF)).
202 !| FLAG |<--|ERROR FLAG; VALUES AND THEIR MEANINGS ARE :
203 !| | |0 NO ERRORS DETECTED
204 !| | |N+K NULL ROW IN A -- ROW = K
205 !| | |2N+K DUPLICATE ENTRY IN A -- ROW = K
206 !| | |3N+K INSUFFICIENT STORAGE IN NSF -- ROW = K
207 !| | |4N+1 INSUFFICIENT STORAGE IN NNF
208 !| | |5N+K NULL PIVOT -- ROW = K
209 !| | |6N+K INSUFFICIENT STORAGE IN NSF -- ROW = K
210 !| | |7N+1 INSUFFICIENT STORAGE IN NNF
211 !| | |8N+K ZERO PIVOT -- ROW = K
212 !| | |10N+1 INSUFFICIENT STORAGE IN NDRV
213 !| | |11N+1 ILLEGAL PATH SPECIFICATION
214 !| IA |-->|POINTERS TO DELIMIT THE ROWS IN A ; SIZE = N+1
215 !| IC |-->|INVERSE OF THE ORDERING OF THE COLUMNS OF MATRIX
216 !| ISP |-->|INTEGER WORKING STORAGE
217 !| JA |-->|COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF A
218 !| N |-->|NUMBER OF VARIABLES/EQUATIONS
219 !| NSP |-->|DECLARED DIMENSION OF RSP : GENERALLY MUST
220 !| | |BE LARGER THAN 5N+3 + 2K (WHERE K = (NUMBER OF
221 !| | |NONZERO ENTRIES IN THE MATRIX M).
222 !| PATH |-->|PATH SPECIFICATION; VALUES AND THEIR MEANINGS ARE:
223 !| | |1 PERFORM NSF AND NNF.
224 !| | |2 PERFORM NNF ONLY (NSF IS ASSUMED TO HAVE BEEN
225 !| | | DONE IN A MANNER COMPATIBLE WITH THE STORAGE
226 !| | | ALLOCATION USED IN THE DRIVER).
227 !| | |3 PERFORM NNS ONLY (NSF AND NNF ARE ASSUMED TO
228 !| | | HAVE BEEN DONE IN A MANNER COMPATIBLE WITH THE
229 !| | | STORAGE ALLOCATION USED IN THE DRIVER).
230 !| | |4 PERFORM NNT ONLY (NSF AND NNF ARE ASSUMED TO
231 !| | | HAVE BEEN DONE IN A MANNER COMPATIBLE WITH THE
232 !| | | STORAGE ALLOCATION USED IN THE DRIVER).
233 !| | |5 PERFORM NSF ONLY.
234 !| R |-->|ORDERING OF THE ROWS OF MATRIX
235 !| RSP |-->|REAL WORKING STORAGE (SIZE = NSP)
236 !| Z |<--|SOLUTION X
237 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
238 !
239  USE bief, ex_sd_cdrv => sd_cdrv
240 !
242  IMPLICIT NONE
243 !
244 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
245 !
246  INTEGER, INTENT(IN) :: N,NSP,PATH,C(n)
247  INTEGER, INTENT(INOUT) :: R(n),FLAG,IA(*),JA(*),ISP(*),ESP,IC(*)
248  DOUBLE PRECISION, INTENT(IN) :: B(n)
249  DOUBLE PRECISION, INTENT(INOUT) :: A(*),Z(n),RSP(nsp)
250 !
251 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
252 !
253  INTEGER D, U, Q, ROW, TMP, AR, UMAX
254  INTEGER IL,IJL,IU,IJU,IRL,JRL,JL,KMAX,JLMAX,IRA,JRA,IRAC,IRU,JRU
255  INTEGER JUTMP , JUMAX,JU
256  INTEGER I,J ,L,LMAX , LRATIO
257 
258 ! SET LRATIO EQUAL TO THE RATIO BETWEEN THE LENGTH OF FLOATING POINT
259 ! AND INTEGER ARRAY DATA. E. G., LRATIO = 1 FOR (REAL, INTEGER),
260 ! LRATIO = 2 FOR (DOUBLE PRECISION, INTEGER)
261 !
262 ! DATA LRATIO/2/
263  lratio=1
264  IF (path.LT.1 .OR. 5.LT.path) GO TO 111
265 !******INITIALIZE AND DIVIDE UP TEMPORARY STORAGE *******************
266  il = 1
267  ijl = il + (n+1)
268  iu = ijl + n
269  iju = iu + (n+1)
270  irl = iju + n
271  jrl = irl + n
272  jl = jrl + n
273 !
274 ! ****** REORDER A IF NECESSARY, CALL NSFC IF FLAG IS SET ***********
275  IF ((path-1) * (path-5) .NE. 0) GO TO 5
276  kmax = (lratio*nsp + 1 - jl) - (n+1) - 5*n
277  jlmax = kmax/2
278  q = jl + jlmax
279  ira = q + (n+1)
280  jra = ira + n
281  irac = jra + n
282  iru = irac + n
283  jru = iru + n
284  jutmp = jru + n
285  jumax = lratio*nsp + 1 - jutmp
286  esp = kmax/lratio
287  IF (jlmax.LE.0 .OR. jumax.LE.0) GO TO 110
288  DO 1 i=1,n
289  IF (c(i).NE.i) GO TO 2
290  1 CONTINUE
291  GO TO 3
292  2 ar = nsp + 1 - n
293 !CP WRITE(LU,*) 'CCP : ON APPELLE SD_NROC_CP'
294  CALL sd_nroc
295  & (n, ic, ia,ja,a, isp(il), rsp(ar), isp(iu), flag)
296  IF (flag.NE.0) GO TO 100
297 !
298 !CP WRITE(LU,*) 'CCP : ON APPELLE SD_NSFC_CP'
299  3 CALL sd_nsfc
300  & (n, r, ic, ia,ja,
301  & jlmax, isp(il), isp(jl), isp(ijl),
302  & jumax, isp(iu), isp(jutmp), isp(iju),
303  & isp(q), isp(ira), isp(jra), isp(irac),
304  & isp(irl), isp(jrl), isp(iru), isp(jru), flag)
305  IF(flag .NE. 0) GO TO 100
306 ! ****** MOVE JU NEXT TO JL *****************************************
307  jlmax = isp(ijl+n-1)
308  ju = jl + jlmax
309  jumax = isp(iju+n-1)
310  IF (jumax.LE.0) GO TO 5
311  DO 4 j=1,jumax
312  4 isp(ju+j-1) = isp(jutmp+j-1)
313 !
314 ! ****** CALL REMAINING SUBROUTINES *********************************
315  5 jlmax = isp(ijl+n-1)
316  ju = jl + jlmax
317  jumax = isp(iju+n-1)
318  l = (ju + jumax - 2 + lratio) / lratio + 1
319  lmax = isp(il+n) - 1
320  d = l + lmax
321  u = d + n
322  row = nsp + 1 - n
323  tmp = row - n
324  umax = tmp - u
325  esp = umax - (isp(iu+n) - 1)
326 !
327  IF ((path-1) * (path-2) .NE. 0) GO TO 6
328  IF (umax.LT.0) GO TO 110
329 !CP WRITE(LU,*) 'CCP : ON APPELLE SD_NNFC_CP'
330  CALL sd_nnfc
331  & (n, r, c, ic, ia, ja, a, z, b,
332  & lmax, isp(il), isp(jl), isp(ijl), rsp(l), rsp(d),
333  & umax, isp(iu), isp(ju), isp(iju), rsp(u),
334  & rsp(row), rsp(tmp), isp(irl), isp(jrl), flag)
335 !CP WRITE(LU,*) 'FLAG=', FLAG
336  IF(flag .NE. 0) GO TO 100
337 !
338  6 IF ((path-3) .NE. 0) GO TO 7
339 !CP WRITE(LU,*) 'CCP : ON APPELLE SD_NNSC_CP'
340  CALL sd_nnsc
341  & (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l),
342  & rsp(d), isp(iu), isp(ju), isp(iju), rsp(u),
343  & z, b, rsp(tmp))
344 !
345  7 IF ((path-4) .NE. 0) GO TO 8
346 !CP WRITE(LU,*) 'CCP : ON APPELLE SD_NNTC_CP'
347  CALL sd_nntc
348  & (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l),
349  & rsp(d), isp(iu), isp(ju), isp(iju), rsp(u),
350  & z, b, rsp(tmp))
351 
352  8 RETURN
353 !
354 ! ** ERROR.. ERROR DETECTED IN NROC, NSFC, NNFC, OR NNSC
355  100 RETURN
356 ! ** ERROR.. INSUFFICIENT STORAGE
357  110 flag = 10*n + 1
358  RETURN
359 ! ** ERROR.. ILLEGAL PATH SPECIFICATION
360  111 flag = 11*n + 1
361  RETURN
362  END
subroutine sd_nnfc(N, R, C, IC, IA, JA, A, Z, B, LMAX, IL, JL, IJL, L, D, UMAX, IU, JU, IJU, U, ROW, TMP, IRL, JRL, FLAG)
Definition: sd_nnfc.f:8
subroutine sd_nnsc(N, R, C, IL, JL, IJL, L, D, IU, JU, IJU, U, Z, B, TMP)
Definition: sd_nnsc.f:7
subroutine sd_nntc(N, R, C, IL, JL, IJL, L, D, IU, JU, IJU, U, Z, B, TMP)
Definition: sd_nntc.f:7
subroutine sd_cdrv(N, R, C, IC, IA, JA, A, B, Z, NSP, ISP, RSP, ESP, PATH, FLAG)
Definition: sd_cdrv.f:7
subroutine sd_nsfc(N, R, IC, IA, JA, JLMAX, IL, JL, IJL, JUMAX, IU, JU, IJU, Q, IRA, JRA, IRAC, IRL, JRL, IRU, JRU, FLAG)
Definition: sd_nsfc.f:8
subroutine sd_nroc(N, IC, IA, JA, A, JAR, AR, P, FLAG)
Definition: sd_nroc.f:7
Definition: bief.f:3