LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cher2k.f
Go to the documentation of this file.
00001 *> \brief \b CHER2K
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE CHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
00012 * 
00013 *       .. Scalar Arguments ..
00014 *       COMPLEX ALPHA
00015 *       REAL BETA
00016 *       INTEGER K,LDA,LDB,LDC,N
00017 *       CHARACTER TRANS,UPLO
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> CHER2K  performs one of the hermitian rank 2k operations
00030 *>
00031 *>    C := alpha*A*B**H + conjg( alpha )*B*A**H + beta*C,
00032 *>
00033 *> or
00034 *>
00035 *>    C := alpha*A**H*B + conjg( alpha )*B**H*A + beta*C,
00036 *>
00037 *> where  alpha and beta  are scalars with  beta  real,  C is an  n by n
00038 *> hermitian matrix and  A and B  are  n by k matrices in the first case
00039 *> and  k by n  matrices in the second case.
00040 *> \endverbatim
00041 *
00042 *  Arguments:
00043 *  ==========
00044 *
00045 *> \param[in] UPLO
00046 *> \verbatim
00047 *>          UPLO is CHARACTER*1
00048 *>           On  entry,   UPLO  specifies  whether  the  upper  or  lower
00049 *>           triangular  part  of the  array  C  is to be  referenced  as
00050 *>           follows:
00051 *>
00052 *>              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
00053 *>                                  is to be referenced.
00054 *>
00055 *>              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
00056 *>                                  is to be referenced.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] TRANS
00060 *> \verbatim
00061 *>          TRANS is CHARACTER*1
00062 *>           On entry,  TRANS  specifies the operation to be performed as
00063 *>           follows:
00064 *>
00065 *>              TRANS = 'N' or 'n'    C := alpha*A*B**H          +
00066 *>                                         conjg( alpha )*B*A**H +
00067 *>                                         beta*C.
00068 *>
00069 *>              TRANS = 'C' or 'c'    C := alpha*A**H*B          +
00070 *>                                         conjg( alpha )*B**H*A +
00071 *>                                         beta*C.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] N
00075 *> \verbatim
00076 *>          N is INTEGER
00077 *>           On entry,  N specifies the order of the matrix C.  N must be
00078 *>           at least zero.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] K
00082 *> \verbatim
00083 *>          K is INTEGER
00084 *>           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
00085 *>           of  columns  of the  matrices  A and B,  and on  entry  with
00086 *>           TRANS = 'C' or 'c',  K  specifies  the number of rows of the
00087 *>           matrices  A and B.  K must be at least zero.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] ALPHA
00091 *> \verbatim
00092 *>          ALPHA is COMPLEX
00093 *>           On entry, ALPHA specifies the scalar alpha.
00094 *> \endverbatim
00095 *>
00096 *> \param[in] A
00097 *> \verbatim
00098 *>          A is COMPLEX array of DIMENSION ( LDA, ka ), where ka is
00099 *>           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
00100 *>           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
00101 *>           part of the array  A  must contain the matrix  A,  otherwise
00102 *>           the leading  k by n  part of the array  A  must contain  the
00103 *>           matrix A.
00104 *> \endverbatim
00105 *>
00106 *> \param[in] LDA
00107 *> \verbatim
00108 *>          LDA is INTEGER
00109 *>           On entry, LDA specifies the first dimension of A as declared
00110 *>           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
00111 *>           then  LDA must be at least  max( 1, n ), otherwise  LDA must
00112 *>           be at least  max( 1, k ).
00113 *> \endverbatim
00114 *>
00115 *> \param[in] B
00116 *> \verbatim
00117 *>          B is COMPLEX array of DIMENSION ( LDB, kb ), where kb is
00118 *>           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
00119 *>           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
00120 *>           part of the array  B  must contain the matrix  B,  otherwise
00121 *>           the leading  k by n  part of the array  B  must contain  the
00122 *>           matrix B.
00123 *> \endverbatim
00124 *>
00125 *> \param[in] LDB
00126 *> \verbatim
00127 *>          LDB is INTEGER
00128 *>           On entry, LDB specifies the first dimension of B as declared
00129 *>           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
00130 *>           then  LDB must be at least  max( 1, n ), otherwise  LDB must
00131 *>           be at least  max( 1, k ).
00132 *> \endverbatim
00133 *>
00134 *> \param[in] BETA
00135 *> \verbatim
00136 *>          BETA is REAL
00137 *>           On entry, BETA specifies the scalar beta.
00138 *> \endverbatim
00139 *>
00140 *> \param[in,out] C
00141 *> \verbatim
00142 *>          C is COMPLEX array of DIMENSION ( LDC, n ).
00143 *>           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
00144 *>           upper triangular part of the array C must contain the upper
00145 *>           triangular part  of the  hermitian matrix  and the strictly
00146 *>           lower triangular part of C is not referenced.  On exit, the
00147 *>           upper triangular part of the array  C is overwritten by the
00148 *>           upper triangular part of the updated matrix.
00149 *>           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
00150 *>           lower triangular part of the array C must contain the lower
00151 *>           triangular part  of the  hermitian matrix  and the strictly
00152 *>           upper triangular part of C is not referenced.  On exit, the
00153 *>           lower triangular part of the array  C is overwritten by the
00154 *>           lower triangular part of the updated matrix.
00155 *>           Note that the imaginary parts of the diagonal elements need
00156 *>           not be set,  they are assumed to be zero,  and on exit they
00157 *>           are set to zero.
00158 *> \endverbatim
00159 *>
00160 *> \param[in] LDC
00161 *> \verbatim
00162 *>          LDC is INTEGER
00163 *>           On entry, LDC specifies the first dimension of C as declared
00164 *>           in  the  calling  (sub)  program.   LDC  must  be  at  least
00165 *>           max( 1, n ).
00166 *> \endverbatim
00167 *
00168 *  Authors:
00169 *  ========
00170 *
00171 *> \author Univ. of Tennessee 
00172 *> \author Univ. of California Berkeley 
00173 *> \author Univ. of Colorado Denver 
00174 *> \author NAG Ltd. 
00175 *
00176 *> \date November 2011
00177 *
00178 *> \ingroup complex_blas_level3
00179 *
00180 *> \par Further Details:
00181 *  =====================
00182 *>
00183 *> \verbatim
00184 *>
00185 *>  Level 3 Blas routine.
00186 *>
00187 *>  -- Written on 8-February-1989.
00188 *>     Jack Dongarra, Argonne National Laboratory.
00189 *>     Iain Duff, AERE Harwell.
00190 *>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
00191 *>     Sven Hammarling, Numerical Algorithms Group Ltd.
00192 *>
00193 *>  -- Modified 8-Nov-93 to set C(J,J) to REAL( C(J,J) ) when BETA = 1.
00194 *>     Ed Anderson, Cray Research Inc.
00195 *> \endverbatim
00196 *>
00197 *  =====================================================================
00198       SUBROUTINE CHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
00199 *
00200 *  -- Reference BLAS level3 routine (version 3.4.0) --
00201 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
00202 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00203 *     November 2011
00204 *
00205 *     .. Scalar Arguments ..
00206       COMPLEX ALPHA
00207       REAL BETA
00208       INTEGER K,LDA,LDB,LDC,N
00209       CHARACTER TRANS,UPLO
00210 *     ..
00211 *     .. Array Arguments ..
00212       COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
00213 *     ..
00214 *
00215 *  =====================================================================
00216 *
00217 *     .. External Functions ..
00218       LOGICAL LSAME
00219       EXTERNAL LSAME
00220 *     ..
00221 *     .. External Subroutines ..
00222       EXTERNAL XERBLA
00223 *     ..
00224 *     .. Intrinsic Functions ..
00225       INTRINSIC CONJG,MAX,REAL
00226 *     ..
00227 *     .. Local Scalars ..
00228       COMPLEX TEMP1,TEMP2
00229       INTEGER I,INFO,J,L,NROWA
00230       LOGICAL UPPER
00231 *     ..
00232 *     .. Parameters ..
00233       REAL ONE
00234       PARAMETER (ONE=1.0E+0)
00235       COMPLEX ZERO
00236       PARAMETER (ZERO= (0.0E+0,0.0E+0))
00237 *     ..
00238 *
00239 *     Test the input parameters.
00240 *
00241       IF (LSAME(TRANS,'N')) THEN
00242           NROWA = N
00243       ELSE
00244           NROWA = K
00245       END IF
00246       UPPER = LSAME(UPLO,'U')
00247 *
00248       INFO = 0
00249       IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
00250           INFO = 1
00251       ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
00252      +         (.NOT.LSAME(TRANS,'C'))) THEN
00253           INFO = 2
00254       ELSE IF (N.LT.0) THEN
00255           INFO = 3
00256       ELSE IF (K.LT.0) THEN
00257           INFO = 4
00258       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
00259           INFO = 7
00260       ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
00261           INFO = 9
00262       ELSE IF (LDC.LT.MAX(1,N)) THEN
00263           INFO = 12
00264       END IF
00265       IF (INFO.NE.0) THEN
00266           CALL XERBLA('CHER2K',INFO)
00267           RETURN
00268       END IF
00269 *
00270 *     Quick return if possible.
00271 *
00272       IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
00273      +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
00274 *
00275 *     And when  alpha.eq.zero.
00276 *
00277       IF (ALPHA.EQ.ZERO) THEN
00278           IF (UPPER) THEN
00279               IF (BETA.EQ.REAL(ZERO)) THEN
00280                   DO 20 J = 1,N
00281                       DO 10 I = 1,J
00282                           C(I,J) = ZERO
00283    10                 CONTINUE
00284    20             CONTINUE
00285               ELSE
00286                   DO 40 J = 1,N
00287                       DO 30 I = 1,J - 1
00288                           C(I,J) = BETA*C(I,J)
00289    30                 CONTINUE
00290                       C(J,J) = BETA*REAL(C(J,J))
00291    40             CONTINUE
00292               END IF
00293           ELSE
00294               IF (BETA.EQ.REAL(ZERO)) THEN
00295                   DO 60 J = 1,N
00296                       DO 50 I = J,N
00297                           C(I,J) = ZERO
00298    50                 CONTINUE
00299    60             CONTINUE
00300               ELSE
00301                   DO 80 J = 1,N
00302                       C(J,J) = BETA*REAL(C(J,J))
00303                       DO 70 I = J + 1,N
00304                           C(I,J) = BETA*C(I,J)
00305    70                 CONTINUE
00306    80             CONTINUE
00307               END IF
00308           END IF
00309           RETURN
00310       END IF
00311 *
00312 *     Start the operations.
00313 *
00314       IF (LSAME(TRANS,'N')) THEN
00315 *
00316 *        Form  C := alpha*A*B**H + conjg( alpha )*B*A**H +
00317 *                   C.
00318 *
00319           IF (UPPER) THEN
00320               DO 130 J = 1,N
00321                   IF (BETA.EQ.REAL(ZERO)) THEN
00322                       DO 90 I = 1,J
00323                           C(I,J) = ZERO
00324    90                 CONTINUE
00325                   ELSE IF (BETA.NE.ONE) THEN
00326                       DO 100 I = 1,J - 1
00327                           C(I,J) = BETA*C(I,J)
00328   100                 CONTINUE
00329                       C(J,J) = BETA*REAL(C(J,J))
00330                   ELSE
00331                       C(J,J) = REAL(C(J,J))
00332                   END IF
00333                   DO 120 L = 1,K
00334                       IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
00335                           TEMP1 = ALPHA*CONJG(B(J,L))
00336                           TEMP2 = CONJG(ALPHA*A(J,L))
00337                           DO 110 I = 1,J - 1
00338                               C(I,J) = C(I,J) + A(I,L)*TEMP1 +
00339      +                                 B(I,L)*TEMP2
00340   110                     CONTINUE
00341                           C(J,J) = REAL(C(J,J)) +
00342      +                             REAL(A(J,L)*TEMP1+B(J,L)*TEMP2)
00343                       END IF
00344   120             CONTINUE
00345   130         CONTINUE
00346           ELSE
00347               DO 180 J = 1,N
00348                   IF (BETA.EQ.REAL(ZERO)) THEN
00349                       DO 140 I = J,N
00350                           C(I,J) = ZERO
00351   140                 CONTINUE
00352                   ELSE IF (BETA.NE.ONE) THEN
00353                       DO 150 I = J + 1,N
00354                           C(I,J) = BETA*C(I,J)
00355   150                 CONTINUE
00356                       C(J,J) = BETA*REAL(C(J,J))
00357                   ELSE
00358                       C(J,J) = REAL(C(J,J))
00359                   END IF
00360                   DO 170 L = 1,K
00361                       IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
00362                           TEMP1 = ALPHA*CONJG(B(J,L))
00363                           TEMP2 = CONJG(ALPHA*A(J,L))
00364                           DO 160 I = J + 1,N
00365                               C(I,J) = C(I,J) + A(I,L)*TEMP1 +
00366      +                                 B(I,L)*TEMP2
00367   160                     CONTINUE
00368                           C(J,J) = REAL(C(J,J)) +
00369      +                             REAL(A(J,L)*TEMP1+B(J,L)*TEMP2)
00370                       END IF
00371   170             CONTINUE
00372   180         CONTINUE
00373           END IF
00374       ELSE
00375 *
00376 *        Form  C := alpha*A**H*B + conjg( alpha )*B**H*A +
00377 *                   C.
00378 *
00379           IF (UPPER) THEN
00380               DO 210 J = 1,N
00381                   DO 200 I = 1,J
00382                       TEMP1 = ZERO
00383                       TEMP2 = ZERO
00384                       DO 190 L = 1,K
00385                           TEMP1 = TEMP1 + CONJG(A(L,I))*B(L,J)
00386                           TEMP2 = TEMP2 + CONJG(B(L,I))*A(L,J)
00387   190                 CONTINUE
00388                       IF (I.EQ.J) THEN
00389                           IF (BETA.EQ.REAL(ZERO)) THEN
00390                               C(J,J) = REAL(ALPHA*TEMP1+
00391      +                                 CONJG(ALPHA)*TEMP2)
00392                           ELSE
00393                               C(J,J) = BETA*REAL(C(J,J)) +
00394      +                                 REAL(ALPHA*TEMP1+
00395      +                                 CONJG(ALPHA)*TEMP2)
00396                           END IF
00397                       ELSE
00398                           IF (BETA.EQ.REAL(ZERO)) THEN
00399                               C(I,J) = ALPHA*TEMP1 + CONJG(ALPHA)*TEMP2
00400                           ELSE
00401                               C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
00402      +                                 CONJG(ALPHA)*TEMP2
00403                           END IF
00404                       END IF
00405   200             CONTINUE
00406   210         CONTINUE
00407           ELSE
00408               DO 240 J = 1,N
00409                   DO 230 I = J,N
00410                       TEMP1 = ZERO
00411                       TEMP2 = ZERO
00412                       DO 220 L = 1,K
00413                           TEMP1 = TEMP1 + CONJG(A(L,I))*B(L,J)
00414                           TEMP2 = TEMP2 + CONJG(B(L,I))*A(L,J)
00415   220                 CONTINUE
00416                       IF (I.EQ.J) THEN
00417                           IF (BETA.EQ.REAL(ZERO)) THEN
00418                               C(J,J) = REAL(ALPHA*TEMP1+
00419      +                                 CONJG(ALPHA)*TEMP2)
00420                           ELSE
00421                               C(J,J) = BETA*REAL(C(J,J)) +
00422      +                                 REAL(ALPHA*TEMP1+
00423      +                                 CONJG(ALPHA)*TEMP2)
00424                           END IF
00425                       ELSE
00426                           IF (BETA.EQ.REAL(ZERO)) THEN
00427                               C(I,J) = ALPHA*TEMP1 + CONJG(ALPHA)*TEMP2
00428                           ELSE
00429                               C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
00430      +                                 CONJG(ALPHA)*TEMP2
00431                           END IF
00432                       END IF
00433   230             CONTINUE
00434   240         CONTINUE
00435           END IF
00436       END IF
00437 *
00438       RETURN
00439 *
00440 *     End of CHER2K.
00441 *
00442       END
 All Files Functions