The TELEMAC-MASCARET system  trunk
sd_nsfc.f
Go to the documentation of this file.
1 ! ******************
2  SUBROUTINE sd_nsfc
3 ! ******************
4 !
5  &(n,r,ic,ia,ja,jlmax,il,jl,ijl,jumax,iu,ju,iju,q,ira,jra,irac,
6  & irl,jrl,iru,jru,flag)
7 !
8 !***********************************************************************
9 ! BIEF V6P3 30/06/2013
10 !***********************************************************************
11 !
12 !brief SYMBOLIC LDU-FACTORIZATION OF NONSYMMETRIC SPARSE MATRIX
13 !+ (COMPRESSED POINTER STORAGE)
14 !code
15 !+ INPUT VARIABLES.. N, R, IC, IA, JA, JLMAX, JUMAX.
16 !+ OUTPUT VARIABLES.. IL, JL, IJL, IU, JU, IJU, FLAG.
17 !+
18 !+ PARAMETERS USED INTERNALLY..
19 !+ NIA \ Q - SUPPOSE M* IS THE RESULT OF REORDERING M. IF
20 !+ \ PROCESSING OF THE ITH ROW OF M* (HENCE THE ITH
21 !+ \ ROW OF U) IS BEING DONE, Q(J) IS INITIALLY
22 !+ \ NONZERO IF M*(I,J) IS NONZERO (J.GE.I). SINCE
23 !+ \ VALUES NEED NOT BE STORED, EACH ENTRY POINTS TO THE
24 !+ \ NEXT NONZERO AND Q(N+1) POINTS TO THE FIRST. N+1
25 !+ \ INDICATES THE END OF THE LIST. FOR EXAMPLE, IF N=9
26 !+ \ AND THE 5TH ROW OF M* IS
27 !+ \ 0 X X 0 X 0 0 X 0
28 !+ \ THEN Q WILL INITIALLY BE
29 !+ \ A A A A 8 A A 10 5 (A - ARBITRARY).
30 !+ \ AS THE ALGORITHM PROCEEDS, OTHER ELEMENTS OF Q
31 !+ \ ARE INSERTED IN THE LIST BECAUSE OF FILLIN.
32 !+ \ Q IS USED IN AN ANALOGOUS MANNER TO COMPUTE THE
33 !+ \ ITH COLUMN OF L.
34 !+ \ SIZE = N+1.
35 !+ NIA \ IRA, - VECTORS USED TO FIND THE COLUMNS OF M. AT THE KTH
36 !+ NIA \ JRA, STEP OF THE FACTORIZATION, IRAC(K) POINTS TO THE
37 !+ NIA \ IRAC HEAD OF A LINKED LIST IN JRA OF ROW INDICES I
38 !+ \ SUCH THAT I .GE. K AND M(I,K) IS NONZERO. ZERO
39 !+ \ INDICATES THE END OF THE LIST. IRA(I) (I.GE.K)
40 !+ \ POINTS TO THE SMALLEST J SUCH THAT J .GE. K AND
41 !+ \ M(I,J) IS NONZERO.
42 !+ \ SIZE OF EACH = N.
43 !+ NIA \ IRL, - VECTORS USED TO FIND THE ROWS OF L. AT THE KTH STEP
44 !+ NIA \ JRL OF THE FACTORIZATION, JRL(K) POINTS TO THE HEAD
45 !+ \ OF A LINKED LIST IN JRL OF COLUMN INDICES J
46 !+ \ SUCH J .LT. K AND L(K,J) IS NONZERO. ZERO
47 !+ \ INDICATES THE END OF THE LIST. IRL(J) (J.LT.K)
48 !+ \ POINTS TO THE SMALLEST I SUCH THAT I .GE. K AND
49 !+ \ L(I,J) IS NONZERO.
50 !+ \ SIZE OF EACH = N.
51 !+ NIA \ IRU, - VECTORS USED IN A MANNER ANALOGOUS TO IRL AND JRL
52 !+ NIA \ JRU TO FIND THE COLUMNS OF U.
53 !+ \ SIZE OF EACH = N.
54 !+
55 !+ INTERNAL VARIABLES..
56 !+ JLPTR - POINTS TO THE LAST POSITION USED IN JL.
57 !+ JUPTR - POINTS TO THE LAST POSITION USED IN JU.
58 !+ JMIN,JMAX - ARE THE INDICES IN A OR U OF THE FIRST AND LAST
59 !+ ELEMENTS TO BE EXAMINED IN A GIVEN ROW.
60 !+ FOR EXAMPLE, JMIN=IA(K), JMAX=IA(K+1)-1.
61 !+
62 !note IMPORTANT : INSPIRED FROM PACKAGE CMLIB3 - YALE UNIVERSITE-YSMP
63 ! DON'T HESITATE TO CHANGE IN/OUTPUT VARIABLES COMMENTS
64 ! FOR CLARITY
65 !
66 !history C. PEYRARD (LNHE)
67 !+ 30/01/13
68 !+ V6P3
69 !+
70 !
71 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
72 !| FLAG |<--| INDICATOR ERROR :
73 !| | |= N + R(K): NULL ROW IN A
74 !| | |= 2*N + R(K):DUPLICATE ENTRY IN A
75 !| | |= 3*N + K:INSUFFICIENT STORAGE FOR JL
76 !| | |= 6*N + K:INSUFFICIENT STORAGE FOR JU
77 !| | |= 5*N + K: ZERO PIVOT
78 !| IA, JA |-->| STRUCTURE OF A NONSYMMETRICAL MATRIX
79 !| IC |-->| INVERSE OF THE ORDERING OF THE COLUMNS OF MATRIX
80 !| IL, JL |<--| STRUCTURE OF LOWER FACTORISED TRIANGULAR MATRIX
81 !| IU, JU |<--| STRUCTURE OF UPPER FACTORISED TRIANGULAR MATRIX
82 !| IJU,IJL |-->| USED TO COMPRESS STORAGE OF JU and JL
83 !| IM |---| INTEGER ONE-DIMENSIONAL WORK ARRAY;SIZE = N
84 !| JLMAX |-->| PREVISIONAL MAXIMUM DIMENSION OF JL
85 !| JUMAX |-->| PREVISIONAL MAXIMUM DIMENSION OF JU
86 !| N |-->| RANK OF MATRIX
87 !| Q |---| INTEGER ONE-DIMENSIONAL WORK ARRAY;SIZE = N+1
88 !| R |-->| ORDERING OF THE ROWS OF MATRIX
89 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
90 !
91  USE declarations_special
92  IMPLICIT NONE
93 !
94 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
95 !
96  INTEGER IA(*), JA(*), IRA(*), JRA(*), IL(*), JL(*), IJL(*)
97  INTEGER IU(*), JU(*), IJU(*), IRL(*), JRL(*), IRU(*), JRU(*)
98  INTEGER R(*), IC(*), Q(*), IRAC(*),FLAG
99 !
100 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
101 !
102  INTEGER CEND, QM, REND, RK, VJ
103  INTEGER NP1,N,JLMIN,JLPTR,JUMIN,JUPTR,K,IAK,JAIAK,LUK,M
104  INTEGER LASTID,LASTI,I,JMIN,JMAX,LONG,JTMP,J,I1,IRLL,JLMAX
105  INTEGER JUMAX,IRAI,JAIRAI,IRUL
106 !
107 ! ****** INITIALIZE POINTERS ****************************************
108  np1 = n + 1
109  jlmin = 1
110  jlptr = 0
111  il(1) = 1
112  jumin = 1
113  juptr = 0
114  iu(1) = 1
115  DO 1 k=1,n
116  irac(k) = 0
117  jra(k) = 0
118  jrl(k) = 0
119  1 jru(k) = 0
120 ! ****** INITIALIZE COLUMN POINTERS FOR A ***************************
121  DO 2 k=1,n
122  rk = r(k)
123  iak = ia(rk)
124  IF (iak .GE. ia(rk+1)) GO TO 101
125  jaiak = ic(ja(iak))
126  IF (jaiak .GT. k) GO TO 105
127  jra(k) = irac(jaiak)
128  irac(jaiak) = k
129  2 ira(k) = iak
130 !
131 ! ****** FOR EACH COLUMN OF L AND ROW OF U **************************
132  DO 41 k=1,n
133 !
134 ! ****** INITIALIZE Q FOR COMPUTING KTH COLUMN OF L *****************
135  q(np1) = np1
136  luk = -1
137 ! ****** BY FILLING IN KTH COLUMN OF A ******************************
138  vj = irac(k)
139  IF (vj .EQ. 0) GO TO 5
140  3 qm = np1
141  4 m = qm
142  qm = q(m)
143  IF (qm .LT. vj) GO TO 4
144  IF (qm .EQ. vj) GO TO 102
145  luk = luk + 1
146  q(m) = vj
147  q(vj) = qm
148  vj = jra(vj)
149  IF (vj .NE. 0) GO TO 3
150 ! ****** LINK THROUGH JRU *******************************************
151  5 lastid = 0
152  lasti = 0
153  ijl(k) = jlptr
154  i = k
155  6 i = jru(i)
156  IF (i .EQ. 0) GO TO 10
157  qm = np1
158  jmin = irl(i)
159  jmax = ijl(i) + il(i+1) - il(i) - 1
160  long = jmax - jmin
161  IF (long .LT. 0) GO TO 6
162  jtmp = jl(jmin)
163  IF (jtmp .NE. k) long = long + 1
164  IF (jtmp .EQ. k) r(i) = -r(i)
165  IF (lastid .GE. long) GO TO 7
166  lasti = i
167  lastid = long
168 ! ****** AND MERGE THE CORRESPONDING COLUMNS INTO THE KTH COLUMN ****
169  7 DO 9 j=jmin,jmax
170  vj = jl(j)
171  8 m = qm
172  qm = q(m)
173  IF (qm .LT. vj) GO TO 8
174  IF (qm .EQ. vj) GO TO 9
175  luk = luk + 1
176  q(m) = vj
177  q(vj) = qm
178  qm = vj
179  9 CONTINUE
180  GO TO 6
181 ! ****** LASTI IS THE LONGEST COLUMN MERGED INTO THE KTH ************
182 ! ****** SEE IF IT EQUALS THE ENTIRE KTH COLUMN *********************
183  10 qm = q(np1)
184  IF (qm .NE. k) GO TO 105
185  IF (luk .EQ. 0) GO TO 17
186  IF (lastid .NE. luk) GO TO 11
187 ! ****** IF SO, JL CAN BE COMPRESSED ********************************
188  irll = irl(lasti)
189  ijl(k) = irll + 1
190  IF (jl(irll) .NE. k) ijl(k) = ijl(k) - 1
191  GO TO 17
192 ! ****** IF NOT, SEE IF KTH COLUMN CAN OVERLAP THE PREVIOUS ONE *****
193  11 IF (jlmin .GT. jlptr) GO TO 15
194  qm = q(qm)
195  DO 12 j=jlmin,jlptr
196  IF (jl(j) < qm ) GOTO 12
197  IF (jl(j) == qm ) GOTO 13
198  IF (jl(j) > qm ) GOTO 15
199  12 CONTINUE
200  GO TO 15
201  13 ijl(k) = j
202  DO 14 i=j,jlptr
203  IF (jl(i) .NE. qm) GO TO 15
204  qm = q(qm)
205  IF (qm .GT. n) GO TO 17
206  14 CONTINUE
207  jlptr = j - 1
208 ! ****** MOVE COLUMN INDICES FROM Q TO JL, UPDATE VECTORS ***********
209  15 jlmin = jlptr + 1
210  ijl(k) = jlmin
211  IF (luk .EQ. 0) GO TO 17
212  jlptr = jlptr + luk
213  IF (jlptr .GT. jlmax) GO TO 103
214  qm = q(np1)
215  DO 16 j=jlmin,jlptr
216  qm = q(qm)
217  16 jl(j) = qm
218  17 irl(k) = ijl(k)
219  il(k+1) = il(k) + luk
220 !
221 ! ****** INITIALIZE Q FOR COMPUTING KTH ROW OF U ********************
222  q(np1) = np1
223  luk = -1
224 ! ****** BY FILLING IN KTH ROW OF REORDERED A ***********************
225  rk = r(k)
226  jmin = ira(k)
227  jmax = ia(rk+1) - 1
228  IF (jmin .GT. jmax) GO TO 20
229  DO 19 j=jmin,jmax
230  vj = ic(ja(j))
231  qm = np1
232  18 m = qm
233  qm = q(m)
234  IF (qm .LT. vj) GO TO 18
235  IF (qm .EQ. vj) GO TO 102
236  luk = luk + 1
237  q(m) = vj
238  q(vj) = qm
239  19 CONTINUE
240 ! ****** LINK THROUGH JRL, ******************************************
241  20 lastid = 0
242  lasti = 0
243  iju(k) = juptr
244  i = k
245  i1 = jrl(k)
246  21 i = i1
247  IF (i .EQ. 0) GO TO 26
248  i1 = jrl(i)
249  qm = np1
250  jmin = iru(i)
251  jmax = iju(i) + iu(i+1) - iu(i) - 1
252  long = jmax - jmin
253  IF (long .LT. 0) GO TO 21
254  jtmp = ju(jmin)
255  IF (jtmp .EQ. k) GO TO 22
256 ! ****** UPDATE IRL AND JRL, *****************************************
257  long = long + 1
258  cend = ijl(i) + il(i+1) - il(i)
259  irl(i) = irl(i) + 1
260  IF (irl(i) .GE. cend) GO TO 22
261  j = jl(irl(i))
262  jrl(i) = jrl(j)
263  jrl(j) = i
264  22 IF (lastid .GE. long) GO TO 23
265  lasti = i
266  lastid = long
267 ! ****** AND MERGE THE CORRESPONDING ROWS INTO THE KTH ROW **********
268  23 DO 25 j=jmin,jmax
269  vj = ju(j)
270  24 m = qm
271  qm = q(m)
272  IF (qm .LT. vj) GO TO 24
273  IF (qm .EQ. vj) GO TO 25
274  luk = luk + 1
275  q(m) = vj
276  q(vj) = qm
277  qm = vj
278  25 CONTINUE
279  GO TO 21
280 ! ****** UPDATE JRL(K) AND IRL(K) ***********************************
281  26 IF (il(k+1) .LE. il(k)) GO TO 27
282  j = jl(irl(k))
283  jrl(k) = jrl(j)
284  jrl(j) = k
285 ! ****** LASTI IS THE LONGEST ROW MERGED INTO THE KTH ***************
286 ! ****** SEE IF IT EQUALS THE ENTIRE KTH ROW ************************
287  27 qm = q(np1)
288  IF (qm .NE. k) GO TO 105
289  IF (luk .EQ. 0) GO TO 34
290  IF (lastid .NE. luk) GO TO 28
291 ! ****** IF SO, JU CAN BE COMPRESSED ********************************
292  irul = iru(lasti)
293  iju(k) = irul + 1
294  IF (ju(irul) .NE. k) iju(k) = iju(k) - 1
295  GO TO 34
296 ! ****** IF NOT, SEE IF KTH ROW CAN OVERLAP THE PREVIOUS ONE ********
297  28 IF (jumin .GT. juptr) GO TO 32
298  qm = q(qm)
299  DO 29 j=jumin,juptr
300  IF (ju(j) < qm) GOTO 29
301  IF (ju(j) == qm) GOTO 30
302  IF (ju(j) > qm) GOTO 32
303  29 CONTINUE
304  GO TO 32
305  30 iju(k) = j
306  DO 31 i=j,juptr
307  IF (ju(i) .NE. qm) GO TO 32
308  qm = q(qm)
309  IF (qm .GT. n) GO TO 34
310  31 CONTINUE
311  juptr = j - 1
312 ! ****** MOVE ROW INDICES FROM Q TO JU, UPDATE VECTORS **************
313  32 jumin = juptr + 1
314  iju(k) = jumin
315  IF (luk .EQ. 0) GO TO 34
316  juptr = juptr + luk
317  IF (juptr .GT. jumax) GO TO 106
318  qm = q(np1)
319  DO 33 j=jumin,juptr
320  qm = q(qm)
321  33 ju(j) = qm
322  34 iru(k) = iju(k)
323  iu(k+1) = iu(k) + luk
324 !
325 ! ****** UPDATE IRU, JRU ********************************************
326  i = k
327  35 i1 = jru(i)
328  IF (r(i) .LT. 0) GO TO 36
329  rend = iju(i) + iu(i+1) - iu(i)
330  IF (iru(i) .GE. rend) GO TO 37
331  j = ju(iru(i))
332  jru(i) = jru(j)
333  jru(j) = i
334  GO TO 37
335  36 r(i) = -r(i)
336  37 i = i1
337  IF (i .EQ. 0) GO TO 38
338  iru(i) = iru(i) + 1
339  GO TO 35
340 !
341 ! ****** UPDATE IRA, JRA, IRAC **************************************
342  38 i = irac(k)
343  IF (i .EQ. 0) GO TO 41
344  39 i1 = jra(i)
345  ira(i) = ira(i) + 1
346  IF (ira(i) .GE. ia(r(i)+1)) GO TO 40
347  irai = ira(i)
348  jairai = ic(ja(irai))
349  IF (jairai .GT. i) GO TO 40
350  jra(i) = irac(jairai)
351  irac(jairai) = i
352  40 i = i1
353  IF (i .NE. 0) GO TO 39
354  41 CONTINUE
355 !
356  ijl(n) = jlptr
357  iju(n) = juptr
358  flag = 0
359  RETURN
360 !
361 ! ** ERROR.. NULL ROW IN A
362  101 flag = n + rk
363  RETURN
364 ! ** ERROR.. DUPLICATE ENTRY IN A
365  102 flag = 2*n + rk
366  RETURN
367 ! ** ERROR.. INSUFFICIENT STORAGE FOR JL
368  103 flag = 3*n + k
369  RETURN
370 ! ** ERROR.. NULL PIVOT
371  105 flag = 5*n + k
372  RETURN
373 ! ** ERROR.. INSUFFICIENT STORAGE FOR JU
374  106 flag = 6*n + k
375  RETURN
376  END
declarations_special
Definition: declarations_special.F:3
q
double precision function q(I)
Definition: q.f:7
sd_nsfc
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
Generated on Tue Jan 5 2021 14:37:29 for The TELEMAC-MASCARET system by doxygen 1.8.13