LAPACK: BLAS/SRC/cher2k.f Source File

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