OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dsfrk.f
Go to the documentation of this file.
1*> \brief \b DSFRK performs a symmetric rank-k operation for matrix in RFP format.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSFRK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsfrk.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsfrk.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsfrk.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
22* C )
23*
24* .. Scalar Arguments ..
25* DOUBLE PRECISION ALPHA, BETA
26* INTEGER K, LDA, N
27* CHARACTER TRANS, TRANSR, UPLO
28* ..
29* .. Array Arguments ..
30* DOUBLE PRECISION A( LDA, * ), C( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> Level 3 BLAS like routine for C in RFP Format.
40*>
41*> DSFRK performs one of the symmetric rank--k operations
42*>
43*> C := alpha*A*A**T + beta*C,
44*>
45*> or
46*>
47*> C := alpha*A**T*A + beta*C,
48*>
49*> where alpha and beta are real scalars, C is an n--by--n symmetric
50*> matrix and A is an n--by--k matrix in the first case and a k--by--n
51*> matrix in the second case.
52*> \endverbatim
53*
54* Arguments:
55* ==========
56*
57*> \param[in] TRANSR
58*> \verbatim
59*> TRANSR is CHARACTER*1
60*> = 'N': The Normal Form of RFP A is stored;
61*> = 'T': The Transpose Form of RFP A is stored.
62*> \endverbatim
63*>
64*> \param[in] UPLO
65*> \verbatim
66*> UPLO is CHARACTER*1
67*> On entry, UPLO specifies whether the upper or lower
68*> triangular part of the array C is to be referenced as
69*> follows:
70*>
71*> UPLO = 'U' or 'u' Only the upper triangular part of C
72*> is to be referenced.
73*>
74*> UPLO = 'L' or 'l' Only the lower triangular part of C
75*> is to be referenced.
76*>
77*> Unchanged on exit.
78*> \endverbatim
79*>
80*> \param[in] TRANS
81*> \verbatim
82*> TRANS is CHARACTER*1
83*> On entry, TRANS specifies the operation to be performed as
84*> follows:
85*>
86*> TRANS = 'N' or 'n' C := alpha*A*A**T + beta*C.
87*>
88*> TRANS = 'T' or 't' C := alpha*A**T*A + beta*C.
89*>
90*> Unchanged on exit.
91*> \endverbatim
92*>
93*> \param[in] N
94*> \verbatim
95*> N is INTEGER
96*> On entry, N specifies the order of the matrix C. N must be
97*> at least zero.
98*> Unchanged on exit.
99*> \endverbatim
100*>
101*> \param[in] K
102*> \verbatim
103*> K is INTEGER
104*> On entry with TRANS = 'N' or 'n', K specifies the number
105*> of columns of the matrix A, and on entry with TRANS = 'T'
106*> or 't', K specifies the number of rows of the matrix A. K
107*> must be at least zero.
108*> Unchanged on exit.
109*> \endverbatim
110*>
111*> \param[in] ALPHA
112*> \verbatim
113*> ALPHA is DOUBLE PRECISION
114*> On entry, ALPHA specifies the scalar alpha.
115*> Unchanged on exit.
116*> \endverbatim
117*>
118*> \param[in] A
119*> \verbatim
120*> A is DOUBLE PRECISION array, dimension (LDA,ka)
121*> where KA
122*> is K when TRANS = 'N' or 'n', and is N otherwise. Before
123*> entry with TRANS = 'N' or 'n', the leading N--by--K part of
124*> the array A must contain the matrix A, otherwise the leading
125*> K--by--N part of the array A must contain the matrix A.
126*> Unchanged on exit.
127*> \endverbatim
128*>
129*> \param[in] LDA
130*> \verbatim
131*> LDA is INTEGER
132*> On entry, LDA specifies the first dimension of A as declared
133*> in the calling (sub) program. When TRANS = 'N' or 'n'
134*> then LDA must be at least max( 1, n ), otherwise LDA must
135*> be at least max( 1, k ).
136*> Unchanged on exit.
137*> \endverbatim
138*>
139*> \param[in] BETA
140*> \verbatim
141*> BETA is DOUBLE PRECISION
142*> On entry, BETA specifies the scalar beta.
143*> Unchanged on exit.
144*> \endverbatim
145*>
146*> \param[in,out] C
147*> \verbatim
148*> C is DOUBLE PRECISION array, dimension (NT)
149*> NT = N*(N+1)/2. On entry, the symmetric matrix C in RFP
150*> Format. RFP Format is described by TRANSR, UPLO and N.
151*> \endverbatim
152*
153* Authors:
154* ========
155*
156*> \author Univ. of Tennessee
157*> \author Univ. of California Berkeley
158*> \author Univ. of Colorado Denver
159*> \author NAG Ltd.
160*
161*> \ingroup doubleOTHERcomputational
162*
163* =====================================================================
164 SUBROUTINE dsfrk( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
165 $ C )
166*
167* -- LAPACK computational routine --
168* -- LAPACK is a software package provided by Univ. of Tennessee, --
169* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
170*
171* .. Scalar Arguments ..
172 DOUBLE PRECISION ALPHA, BETA
173 INTEGER K, LDA, N
174 CHARACTER TRANS, TRANSR, UPLO
175* ..
176* .. Array Arguments ..
177 DOUBLE PRECISION A( LDA, * ), C( * )
178* ..
179*
180* =====================================================================
181*
182* ..
183* .. Parameters ..
184 DOUBLE PRECISION ONE, ZERO
185 parameter( one = 1.0d+0, zero = 0.0d+0 )
186* ..
187* .. Local Scalars ..
188 LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS
189 INTEGER INFO, NROWA, J, NK, N1, N2
190* ..
191* .. External Functions ..
192 LOGICAL LSAME
193 EXTERNAL lsame
194* ..
195* .. External Subroutines ..
196 EXTERNAL xerbla, dgemm, dsyrk
197* ..
198* .. Intrinsic Functions ..
199 INTRINSIC max
200* ..
201* .. Executable Statements ..
202*
203* Test the input parameters.
204*
205 info = 0
206 normaltransr = lsame( transr, 'N' )
207 lower = lsame( uplo, 'L' )
208 notrans = lsame( trans, 'N' )
209*
210 IF( notrans ) THEN
211 nrowa = n
212 ELSE
213 nrowa = k
214 END IF
215*
216 IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'T' ) ) THEN
217 info = -1
218 ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
219 info = -2
220 ELSE IF( .NOT.notrans .AND. .NOT.lsame( trans, 'T' ) ) THEN
221 info = -3
222 ELSE IF( n.LT.0 ) THEN
223 info = -4
224 ELSE IF( k.LT.0 ) THEN
225 info = -5
226 ELSE IF( lda.LT.max( 1, nrowa ) ) THEN
227 info = -8
228 END IF
229 IF( info.NE.0 ) THEN
230 CALL xerbla( 'DSFRK ', -info )
231 RETURN
232 END IF
233*
234* Quick return if possible.
235*
236* The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
237* done (it is in DSYRK for example) and left in the general case.
238*
239 IF( ( n.EQ.0 ) .OR. ( ( ( alpha.EQ.zero ) .OR. ( k.EQ.0 ) ) .AND.
240 $ ( beta.EQ.one ) ) )RETURN
241*
242 IF( ( alpha.EQ.zero ) .AND. ( beta.EQ.zero ) ) THEN
243 DO j = 1, ( ( n*( n+1 ) ) / 2 )
244 c( j ) = zero
245 END DO
246 RETURN
247 END IF
248*
249* C is N-by-N.
250* If N is odd, set NISODD = .TRUE., and N1 and N2.
251* If N is even, NISODD = .FALSE., and NK.
252*
253 IF( mod( n, 2 ).EQ.0 ) THEN
254 nisodd = .false.
255 nk = n / 2
256 ELSE
257 nisodd = .true.
258 IF( lower ) THEN
259 n2 = n / 2
260 n1 = n - n2
261 ELSE
262 n1 = n / 2
263 n2 = n - n1
264 END IF
265 END IF
266*
267 IF( nisodd ) THEN
268*
269* N is odd
270*
271 IF( normaltransr ) THEN
272*
273* N is odd and TRANSR = 'N'
274*
275 IF( lower ) THEN
276*
277* N is odd, TRANSR = 'N', and UPLO = 'L'
278*
279 IF( notrans ) THEN
280*
281* N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
282*
283 CALL dsyrk( 'L', 'N', n1, k, alpha, a( 1, 1 ), lda,
284 $ beta, c( 1 ), n )
285 CALL dsyrk( 'U', 'N', n2, k, alpha, a( n1+1, 1 ), lda,
286 $ beta, c( n+1 ), n )
287 CALL dgemm( 'N', 'T', n2, n1, k, alpha, a( n1+1, 1 ),
288 $ lda, a( 1, 1 ), lda, beta, c( n1+1 ), n )
289*
290 ELSE
291*
292* N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
293*
294 CALL dsyrk( 'L', 'T', n1, k, alpha, a( 1, 1 ), lda,
295 $ beta, c( 1 ), n )
296 CALL dsyrk( 'U', 'T', n2, k, alpha, a( 1, n1+1 ), lda,
297 $ beta, c( n+1 ), n )
298 CALL dgemm( 'T', 'N', n2, n1, k, alpha, a( 1, n1+1 ),
299 $ lda, a( 1, 1 ), lda, beta, c( n1+1 ), n )
300*
301 END IF
302*
303 ELSE
304*
305* N is odd, TRANSR = 'N', and UPLO = 'U'
306*
307 IF( notrans ) THEN
308*
309* N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
310*
311 CALL dsyrk( 'L', 'N', n1, k, alpha, a( 1, 1 ), lda,
312 $ beta, c( n2+1 ), n )
313 CALL dsyrk( 'U', 'N', n2, k, alpha, a( n2, 1 ), lda,
314 $ beta, c( n1+1 ), n )
315 CALL dgemm( 'N', 'T', n1, n2, k, alpha, a( 1, 1 ),
316 $ lda, a( n2, 1 ), lda, beta, c( 1 ), n )
317*
318 ELSE
319*
320* N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
321*
322 CALL dsyrk( 'L', 'T', n1, k, alpha, a( 1, 1 ), lda,
323 $ beta, c( n2+1 ), n )
324 CALL dsyrk( 'U', 'T', n2, k, alpha, a( 1, n2 ), lda,
325 $ beta, c( n1+1 ), n )
326 CALL dgemm( 'T', 'N', n1, n2, k, alpha, a( 1, 1 ),
327 $ lda, a( 1, n2 ), lda, beta, c( 1 ), n )
328*
329 END IF
330*
331 END IF
332*
333 ELSE
334*
335* N is odd, and TRANSR = 'T'
336*
337 IF( lower ) THEN
338*
339* N is odd, TRANSR = 'T', and UPLO = 'L'
340*
341 IF( notrans ) THEN
342*
343* N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
344*
345 CALL dsyrk( 'U', 'N', n1, k, alpha, a( 1, 1 ), lda,
346 $ beta, c( 1 ), n1 )
347 CALL dsyrk( 'L', 'N', n2, k, alpha, a( n1+1, 1 ), lda,
348 $ beta, c( 2 ), n1 )
349 CALL dgemm( 'N', 'T', n1, n2, k, alpha, a( 1, 1 ),
350 $ lda, a( n1+1, 1 ), lda, beta,
351 $ c( n1*n1+1 ), n1 )
352*
353 ELSE
354*
355* N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
356*
357 CALL dsyrk( 'U', 'T', n1, k, alpha, a( 1, 1 ), lda,
358 $ beta, c( 1 ), n1 )
359 CALL dsyrk( 'L', 'T', n2, k, alpha, a( 1, n1+1 ), lda,
360 $ beta, c( 2 ), n1 )
361 CALL dgemm( 'T', 'n', N1, N2, K, ALPHA, A( 1, 1 ),
362 $ LDA, A( 1, N1+1 ), LDA, BETA,
363 $ C( N1*N1+1 ), N1 )
364*
365 END IF
366*
367 ELSE
368*
369* N is odd, TRANSR = 'T', and UPLO = 'U'
370*
371 IF( NOTRANS ) THEN
372*
373* N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
374*
375 CALL DSYRK( 'u', 'n', N1, K, ALPHA, A( 1, 1 ), LDA,
376 $ BETA, C( N2*N2+1 ), N2 )
377 CALL DSYRK( 'l', 'n', N2, K, ALPHA, A( N1+1, 1 ), LDA,
378 $ BETA, C( N1*N2+1 ), N2 )
379 CALL DGEMM( 'n', 't', N2, N1, K, ALPHA, A( N1+1, 1 ),
380 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
381*
382 ELSE
383*
384* N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
385*
386 CALL DSYRK( 'u', 't', N1, K, ALPHA, A( 1, 1 ), LDA,
387 $ BETA, C( N2*N2+1 ), N2 )
388 CALL DSYRK( 'l', 't', N2, K, ALPHA, A( 1, N1+1 ), LDA,
389 $ BETA, C( N1*N2+1 ), N2 )
390 CALL DGEMM( 't', 'n', N2, N1, K, ALPHA, A( 1, N1+1 ),
391 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
392*
393 END IF
394*
395 END IF
396*
397 END IF
398*
399 ELSE
400*
401* N is even
402*
403 IF( NORMALTRANSR ) THEN
404*
405* N is even and TRANSR = 'N'
406*
407 IF( LOWER ) THEN
408*
409* N is even, TRANSR = 'N', and UPLO = 'L'
410*
411 IF( NOTRANS ) THEN
412*
413* N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
414*
415 CALL DSYRK( 'l', 'n', NK, K, ALPHA, A( 1, 1 ), LDA,
416 $ BETA, C( 2 ), N+1 )
417 CALL DSYRK( 'u', 'n', NK, K, ALPHA, A( NK+1, 1 ), LDA,
418 $ BETA, C( 1 ), N+1 )
419 CALL DGEMM( 'n', 't', NK, NK, K, ALPHA, A( NK+1, 1 ),
420 $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
421 $ N+1 )
422*
423 ELSE
424*
425* N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
426*
427 CALL DSYRK( 'l', 't', NK, K, ALPHA, A( 1, 1 ), LDA,
428 $ BETA, C( 2 ), N+1 )
429 CALL DSYRK( 'u', 't', NK, K, ALPHA, A( 1, NK+1 ), LDA,
430 $ BETA, C( 1 ), N+1 )
431 CALL DGEMM( 't', 'n', NK, NK, K, ALPHA, A( 1, NK+1 ),
432 $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
433 $ N+1 )
434*
435 END IF
436*
437 ELSE
438*
439* N is even, TRANSR = 'N', and UPLO = 'U'
440*
441 IF( NOTRANS ) THEN
442*
443* N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
444*
445 CALL DSYRK( 'l', 'n', NK, K, ALPHA, A( 1, 1 ), LDA,
446 $ BETA, C( NK+2 ), N+1 )
447 CALL DSYRK( 'u', 'n', NK, K, ALPHA, A( NK+1, 1 ), LDA,
448 $ BETA, C( NK+1 ), N+1 )
449 CALL DGEMM( 'n', 't', NK, NK, K, ALPHA, A( 1, 1 ),
450 $ LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ),
451 $ N+1 )
452*
453 ELSE
454*
455* N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
456*
457 CALL DSYRK( 'l', 't', NK, K, ALPHA, A( 1, 1 ), LDA,
458 $ BETA, C( NK+2 ), N+1 )
459 CALL DSYRK( 'u', 't', NK, K, ALPHA, A( 1, NK+1 ), LDA,
460 $ BETA, C( NK+1 ), N+1 )
461 CALL DGEMM( 't', 'n', NK, NK, K, ALPHA, A( 1, 1 ),
462 $ LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ),
463 $ N+1 )
464*
465 END IF
466*
467 END IF
468*
469 ELSE
470*
471* N is even, and TRANSR = 'T'
472*
473 IF( LOWER ) THEN
474*
475* N is even, TRANSR = 'T', and UPLO = 'L'
476*
477 IF( NOTRANS ) THEN
478*
479* N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
480*
481 CALL DSYRK( 'u', 'n', NK, K, ALPHA, A( 1, 1 ), LDA,
482 $ BETA, C( NK+1 ), NK )
483 CALL DSYRK( 'l', 'n', NK, K, ALPHA, A( NK+1, 1 ), LDA,
484 $ BETA, C( 1 ), NK )
485 CALL DGEMM( 'n', 't', NK, NK, K, ALPHA, A( 1, 1 ),
486 $ LDA, A( NK+1, 1 ), LDA, BETA,
487 $ C( ( ( NK+1 )*NK )+1 ), NK )
488*
489 ELSE
490*
491* N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
492*
493 CALL DSYRK( 'u', 't', NK, K, ALPHA, A( 1, 1 ), LDA,
494 $ BETA, C( NK+1 ), NK )
495 CALL DSYRK( 'l', 't', NK, K, ALPHA, A( 1, NK+1 ), LDA,
496 $ BETA, C( 1 ), NK )
497 CALL DGEMM( 't', 'n', NK, NK, K, ALPHA, A( 1, 1 ),
498 $ LDA, A( 1, NK+1 ), LDA, BETA,
499 $ C( ( ( NK+1 )*NK )+1 ), NK )
500*
501 END IF
502*
503 ELSE
504*
505* N is even, TRANSR = 'T', and UPLO = 'U'
506*
507 IF( NOTRANS ) THEN
508*
509* N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
510*
511 CALL DSYRK( 'u', 'n', NK, K, ALPHA, A( 1, 1 ), LDA,
512 $ BETA, C( NK*( NK+1 )+1 ), NK )
513 CALL DSYRK( 'l', 'n', NK, K, ALPHA, A( NK+1, 1 ), LDA,
514 $ BETA, C( NK*NK+1 ), NK )
515 CALL DGEMM( 'n', 't', NK, NK, K, ALPHA, A( NK+1, 1 ),
516 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
517*
518 ELSE
519*
520* N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
521*
522 CALL DSYRK( 'u', 't', NK, K, ALPHA, A( 1, 1 ), LDA,
523 $ BETA, C( NK*( NK+1 )+1 ), NK )
524 CALL DSYRK( 'l', 't', NK, K, ALPHA, A( 1, NK+1 ), LDA,
525 $ BETA, C( NK*NK+1 ), NK )
526 CALL DGEMM( 't', 'n', NK, NK, K, ALPHA, A( 1, NK+1 ),
527 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
528*
529 END IF
530*
531 END IF
532*
533 END IF
534*
535 END IF
536*
537 RETURN
538*
539* End of DSFRK
540*
541 END
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
dsfrk
subroutine dsfrk(transr, uplo, trans, n, k, alpha, a, lda, beta, c)
DSFRK performs a symmetric rank-k operation for matrix in RFP format.
Definition dsfrk.f:166
dsyrk
subroutine dsyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
DSYRK
Definition dsyrk.f:169
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187
max
#define max(a, b)
Definition macros.h:21