LAPACK 3.3.1
Linear Algebra PACKage

ssfrk.f

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