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
124
125
126
127
128
129
130
131
132
133
134
135
|
SUBROUTINE MPYA3D (AA,BB,NROW,BAND,CC)
C
C WITH ENTRY MPYA3S (A,B,NROW,BAND,C)
C
C WAS NAMED DATBAD/DATBAS IN UAI CODE
C
C THESE ROUTINES PERFORM TRIPLE MATRIX MULTIPLY OF THE FORM
C
C T
C C = C + A * B * A
C
C ON TWO INCOMING ROW-LOADED MATRICES A AND B, AND ADD THEM TO
C MATRIX C
C
C THE INCOMING MATRICES MUST BE SQUARE (AND OBVIOUSLY OF THE SAME
C SIZE, NROW.) AND
C SYMMETRICAL (SINCE WE OPERATE ONLY ON LOWER TRIANGULAR MATRICES)
C
C MATRIX A CAN BE A PSUEDO-DIAGONAL MATRIX, I.E. A MATRIX HAVING
C SQUARE PARTITIONS OF NON-ZERO TERMS ALONG ITS DIAGONAL.
C THESE PARTITIONS ARE OF THE SIZE BAND X BAND.
C NOTE THAT NROW MUST BE AN INTEGER MULTIPLE OF BAND.
C
C THIS ALGORITHM IS SUITABLE FOR TRIPLE MULTIPLIES INVOLVING GLOBAL
C TRANSFORMATIONS.
C
C
INTEGER BAND
REAL A(1) ,B(1) ,C(1)
DOUBLE PRECISION AA(1),BB(1),CC(1),DD
C
C
C DOUBLE PRECISION VERSION
C
II = 0
DO 50 IB = 1,NROW
IA1 = ((IB-1)/BAND+1)*BAND
C
DO 40 ID1 = 1,NROW,BAND
ID2 = ID1 + BAND - 1
IF (ID1 .GT. IA1) GO TO 50
C
ID11N = (ID1-1)*NROW
DO 30 ID = ID1,ID2
JJ = ID11N
DD = 0.0D0
C
DO 10 IC = ID1,ID2
IBIC = II + IC
ICID = JJ + ID
IF (AA(ICID) .EQ. 0.0D0) GO TO 10
DD = DD + BB(IBIC)*AA(ICID)
10 JJ = JJ + NROW
C
IF (DD .EQ. 0.0D0) GO TO 30
KK = (ID-1)*NROW
C
DO 20 IA = ID,IA1
IBIA = II + IA
IF (AA(IBIA) .EQ. 0.0D0) GO TO 20
IAID = KK + ID
CC(IAID) = CC(IAID) + DD*AA(IBIA)
20 KK = KK + NROW
C
30 CONTINUE
40 CONTINUE
50 II = II + NROW
C
C COPY THE LOWER TRIANGLE TO THE UPPER
C
KK = NROW - 1
II = 0
DO 70 I = 1,KK
IB = I + 1
JJ = I*NROW
DO 60 J = IB,NROW
CC(II+J) = CC(JJ+I)
60 JJ = JJ + NROW
70 II = II + NROW
C
RETURN
C
C
ENTRY MPYA3S (A,B,NROW,BAND,C)
C ==============================
C
C SINGLE PRECISION VERSION
C
II = 0
DO 150 IB = 1,NROW
IA1 = ((IB-1)/BAND+1)*BAND
C
DO 140 ID1 = 1,NROW,BAND
ID2 = ID1 + BAND - 1
IF (ID1 .GT. IA1) GO TO 150
C
ID11N = (ID1-1)*NROW
DO 130 ID = ID1,ID2
JJ = ID11N
DD = 0.0D0
C
DO 110 IC = ID1,ID2
IBIC = II + IC
ICID = JJ + ID
IF (A(ICID) .EQ. 0.0) GO TO 110
DD = DD + DBLE(B(IBIC))*DBLE(A(ICID))
110 JJ = JJ + NROW
IF (DD .EQ. 0.0D0) GO TO 130
KK = (ID-1)*NROW
C
DO 120 IA = ID,IA1
IBIA = II + IA
IF (A(IBIA) .EQ. 0.0) GO TO 120
IAID = KK + ID
C(IAID) = SNGL(DBLE(C(IAID)) + DD*DBLE(A(IBIA)))
120 KK = KK + NROW
C
130 CONTINUE
140 CONTINUE
150 II = II + NROW
C
C COPY THE LOWER TRIANGLE TO THE UPPER
C
KK = NROW - 1
II = 0
DO 170 I = 1,KK
IB = I + 1
JJ = I*NROW
DO 160 J = IB,NROW
C(II+J) = C(JJ+I)
160 JJ = JJ + NROW
170 II = II + NROW
C
RETURN
END
|
