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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
|
SUBROUTINE PBDMATADD( ICONTXT, MODE, M, N, ALPHA, A, LDA, BETA, B,
$ LDB )
*
* -- PB-BLAS routine (version 2.1) --
* University of Tennessee, Knoxville, Oak Ridge National Laboratory.
* April 28, 1996
*
* .. Scalar Arguments ..
CHARACTER*1 MODE
INTEGER ICONTXT, LDA, LDB, M, N
DOUBLE PRECISION ALPHA, BETA
* ..
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
* Purpose
* =======
*
* PBDMATADD performs the matrix add operation B := alpha*A + beta*B,
* where alpha and beta are scalars, and A and B are m-by-n
* upper/lower trapezoidal matrices, or rectangular matrices.
*
* Arguments
* =========
*
* ICONTXT (input) INTEGER
* ICONTXT is the BLACS mechanism for partitioning communication
* space. A defining property of a context is that a message in
* a context cannot be sent or received in another context. The
* BLACS context includes the definition of a grid, and each
* process' coordinates in it.
*
* MODE (input) CHARACTER*1
* Specifies the part of the matrix A, or (conjugate) transposed
* matrix A to be added to the matrix B,
* = 'U': Upper triangular part
* up(B) = alpha*up(A) + beta*up(B)
* = 'L': Lower triangular part
* lo(B) = alpha*lo(A) + beta*lo(B)
* = 'T': Transposed matrix A
* B = alpha*A**T + beta*B
* = 'C': Conjugate transposed matrix A
* B = alpha*A**H + beta*B
* Otherwise: B = alpha*A + beta*B
* if M = LDA = LDB: use one BLAS loop
* if MODE = 'V' : columnwise copy using BLAS if possible
* else : use double loops
*
* M (input) INTEGER
* M specifies the number of columns of the matrix A if
* MODE != 'T'/'C', and it specifies the number of rows of the
* matrix A otherwise. It also specifies the number of rows of
* the matrix B. M >= 0.
*
* N (input) INTEGER
* N specifies the number of rows of the matrix A if
* MODE != 'T'/'C', and it specifies the number of columns of
* the matrix A otherwise. It also specifies the number of
* columns of the matrix B. N >= 0.
*
* ALPHA (input) DOUBLE PRECISION
* ALPHA specifies the scalar alpha.
*
* A (input) DOUBLE PRECISION array, dimension (LDA,N)
* The m by n matrix A if MODE != 'T'/'C'.
* If MODE = 'U', only the upper triangle or trapezoid is
* accessed; if MODE = 'L', only the lower triangle or
* trapezoid is accessed. Otherwise all m-by-n data matrix
* is accessed.
* And the n by m matrix A if MODE = 'T'/'C'.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,M) if
* MODE != 'T'/'C'. And LDA >= max(1,N) if MODE = 'T'/'C'.
*
* BETA (input) DOUBLE PRECISION
* BETA specifies the scalar beta.
*
* B (input) DOUBLE PRECISION array, dimension (LDB,N)
* On exit, B = alpha*A + beta*B
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,M).
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0)
* ..
* .. Local Scalars ..
INTEGER I, J
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL DSCAL, DCOPY, DAXPY
* ..
* .. Intrinsic Functions ..
INTRINSIC MIN
* ..
* .. Executable Statements ..
*
IF( M.LE.0 .OR. N.LE.0 .OR. ( ALPHA.EQ.ZERO.AND.BETA.EQ.ONE ) )
$ RETURN
*
* A is upper triangular or upper trapezoidal,
*
IF( LSAME( MODE, 'U' ) ) THEN
IF( ALPHA.EQ.ZERO ) THEN
IF( BETA.EQ.ZERO ) THEN
DO 20 J = 1, N
DO 10 I = 1, MIN( J, M )
B( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40 J = 1, N
DO 30 I = 1, MIN( J, M )
B( I, J ) = BETA * B( I, J )
30 CONTINUE
40 CONTINUE
END IF
*
ELSE IF( ALPHA.EQ.ONE ) THEN
IF( BETA.EQ.ZERO ) THEN
DO 60 J = 1, N
DO 50 I = 1, MIN( J, M )
B( I, J ) = A( I, J )
50 CONTINUE
60 CONTINUE
ELSE IF( BETA.EQ.ONE ) THEN
DO 80 J = 1, N
DO 70 I = 1, MIN( J, M )
B( I, J ) = A( I, J ) + B( I, J )
70 CONTINUE
80 CONTINUE
ELSE
DO 100 J = 1, N
DO 90 I = 1, MIN( J, M )
B( I, J ) = A( I, J ) + BETA * B( I, J )
90 CONTINUE
100 CONTINUE
END IF
*
ELSE
IF( BETA.EQ.ZERO ) THEN
DO 120 J = 1, N
DO 110 I = 1, MIN( J, M )
B( I, J ) = ALPHA * A( I, J )
110 CONTINUE
120 CONTINUE
ELSE IF( BETA.EQ.ONE ) THEN
DO 140 J = 1, N
DO 130 I = 1, MIN( J, M )
B( I, J ) = ALPHA * A( I, J ) + B( I, J )
130 CONTINUE
140 CONTINUE
ELSE
DO 160 J = 1, N
DO 150 I = 1, MIN( J, M )
B( I, J ) = ALPHA * A( I, J ) + BETA * B( I, J )
150 CONTINUE
160 CONTINUE
END IF
END IF
*
* A is lower triangular or upper trapezoidal,
*
ELSE IF( LSAME( MODE, 'L' ) ) THEN
IF( ALPHA.EQ.ZERO ) THEN
IF( BETA.EQ.ZERO ) THEN
DO 180 J = 1, N
DO 170 I = J, M
B( I, J ) = ZERO
170 CONTINUE
180 CONTINUE
ELSE
DO 200 J = 1, N
DO 190 I = J, M
B( I, J ) = BETA * B( I, J )
190 CONTINUE
200 CONTINUE
END IF
*
ELSE IF( ALPHA.EQ.ONE ) THEN
IF( BETA.EQ.ZERO ) THEN
DO 220 J = 1, N
DO 210 I = J, M
B( I, J ) = A( I, J )
210 CONTINUE
220 CONTINUE
ELSE IF( BETA.EQ.ONE ) THEN
DO 240 J = 1, N
DO 230 I = J, M
B( I, J ) = A( I, J ) + B( I, J )
230 CONTINUE
240 CONTINUE
ELSE
DO 260 J = 1, N
DO 250 I = J, M
B( I, J ) = A( I, J ) + BETA * B( I, J )
250 CONTINUE
260 CONTINUE
END IF
*
ELSE
IF( BETA.EQ.ZERO ) THEN
DO 280 J = 1, N
DO 270 I = J, M
B( I, J ) = ALPHA * A( I, J )
270 CONTINUE
280 CONTINUE
ELSE IF( BETA.EQ.ONE ) THEN
DO 300 J = 1, N
DO 290 I = J, M
B( I, J ) = ALPHA * A( I, J ) + B( I, J )
290 CONTINUE
300 CONTINUE
ELSE
DO 320 J = 1, N
DO 310 I = J, M
B( I, J ) = ALPHA * A( I, J ) + BETA * B( I, J )
310 CONTINUE
320 CONTINUE
END IF
END IF
*
* If MODE = 'Transpose'/'Conjugate'
*
ELSE IF( LSAME( MODE, 'T' ) .OR. LSAME( MODE, 'C' ) ) THEN
IF( ALPHA.EQ.ZERO ) THEN
IF( BETA.EQ.ZERO ) THEN
DO 340 J = 1, N
DO 330 I = 1, M
B( I, J ) = ZERO
330 CONTINUE
340 CONTINUE
ELSE
DO 360 J = 1, N
DO 350 I = 1, M
B( I, J ) = BETA * B( I, J )
350 CONTINUE
360 CONTINUE
END IF
*
ELSE IF( ALPHA.EQ.ONE ) THEN
IF( BETA.EQ.ZERO ) THEN
DO 380 J = 1, N
DO 370 I = 1, M
B( I, J ) = A( J, I )
370 CONTINUE
380 CONTINUE
ELSE IF( BETA.EQ.ONE ) THEN
DO 400 J = 1, N
DO 390 I = 1, M
B( I, J ) = A( J, I ) + B( I, J )
390 CONTINUE
400 CONTINUE
ELSE
DO 420 J = 1, N
DO 410 I = 1, M
B( I, J ) = A( J, I ) + BETA * B( I, J )
410 CONTINUE
420 CONTINUE
END IF
*
ELSE
IF( BETA.EQ.ZERO ) THEN
DO 440 J = 1, N
DO 430 I = 1, M
B( I, J ) = ALPHA * A( J, I )
430 CONTINUE
440 CONTINUE
ELSE IF( BETA.EQ.ONE ) THEN
DO 460 J = 1, N
DO 450 I = 1, M
B( I, J ) = ALPHA * A( J, I ) + B( I, J )
450 CONTINUE
460 CONTINUE
ELSE
DO 480 J = 1, N
DO 470 I = 1, M
B( I, J ) = ALPHA * A( J, I ) + BETA * B( I, J )
470 CONTINUE
480 CONTINUE
END IF
END IF
*
* Other cases (for genral matrix additions)
*
ELSE
IF( ALPHA.EQ.ZERO ) THEN
IF( BETA.EQ.ZERO ) THEN
DO 500 J = 1, N
DO 490 I = 1, M
B( I, J ) = ZERO
490 CONTINUE
500 CONTINUE
*
ELSE
IF( M.EQ.LDB ) THEN
CALL DSCAL( M*N, BETA, B( 1, 1 ), 1 )
ELSE IF( LSAME( MODE, 'V' ) ) THEN
DO 510 J = 1, N
CALL DSCAL( M, BETA, B( 1, J ), 1 )
510 CONTINUE
ELSE
DO 530 J = 1, N
DO 520 I = 1, M
B( I, J ) = BETA * B( I, J )
520 CONTINUE
530 CONTINUE
END IF
END IF
*
ELSE IF( ALPHA.EQ.ONE ) THEN
IF( BETA.EQ.ZERO ) THEN
IF( M.EQ.LDA .AND. M.EQ.LDB ) THEN
CALL DCOPY( M*N, A( 1, 1 ), 1, B( 1, 1 ), 1 )
ELSE IF( LSAME( MODE, 'V' ) ) THEN
DO 540 J = 1, N
CALL DCOPY( M, A( 1, J ), 1, B( 1, J ), 1 )
540 CONTINUE
ELSE
DO 560 J = 1, N
DO 550 I = 1, M
B( I, J ) = A( I, J )
550 CONTINUE
560 CONTINUE
END IF
*
ELSE IF( BETA.EQ.ONE ) THEN
DO 580 J = 1, N
DO 570 I = 1, M
B( I, J ) = A( I, J ) + B( I, J )
570 CONTINUE
580 CONTINUE
*
ELSE
DO 600 J = 1, N
DO 590 I = 1, M
B( I, J ) = A( I, J ) + BETA * B( I, J )
590 CONTINUE
600 CONTINUE
END IF
*
ELSE
IF( BETA.EQ.ZERO ) THEN
DO 620 J = 1, N
DO 610 I = 1, M
B( I, J ) = ALPHA * A( I, J )
610 CONTINUE
620 CONTINUE
*
ELSE IF( BETA.EQ.ONE ) THEN
IF( M.EQ.LDA .AND. M.EQ.LDB ) THEN
CALL DAXPY( M*N, ALPHA, A( 1, 1 ), 1, B( 1, 1 ), 1 )
ELSE IF( LSAME( MODE, 'V' ) ) THEN
DO 630 J = 1, N
CALL DAXPY( M, ALPHA, A( 1, J ), 1, B( 1, J ), 1 )
630 CONTINUE
ELSE
DO 650 J = 1, N
DO 640 I = 1, M
B( I, J ) = ALPHA * A( I, J ) + B( I, J )
640 CONTINUE
650 CONTINUE
END IF
*
ELSE
DO 670 J = 1, N
DO 660 I = 1, M
B( I, J ) = ALPHA * A( I, J ) + BETA * B( I, J )
660 CONTINUE
670 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of PBDMATADD
*
END
|
