Annotation of rpl/lapack/lapack/dsfrk.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
! 2: + C )
! 3: *
! 4: * -- LAPACK routine (version 3.2.2) --
! 5: *
! 6: * -- Contributed by Julien Langou of the Univ. of Colorado Denver --
! 7: * -- June 2010 --
! 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' + beta*C,
! 30: *
! 31: * or
! 32: *
! 33: * C := alpha*A'*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
! 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
! 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
! 60: * On entry, TRANS specifies the operation to be performed as
! 61: * follows:
! 62: *
! 63: * TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
! 64: *
! 65: * TRANS = 'T' or 't' C := alpha*A'*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: * Arguments
! 110: * ==========
! 111: *
! 112: * ..
! 113: * .. Parameters ..
! 114: DOUBLE PRECISION ONE, ZERO
! 115: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! 116: * ..
! 117: * .. Local Scalars ..
! 118: LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS
! 119: INTEGER INFO, NROWA, J, NK, N1, N2
! 120: * ..
! 121: * .. External Functions ..
! 122: LOGICAL LSAME
! 123: EXTERNAL LSAME
! 124: * ..
! 125: * .. External Subroutines ..
! 126: EXTERNAL XERBLA, DGEMM, DSYRK
! 127: * ..
! 128: * .. Intrinsic Functions ..
! 129: INTRINSIC MAX
! 130: * ..
! 131: * .. Executable Statements ..
! 132: *
! 133: * Test the input parameters.
! 134: *
! 135: INFO = 0
! 136: NORMALTRANSR = LSAME( TRANSR, 'N' )
! 137: LOWER = LSAME( UPLO, 'L' )
! 138: NOTRANS = LSAME( TRANS, 'N' )
! 139: *
! 140: IF( NOTRANS ) THEN
! 141: NROWA = N
! 142: ELSE
! 143: NROWA = K
! 144: END IF
! 145: *
! 146: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
! 147: INFO = -1
! 148: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
! 149: INFO = -2
! 150: ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
! 151: INFO = -3
! 152: ELSE IF( N.LT.0 ) THEN
! 153: INFO = -4
! 154: ELSE IF( K.LT.0 ) THEN
! 155: INFO = -5
! 156: ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
! 157: INFO = -8
! 158: END IF
! 159: IF( INFO.NE.0 ) THEN
! 160: CALL XERBLA( 'DSFRK ', -INFO )
! 161: RETURN
! 162: END IF
! 163: *
! 164: * Quick return if possible.
! 165: *
! 166: * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
! 167: * done (it is in DSYRK for example) and left in the general case.
! 168: *
! 169: IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
! 170: + ( BETA.EQ.ONE ) ) )RETURN
! 171: *
! 172: IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
! 173: DO J = 1, ( ( N*( N+1 ) ) / 2 )
! 174: C( J ) = ZERO
! 175: END DO
! 176: RETURN
! 177: END IF
! 178: *
! 179: * C is N-by-N.
! 180: * If N is odd, set NISODD = .TRUE., and N1 and N2.
! 181: * If N is even, NISODD = .FALSE., and NK.
! 182: *
! 183: IF( MOD( N, 2 ).EQ.0 ) THEN
! 184: NISODD = .FALSE.
! 185: NK = N / 2
! 186: ELSE
! 187: NISODD = .TRUE.
! 188: IF( LOWER ) THEN
! 189: N2 = N / 2
! 190: N1 = N - N2
! 191: ELSE
! 192: N1 = N / 2
! 193: N2 = N - N1
! 194: END IF
! 195: END IF
! 196: *
! 197: IF( NISODD ) THEN
! 198: *
! 199: * N is odd
! 200: *
! 201: IF( NORMALTRANSR ) THEN
! 202: *
! 203: * N is odd and TRANSR = 'N'
! 204: *
! 205: IF( LOWER ) THEN
! 206: *
! 207: * N is odd, TRANSR = 'N', and UPLO = 'L'
! 208: *
! 209: IF( NOTRANS ) THEN
! 210: *
! 211: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
! 212: *
! 213: CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
! 214: + BETA, C( 1 ), N )
! 215: CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
! 216: + BETA, C( N+1 ), N )
! 217: CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
! 218: + LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
! 219: *
! 220: ELSE
! 221: *
! 222: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
! 223: *
! 224: CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
! 225: + BETA, C( 1 ), N )
! 226: CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
! 227: + BETA, C( N+1 ), N )
! 228: CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
! 229: + LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
! 230: *
! 231: END IF
! 232: *
! 233: ELSE
! 234: *
! 235: * N is odd, TRANSR = 'N', and UPLO = 'U'
! 236: *
! 237: IF( NOTRANS ) THEN
! 238: *
! 239: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
! 240: *
! 241: CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
! 242: + BETA, C( N2+1 ), N )
! 243: CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
! 244: + BETA, C( N1+1 ), N )
! 245: CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
! 246: + LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N )
! 247: *
! 248: ELSE
! 249: *
! 250: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
! 251: *
! 252: CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
! 253: + BETA, C( N2+1 ), N )
! 254: CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA,
! 255: + BETA, C( N1+1 ), N )
! 256: CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
! 257: + LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N )
! 258: *
! 259: END IF
! 260: *
! 261: END IF
! 262: *
! 263: ELSE
! 264: *
! 265: * N is odd, and TRANSR = 'T'
! 266: *
! 267: IF( LOWER ) THEN
! 268: *
! 269: * N is odd, TRANSR = 'T', and UPLO = 'L'
! 270: *
! 271: IF( NOTRANS ) THEN
! 272: *
! 273: * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
! 274: *
! 275: CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
! 276: + BETA, C( 1 ), N1 )
! 277: CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
! 278: + BETA, C( 2 ), N1 )
! 279: CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
! 280: + LDA, A( N1+1, 1 ), LDA, BETA,
! 281: + C( N1*N1+1 ), N1 )
! 282: *
! 283: ELSE
! 284: *
! 285: * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
! 286: *
! 287: CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
! 288: + BETA, C( 1 ), N1 )
! 289: CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
! 290: + BETA, C( 2 ), N1 )
! 291: CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
! 292: + LDA, A( 1, N1+1 ), LDA, BETA,
! 293: + C( N1*N1+1 ), N1 )
! 294: *
! 295: END IF
! 296: *
! 297: ELSE
! 298: *
! 299: * N is odd, TRANSR = 'T', and UPLO = 'U'
! 300: *
! 301: IF( NOTRANS ) THEN
! 302: *
! 303: * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
! 304: *
! 305: CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
! 306: + BETA, C( N2*N2+1 ), N2 )
! 307: CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
! 308: + BETA, C( N1*N2+1 ), N2 )
! 309: CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
! 310: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
! 311: *
! 312: ELSE
! 313: *
! 314: * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
! 315: *
! 316: CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
! 317: + BETA, C( N2*N2+1 ), N2 )
! 318: CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
! 319: + BETA, C( N1*N2+1 ), N2 )
! 320: CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
! 321: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
! 322: *
! 323: END IF
! 324: *
! 325: END IF
! 326: *
! 327: END IF
! 328: *
! 329: ELSE
! 330: *
! 331: * N is even
! 332: *
! 333: IF( NORMALTRANSR ) THEN
! 334: *
! 335: * N is even and TRANSR = 'N'
! 336: *
! 337: IF( LOWER ) THEN
! 338: *
! 339: * N is even, TRANSR = 'N', and UPLO = 'L'
! 340: *
! 341: IF( NOTRANS ) THEN
! 342: *
! 343: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
! 344: *
! 345: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
! 346: + BETA, C( 2 ), N+1 )
! 347: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
! 348: + BETA, C( 1 ), N+1 )
! 349: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
! 350: + LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
! 351: + N+1 )
! 352: *
! 353: ELSE
! 354: *
! 355: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
! 356: *
! 357: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
! 358: + BETA, C( 2 ), N+1 )
! 359: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
! 360: + BETA, C( 1 ), N+1 )
! 361: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
! 362: + LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
! 363: + N+1 )
! 364: *
! 365: END IF
! 366: *
! 367: ELSE
! 368: *
! 369: * N is even, TRANSR = 'N', and UPLO = 'U'
! 370: *
! 371: IF( NOTRANS ) THEN
! 372: *
! 373: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
! 374: *
! 375: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
! 376: + BETA, C( NK+2 ), N+1 )
! 377: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
! 378: + BETA, C( NK+1 ), N+1 )
! 379: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
! 380: + LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ),
! 381: + N+1 )
! 382: *
! 383: ELSE
! 384: *
! 385: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
! 386: *
! 387: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
! 388: + BETA, C( NK+2 ), N+1 )
! 389: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
! 390: + BETA, C( NK+1 ), N+1 )
! 391: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
! 392: + LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ),
! 393: + N+1 )
! 394: *
! 395: END IF
! 396: *
! 397: END IF
! 398: *
! 399: ELSE
! 400: *
! 401: * N is even, and TRANSR = 'T'
! 402: *
! 403: IF( LOWER ) THEN
! 404: *
! 405: * N is even, TRANSR = 'T', and UPLO = 'L'
! 406: *
! 407: IF( NOTRANS ) THEN
! 408: *
! 409: * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
! 410: *
! 411: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
! 412: + BETA, C( NK+1 ), NK )
! 413: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
! 414: + BETA, C( 1 ), NK )
! 415: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
! 416: + LDA, A( NK+1, 1 ), LDA, BETA,
! 417: + C( ( ( NK+1 )*NK )+1 ), NK )
! 418: *
! 419: ELSE
! 420: *
! 421: * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
! 422: *
! 423: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
! 424: + BETA, C( NK+1 ), NK )
! 425: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
! 426: + BETA, C( 1 ), NK )
! 427: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
! 428: + LDA, A( 1, NK+1 ), LDA, BETA,
! 429: + C( ( ( NK+1 )*NK )+1 ), NK )
! 430: *
! 431: END IF
! 432: *
! 433: ELSE
! 434: *
! 435: * N is even, TRANSR = 'T', and UPLO = 'U'
! 436: *
! 437: IF( NOTRANS ) THEN
! 438: *
! 439: * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
! 440: *
! 441: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
! 442: + BETA, C( NK*( NK+1 )+1 ), NK )
! 443: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
! 444: + BETA, C( NK*NK+1 ), NK )
! 445: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
! 446: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
! 447: *
! 448: ELSE
! 449: *
! 450: * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
! 451: *
! 452: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
! 453: + BETA, C( NK*( NK+1 )+1 ), NK )
! 454: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
! 455: + BETA, C( NK*NK+1 ), NK )
! 456: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
! 457: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
! 458: *
! 459: END IF
! 460: *
! 461: END IF
! 462: *
! 463: END IF
! 464: *
! 465: END IF
! 466: *
! 467: RETURN
! 468: *
! 469: * End of DSFRK
! 470: *
! 471: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsfrk.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack