LAPACK: BLAS/SRC/zsyr2k.f Source File

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