LAPACK 3.3.1
Linear Algebra PACKage

chfrk.f

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