LAPACK: BLAS/SRC/zherk.f Source File

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