BLAS/SRC/ssyrk.f

   1 *> \brief \b SSYRK
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE SSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
  12 *
  13 *       .. Scalar Arguments ..
  14 *       REAL ALPHA,BETA
  15 *       INTEGER K,LDA,LDC,N
  16 *       CHARACTER TRANS,UPLO
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       REAL A(LDA,*),C(LDC,*)
  20 *       ..
  21 *
  22 *
  23 *> \par Purpose:
  24 *  =============
  25 *>
  26 *> \verbatim
  27 *>
  28 *> SSYRK  performs one of the symmetric rank k operations
  29 *>
  30 *>    C := alpha*A*A**T + beta*C,
  31 *>
  32 *> or
  33 *>
  34 *>    C := alpha*A**T*A + beta*C,
  35 *>
  36 *> where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
  37 *> and  A  is an  n by k  matrix in the first case and a  k by n  matrix
  38 *> in the second case.
  39 *> \endverbatim
  40 *
  41 *  Arguments:
  42 *  ==========
  43 *
  44 *> \param[in] UPLO
  45 *> \verbatim
  46 *>          UPLO is CHARACTER*1
  47 *>           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  48 *>           triangular  part  of the  array  C  is to be  referenced  as
  49 *>           follows:
  50 *>
  51 *>              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
  52 *>                                  is to be referenced.
  53 *>
  54 *>              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
  55 *>                                  is to be referenced.
  56 *> \endverbatim
  57 *>
  58 *> \param[in] TRANS
  59 *> \verbatim
  60 *>          TRANS is CHARACTER*1
  61 *>           On entry,  TRANS  specifies the operation to be performed as
  62 *>           follows:
  63 *>
  64 *>              TRANS = 'N' or 'n'   C := alpha*A*A**T + beta*C.
  65 *>
  66 *>              TRANS = 'T' or 't'   C := alpha*A**T*A + beta*C.
  67 *>
  68 *>              TRANS = 'C' or 'c'   C := alpha*A**T*A + beta*C.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] N
  72 *> \verbatim
  73 *>          N is INTEGER
  74 *>           On entry,  N specifies the order of the matrix C.  N must be
  75 *>           at least zero.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] K
  79 *> \verbatim
  80 *>          K is INTEGER
  81 *>           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
  82 *>           of  columns   of  the   matrix   A,   and  on   entry   with
  83 *>           TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
  84 *>           of rows of the matrix  A.  K must be at least zero.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] ALPHA
  88 *> \verbatim
  89 *>          ALPHA is REAL
  90 *>           On entry, ALPHA specifies the scalar alpha.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] A
  94 *> \verbatim
  95 *>          A is REAL array of DIMENSION ( LDA, ka ), where ka is
  96 *>           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  97 *>           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  98 *>           part of the array  A  must contain the matrix  A,  otherwise
  99 *>           the leading  k by n  part of the array  A  must contain  the
 100 *>           matrix A.
 101 *> \endverbatim
 102 *>
 103 *> \param[in] LDA
 104 *> \verbatim
 105 *>          LDA is INTEGER
 106 *>           On entry, LDA specifies the first dimension of A as declared
 107 *>           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
 108 *>           then  LDA must be at least  max( 1, n ), otherwise  LDA must
 109 *>           be at least  max( 1, k ).
 110 *> \endverbatim
 111 *>
 112 *> \param[in] BETA
 113 *> \verbatim
 114 *>          BETA is REAL
 115 *>           On entry, BETA specifies the scalar beta.
 116 *> \endverbatim
 117 *>
 118 *> \param[in,out] C
 119 *> \verbatim
 120 *>          C is REAL array of DIMENSION ( LDC, n ).
 121 *>           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
 122 *>           upper triangular part of the array C must contain the upper
 123 *>           triangular part  of the  symmetric matrix  and the strictly
 124 *>           lower triangular part of C is not referenced.  On exit, the
 125 *>           upper triangular part of the array  C is overwritten by the
 126 *>           upper triangular part of the updated matrix.
 127 *>           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
 128 *>           lower triangular part of the array C must contain the lower
 129 *>           triangular part  of the  symmetric matrix  and the strictly
 130 *>           upper triangular part of C is not referenced.  On exit, the
 131 *>           lower triangular part of the array  C is overwritten by the
 132 *>           lower triangular part of the updated matrix.
 133 *> \endverbatim
 134 *>
 135 *> \param[in] LDC
 136 *> \verbatim
 137 *>          LDC is INTEGER
 138 *>           On entry, LDC specifies the first dimension of C as declared
 139 *>           in  the  calling  (sub)  program.   LDC  must  be  at  least
 140 *>           max( 1, n ).
 141 *> \endverbatim
 142 *
 143 *  Authors:
 144 *  ========
 145 *
 146 *> \author Univ. of Tennessee
 147 *> \author Univ. of California Berkeley
 148 *> \author Univ. of Colorado Denver
 149 *> \author NAG Ltd.
 150 *
 151 *> \date November 2011
 152 *
 153 *> \ingroup single_blas_level3
 154 *
 155 *> \par Further Details:
 156 *  =====================
 157 *>
 158 *> \verbatim
 159 *>
 160 *>  Level 3 Blas routine.
 161 *>
 162 *>  -- Written on 8-February-1989.
 163 *>     Jack Dongarra, Argonne National Laboratory.
 164 *>     Iain Duff, AERE Harwell.
 165 *>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
 166 *>     Sven Hammarling, Numerical Algorithms Group Ltd.
 167 *> \endverbatim
 168 *>
 169 *  =====================================================================
 170       SUBROUTINE SSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
 171 *
 172 *  -- Reference BLAS level3 routine (version 3.4.0) --
 173 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 174 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 175 *     November 2011
 176 *
 177 *     .. Scalar Arguments ..
 178       REAL ALPHA,BETA
 179       INTEGER K,LDA,LDC,N
 180       CHARACTER TRANS,UPLO
 181 *     ..
 182 *     .. Array Arguments ..
 183       REAL A(LDA,*),C(LDC,*)
 184 *     ..
 185 *
 186 *  =====================================================================
 187 *
 188 *     .. External Functions ..
 189       LOGICAL LSAME
 190       EXTERNAL LSAME
 191 *     ..
 192 *     .. External Subroutines ..
 193       EXTERNAL XERBLA
 194 *     ..
 195 *     .. Intrinsic Functions ..
 196       INTRINSIC MAX
 197 *     ..
 198 *     .. Local Scalars ..
 199       REAL TEMP
 200       INTEGER I,INFO,J,L,NROWA
 201       LOGICAL UPPER
 202 *     ..
 203 *     .. Parameters ..
 204       REAL ONE,ZERO
 205       PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
 206 *     ..
 207 *
 208 *     Test the input parameters.
 209 *
 210       IF (LSAME(TRANS,'N')) THEN
 211           NROWA = N
 212       ELSE
 213           NROWA = K
 214       END IF
 215       UPPER = LSAME(UPLO,'U')
 216 *
 217       INFO = 0
 218       IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
 219           INFO = 1
 220       ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
 221      +         (.NOT.LSAME(TRANS,'T')) .AND.
 222      +         (.NOT.LSAME(TRANS,'C'))) THEN
 223           INFO = 2
 224       ELSE IF (N.LT.0) THEN
 225           INFO = 3
 226       ELSE IF (K.LT.0) THEN
 227           INFO = 4
 228       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
 229           INFO = 7
 230       ELSE IF (LDC.LT.MAX(1,N)) THEN
 231           INFO = 10
 232       END IF
 233       IF (INFO.NE.0) THEN
 234           CALL XERBLA('SSYRK ',INFO)
 235           RETURN
 236       END IF
 237 *
 238 *     Quick return if possible.
 239 *
 240       IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
 241      +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
 242 *
 243 *     And when  alpha.eq.zero.
 244 *
 245       IF (ALPHA.EQ.ZERO) THEN
 246           IF (UPPER) THEN
 247               IF (BETA.EQ.ZERO) THEN
 248                   DO 20 J = 1,N
 249                       DO 10 I = 1,J
 250                           C(I,J) = ZERO
 251    10                 CONTINUE
 252    20             CONTINUE
 253               ELSE
 254                   DO 40 J = 1,N
 255                       DO 30 I = 1,J
 256                           C(I,J) = BETA*C(I,J)
 257    30                 CONTINUE
 258    40             CONTINUE
 259               END IF
 260           ELSE
 261               IF (BETA.EQ.ZERO) THEN
 262                   DO 60 J = 1,N
 263                       DO 50 I = J,N
 264                           C(I,J) = ZERO
 265    50                 CONTINUE
 266    60             CONTINUE
 267               ELSE
 268                   DO 80 J = 1,N
 269                       DO 70 I = J,N
 270                           C(I,J) = BETA*C(I,J)
 271    70                 CONTINUE
 272    80             CONTINUE
 273               END IF
 274           END IF
 275           RETURN
 276       END IF
 277 *
 278 *     Start the operations.
 279 *
 280       IF (LSAME(TRANS,'N')) THEN
 281 *
 282 *        Form  C := alpha*A*A**T + beta*C.
 283 *
 284           IF (UPPER) THEN
 285               DO 130 J = 1,N
 286                   IF (BETA.EQ.ZERO) THEN
 287                       DO 90 I = 1,J
 288                           C(I,J) = ZERO
 289    90                 CONTINUE
 290                   ELSE IF (BETA.NE.ONE) THEN
 291                       DO 100 I = 1,J
 292                           C(I,J) = BETA*C(I,J)
 293   100                 CONTINUE
 294                   END IF
 295                   DO 120 L = 1,K
 296                       IF (A(J,L).NE.ZERO) THEN
 297                           TEMP = ALPHA*A(J,L)
 298                           DO 110 I = 1,J
 299                               C(I,J) = C(I,J) + TEMP*A(I,L)
 300   110                     CONTINUE
 301                       END IF
 302   120             CONTINUE
 303   130         CONTINUE
 304           ELSE
 305               DO 180 J = 1,N
 306                   IF (BETA.EQ.ZERO) THEN
 307                       DO 140 I = J,N
 308                           C(I,J) = ZERO
 309   140                 CONTINUE
 310                   ELSE IF (BETA.NE.ONE) THEN
 311                       DO 150 I = J,N
 312                           C(I,J) = BETA*C(I,J)
 313   150                 CONTINUE
 314                   END IF
 315                   DO 170 L = 1,K
 316                       IF (A(J,L).NE.ZERO) THEN
 317                           TEMP = ALPHA*A(J,L)
 318                           DO 160 I = J,N
 319                               C(I,J) = C(I,J) + TEMP*A(I,L)
 320   160                     CONTINUE
 321                       END IF
 322   170             CONTINUE
 323   180         CONTINUE
 324           END IF
 325       ELSE
 326 *
 327 *        Form  C := alpha*A**T*A + beta*C.
 328 *
 329           IF (UPPER) THEN
 330               DO 210 J = 1,N
 331                   DO 200 I = 1,J
 332                       TEMP = ZERO
 333                       DO 190 L = 1,K
 334                           TEMP = TEMP + A(L,I)*A(L,J)
 335   190                 CONTINUE
 336                       IF (BETA.EQ.ZERO) THEN
 337                           C(I,J) = ALPHA*TEMP
 338                       ELSE
 339                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
 340                       END IF
 341   200             CONTINUE
 342   210         CONTINUE
 343           ELSE
 344               DO 240 J = 1,N
 345                   DO 230 I = J,N
 346                       TEMP = ZERO
 347                       DO 220 L = 1,K
 348                           TEMP = TEMP + A(L,I)*A(L,J)
 349   220                 CONTINUE
 350                       IF (BETA.EQ.ZERO) THEN
 351                           C(I,J) = ALPHA*TEMP
 352                       ELSE
 353                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
 354                       END IF
 355   230             CONTINUE
 356   240         CONTINUE
 357           END IF
 358       END IF
 359 *
 360       RETURN
 361 *
 362 *     End of SSYRK .
 363 *
 364       END