LAPACK: BLAS/SRC/dsyrk.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSYRK(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 PRECISION A(LDA,*),C(LDC,*)
00009 *     ..
00010 *
00011 *  Purpose
00012 *  =======
00013 *
00014 *  DSYRK  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 *              TRANS = 'C' or 'c'   C := alpha*A'*A + 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   matrix   A,   and  on   entry   with
00062 *           TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
00063 *           of rows of the matrix  A.  K must be at least zero.
00064 *           Unchanged on exit.
00065 *
00066 *  ALPHA  - DOUBLE PRECISION.
00067 *           On entry, ALPHA specifies the scalar alpha.
00068 *           Unchanged on exit.
00069 *
00070 *  A      - DOUBLE PRECISION 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 *  BETA   - DOUBLE PRECISION.
00086 *           On entry, BETA specifies the scalar beta.
00087 *           Unchanged on exit.
00088 *
00089 *  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
00090 *           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
00091 *           upper triangular part of the array C must contain the upper
00092 *           triangular part  of the  symmetric matrix  and the strictly
00093 *           lower triangular part of C is not referenced.  On exit, the
00094 *           upper triangular part of the array  C is overwritten by the
00095 *           upper triangular part of the updated matrix.
00096 *           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
00097 *           lower triangular part of the array C must contain the lower
00098 *           triangular part  of the  symmetric matrix  and the strictly
00099 *           upper triangular part of C is not referenced.  On exit, the
00100 *           lower triangular part of the array  C is overwritten by the
00101 *           lower triangular part of the updated matrix.
00102 *
00103 *  LDC    - INTEGER.
00104 *           On entry, LDC specifies the first dimension of C as declared
00105 *           in  the  calling  (sub)  program.   LDC  must  be  at  least
00106 *           max( 1, n ).
00107 *           Unchanged on exit.
00108 *
00109 *  Further Details
00110 *  ===============
00111 *
00112 *  Level 3 Blas routine.
00113 *
00114 *  -- Written on 8-February-1989.
00115 *     Jack Dongarra, Argonne National Laboratory.
00116 *     Iain Duff, AERE Harwell.
00117 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
00118 *     Sven Hammarling, Numerical Algorithms Group Ltd.
00119 *
00120 *  =====================================================================
00121 *
00122 *     .. External Functions ..
00123       LOGICAL LSAME
00124       EXTERNAL LSAME
00125 *     ..
00126 *     .. External Subroutines ..
00127       EXTERNAL XERBLA
00128 *     ..
00129 *     .. Intrinsic Functions ..
00130       INTRINSIC MAX
00131 *     ..
00132 *     .. Local Scalars ..
00133       DOUBLE PRECISION TEMP
00134       INTEGER I,INFO,J,L,NROWA
00135       LOGICAL UPPER
00136 *     ..
00137 *     .. Parameters ..
00138       DOUBLE PRECISION ONE,ZERO
00139       PARAMETER (ONE=1.0D+0,ZERO=0.0D+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')) .AND.
00156      +         (.NOT.LSAME(TRANS,'C'))) THEN
00157           INFO = 2
00158       ELSE IF (N.LT.0) THEN
00159           INFO = 3
00160       ELSE IF (K.LT.0) THEN
00161           INFO = 4
00162       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
00163           INFO = 7
00164       ELSE IF (LDC.LT.MAX(1,N)) THEN
00165           INFO = 10
00166       END IF
00167       IF (INFO.NE.0) THEN
00168           CALL XERBLA('DSYRK ',INFO)
00169           RETURN
00170       END IF
00171 *
00172 *     Quick return if possible.
00173 *
00174       IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
00175      +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
00176 *
00177 *     And when  alpha.eq.zero.
00178 *
00179       IF (ALPHA.EQ.ZERO) THEN
00180           IF (UPPER) THEN
00181               IF (BETA.EQ.ZERO) THEN
00182                   DO 20 J = 1,N
00183                       DO 10 I = 1,J
00184                           C(I,J) = ZERO
00185    10                 CONTINUE
00186    20             CONTINUE
00187               ELSE
00188                   DO 40 J = 1,N
00189                       DO 30 I = 1,J
00190                           C(I,J) = BETA*C(I,J)
00191    30                 CONTINUE
00192    40             CONTINUE
00193               END IF
00194           ELSE
00195               IF (BETA.EQ.ZERO) THEN
00196                   DO 60 J = 1,N
00197                       DO 50 I = J,N
00198                           C(I,J) = ZERO
00199    50                 CONTINUE
00200    60             CONTINUE
00201               ELSE
00202                   DO 80 J = 1,N
00203                       DO 70 I = J,N
00204                           C(I,J) = BETA*C(I,J)
00205    70                 CONTINUE
00206    80             CONTINUE
00207               END IF
00208           END IF
00209           RETURN
00210       END IF
00211 *
00212 *     Start the operations.
00213 *
00214       IF (LSAME(TRANS,'N')) THEN
00215 *
00216 *        Form  C := alpha*A*A' + beta*C.
00217 *
00218           IF (UPPER) THEN
00219               DO 130 J = 1,N
00220                   IF (BETA.EQ.ZERO) THEN
00221                       DO 90 I = 1,J
00222                           C(I,J) = ZERO
00223    90                 CONTINUE
00224                   ELSE IF (BETA.NE.ONE) THEN
00225                       DO 100 I = 1,J
00226                           C(I,J) = BETA*C(I,J)
00227   100                 CONTINUE
00228                   END IF
00229                   DO 120 L = 1,K
00230                       IF (A(J,L).NE.ZERO) THEN
00231                           TEMP = ALPHA*A(J,L)
00232                           DO 110 I = 1,J
00233                               C(I,J) = C(I,J) + TEMP*A(I,L)
00234   110                     CONTINUE
00235                       END IF
00236   120             CONTINUE
00237   130         CONTINUE
00238           ELSE
00239               DO 180 J = 1,N
00240                   IF (BETA.EQ.ZERO) THEN
00241                       DO 140 I = J,N
00242                           C(I,J) = ZERO
00243   140                 CONTINUE
00244                   ELSE IF (BETA.NE.ONE) THEN
00245                       DO 150 I = J,N
00246                           C(I,J) = BETA*C(I,J)
00247   150                 CONTINUE
00248                   END IF
00249                   DO 170 L = 1,K
00250                       IF (A(J,L).NE.ZERO) THEN
00251                           TEMP = ALPHA*A(J,L)
00252                           DO 160 I = J,N
00253                               C(I,J) = C(I,J) + TEMP*A(I,L)
00254   160                     CONTINUE
00255                       END IF
00256   170             CONTINUE
00257   180         CONTINUE
00258           END IF
00259       ELSE
00260 *
00261 *        Form  C := alpha*A'*A + beta*C.
00262 *
00263           IF (UPPER) THEN
00264               DO 210 J = 1,N
00265                   DO 200 I = 1,J
00266                       TEMP = ZERO
00267                       DO 190 L = 1,K
00268                           TEMP = TEMP + A(L,I)*A(L,J)
00269   190                 CONTINUE
00270                       IF (BETA.EQ.ZERO) THEN
00271                           C(I,J) = ALPHA*TEMP
00272                       ELSE
00273                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
00274                       END IF
00275   200             CONTINUE
00276   210         CONTINUE
00277           ELSE
00278               DO 240 J = 1,N
00279                   DO 230 I = J,N
00280                       TEMP = ZERO
00281                       DO 220 L = 1,K
00282                           TEMP = TEMP + A(L,I)*A(L,J)
00283   220                 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   230             CONTINUE
00290   240         CONTINUE
00291           END IF
00292       END IF
00293 *
00294       RETURN
00295 *
00296 *     End of DSYRK .
00297 *
00298       END