1 SUBROUTINE DSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
2 $ C )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 *
6 * -- Contributed by Julien Langou of the Univ. of Colorado Denver --
7 * -- April 2011 --
8 *
9 * -- LAPACK is a software package provided by Univ. of Tennessee, --
10 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
11 *
12 * ..
13 * .. Scalar Arguments ..
14 DOUBLE PRECISION ALPHA, BETA
15 INTEGER K, LDA, N
16 CHARACTER TRANS, TRANSR, UPLO
17 * ..
18 * .. Array Arguments ..
19 DOUBLE PRECISION A( LDA, * ), C( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * Level 3 BLAS like routine for C in RFP Format.
26 *
27 * DSFRK performs one of the symmetric rank--k operations
28 *
29 * C := alpha*A*A**T + beta*C,
30 *
31 * or
32 *
33 * C := alpha*A**T*A + beta*C,
34 *
35 * where alpha and beta are real scalars, C is an n--by--n symmetric
36 * matrix and A is an n--by--k matrix in the first case and a k--by--n
37 * matrix in the second case.
38 *
39 * Arguments
40 * ==========
41 *
42 * TRANSR (input) CHARACTER*1
43 * = 'N': The Normal Form of RFP A is stored;
44 * = 'T': The Transpose Form of RFP A is stored.
45 *
46 * UPLO (input) CHARACTER*1
47 * On entry, UPLO specifies whether the upper or lower
48 * triangular part of the array C is to be referenced as
49 * follows:
50 *
51 * UPLO = 'U' or 'u' Only the upper triangular part of C
52 * is to be referenced.
53 *
54 * UPLO = 'L' or 'l' Only the lower triangular part of C
55 * is to be referenced.
56 *
57 * Unchanged on exit.
58 *
59 * TRANS (input) CHARACTER*1
60 * On entry, TRANS specifies the operation to be performed as
61 * follows:
62 *
63 * TRANS = 'N' or 'n' C := alpha*A*A**T + beta*C.
64 *
65 * TRANS = 'T' or 't' C := alpha*A**T*A + beta*C.
66 *
67 * Unchanged on exit.
68 *
69 * N (input) INTEGER
70 * On entry, N specifies the order of the matrix C. N must be
71 * at least zero.
72 * Unchanged on exit.
73 *
74 * K (input) INTEGER
75 * On entry with TRANS = 'N' or 'n', K specifies the number
76 * of columns of the matrix A, and on entry with TRANS = 'T'
77 * or 't', K specifies the number of rows of the matrix A. K
78 * must be at least zero.
79 * Unchanged on exit.
80 *
81 * ALPHA (input) DOUBLE PRECISION
82 * On entry, ALPHA specifies the scalar alpha.
83 * Unchanged on exit.
84 *
85 * A (input) DOUBLE PRECISION array, dimension (LDA,ka)
86 * where KA
87 * is K when TRANS = 'N' or 'n', and is N otherwise. Before
88 * entry with TRANS = 'N' or 'n', the leading N--by--K part of
89 * the array A must contain the matrix A, otherwise the leading
90 * K--by--N part of the array A must contain the matrix A.
91 * Unchanged on exit.
92 *
93 * LDA (input) INTEGER
94 * On entry, LDA specifies the first dimension of A as declared
95 * in the calling (sub) program. When TRANS = 'N' or 'n'
96 * then LDA must be at least max( 1, n ), otherwise LDA must
97 * be at least max( 1, k ).
98 * Unchanged on exit.
99 *
100 * BETA (input) DOUBLE PRECISION
101 * On entry, BETA specifies the scalar beta.
102 * Unchanged on exit.
103 *
104 *
105 * C (input/output) DOUBLE PRECISION array, dimension (NT)
106 * NT = N*(N+1)/2. On entry, the symmetric matrix C in RFP
107 * Format. RFP Format is described by TRANSR, UPLO and N.
108 *
109 * =====================================================================
110 *
111 * ..
112 * .. Parameters ..
113 DOUBLE PRECISION ONE, ZERO
114 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
115 * ..
116 * .. Local Scalars ..
117 LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS
118 INTEGER INFO, NROWA, J, NK, N1, N2
119 * ..
120 * .. External Functions ..
121 LOGICAL LSAME
122 EXTERNAL LSAME
123 * ..
124 * .. External Subroutines ..
125 EXTERNAL XERBLA, DGEMM, DSYRK
126 * ..
127 * .. Intrinsic Functions ..
128 INTRINSIC MAX
129 * ..
130 * .. Executable Statements ..
131 *
132 * Test the input parameters.
133 *
134 INFO = 0
135 NORMALTRANSR = LSAME( TRANSR, 'N' )
136 LOWER = LSAME( UPLO, 'L' )
137 NOTRANS = LSAME( TRANS, 'N' )
138 *
139 IF( NOTRANS ) THEN
140 NROWA = N
141 ELSE
142 NROWA = K
143 END IF
144 *
145 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
146 INFO = -1
147 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
148 INFO = -2
149 ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
150 INFO = -3
151 ELSE IF( N.LT.0 ) THEN
152 INFO = -4
153 ELSE IF( K.LT.0 ) THEN
154 INFO = -5
155 ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
156 INFO = -8
157 END IF
158 IF( INFO.NE.0 ) THEN
159 CALL XERBLA( 'DSFRK ', -INFO )
160 RETURN
161 END IF
162 *
163 * Quick return if possible.
164 *
165 * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
166 * done (it is in DSYRK for example) and left in the general case.
167 *
168 IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
169 $ ( BETA.EQ.ONE ) ) )RETURN
170 *
171 IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
172 DO J = 1, ( ( N*( N+1 ) ) / 2 )
173 C( J ) = ZERO
174 END DO
175 RETURN
176 END IF
177 *
178 * C is N-by-N.
179 * If N is odd, set NISODD = .TRUE., and N1 and N2.
180 * If N is even, NISODD = .FALSE., and NK.
181 *
182 IF( MOD( N, 2 ).EQ.0 ) THEN
183 NISODD = .FALSE.
184 NK = N / 2
185 ELSE
186 NISODD = .TRUE.
187 IF( LOWER ) THEN
188 N2 = N / 2
189 N1 = N - N2
190 ELSE
191 N1 = N / 2
192 N2 = N - N1
193 END IF
194 END IF
195 *
196 IF( NISODD ) THEN
197 *
198 * N is odd
199 *
200 IF( NORMALTRANSR ) THEN
201 *
202 * N is odd and TRANSR = 'N'
203 *
204 IF( LOWER ) THEN
205 *
206 * N is odd, TRANSR = 'N', and UPLO = 'L'
207 *
208 IF( NOTRANS ) THEN
209 *
210 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
211 *
212 CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
213 $ BETA, C( 1 ), N )
214 CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
215 $ BETA, C( N+1 ), N )
216 CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
217 $ LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
218 *
219 ELSE
220 *
221 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
222 *
223 CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
224 $ BETA, C( 1 ), N )
225 CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
226 $ BETA, C( N+1 ), N )
227 CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
228 $ LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
229 *
230 END IF
231 *
232 ELSE
233 *
234 * N is odd, TRANSR = 'N', and UPLO = 'U'
235 *
236 IF( NOTRANS ) THEN
237 *
238 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
239 *
240 CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
241 $ BETA, C( N2+1 ), N )
242 CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
243 $ BETA, C( N1+1 ), N )
244 CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
245 $ LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N )
246 *
247 ELSE
248 *
249 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
250 *
251 CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
252 $ BETA, C( N2+1 ), N )
253 CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA,
254 $ BETA, C( N1+1 ), N )
255 CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
256 $ LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N )
257 *
258 END IF
259 *
260 END IF
261 *
262 ELSE
263 *
264 * N is odd, and TRANSR = 'T'
265 *
266 IF( LOWER ) THEN
267 *
268 * N is odd, TRANSR = 'T', and UPLO = 'L'
269 *
270 IF( NOTRANS ) THEN
271 *
272 * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
273 *
274 CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
275 $ BETA, C( 1 ), N1 )
276 CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
277 $ BETA, C( 2 ), N1 )
278 CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
279 $ LDA, A( N1+1, 1 ), LDA, BETA,
280 $ C( N1*N1+1 ), N1 )
281 *
282 ELSE
283 *
284 * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
285 *
286 CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
287 $ BETA, C( 1 ), N1 )
288 CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
289 $ BETA, C( 2 ), N1 )
290 CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
291 $ LDA, A( 1, N1+1 ), LDA, BETA,
292 $ C( N1*N1+1 ), N1 )
293 *
294 END IF
295 *
296 ELSE
297 *
298 * N is odd, TRANSR = 'T', and UPLO = 'U'
299 *
300 IF( NOTRANS ) THEN
301 *
302 * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
303 *
304 CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
305 $ BETA, C( N2*N2+1 ), N2 )
306 CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
307 $ BETA, C( N1*N2+1 ), N2 )
308 CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
309 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
310 *
311 ELSE
312 *
313 * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
314 *
315 CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
316 $ BETA, C( N2*N2+1 ), N2 )
317 CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
318 $ BETA, C( N1*N2+1 ), N2 )
319 CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
320 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
321 *
322 END IF
323 *
324 END IF
325 *
326 END IF
327 *
328 ELSE
329 *
330 * N is even
331 *
332 IF( NORMALTRANSR ) THEN
333 *
334 * N is even and TRANSR = 'N'
335 *
336 IF( LOWER ) THEN
337 *
338 * N is even, TRANSR = 'N', and UPLO = 'L'
339 *
340 IF( NOTRANS ) THEN
341 *
342 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
343 *
344 CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
345 $ BETA, C( 2 ), N+1 )
346 CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
347 $ BETA, C( 1 ), N+1 )
348 CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
349 $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
350 $ N+1 )
351 *
352 ELSE
353 *
354 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
355 *
356 CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
357 $ BETA, C( 2 ), N+1 )
358 CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
359 $ BETA, C( 1 ), N+1 )
360 CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
361 $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
362 $ N+1 )
363 *
364 END IF
365 *
366 ELSE
367 *
368 * N is even, TRANSR = 'N', and UPLO = 'U'
369 *
370 IF( NOTRANS ) THEN
371 *
372 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
373 *
374 CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
375 $ BETA, C( NK+2 ), N+1 )
376 CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
377 $ BETA, C( NK+1 ), N+1 )
378 CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
379 $ LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ),
380 $ N+1 )
381 *
382 ELSE
383 *
384 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
385 *
386 CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
387 $ BETA, C( NK+2 ), N+1 )
388 CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
389 $ BETA, C( NK+1 ), N+1 )
390 CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
391 $ LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ),
392 $ N+1 )
393 *
394 END IF
395 *
396 END IF
397 *
398 ELSE
399 *
400 * N is even, and TRANSR = 'T'
401 *
402 IF( LOWER ) THEN
403 *
404 * N is even, TRANSR = 'T', and UPLO = 'L'
405 *
406 IF( NOTRANS ) THEN
407 *
408 * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
409 *
410 CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
411 $ BETA, C( NK+1 ), NK )
412 CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
413 $ BETA, C( 1 ), NK )
414 CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
415 $ LDA, A( NK+1, 1 ), LDA, BETA,
416 $ C( ( ( NK+1 )*NK )+1 ), NK )
417 *
418 ELSE
419 *
420 * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
421 *
422 CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
423 $ BETA, C( NK+1 ), NK )
424 CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
425 $ BETA, C( 1 ), NK )
426 CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
427 $ LDA, A( 1, NK+1 ), LDA, BETA,
428 $ C( ( ( NK+1 )*NK )+1 ), NK )
429 *
430 END IF
431 *
432 ELSE
433 *
434 * N is even, TRANSR = 'T', and UPLO = 'U'
435 *
436 IF( NOTRANS ) THEN
437 *
438 * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
439 *
440 CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
441 $ BETA, C( NK*( NK+1 )+1 ), NK )
442 CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
443 $ BETA, C( NK*NK+1 ), NK )
444 CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
445 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
446 *
447 ELSE
448 *
449 * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
450 *
451 CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
452 $ BETA, C( NK*( NK+1 )+1 ), NK )
453 CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
454 $ BETA, C( NK*NK+1 ), NK )
455 CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
456 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
457 *
458 END IF
459 *
460 END IF
461 *
462 END IF
463 *
464 END IF
465 *
466 RETURN
467 *
468 * End of DSFRK
469 *
470 END
2 $ C )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 *
6 * -- Contributed by Julien Langou of the Univ. of Colorado Denver --
7 * -- April 2011 --
8 *
9 * -- LAPACK is a software package provided by Univ. of Tennessee, --
10 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
11 *
12 * ..
13 * .. Scalar Arguments ..
14 DOUBLE PRECISION ALPHA, BETA
15 INTEGER K, LDA, N
16 CHARACTER TRANS, TRANSR, UPLO
17 * ..
18 * .. Array Arguments ..
19 DOUBLE PRECISION A( LDA, * ), C( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * Level 3 BLAS like routine for C in RFP Format.
26 *
27 * DSFRK performs one of the symmetric rank--k operations
28 *
29 * C := alpha*A*A**T + beta*C,
30 *
31 * or
32 *
33 * C := alpha*A**T*A + beta*C,
34 *
35 * where alpha and beta are real scalars, C is an n--by--n symmetric
36 * matrix and A is an n--by--k matrix in the first case and a k--by--n
37 * matrix in the second case.
38 *
39 * Arguments
40 * ==========
41 *
42 * TRANSR (input) CHARACTER*1
43 * = 'N': The Normal Form of RFP A is stored;
44 * = 'T': The Transpose Form of RFP A is stored.
45 *
46 * UPLO (input) CHARACTER*1
47 * On entry, UPLO specifies whether the upper or lower
48 * triangular part of the array C is to be referenced as
49 * follows:
50 *
51 * UPLO = 'U' or 'u' Only the upper triangular part of C
52 * is to be referenced.
53 *
54 * UPLO = 'L' or 'l' Only the lower triangular part of C
55 * is to be referenced.
56 *
57 * Unchanged on exit.
58 *
59 * TRANS (input) CHARACTER*1
60 * On entry, TRANS specifies the operation to be performed as
61 * follows:
62 *
63 * TRANS = 'N' or 'n' C := alpha*A*A**T + beta*C.
64 *
65 * TRANS = 'T' or 't' C := alpha*A**T*A + beta*C.
66 *
67 * Unchanged on exit.
68 *
69 * N (input) INTEGER
70 * On entry, N specifies the order of the matrix C. N must be
71 * at least zero.
72 * Unchanged on exit.
73 *
74 * K (input) INTEGER
75 * On entry with TRANS = 'N' or 'n', K specifies the number
76 * of columns of the matrix A, and on entry with TRANS = 'T'
77 * or 't', K specifies the number of rows of the matrix A. K
78 * must be at least zero.
79 * Unchanged on exit.
80 *
81 * ALPHA (input) DOUBLE PRECISION
82 * On entry, ALPHA specifies the scalar alpha.
83 * Unchanged on exit.
84 *
85 * A (input) DOUBLE PRECISION array, dimension (LDA,ka)
86 * where KA
87 * is K when TRANS = 'N' or 'n', and is N otherwise. Before
88 * entry with TRANS = 'N' or 'n', the leading N--by--K part of
89 * the array A must contain the matrix A, otherwise the leading
90 * K--by--N part of the array A must contain the matrix A.
91 * Unchanged on exit.
92 *
93 * LDA (input) INTEGER
94 * On entry, LDA specifies the first dimension of A as declared
95 * in the calling (sub) program. When TRANS = 'N' or 'n'
96 * then LDA must be at least max( 1, n ), otherwise LDA must
97 * be at least max( 1, k ).
98 * Unchanged on exit.
99 *
100 * BETA (input) DOUBLE PRECISION
101 * On entry, BETA specifies the scalar beta.
102 * Unchanged on exit.
103 *
104 *
105 * C (input/output) DOUBLE PRECISION array, dimension (NT)
106 * NT = N*(N+1)/2. On entry, the symmetric matrix C in RFP
107 * Format. RFP Format is described by TRANSR, UPLO and N.
108 *
109 * =====================================================================
110 *
111 * ..
112 * .. Parameters ..
113 DOUBLE PRECISION ONE, ZERO
114 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
115 * ..
116 * .. Local Scalars ..
117 LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS
118 INTEGER INFO, NROWA, J, NK, N1, N2
119 * ..
120 * .. External Functions ..
121 LOGICAL LSAME
122 EXTERNAL LSAME
123 * ..
124 * .. External Subroutines ..
125 EXTERNAL XERBLA, DGEMM, DSYRK
126 * ..
127 * .. Intrinsic Functions ..
128 INTRINSIC MAX
129 * ..
130 * .. Executable Statements ..
131 *
132 * Test the input parameters.
133 *
134 INFO = 0
135 NORMALTRANSR = LSAME( TRANSR, 'N' )
136 LOWER = LSAME( UPLO, 'L' )
137 NOTRANS = LSAME( TRANS, 'N' )
138 *
139 IF( NOTRANS ) THEN
140 NROWA = N
141 ELSE
142 NROWA = K
143 END IF
144 *
145 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
146 INFO = -1
147 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
148 INFO = -2
149 ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
150 INFO = -3
151 ELSE IF( N.LT.0 ) THEN
152 INFO = -4
153 ELSE IF( K.LT.0 ) THEN
154 INFO = -5
155 ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
156 INFO = -8
157 END IF
158 IF( INFO.NE.0 ) THEN
159 CALL XERBLA( 'DSFRK ', -INFO )
160 RETURN
161 END IF
162 *
163 * Quick return if possible.
164 *
165 * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
166 * done (it is in DSYRK for example) and left in the general case.
167 *
168 IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
169 $ ( BETA.EQ.ONE ) ) )RETURN
170 *
171 IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
172 DO J = 1, ( ( N*( N+1 ) ) / 2 )
173 C( J ) = ZERO
174 END DO
175 RETURN
176 END IF
177 *
178 * C is N-by-N.
179 * If N is odd, set NISODD = .TRUE., and N1 and N2.
180 * If N is even, NISODD = .FALSE., and NK.
181 *
182 IF( MOD( N, 2 ).EQ.0 ) THEN
183 NISODD = .FALSE.
184 NK = N / 2
185 ELSE
186 NISODD = .TRUE.
187 IF( LOWER ) THEN
188 N2 = N / 2
189 N1 = N - N2
190 ELSE
191 N1 = N / 2
192 N2 = N - N1
193 END IF
194 END IF
195 *
196 IF( NISODD ) THEN
197 *
198 * N is odd
199 *
200 IF( NORMALTRANSR ) THEN
201 *
202 * N is odd and TRANSR = 'N'
203 *
204 IF( LOWER ) THEN
205 *
206 * N is odd, TRANSR = 'N', and UPLO = 'L'
207 *
208 IF( NOTRANS ) THEN
209 *
210 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
211 *
212 CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
213 $ BETA, C( 1 ), N )
214 CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
215 $ BETA, C( N+1 ), N )
216 CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
217 $ LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
218 *
219 ELSE
220 *
221 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
222 *
223 CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
224 $ BETA, C( 1 ), N )
225 CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
226 $ BETA, C( N+1 ), N )
227 CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
228 $ LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
229 *
230 END IF
231 *
232 ELSE
233 *
234 * N is odd, TRANSR = 'N', and UPLO = 'U'
235 *
236 IF( NOTRANS ) THEN
237 *
238 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
239 *
240 CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
241 $ BETA, C( N2+1 ), N )
242 CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
243 $ BETA, C( N1+1 ), N )
244 CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
245 $ LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N )
246 *
247 ELSE
248 *
249 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
250 *
251 CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
252 $ BETA, C( N2+1 ), N )
253 CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA,
254 $ BETA, C( N1+1 ), N )
255 CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
256 $ LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N )
257 *
258 END IF
259 *
260 END IF
261 *
262 ELSE
263 *
264 * N is odd, and TRANSR = 'T'
265 *
266 IF( LOWER ) THEN
267 *
268 * N is odd, TRANSR = 'T', and UPLO = 'L'
269 *
270 IF( NOTRANS ) THEN
271 *
272 * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
273 *
274 CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
275 $ BETA, C( 1 ), N1 )
276 CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
277 $ BETA, C( 2 ), N1 )
278 CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
279 $ LDA, A( N1+1, 1 ), LDA, BETA,
280 $ C( N1*N1+1 ), N1 )
281 *
282 ELSE
283 *
284 * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
285 *
286 CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
287 $ BETA, C( 1 ), N1 )
288 CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
289 $ BETA, C( 2 ), N1 )
290 CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
291 $ LDA, A( 1, N1+1 ), LDA, BETA,
292 $ C( N1*N1+1 ), N1 )
293 *
294 END IF
295 *
296 ELSE
297 *
298 * N is odd, TRANSR = 'T', and UPLO = 'U'
299 *
300 IF( NOTRANS ) THEN
301 *
302 * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
303 *
304 CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
305 $ BETA, C( N2*N2+1 ), N2 )
306 CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
307 $ BETA, C( N1*N2+1 ), N2 )
308 CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
309 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
310 *
311 ELSE
312 *
313 * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
314 *
315 CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
316 $ BETA, C( N2*N2+1 ), N2 )
317 CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
318 $ BETA, C( N1*N2+1 ), N2 )
319 CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
320 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
321 *
322 END IF
323 *
324 END IF
325 *
326 END IF
327 *
328 ELSE
329 *
330 * N is even
331 *
332 IF( NORMALTRANSR ) THEN
333 *
334 * N is even and TRANSR = 'N'
335 *
336 IF( LOWER ) THEN
337 *
338 * N is even, TRANSR = 'N', and UPLO = 'L'
339 *
340 IF( NOTRANS ) THEN
341 *
342 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
343 *
344 CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
345 $ BETA, C( 2 ), N+1 )
346 CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
347 $ BETA, C( 1 ), N+1 )
348 CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
349 $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
350 $ N+1 )
351 *
352 ELSE
353 *
354 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
355 *
356 CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
357 $ BETA, C( 2 ), N+1 )
358 CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
359 $ BETA, C( 1 ), N+1 )
360 CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
361 $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
362 $ N+1 )
363 *
364 END IF
365 *
366 ELSE
367 *
368 * N is even, TRANSR = 'N', and UPLO = 'U'
369 *
370 IF( NOTRANS ) THEN
371 *
372 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
373 *
374 CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
375 $ BETA, C( NK+2 ), N+1 )
376 CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
377 $ BETA, C( NK+1 ), N+1 )
378 CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
379 $ LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ),
380 $ N+1 )
381 *
382 ELSE
383 *
384 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
385 *
386 CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
387 $ BETA, C( NK+2 ), N+1 )
388 CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
389 $ BETA, C( NK+1 ), N+1 )
390 CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
391 $ LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ),
392 $ N+1 )
393 *
394 END IF
395 *
396 END IF
397 *
398 ELSE
399 *
400 * N is even, and TRANSR = 'T'
401 *
402 IF( LOWER ) THEN
403 *
404 * N is even, TRANSR = 'T', and UPLO = 'L'
405 *
406 IF( NOTRANS ) THEN
407 *
408 * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
409 *
410 CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
411 $ BETA, C( NK+1 ), NK )
412 CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
413 $ BETA, C( 1 ), NK )
414 CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
415 $ LDA, A( NK+1, 1 ), LDA, BETA,
416 $ C( ( ( NK+1 )*NK )+1 ), NK )
417 *
418 ELSE
419 *
420 * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
421 *
422 CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
423 $ BETA, C( NK+1 ), NK )
424 CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
425 $ BETA, C( 1 ), NK )
426 CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
427 $ LDA, A( 1, NK+1 ), LDA, BETA,
428 $ C( ( ( NK+1 )*NK )+1 ), NK )
429 *
430 END IF
431 *
432 ELSE
433 *
434 * N is even, TRANSR = 'T', and UPLO = 'U'
435 *
436 IF( NOTRANS ) THEN
437 *
438 * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
439 *
440 CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
441 $ BETA, C( NK*( NK+1 )+1 ), NK )
442 CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
443 $ BETA, C( NK*NK+1 ), NK )
444 CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
445 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
446 *
447 ELSE
448 *
449 * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
450 *
451 CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
452 $ BETA, C( NK*( NK+1 )+1 ), NK )
453 CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
454 $ BETA, C( NK*NK+1 ), NK )
455 CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
456 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
457 *
458 END IF
459 *
460 END IF
461 *
462 END IF
463 *
464 END IF
465 *
466 RETURN
467 *
468 * End of DSFRK
469 *
470 END