1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
SUBROUTINE FBS2 (BLOCK,Y,YN,NWDS)
C
C FBS2 EXECUTES THE FORWARD/BACKWARD PASS FOR FBSF IN RDP
C
INTEGER BLOCK(8), DBL, BUF(2), SUBNAM, BEGN, END
DOUBLE PRECISION Y(1), YN(1), LJJ, L, YJK, SUM, ZERO
CHARACTER UFM*23, UWM*25, UIM*29, SFM*25
COMMON /XMSSG / UFM, UWM, UIM, SFM
COMMON /MACHIN/ MACH
COMMON /SYSTEM/ SYSBUF, NOUT
COMMON /ZZZZZZ/ L(1)
COMMON /FBSX / DBL, N
DATA ZERO / 0.0D+0 /
DATA SUBNAM, BEGN, END / 4HFBS2, 4HBEGN, 4HEND /
C
BUF(1) = SUBNAM
BUF(2) = BEGN
CALL CONMSG (BUF,2,0)
NBRITM = NWDS/2
J = (LOCFX(YN)-LOCFX(Y)+1)/NWDS
LAST = MAX0(J,1)*NBRITM
DO 35 J = 1,N
J1 = J - 1
DO 5 K = J,LAST,NBRITM
IF (Y(K) .NE. ZERO) GO TO 7
5 CONTINUE
CALL SKPREC (BLOCK(1),1)
GO TO 35
C
C MAKE 1ST STRING CALL FOR COLUMN AND SAVE DIAGONAL ELEMENT
C
7 BLOCK(8) = -1
CALL GETSTR (*80,BLOCK)
IF (BLOCK(4) .NE. J) GO TO 80
JSTR = BLOCK(5)
LJJ = 1.0D+0/L(JSTR)
IF (BLOCK(6) .EQ. 1) GO TO 20
NSTR = JSTR + BLOCK(6) - 1
JSTR = JSTR + 1
BLOCK(4) = BLOCK(4) + 1
C
C PROCESS CURRENT STRING IN TRIANGULAR FACTOR AGAINST EACH
C LOAD VECTOR IN CORE -- Y(I,K) = Y(I,K) + L(I,J)*Y(J,K)
C
10 DO 15 K = 1,LAST,NBRITM
YJK = Y(J1+K)
IF (YJK .EQ. ZERO) GO TO 15
IK = BLOCK(4) + K - 1
DO 12 IJ = JSTR,NSTR
Y(IK) = Y(IK) + L(IJ)*YJK
12 IK = IK + 1
15 CONTINUE
C
C GET NEXT STRING IN TRIANGULAR FACTOR
C
20 CALL ENDGET (BLOCK)
CALL GETSTR (*30,BLOCK)
JSTR = BLOCK(5)
NSTR = JSTR + BLOCK(6) - 1
GO TO 10
C
C END-OF-COLUMN ON TRIANGULAR FACTOR -- DIVIDE BY DIAGONAL
C
30 DO 32 K = J,LAST,NBRITM
32 Y(K) = Y(K)*LJJ
C
35 CONTINUE
C
C INITIALIZE FOR BACKWARD PASS BY SKIPPING THE NTH COLUMN
C
IF (N .EQ. 1) GO TO 70
CALL BCKREC (BLOCK)
J = N - 1
C
C GET A STRING IN CURRENT COLUMN. IF THIS STRING INCLUDES DIAGONAL,
C ADJUST STRING TO SKIP IT.
C
40 J1 = J - 1
BLOCK(8) = -1
42 CALL GETSTB (*60,BLOCK)
IF (BLOCK(4)-BLOCK(6) .EQ. J1) BLOCK(6) = BLOCK(6) - 1
IF (BLOCK(6) .EQ. 0) GO TO 58
NTERMS = BLOCK(6)
C
C PROCESS CURRENT STRING IN TRIANGULAR FACTOR AGAINST EACH
C LOAD VECTOR IN CORE -- Y(J,K) = Y(J,K) + L(J,I)*Y(I,K)
C
DO 55 K = 1,LAST,NBRITM
JI = BLOCK(5) + 1
IK = BLOCK(4) + K
SUM = 0.0D+0
DO 53 II = 1,NTERMS
JI = JI - 1
IK = IK - 1
SUM = SUM + L(JI)*Y(IK)
53 CONTINUE
Y(J1+K) = Y(J1+K) + SUM
55 CONTINUE
C
C TERMINATE CURRENT STRING AND GET NEXT STRING
C
58 CALL ENDGTB (BLOCK)
GO TO 42
C
C END-OF-COLUMN -- TEST FOR COMPLETION
C
60 IF (J .EQ. 1) GO TO 70
J = J - 1
GO TO 40
C
70 BUF(2) = END
CALL CONMSG (BUF,2,0)
RETURN
C
C
C FATAL ERROR MESSAGE
C
80 WRITE (NOUT,82) SFM,SUBNAM
82 FORMAT (A25,' 2149, SUBROUTINE ',A4,/5X,'FIRST ELEMENT OF A COLU',
1 'MN OF LOWER TRIANGULAR MATRIX IS NOT THE DIAGONAL ELEMENT')
CALL MESAGE (-61,0,0)
RETURN
END
|
