LAPACK 3.3.1
Linear Algebra PACKage

dsfrk.f

Go to the documentation of this file.
00001       SUBROUTINE DSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
00002      $                  C )
00003 *
00004 *  -- LAPACK routine (version 3.3.1)                                    --
00005 *
00006 *  -- Contributed by Julien Langou of the Univ. of Colorado Denver    --
00007 *  -- April 2011                                                      --
00008 *
00009 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00010 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00011 *
00012 *     ..
00013 *     .. Scalar Arguments ..
00014       DOUBLE PRECISION   ALPHA, BETA
00015       INTEGER            K, LDA, N
00016       CHARACTER          TRANS, TRANSR, UPLO
00017 *     ..
00018 *     .. Array Arguments ..
00019       DOUBLE PRECISION   A( LDA, * ), C( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  Level 3 BLAS like routine for C in RFP Format.
00026 *
00027 *  DSFRK performs one of the symmetric rank--k operations
00028 *
00029 *     C := alpha*A*A**T + beta*C,
00030 *
00031 *  or
00032 *
00033 *     C := alpha*A**T*A + beta*C,
00034 *
00035 *  where alpha and beta are real scalars, C is an n--by--n symmetric
00036 *  matrix and A is an n--by--k matrix in the first case and a k--by--n
00037 *  matrix in the second case.
00038 *
00039 *  Arguments
00040 *  ==========
00041 *
00042 *  TRANSR  (input) CHARACTER*1
00043 *          = 'N':  The Normal Form of RFP A is stored;
00044 *          = 'T':  The Transpose Form of RFP A is stored.
00045 *
00046 *  UPLO    (input) CHARACTER*1
00047 *           On  entry, UPLO specifies whether the upper or lower
00048 *           triangular part of the array C is to be referenced as
00049 *           follows:
00050 *
00051 *              UPLO = 'U' or 'u'   Only the upper triangular part of C
00052 *                                  is to be referenced.
00053 *
00054 *              UPLO = 'L' or 'l'   Only the lower triangular part of C
00055 *                                  is to be referenced.
00056 *
00057 *           Unchanged on exit.
00058 *
00059 *  TRANS   (input) CHARACTER*1
00060 *           On entry, TRANS specifies the operation to be performed as
00061 *           follows:
00062 *
00063 *              TRANS = 'N' or 'n'   C := alpha*A*A**T + beta*C.
00064 *
00065 *              TRANS = 'T' or 't'   C := alpha*A**T*A + beta*C.
00066 *
00067 *           Unchanged on exit.
00068 *
00069 *  N       (input) INTEGER
00070 *           On entry, N specifies the order of the matrix C. N must be
00071 *           at least zero.
00072 *           Unchanged on exit.
00073 *
00074 *  K       (input) INTEGER
00075 *           On entry with TRANS = 'N' or 'n', K specifies the number
00076 *           of  columns of the matrix A, and on entry with TRANS = 'T'
00077 *           or 't', K specifies the number of rows of the matrix A. K
00078 *           must be at least zero.
00079 *           Unchanged on exit.
00080 *
00081 *  ALPHA   (input) DOUBLE PRECISION
00082 *           On entry, ALPHA specifies the scalar alpha.
00083 *           Unchanged on exit.
00084 *
00085 *  A       (input) DOUBLE PRECISION array, dimension (LDA,ka)
00086 *           where KA
00087 *           is K  when TRANS = 'N' or 'n', and is N otherwise. Before
00088 *           entry with TRANS = 'N' or 'n', the leading N--by--K part of
00089 *           the array A must contain the matrix A, otherwise the leading
00090 *           K--by--N part of the array A must contain the matrix A.
00091 *           Unchanged on exit.
00092 *
00093 *  LDA     (input) INTEGER
00094 *           On entry, LDA specifies the first dimension of A as declared
00095 *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
00096 *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
00097 *           be at least  max( 1, k ).
00098 *           Unchanged on exit.
00099 *
00100 *  BETA    (input) DOUBLE PRECISION
00101 *           On entry, BETA specifies the scalar beta.
00102 *           Unchanged on exit.
00103 *
00104 *
00105 *  C       (input/output) DOUBLE PRECISION array, dimension (NT)
00106 *           NT = N*(N+1)/2. On entry, the symmetric matrix C in RFP
00107 *           Format. RFP Format is described by TRANSR, UPLO and N.
00108 *
00109 *  =====================================================================
00110 *
00111 *     ..
00112 *     .. Parameters ..
00113       DOUBLE PRECISION   ONE, ZERO
00114       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00115 *     ..
00116 *     .. Local Scalars ..
00117       LOGICAL            LOWER, NORMALTRANSR, NISODD, NOTRANS
00118       INTEGER            INFO, NROWA, J, NK, N1, N2
00119 *     ..
00120 *     .. External Functions ..
00121       LOGICAL            LSAME
00122       EXTERNAL           LSAME
00123 *     ..
00124 *     .. External Subroutines ..
00125       EXTERNAL           XERBLA, DGEMM, DSYRK
00126 *     ..
00127 *     .. Intrinsic Functions ..
00128       INTRINSIC          MAX
00129 *     ..
00130 *     .. Executable Statements ..
00131 *
00132 *     Test the input parameters.
00133 *
00134       INFO = 0
00135       NORMALTRANSR = LSAME( TRANSR, 'N' )
00136       LOWER = LSAME( UPLO, 'L' )
00137       NOTRANS = LSAME( TRANS, 'N' )
00138 *
00139       IF( NOTRANS ) THEN
00140          NROWA = N
00141       ELSE
00142          NROWA = K
00143       END IF
00144 *
00145       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
00146          INFO = -1
00147       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
00148          INFO = -2
00149       ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
00150          INFO = -3
00151       ELSE IF( N.LT.0 ) THEN
00152          INFO = -4
00153       ELSE IF( K.LT.0 ) THEN
00154          INFO = -5
00155       ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
00156          INFO = -8
00157       END IF
00158       IF( INFO.NE.0 ) THEN
00159          CALL XERBLA( 'DSFRK ', -INFO )
00160          RETURN
00161       END IF
00162 *
00163 *     Quick return if possible.
00164 *
00165 *     The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
00166 *     done (it is in DSYRK for example) and left in the general case.
00167 *
00168       IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
00169      $    ( BETA.EQ.ONE ) ) )RETURN
00170 *
00171       IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
00172          DO J = 1, ( ( N*( N+1 ) ) / 2 )
00173             C( J ) = ZERO
00174          END DO
00175          RETURN
00176       END IF
00177 *
00178 *     C is N-by-N.
00179 *     If N is odd, set NISODD = .TRUE., and N1 and N2.
00180 *     If N is even, NISODD = .FALSE., and NK.
00181 *
00182       IF( MOD( N, 2 ).EQ.0 ) THEN
00183          NISODD = .FALSE.
00184          NK = N / 2
00185       ELSE
00186          NISODD = .TRUE.
00187          IF( LOWER ) THEN
00188             N2 = N / 2
00189             N1 = N - N2
00190          ELSE
00191             N1 = N / 2
00192             N2 = N - N1
00193          END IF
00194       END IF
00195 *
00196       IF( NISODD ) THEN
00197 *
00198 *        N is odd
00199 *
00200          IF( NORMALTRANSR ) THEN
00201 *
00202 *           N is odd and TRANSR = 'N'
00203 *
00204             IF( LOWER ) THEN
00205 *
00206 *              N is odd, TRANSR = 'N', and UPLO = 'L'
00207 *
00208                IF( NOTRANS ) THEN
00209 *
00210 *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
00211 *
00212                   CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
00213      $                        BETA, C( 1 ), N )
00214                   CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
00215      $                        BETA, C( N+1 ), N )
00216                   CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
00217      $                        LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
00218 *
00219                ELSE
00220 *
00221 *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
00222 *
00223                   CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
00224      $                        BETA, C( 1 ), N )
00225                   CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
00226      $                        BETA, C( N+1 ), N )
00227                   CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
00228      $                        LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
00229 *
00230                END IF
00231 *
00232             ELSE
00233 *
00234 *              N is odd, TRANSR = 'N', and UPLO = 'U'
00235 *
00236                IF( NOTRANS ) THEN
00237 *
00238 *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
00239 *
00240                   CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
00241      $                        BETA, C( N2+1 ), N )
00242                   CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
00243      $                        BETA, C( N1+1 ), N )
00244                   CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
00245      $                        LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N )
00246 *
00247                ELSE
00248 *
00249 *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
00250 *
00251                   CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
00252      $                        BETA, C( N2+1 ), N )
00253                   CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA,
00254      $                        BETA, C( N1+1 ), N )
00255                   CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
00256      $                        LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N )
00257 *
00258                END IF
00259 *
00260             END IF
00261 *
00262          ELSE
00263 *
00264 *           N is odd, and TRANSR = 'T'
00265 *
00266             IF( LOWER ) THEN
00267 *
00268 *              N is odd, TRANSR = 'T', and UPLO = 'L'
00269 *
00270                IF( NOTRANS ) THEN
00271 *
00272 *                 N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
00273 *
00274                   CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
00275      $                        BETA, C( 1 ), N1 )
00276                   CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
00277      $                        BETA, C( 2 ), N1 )
00278                   CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
00279      $                        LDA, A( N1+1, 1 ), LDA, BETA,
00280      $                        C( N1*N1+1 ), N1 )
00281 *
00282                ELSE
00283 *
00284 *                 N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
00285 *
00286                   CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
00287      $                        BETA, C( 1 ), N1 )
00288                   CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
00289      $                        BETA, C( 2 ), N1 )
00290                   CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
00291      $                        LDA, A( 1, N1+1 ), LDA, BETA,
00292      $                        C( N1*N1+1 ), N1 )
00293 *
00294                END IF
00295 *
00296             ELSE
00297 *
00298 *              N is odd, TRANSR = 'T', and UPLO = 'U'
00299 *
00300                IF( NOTRANS ) THEN
00301 *
00302 *                 N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
00303 *
00304                   CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
00305      $                        BETA, C( N2*N2+1 ), N2 )
00306                   CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
00307      $                        BETA, C( N1*N2+1 ), N2 )
00308                   CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
00309      $                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
00310 *
00311                ELSE
00312 *
00313 *                 N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
00314 *
00315                   CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
00316      $                        BETA, C( N2*N2+1 ), N2 )
00317                   CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
00318      $                        BETA, C( N1*N2+1 ), N2 )
00319                   CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
00320      $                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
00321 *
00322                END IF
00323 *
00324             END IF
00325 *
00326          END IF
00327 *
00328       ELSE
00329 *
00330 *        N is even
00331 *
00332          IF( NORMALTRANSR ) THEN
00333 *
00334 *           N is even and TRANSR = 'N'
00335 *
00336             IF( LOWER ) THEN
00337 *
00338 *              N is even, TRANSR = 'N', and UPLO = 'L'
00339 *
00340                IF( NOTRANS ) THEN
00341 *
00342 *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
00343 *
00344                   CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
00345      $                        BETA, C( 2 ), N+1 )
00346                   CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
00347      $                        BETA, C( 1 ), N+1 )
00348                   CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
00349      $                        LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
00350      $                        N+1 )
00351 *
00352                ELSE
00353 *
00354 *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
00355 *
00356                   CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
00357      $                        BETA, C( 2 ), N+1 )
00358                   CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
00359      $                        BETA, C( 1 ), N+1 )
00360                   CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
00361      $                        LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
00362      $                        N+1 )
00363 *
00364                END IF
00365 *
00366             ELSE
00367 *
00368 *              N is even, TRANSR = 'N', and UPLO = 'U'
00369 *
00370                IF( NOTRANS ) THEN
00371 *
00372 *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
00373 *
00374                   CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
00375      $                        BETA, C( NK+2 ), N+1 )
00376                   CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
00377      $                        BETA, C( NK+1 ), N+1 )
00378                   CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
00379      $                        LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ),
00380      $                        N+1 )
00381 *
00382                ELSE
00383 *
00384 *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
00385 *
00386                   CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
00387      $                        BETA, C( NK+2 ), N+1 )
00388                   CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
00389      $                        BETA, C( NK+1 ), N+1 )
00390                   CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
00391      $                        LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ),
00392      $                        N+1 )
00393 *
00394                END IF
00395 *
00396             END IF
00397 *
00398          ELSE
00399 *
00400 *           N is even, and TRANSR = 'T'
00401 *
00402             IF( LOWER ) THEN
00403 *
00404 *              N is even, TRANSR = 'T', and UPLO = 'L'
00405 *
00406                IF( NOTRANS ) THEN
00407 *
00408 *                 N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
00409 *
00410                   CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
00411      $                        BETA, C( NK+1 ), NK )
00412                   CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
00413      $                        BETA, C( 1 ), NK )
00414                   CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
00415      $                        LDA, A( NK+1, 1 ), LDA, BETA,
00416      $                        C( ( ( NK+1 )*NK )+1 ), NK )
00417 *
00418                ELSE
00419 *
00420 *                 N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
00421 *
00422                   CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
00423      $                        BETA, C( NK+1 ), NK )
00424                   CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
00425      $                        BETA, C( 1 ), NK )
00426                   CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
00427      $                        LDA, A( 1, NK+1 ), LDA, BETA,
00428      $                        C( ( ( NK+1 )*NK )+1 ), NK )
00429 *
00430                END IF
00431 *
00432             ELSE
00433 *
00434 *              N is even, TRANSR = 'T', and UPLO = 'U'
00435 *
00436                IF( NOTRANS ) THEN
00437 *
00438 *                 N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
00439 *
00440                   CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
00441      $                        BETA, C( NK*( NK+1 )+1 ), NK )
00442                   CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
00443      $                        BETA, C( NK*NK+1 ), NK )
00444                   CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
00445      $                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
00446 *
00447                ELSE
00448 *
00449 *                 N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
00450 *
00451                   CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
00452      $                        BETA, C( NK*( NK+1 )+1 ), NK )
00453                   CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
00454      $                        BETA, C( NK*NK+1 ), NK )
00455                   CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
00456      $                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
00457 *
00458                END IF
00459 *
00460             END IF
00461 *
00462          END IF
00463 *
00464       END IF
00465 *
00466       RETURN
00467 *
00468 *     End of DSFRK
00469 *
00470       END
 All Files Functions