LAPACK: BLAS/SRC/csyrk.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
00002 *     .. Scalar Arguments ..
00003       COMPLEX ALPHA,BETA
00004       INTEGER K,LDA,LDC,N
00005       CHARACTER TRANS,UPLO
00006 *     ..
00007 *     .. Array Arguments ..
00008       COMPLEX A(LDA,*),C(LDC,*)
00009 *     ..
00010 *
00011 *  Purpose
00012 *  =======
00013 *
00014 *  CSYRK  performs one of the symmetric rank k operations
00015 *
00016 *     C := alpha*A*A' + beta*C,
00017 *
00018 *  or
00019 *
00020 *     C := alpha*A'*A + beta*C,
00021 *
00022 *  where  alpha and beta  are scalars,  C is an  n by n symmetric matrix
00023 *  and  A  is an  n by k  matrix in the first case and a  k by n  matrix
00024 *  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*A' + beta*C.
00047 *
00048 *              TRANS = 'T' or 't'   C := alpha*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 = 'T' or 't',  K  specifies  the number of rows of the
00061 *           matrix A.  K must be at least zero.
00062 *           Unchanged on exit.
00063 *
00064 *  ALPHA  - COMPLEX         .
00065 *           On entry, ALPHA specifies the scalar alpha.
00066 *           Unchanged on exit.
00067 *
00068 *  A      - COMPLEX          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   - COMPLEX         .
00084 *           On entry, BETA specifies the scalar beta.
00085 *           Unchanged on exit.
00086 *
00087 *  C      - COMPLEX          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  symmetric 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  symmetric 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 *
00101 *  LDC    - INTEGER.
00102 *           On entry, LDC specifies the first dimension of C as declared
00103 *           in  the  calling  (sub)  program.   LDC  must  be  at  least
00104 *           max( 1, n ).
00105 *           Unchanged on exit.
00106 *
00107 *  Further Details
00108 *  ===============
00109 *
00110 *  Level 3 Blas routine.
00111 *
00112 *  -- Written on 8-February-1989.
00113 *     Jack Dongarra, Argonne National Laboratory.
00114 *     Iain Duff, AERE Harwell.
00115 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
00116 *     Sven Hammarling, Numerical Algorithms Group Ltd.
00117 *
00118 *  =====================================================================
00119 *
00120 *     .. External Functions ..
00121       LOGICAL LSAME
00122       EXTERNAL LSAME
00123 *     ..
00124 *     .. External Subroutines ..
00125       EXTERNAL XERBLA
00126 *     ..
00127 *     .. Intrinsic Functions ..
00128       INTRINSIC MAX
00129 *     ..
00130 *     .. Local Scalars ..
00131       COMPLEX TEMP
00132       INTEGER I,INFO,J,L,NROWA
00133       LOGICAL UPPER
00134 *     ..
00135 *     .. Parameters ..
00136       COMPLEX ONE
00137       PARAMETER (ONE= (1.0E+0,0.0E+0))
00138       COMPLEX ZERO
00139       PARAMETER (ZERO= (0.0E+0,0.0E+0))
00140 *     ..
00141 *
00142 *     Test the input parameters.
00143 *
00144       IF (LSAME(TRANS,'N')) THEN
00145           NROWA = N
00146       ELSE
00147           NROWA = K
00148       END IF
00149       UPPER = LSAME(UPLO,'U')
00150 *
00151       INFO = 0
00152       IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
00153           INFO = 1
00154       ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
00155      +         (.NOT.LSAME(TRANS,'T'))) THEN
00156           INFO = 2
00157       ELSE IF (N.LT.0) THEN
00158           INFO = 3
00159       ELSE IF (K.LT.0) THEN
00160           INFO = 4
00161       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
00162           INFO = 7
00163       ELSE IF (LDC.LT.MAX(1,N)) THEN
00164           INFO = 10
00165       END IF
00166       IF (INFO.NE.0) THEN
00167           CALL XERBLA('CSYRK ',INFO)
00168           RETURN
00169       END IF
00170 *
00171 *     Quick return if possible.
00172 *
00173       IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
00174      +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
00175 *
00176 *     And when  alpha.eq.zero.
00177 *
00178       IF (ALPHA.EQ.ZERO) THEN
00179           IF (UPPER) THEN
00180               IF (BETA.EQ.ZERO) THEN
00181                   DO 20 J = 1,N
00182                       DO 10 I = 1,J
00183                           C(I,J) = ZERO
00184    10                 CONTINUE
00185    20             CONTINUE
00186               ELSE
00187                   DO 40 J = 1,N
00188                       DO 30 I = 1,J
00189                           C(I,J) = BETA*C(I,J)
00190    30                 CONTINUE
00191    40             CONTINUE
00192               END IF
00193           ELSE
00194               IF (BETA.EQ.ZERO) THEN
00195                   DO 60 J = 1,N
00196                       DO 50 I = J,N
00197                           C(I,J) = ZERO
00198    50                 CONTINUE
00199    60             CONTINUE
00200               ELSE
00201                   DO 80 J = 1,N
00202                       DO 70 I = J,N
00203                           C(I,J) = BETA*C(I,J)
00204    70                 CONTINUE
00205    80             CONTINUE
00206               END IF
00207           END IF
00208           RETURN
00209       END IF
00210 *
00211 *     Start the operations.
00212 *
00213       IF (LSAME(TRANS,'N')) THEN
00214 *
00215 *        Form  C := alpha*A*A' + beta*C.
00216 *
00217           IF (UPPER) THEN
00218               DO 130 J = 1,N
00219                   IF (BETA.EQ.ZERO) THEN
00220                       DO 90 I = 1,J
00221                           C(I,J) = ZERO
00222    90                 CONTINUE
00223                   ELSE IF (BETA.NE.ONE) THEN
00224                       DO 100 I = 1,J
00225                           C(I,J) = BETA*C(I,J)
00226   100                 CONTINUE
00227                   END IF
00228                   DO 120 L = 1,K
00229                       IF (A(J,L).NE.ZERO) THEN
00230                           TEMP = ALPHA*A(J,L)
00231                           DO 110 I = 1,J
00232                               C(I,J) = C(I,J) + TEMP*A(I,L)
00233   110                     CONTINUE
00234                       END IF
00235   120             CONTINUE
00236   130         CONTINUE
00237           ELSE
00238               DO 180 J = 1,N
00239                   IF (BETA.EQ.ZERO) THEN
00240                       DO 140 I = J,N
00241                           C(I,J) = ZERO
00242   140                 CONTINUE
00243                   ELSE IF (BETA.NE.ONE) THEN
00244                       DO 150 I = J,N
00245                           C(I,J) = BETA*C(I,J)
00246   150                 CONTINUE
00247                   END IF
00248                   DO 170 L = 1,K
00249                       IF (A(J,L).NE.ZERO) THEN
00250                           TEMP = ALPHA*A(J,L)
00251                           DO 160 I = J,N
00252                               C(I,J) = C(I,J) + TEMP*A(I,L)
00253   160                     CONTINUE
00254                       END IF
00255   170             CONTINUE
00256   180         CONTINUE
00257           END IF
00258       ELSE
00259 *
00260 *        Form  C := alpha*A'*A + beta*C.
00261 *
00262           IF (UPPER) THEN
00263               DO 210 J = 1,N
00264                   DO 200 I = 1,J
00265                       TEMP = ZERO
00266                       DO 190 L = 1,K
00267                           TEMP = TEMP + A(L,I)*A(L,J)
00268   190                 CONTINUE
00269                       IF (BETA.EQ.ZERO) THEN
00270                           C(I,J) = ALPHA*TEMP
00271                       ELSE
00272                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
00273                       END IF
00274   200             CONTINUE
00275   210         CONTINUE
00276           ELSE
00277               DO 240 J = 1,N
00278                   DO 230 I = J,N
00279                       TEMP = ZERO
00280                       DO 220 L = 1,K
00281                           TEMP = TEMP + A(L,I)*A(L,J)
00282   220                 CONTINUE
00283                       IF (BETA.EQ.ZERO) THEN
00284                           C(I,J) = ALPHA*TEMP
00285                       ELSE
00286                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
00287                       END IF
00288   230             CONTINUE
00289   240         CONTINUE
00290           END IF
00291       END IF
00292 *
00293       RETURN
00294 *
00295 *     End of CSYRK .
00296 *
00297       END