LAPACK: SRC/chfrk.f Source File

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