|
1644
|
1 SUBROUTINE SOLSY (WM, IWM, X, TEM) |
|
|
2 CLLL. OPTIMIZE |
|
|
3 INTEGER IWM |
|
|
4 INTEGER IOWND, IOWNS, |
|
|
5 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, |
|
|
6 2 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
|
|
7 INTEGER I, MEBAND, ML, MU |
|
|
8 DOUBLE PRECISION WM, X, TEM |
|
|
9 DOUBLE PRECISION ROWNS, |
|
|
10 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND |
|
|
11 DOUBLE PRECISION DI, HL0, PHL0, R |
|
1645
|
12 DIMENSION WM(*), IWM(*), X(*), TEM(*) |
|
1644
|
13 COMMON /LS0001/ ROWNS(209), |
|
|
14 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, |
|
|
15 3 IOWND(14), IOWNS(6), |
|
|
16 4 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, |
|
|
17 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
|
|
18 C----------------------------------------------------------------------- |
|
|
19 C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM |
|
|
20 C A CHORD ITERATION. IT IS CALLED IF MITER .NE. 0. |
|
|
21 C IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS. |
|
|
22 C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL |
|
|
23 C MATRIX, AND THEN COMPUTES THE SOLUTION. |
|
|
24 C IF MITER IS 4 OR 5, IT CALLS DGBSL. |
|
|
25 C COMMUNICATION WITH SOLSY USES THE FOLLOWING VARIABLES.. |
|
|
26 C WM = REAL WORK SPACE CONTAINING THE INVERSE DIAGONAL MATRIX IF |
|
|
27 C MITER = 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE. |
|
|
28 C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3). |
|
|
29 C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA.. |
|
|
30 C WM(1) = SQRT(UROUND) (NOT USED HERE), |
|
|
31 C WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER = 3. |
|
|
32 C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT |
|
|
33 C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND |
|
|
34 C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5. |
|
|
35 C X = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION VECTOR |
|
|
36 C ON OUTPUT, OF LENGTH N. |
|
|
37 C TEM = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION. |
|
|
38 C IERSL = OUTPUT FLAG (IN COMMON). IERSL = 0 IF NO TROUBLE OCCURRED. |
|
|
39 C IERSL = 1 IF A SINGULAR MATRIX AROSE WITH MITER = 3. |
|
|
40 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N. |
|
|
41 C----------------------------------------------------------------------- |
|
|
42 IERSL = 0 |
|
|
43 GO TO (100, 100, 300, 400, 400), MITER |
|
|
44 100 CALL DGESL (WM(3), N, N, IWM(21), X, 0) |
|
|
45 RETURN |
|
|
46 C |
|
|
47 300 PHL0 = WM(2) |
|
|
48 HL0 = H*EL0 |
|
|
49 WM(2) = HL0 |
|
|
50 IF (HL0 .EQ. PHL0) GO TO 330 |
|
|
51 R = HL0/PHL0 |
|
|
52 DO 320 I = 1,N |
|
|
53 DI = 1.0D0 - R*(1.0D0 - 1.0D0/WM(I+2)) |
|
|
54 IF (DABS(DI) .EQ. 0.0D0) GO TO 390 |
|
|
55 320 WM(I+2) = 1.0D0/DI |
|
|
56 330 DO 340 I = 1,N |
|
|
57 340 X(I) = WM(I+2)*X(I) |
|
|
58 RETURN |
|
|
59 390 IERSL = 1 |
|
|
60 RETURN |
|
|
61 C |
|
|
62 400 ML = IWM(1) |
|
|
63 MU = IWM(2) |
|
|
64 MEBAND = 2*ML + MU + 1 |
|
|
65 CALL DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0) |
|
|
66 RETURN |
|
|
67 C----------------------- END OF SUBROUTINE SOLSY ----------------------- |
|
|
68 END |
