ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
Go to the documentation of this file.
1  SUBROUTINE pbdmatadd( ICONTXT, MODE, M, N, ALPHA, A, LDA, BETA, B,
2  \$ LDB )
3 *
4 * -- PB-BLAS routine (version 2.1) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory.
6 * April 28, 1996
7 *
8 * .. Scalar Arguments ..
9  CHARACTER*1 MODE
10  INTEGER ICONTXT, LDA, LDB, M, N
11  DOUBLE PRECISION ALPHA, BETA
12 * ..
13 * .. Array Arguments ..
14  DOUBLE PRECISION A( LDA, * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * PBDMATADD performs the matrix add operation B := alpha*A + beta*B,
21 * where alpha and beta are scalars, and A and B are m-by-n
22 * upper/lower trapezoidal matrices, or rectangular matrices.
23 *
24 * Arguments
25 * =========
26 *
27 * ICONTXT (input) INTEGER
28 * ICONTXT is the BLACS mechanism for partitioning communication
29 * space. A defining property of a context is that a message in
30 * a context cannot be sent or received in another context. The
31 * BLACS context includes the definition of a grid, and each
32 * process' coordinates in it.
33 *
34 * MODE (input) CHARACTER*1
35 * Specifies the part of the matrix A, or (conjugate) transposed
36 * matrix A to be added to the matrix B,
37 * = 'U': Upper triangular part
38 * up(B) = alpha*up(A) + beta*up(B)
39 * = 'L': Lower triangular part
40 * lo(B) = alpha*lo(A) + beta*lo(B)
41 * = 'T': Transposed matrix A
42 * B = alpha*A**T + beta*B
43 * = 'C': Conjugate transposed matrix A
44 * B = alpha*A**H + beta*B
45 * Otherwise: B = alpha*A + beta*B
46 * if M = LDA = LDB: use one BLAS loop
47 * if MODE = 'V' : columnwise copy using BLAS if possible
48 * else : use double loops
49 *
50 * M (input) INTEGER
51 * M specifies the number of columns of the matrix A if
52 * MODE != 'T'/'C', and it specifies the number of rows of the
53 * matrix A otherwise. It also specifies the number of rows of
54 * the matrix B. M >= 0.
55 *
56 * N (input) INTEGER
57 * N specifies the number of rows of the matrix A if
58 * MODE != 'T'/'C', and it specifies the number of columns of
59 * the matrix A otherwise. It also specifies the number of
60 * columns of the matrix B. N >= 0.
61 *
62 * ALPHA (input) DOUBLE PRECISION
63 * ALPHA specifies the scalar alpha.
64 *
65 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
66 * The m by n matrix A if MODE != 'T'/'C'.
67 * If MODE = 'U', only the upper triangle or trapezoid is
68 * accessed; if MODE = 'L', only the lower triangle or
69 * trapezoid is accessed. Otherwise all m-by-n data matrix
70 * is accessed.
71 * And the n by m matrix A if MODE = 'T'/'C'.
72 *
73 * LDA (input) INTEGER
74 * The leading dimension of the array A. LDA >= max(1,M) if
75 * MODE != 'T'/'C'. And LDA >= max(1,N) if MODE = 'T'/'C'.
76 *
77 * BETA (input) DOUBLE PRECISION
78 * BETA specifies the scalar beta.
79 *
80 * B (input) DOUBLE PRECISION array, dimension (LDB,N)
81 * On exit, B = alpha*A + beta*B
82 *
83 * LDB (input) INTEGER
84 * The leading dimension of the array B. LDB >= max(1,M).
85 *
86 * =====================================================================
87 *
88 * .. Parameters ..
89  DOUBLE PRECISION ZERO, ONE
90  parameter( zero = 0.0d+0, one = 1.0d+0)
91 * ..
92 * .. Local Scalars ..
93  INTEGER I, J
94 * ..
95 * .. External Functions ..
96  LOGICAL LSAME
97  EXTERNAL lsame
98 * ..
99 * .. External Subroutines ..
100  EXTERNAL dscal, dcopy, daxpy
101 * ..
102 * .. Intrinsic Functions ..
103  INTRINSIC min
104 * ..
105 * .. Executable Statements ..
106 *
107  IF( m.LE.0 .OR. n.LE.0 .OR. ( alpha.EQ.zero.AND.beta.EQ.one ) )
108  \$ RETURN
109 *
110 * A is upper triangular or upper trapezoidal,
111 *
112  IF( lsame( mode, 'U' ) ) THEN
113  IF( alpha.EQ.zero ) THEN
114  IF( beta.EQ.zero ) THEN
115  DO 20 j = 1, n
116  DO 10 i = 1, min( j, m )
117  b( i, j ) = zero
118  10 CONTINUE
119  20 CONTINUE
120  ELSE
121  DO 40 j = 1, n
122  DO 30 i = 1, min( j, m )
123  b( i, j ) = beta * b( i, j )
124  30 CONTINUE
125  40 CONTINUE
126  END IF
127 *
128  ELSE IF( alpha.EQ.one ) THEN
129  IF( beta.EQ.zero ) THEN
130  DO 60 j = 1, n
131  DO 50 i = 1, min( j, m )
132  b( i, j ) = a( i, j )
133  50 CONTINUE
134  60 CONTINUE
135  ELSE IF( beta.EQ.one ) THEN
136  DO 80 j = 1, n
137  DO 70 i = 1, min( j, m )
138  b( i, j ) = a( i, j ) + b( i, j )
139  70 CONTINUE
140  80 CONTINUE
141  ELSE
142  DO 100 j = 1, n
143  DO 90 i = 1, min( j, m )
144  b( i, j ) = a( i, j ) + beta * b( i, j )
145  90 CONTINUE
146  100 CONTINUE
147  END IF
148 *
149  ELSE
150  IF( beta.EQ.zero ) THEN
151  DO 120 j = 1, n
152  DO 110 i = 1, min( j, m )
153  b( i, j ) = alpha * a( i, j )
154  110 CONTINUE
155  120 CONTINUE
156  ELSE IF( beta.EQ.one ) THEN
157  DO 140 j = 1, n
158  DO 130 i = 1, min( j, m )
159  b( i, j ) = alpha * a( i, j ) + b( i, j )
160  130 CONTINUE
161  140 CONTINUE
162  ELSE
163  DO 160 j = 1, n
164  DO 150 i = 1, min( j, m )
165  b( i, j ) = alpha * a( i, j ) + beta * b( i, j )
166  150 CONTINUE
167  160 CONTINUE
168  END IF
169  END IF
170 *
171 * A is lower triangular or upper trapezoidal,
172 *
173  ELSE IF( lsame( mode, 'L' ) ) THEN
174  IF( alpha.EQ.zero ) THEN
175  IF( beta.EQ.zero ) THEN
176  DO 180 j = 1, n
177  DO 170 i = j, m
178  b( i, j ) = zero
179  170 CONTINUE
180  180 CONTINUE
181  ELSE
182  DO 200 j = 1, n
183  DO 190 i = j, m
184  b( i, j ) = beta * b( i, j )
185  190 CONTINUE
186  200 CONTINUE
187  END IF
188 *
189  ELSE IF( alpha.EQ.one ) THEN
190  IF( beta.EQ.zero ) THEN
191  DO 220 j = 1, n
192  DO 210 i = j, m
193  b( i, j ) = a( i, j )
194  210 CONTINUE
195  220 CONTINUE
196  ELSE IF( beta.EQ.one ) THEN
197  DO 240 j = 1, n
198  DO 230 i = j, m
199  b( i, j ) = a( i, j ) + b( i, j )
200  230 CONTINUE
201  240 CONTINUE
202  ELSE
203  DO 260 j = 1, n
204  DO 250 i = j, m
205  b( i, j ) = a( i, j ) + beta * b( i, j )
206  250 CONTINUE
207  260 CONTINUE
208  END IF
209 *
210  ELSE
211  IF( beta.EQ.zero ) THEN
212  DO 280 j = 1, n
213  DO 270 i = j, m
214  b( i, j ) = alpha * a( i, j )
215  270 CONTINUE
216  280 CONTINUE
217  ELSE IF( beta.EQ.one ) THEN
218  DO 300 j = 1, n
219  DO 290 i = j, m
220  b( i, j ) = alpha * a( i, j ) + b( i, j )
221  290 CONTINUE
222  300 CONTINUE
223  ELSE
224  DO 320 j = 1, n
225  DO 310 i = j, m
226  b( i, j ) = alpha * a( i, j ) + beta * b( i, j )
227  310 CONTINUE
228  320 CONTINUE
229  END IF
230  END IF
231 *
232 * If MODE = 'Transpose'/'Conjugate'
233 *
234  ELSE IF( lsame( mode, 'T' ) .OR. lsame( mode, 'C' ) ) THEN
235  IF( alpha.EQ.zero ) THEN
236  IF( beta.EQ.zero ) THEN
237  DO 340 j = 1, n
238  DO 330 i = 1, m
239  b( i, j ) = zero
240  330 CONTINUE
241  340 CONTINUE
242  ELSE
243  DO 360 j = 1, n
244  DO 350 i = 1, m
245  b( i, j ) = beta * b( i, j )
246  350 CONTINUE
247  360 CONTINUE
248  END IF
249 *
250  ELSE IF( alpha.EQ.one ) THEN
251  IF( beta.EQ.zero ) THEN
252  DO 380 j = 1, n
253  DO 370 i = 1, m
254  b( i, j ) = a( j, i )
255  370 CONTINUE
256  380 CONTINUE
257  ELSE IF( beta.EQ.one ) THEN
258  DO 400 j = 1, n
259  DO 390 i = 1, m
260  b( i, j ) = a( j, i ) + b( i, j )
261  390 CONTINUE
262  400 CONTINUE
263  ELSE
264  DO 420 j = 1, n
265  DO 410 i = 1, m
266  b( i, j ) = a( j, i ) + beta * b( i, j )
267  410 CONTINUE
268  420 CONTINUE
269  END IF
270 *
271  ELSE
272  IF( beta.EQ.zero ) THEN
273  DO 440 j = 1, n
274  DO 430 i = 1, m
275  b( i, j ) = alpha * a( j, i )
276  430 CONTINUE
277  440 CONTINUE
278  ELSE IF( beta.EQ.one ) THEN
279  DO 460 j = 1, n
280  DO 450 i = 1, m
281  b( i, j ) = alpha * a( j, i ) + b( i, j )
282  450 CONTINUE
283  460 CONTINUE
284  ELSE
285  DO 480 j = 1, n
286  DO 470 i = 1, m
287  b( i, j ) = alpha * a( j, i ) + beta * b( i, j )
288  470 CONTINUE
289  480 CONTINUE
290  END IF
291  END IF
292 *
293 * Other cases (for genral matrix additions)
294 *
295  ELSE
296  IF( alpha.EQ.zero ) THEN
297  IF( beta.EQ.zero ) THEN
298  DO 500 j = 1, n
299  DO 490 i = 1, m
300  b( i, j ) = zero
301  490 CONTINUE
302  500 CONTINUE
303 *
304  ELSE
305  IF( m.EQ.ldb ) THEN
306  CALL dscal( m*n, beta, b( 1, 1 ), 1 )
307  ELSE IF( lsame( mode, 'V' ) ) THEN
308  DO 510 j = 1, n
309  CALL dscal( m, beta, b( 1, j ), 1 )
310  510 CONTINUE
311  ELSE
312  DO 530 j = 1, n
313  DO 520 i = 1, m
314  b( i, j ) = beta * b( i, j )
315  520 CONTINUE
316  530 CONTINUE
317  END IF
318  END IF
319 *
320  ELSE IF( alpha.EQ.one ) THEN
321  IF( beta.EQ.zero ) THEN
322  IF( m.EQ.lda .AND. m.EQ.ldb ) THEN
323  CALL dcopy( m*n, a( 1, 1 ), 1, b( 1, 1 ), 1 )
324  ELSE IF( lsame( mode, 'V' ) ) THEN
325  DO 540 j = 1, n
326  CALL dcopy( m, a( 1, j ), 1, b( 1, j ), 1 )
327  540 CONTINUE
328  ELSE
329  DO 560 j = 1, n
330  DO 550 i = 1, m
331  b( i, j ) = a( i, j )
332  550 CONTINUE
333  560 CONTINUE
334  END IF
335 *
336  ELSE IF( beta.EQ.one ) THEN
337  DO 580 j = 1, n
338  DO 570 i = 1, m
339  b( i, j ) = a( i, j ) + b( i, j )
340  570 CONTINUE
341  580 CONTINUE
342 *
343  ELSE
344  DO 600 j = 1, n
345  DO 590 i = 1, m
346  b( i, j ) = a( i, j ) + beta * b( i, j )
347  590 CONTINUE
348  600 CONTINUE
349  END IF
350 *
351  ELSE
352  IF( beta.EQ.zero ) THEN
353  DO 620 j = 1, n
354  DO 610 i = 1, m
355  b( i, j ) = alpha * a( i, j )
356  610 CONTINUE
357  620 CONTINUE
358 *
359  ELSE IF( beta.EQ.one ) THEN
360  IF( m.EQ.lda .AND. m.EQ.ldb ) THEN
361  CALL daxpy( m*n, alpha, a( 1, 1 ), 1, b( 1, 1 ), 1 )
362  ELSE IF( lsame( mode, 'V' ) ) THEN
363  DO 630 j = 1, n
364  CALL daxpy( m, alpha, a( 1, j ), 1, b( 1, j ), 1 )
365  630 CONTINUE
366  ELSE
367  DO 650 j = 1, n
368  DO 640 i = 1, m
369  b( i, j ) = alpha * a( i, j ) + b( i, j )
370  640 CONTINUE
371  650 CONTINUE
372  END IF
373 *
374  ELSE
375  DO 670 j = 1, n
376  DO 660 i = 1, m
377  b( i, j ) = alpha * a( i, j ) + beta * b( i, j )
378  660 CONTINUE
379  670 CONTINUE
380  END IF
381  END IF
382  END IF
383 *
384  RETURN
385 *