LAPACK: SRC/dsfrk.f Source File

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.0)                                    --
00005 *
00006 *  -- Contributed by Julien Langou of the Univ. of Colorado Denver    --
00007 *     November 2010
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' + beta*C,
00030 *
00031 *  or
00032 *
00033 *     C := alpha*A'*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' + beta*C.
00064 *
00065 *              TRANS = 'T' or 't'   C := alpha*A'*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 *  Arguments
00110 *  ==========
00111 *
00112 *     ..
00113 *     .. Parameters ..
00114       DOUBLE PRECISION   ONE, ZERO
00115       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00116 *     ..
00117 *     .. Local Scalars ..
00118       LOGICAL            LOWER, NORMALTRANSR, NISODD, NOTRANS
00119       INTEGER            INFO, NROWA, J, NK, N1, N2
00120 *     ..
00121 *     .. External Functions ..
00122       LOGICAL            LSAME
00123       EXTERNAL           LSAME
00124 *     ..
00125 *     .. External Subroutines ..
00126       EXTERNAL           XERBLA, DGEMM, DSYRK
00127 *     ..
00128 *     .. Intrinsic Functions ..
00129       INTRINSIC          MAX
00130 *     ..
00131 *     .. Executable Statements ..
00132 *
00133 *     Test the input parameters.
00134 *
00135       INFO = 0
00136       NORMALTRANSR = LSAME( TRANSR, 'N' )
00137       LOWER = LSAME( UPLO, 'L' )
00138       NOTRANS = LSAME( TRANS, 'N' )
00139 *
00140       IF( NOTRANS ) THEN
00141          NROWA = N
00142       ELSE
00143          NROWA = K
00144       END IF
00145 *
00146       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
00147          INFO = -1
00148       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
00149          INFO = -2
00150       ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
00151          INFO = -3
00152       ELSE IF( N.LT.0 ) THEN
00153          INFO = -4
00154       ELSE IF( K.LT.0 ) THEN
00155          INFO = -5
00156       ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
00157          INFO = -8
00158       END IF
00159       IF( INFO.NE.0 ) THEN
00160          CALL XERBLA( 'DSFRK ', -INFO )
00161          RETURN
00162       END IF
00163 *
00164 *     Quick return if possible.
00165 *
00166 *     The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
00167 *     done (it is in DSYRK for example) and left in the general case.
00168 *
00169       IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
00170      +    ( BETA.EQ.ONE ) ) )RETURN
00171 *
00172       IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
00173          DO J = 1, ( ( N*( N+1 ) ) / 2 )
00174             C( J ) = ZERO
00175          END DO
00176          RETURN
00177       END IF
00178 *
00179 *     C is N-by-N.
00180 *     If N is odd, set NISODD = .TRUE., and N1 and N2.
00181 *     If N is even, NISODD = .FALSE., and NK.
00182 *
00183       IF( MOD( N, 2 ).EQ.0 ) THEN
00184          NISODD = .FALSE.
00185          NK = N / 2
00186       ELSE
00187          NISODD = .TRUE.
00188          IF( LOWER ) THEN
00189             N2 = N / 2
00190             N1 = N - N2
00191          ELSE
00192             N1 = N / 2
00193             N2 = N - N1
00194          END IF
00195       END IF
00196 *
00197       IF( NISODD ) THEN
00198 *
00199 *        N is odd
00200 *
00201          IF( NORMALTRANSR ) THEN
00202 *
00203 *           N is odd and TRANSR = 'N'
00204 *
00205             IF( LOWER ) THEN
00206 *
00207 *              N is odd, TRANSR = 'N', and UPLO = 'L'
00208 *
00209                IF( NOTRANS ) THEN
00210 *
00211 *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
00212 *
00213                   CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
00214      +                        BETA, C( 1 ), N )
00215                   CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
00216      +                        BETA, C( N+1 ), N )
00217                   CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
00218      +                        LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
00219 *
00220                ELSE
00221 *
00222 *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
00223 *
00224                   CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
00225      +                        BETA, C( 1 ), N )
00226                   CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
00227      +                        BETA, C( N+1 ), N )
00228                   CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
00229      +                        LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
00230 *
00231                END IF
00232 *
00233             ELSE
00234 *
00235 *              N is odd, TRANSR = 'N', and UPLO = 'U'
00236 *
00237                IF( NOTRANS ) THEN
00238 *
00239 *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
00240 *
00241                   CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
00242      +                        BETA, C( N2+1 ), N )
00243                   CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
00244      +                        BETA, C( N1+1 ), N )
00245                   CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
00246      +                        LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N )
00247 *
00248                ELSE
00249 *
00250 *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
00251 *
00252                   CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
00253      +                        BETA, C( N2+1 ), N )
00254                   CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA,
00255      +                        BETA, C( N1+1 ), N )
00256                   CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
00257      +                        LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N )
00258 *
00259                END IF
00260 *
00261             END IF
00262 *
00263          ELSE
00264 *
00265 *           N is odd, and TRANSR = 'T'
00266 *
00267             IF( LOWER ) THEN
00268 *
00269 *              N is odd, TRANSR = 'T', and UPLO = 'L'
00270 *
00271                IF( NOTRANS ) THEN
00272 *
00273 *                 N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
00274 *
00275                   CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
00276      +                        BETA, C( 1 ), N1 )
00277                   CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
00278      +                        BETA, C( 2 ), N1 )
00279                   CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
00280      +                        LDA, A( N1+1, 1 ), LDA, BETA,
00281      +                        C( N1*N1+1 ), N1 )
00282 *
00283                ELSE
00284 *
00285 *                 N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
00286 *
00287                   CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
00288      +                        BETA, C( 1 ), N1 )
00289                   CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
00290      +                        BETA, C( 2 ), N1 )
00291                   CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
00292      +                        LDA, A( 1, N1+1 ), LDA, BETA,
00293      +                        C( N1*N1+1 ), N1 )
00294 *
00295                END IF
00296 *
00297             ELSE
00298 *
00299 *              N is odd, TRANSR = 'T', and UPLO = 'U'
00300 *
00301                IF( NOTRANS ) THEN
00302 *
00303 *                 N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
00304 *
00305                   CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
00306      +                        BETA, C( N2*N2+1 ), N2 )
00307                   CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
00308      +                        BETA, C( N1*N2+1 ), N2 )
00309                   CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
00310      +                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
00311 *
00312                ELSE
00313 *
00314 *                 N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
00315 *
00316                   CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
00317      +                        BETA, C( N2*N2+1 ), N2 )
00318                   CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
00319      +                        BETA, C( N1*N2+1 ), N2 )
00320                   CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
00321      +                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
00322 *
00323                END IF
00324 *
00325             END IF
00326 *
00327          END IF
00328 *
00329       ELSE
00330 *
00331 *        N is even
00332 *
00333          IF( NORMALTRANSR ) THEN
00334 *
00335 *           N is even and TRANSR = 'N'
00336 *
00337             IF( LOWER ) THEN
00338 *
00339 *              N is even, TRANSR = 'N', and UPLO = 'L'
00340 *
00341                IF( NOTRANS ) THEN
00342 *
00343 *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
00344 *
00345                   CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
00346      +                        BETA, C( 2 ), N+1 )
00347                   CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
00348      +                        BETA, C( 1 ), N+1 )
00349                   CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
00350      +                        LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
00351      +                        N+1 )
00352 *
00353                ELSE
00354 *
00355 *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
00356 *
00357                   CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
00358      +                        BETA, C( 2 ), N+1 )
00359                   CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
00360      +                        BETA, C( 1 ), N+1 )
00361                   CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
00362      +                        LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
00363      +                        N+1 )
00364 *
00365                END IF
00366 *
00367             ELSE
00368 *
00369 *              N is even, TRANSR = 'N', and UPLO = 'U'
00370 *
00371                IF( NOTRANS ) THEN
00372 *
00373 *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
00374 *
00375                   CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
00376      +                        BETA, C( NK+2 ), N+1 )
00377                   CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
00378      +                        BETA, C( NK+1 ), N+1 )
00379                   CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
00380      +                        LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ),
00381      +                        N+1 )
00382 *
00383                ELSE
00384 *
00385 *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
00386 *
00387                   CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
00388      +                        BETA, C( NK+2 ), N+1 )
00389                   CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
00390      +                        BETA, C( NK+1 ), N+1 )
00391                   CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
00392      +                        LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ),
00393      +                        N+1 )
00394 *
00395                END IF
00396 *
00397             END IF
00398 *
00399          ELSE
00400 *
00401 *           N is even, and TRANSR = 'T'
00402 *
00403             IF( LOWER ) THEN
00404 *
00405 *              N is even, TRANSR = 'T', and UPLO = 'L'
00406 *
00407                IF( NOTRANS ) THEN
00408 *
00409 *                 N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
00410 *
00411                   CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
00412      +                        BETA, C( NK+1 ), NK )
00413                   CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
00414      +                        BETA, C( 1 ), NK )
00415                   CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
00416      +                        LDA, A( NK+1, 1 ), LDA, BETA,
00417      +                        C( ( ( NK+1 )*NK )+1 ), NK )
00418 *
00419                ELSE
00420 *
00421 *                 N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
00422 *
00423                   CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
00424      +                        BETA, C( NK+1 ), NK )
00425                   CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
00426      +                        BETA, C( 1 ), NK )
00427                   CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
00428      +                        LDA, A( 1, NK+1 ), LDA, BETA,
00429      +                        C( ( ( NK+1 )*NK )+1 ), NK )
00430 *
00431                END IF
00432 *
00433             ELSE
00434 *
00435 *              N is even, TRANSR = 'T', and UPLO = 'U'
00436 *
00437                IF( NOTRANS ) THEN
00438 *
00439 *                 N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
00440 *
00441                   CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
00442      +                        BETA, C( NK*( NK+1 )+1 ), NK )
00443                   CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
00444      +                        BETA, C( NK*NK+1 ), NK )
00445                   CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
00446      +                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
00447 *
00448                ELSE
00449 *
00450 *                 N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
00451 *
00452                   CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
00453      +                        BETA, C( NK*( NK+1 )+1 ), NK )
00454                   CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
00455      +                        BETA, C( NK*NK+1 ), NK )
00456                   CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
00457      +                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
00458 *
00459                END IF
00460 *
00461             END IF
00462 *
00463          END IF
00464 *
00465       END IF
00466 *
00467       RETURN
00468 *
00469 *     End of DSFRK
00470 *
00471       END