LAPACK: BLAS/SRC/dsyr2k.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
00002 *     .. Scalar Arguments ..
00003       DOUBLE PRECISION ALPHA,BETA
00004       INTEGER K,LDA,LDB,LDC,N
00005       CHARACTER TRANS,UPLO
00006 *     ..
00007 *     .. Array Arguments ..
00008       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
00009 *     ..
00010 *
00011 *  Purpose
00012 *  =======
00013 *
00014 *  DSYR2K  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 *              TRANS = 'C' or 'c'   C := alpha*A'*B + alpha*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 = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
00066 *           of rows of the matrices  A and B.  K must be at least  zero.
00067 *           Unchanged on exit.
00068 *
00069 *  ALPHA  - DOUBLE PRECISION.
00070 *           On entry, ALPHA specifies the scalar alpha.
00071 *           Unchanged on exit.
00072 *
00073 *  A      - DOUBLE PRECISION 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      - DOUBLE PRECISION 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   - DOUBLE PRECISION.
00104 *           On entry, BETA specifies the scalar beta.
00105 *           Unchanged on exit.
00106 *
00107 *  C      - DOUBLE PRECISION 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  symmetric 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  symmetric 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 *
00121 *  LDC    - INTEGER.
00122 *           On entry, LDC specifies the first dimension of C as declared
00123 *           in  the  calling  (sub)  program.   LDC  must  be  at  least
00124 *           max( 1, n ).
00125 *           Unchanged on exit.
00126 *
00127 *  Further Details
00128 *  ===============
00129 *
00130 *  Level 3 Blas routine.
00131 *
00132 *
00133 *  -- Written on 8-February-1989.
00134 *     Jack Dongarra, Argonne National Laboratory.
00135 *     Iain Duff, AERE Harwell.
00136 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
00137 *     Sven Hammarling, Numerical Algorithms Group Ltd.
00138 *
00139 *  =====================================================================
00140 *
00141 *     .. External Functions ..
00142       LOGICAL LSAME
00143       EXTERNAL LSAME
00144 *     ..
00145 *     .. External Subroutines ..
00146       EXTERNAL XERBLA
00147 *     ..
00148 *     .. Intrinsic Functions ..
00149       INTRINSIC MAX
00150 *     ..
00151 *     .. Local Scalars ..
00152       DOUBLE PRECISION TEMP1,TEMP2
00153       INTEGER I,INFO,J,L,NROWA
00154       LOGICAL UPPER
00155 *     ..
00156 *     .. Parameters ..
00157       DOUBLE PRECISION ONE,ZERO
00158       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
00159 *     ..
00160 *
00161 *     Test the input parameters.
00162 *
00163       IF (LSAME(TRANS,'N')) THEN
00164           NROWA = N
00165       ELSE
00166           NROWA = K
00167       END IF
00168       UPPER = LSAME(UPLO,'U')
00169 *
00170       INFO = 0
00171       IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
00172           INFO = 1
00173       ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
00174      +         (.NOT.LSAME(TRANS,'T')) .AND.
00175      +         (.NOT.LSAME(TRANS,'C'))) THEN
00176           INFO = 2
00177       ELSE IF (N.LT.0) THEN
00178           INFO = 3
00179       ELSE IF (K.LT.0) THEN
00180           INFO = 4
00181       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
00182           INFO = 7
00183       ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
00184           INFO = 9
00185       ELSE IF (LDC.LT.MAX(1,N)) THEN
00186           INFO = 12
00187       END IF
00188       IF (INFO.NE.0) THEN
00189           CALL XERBLA('DSYR2K',INFO)
00190           RETURN
00191       END IF
00192 *
00193 *     Quick return if possible.
00194 *
00195       IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
00196      +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
00197 *
00198 *     And when  alpha.eq.zero.
00199 *
00200       IF (ALPHA.EQ.ZERO) THEN
00201           IF (UPPER) THEN
00202               IF (BETA.EQ.ZERO) THEN
00203                   DO 20 J = 1,N
00204                       DO 10 I = 1,J
00205                           C(I,J) = ZERO
00206    10                 CONTINUE
00207    20             CONTINUE
00208               ELSE
00209                   DO 40 J = 1,N
00210                       DO 30 I = 1,J
00211                           C(I,J) = BETA*C(I,J)
00212    30                 CONTINUE
00213    40             CONTINUE
00214               END IF
00215           ELSE
00216               IF (BETA.EQ.ZERO) THEN
00217                   DO 60 J = 1,N
00218                       DO 50 I = J,N
00219                           C(I,J) = ZERO
00220    50                 CONTINUE
00221    60             CONTINUE
00222               ELSE
00223                   DO 80 J = 1,N
00224                       DO 70 I = J,N
00225                           C(I,J) = BETA*C(I,J)
00226    70                 CONTINUE
00227    80             CONTINUE
00228               END IF
00229           END IF
00230           RETURN
00231       END IF
00232 *
00233 *     Start the operations.
00234 *
00235       IF (LSAME(TRANS,'N')) THEN
00236 *
00237 *        Form  C := alpha*A*B' + alpha*B*A' + C.
00238 *
00239           IF (UPPER) THEN
00240               DO 130 J = 1,N
00241                   IF (BETA.EQ.ZERO) THEN
00242                       DO 90 I = 1,J
00243                           C(I,J) = ZERO
00244    90                 CONTINUE
00245                   ELSE IF (BETA.NE.ONE) THEN
00246                       DO 100 I = 1,J
00247                           C(I,J) = BETA*C(I,J)
00248   100                 CONTINUE
00249                   END IF
00250                   DO 120 L = 1,K
00251                       IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
00252                           TEMP1 = ALPHA*B(J,L)
00253                           TEMP2 = ALPHA*A(J,L)
00254                           DO 110 I = 1,J
00255                               C(I,J) = C(I,J) + A(I,L)*TEMP1 +
00256      +                                 B(I,L)*TEMP2
00257   110                     CONTINUE
00258                       END IF
00259   120             CONTINUE
00260   130         CONTINUE
00261           ELSE
00262               DO 180 J = 1,N
00263                   IF (BETA.EQ.ZERO) THEN
00264                       DO 140 I = J,N
00265                           C(I,J) = ZERO
00266   140                 CONTINUE
00267                   ELSE IF (BETA.NE.ONE) THEN
00268                       DO 150 I = J,N
00269                           C(I,J) = BETA*C(I,J)
00270   150                 CONTINUE
00271                   END IF
00272                   DO 170 L = 1,K
00273                       IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
00274                           TEMP1 = ALPHA*B(J,L)
00275                           TEMP2 = ALPHA*A(J,L)
00276                           DO 160 I = J,N
00277                               C(I,J) = C(I,J) + A(I,L)*TEMP1 +
00278      +                                 B(I,L)*TEMP2
00279   160                     CONTINUE
00280                       END IF
00281   170             CONTINUE
00282   180         CONTINUE
00283           END IF
00284       ELSE
00285 *
00286 *        Form  C := alpha*A'*B + alpha*B'*A + C.
00287 *
00288           IF (UPPER) THEN
00289               DO 210 J = 1,N
00290                   DO 200 I = 1,J
00291                       TEMP1 = ZERO
00292                       TEMP2 = ZERO
00293                       DO 190 L = 1,K
00294                           TEMP1 = TEMP1 + A(L,I)*B(L,J)
00295                           TEMP2 = TEMP2 + B(L,I)*A(L,J)
00296   190                 CONTINUE
00297                       IF (BETA.EQ.ZERO) THEN
00298                           C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
00299                       ELSE
00300                           C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
00301      +                             ALPHA*TEMP2
00302                       END IF
00303   200             CONTINUE
00304   210         CONTINUE
00305           ELSE
00306               DO 240 J = 1,N
00307                   DO 230 I = J,N
00308                       TEMP1 = ZERO
00309                       TEMP2 = ZERO
00310                       DO 220 L = 1,K
00311                           TEMP1 = TEMP1 + A(L,I)*B(L,J)
00312                           TEMP2 = TEMP2 + B(L,I)*A(L,J)
00313   220                 CONTINUE
00314                       IF (BETA.EQ.ZERO) THEN
00315                           C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
00316                       ELSE
00317                           C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
00318      +                             ALPHA*TEMP2
00319                       END IF
00320   230             CONTINUE
00321   240         CONTINUE
00322           END IF
00323       END IF
00324 *
00325       RETURN
00326 *
00327 *     End of DSYR2K.
00328 *
00329       END