| 1 | SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC) |
|---|
| 2 | * .. Scalar Arguments .. |
|---|
| 3 | DOUBLE PRECISION ALPHA,BETA |
|---|
| 4 | INTEGER LDA,LDB,LDC,M,N |
|---|
| 5 | CHARACTER SIDE,UPLO |
|---|
| 6 | * .. |
|---|
| 7 | * .. Array Arguments .. |
|---|
| 8 | DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) |
|---|
| 9 | * .. |
|---|
| 10 | * |
|---|
| 11 | * Purpose |
|---|
| 12 | * ======= |
|---|
| 13 | * |
|---|
| 14 | * DSYMM performs one of the matrix-matrix operations |
|---|
| 15 | * |
|---|
| 16 | * C := alpha*A*B + beta*C, |
|---|
| 17 | * |
|---|
| 18 | * or |
|---|
| 19 | * |
|---|
| 20 | * C := alpha*B*A + beta*C, |
|---|
| 21 | * |
|---|
| 22 | * where alpha and beta are scalars, A is a symmetric matrix and B and |
|---|
| 23 | * C are m by n matrices. |
|---|
| 24 | * |
|---|
| 25 | * Arguments |
|---|
| 26 | * ========== |
|---|
| 27 | * |
|---|
| 28 | * SIDE - CHARACTER*1. |
|---|
| 29 | * On entry, SIDE specifies whether the symmetric matrix A |
|---|
| 30 | * appears on the left or right in the operation as follows: |
|---|
| 31 | * |
|---|
| 32 | * SIDE = 'L' or 'l' C := alpha*A*B + beta*C, |
|---|
| 33 | * |
|---|
| 34 | * SIDE = 'R' or 'r' C := alpha*B*A + beta*C, |
|---|
| 35 | * |
|---|
| 36 | * Unchanged on exit. |
|---|
| 37 | * |
|---|
| 38 | * UPLO - CHARACTER*1. |
|---|
| 39 | * On entry, UPLO specifies whether the upper or lower |
|---|
| 40 | * triangular part of the symmetric matrix A is to be |
|---|
| 41 | * referenced as follows: |
|---|
| 42 | * |
|---|
| 43 | * UPLO = 'U' or 'u' Only the upper triangular part of the |
|---|
| 44 | * symmetric matrix is to be referenced. |
|---|
| 45 | * |
|---|
| 46 | * UPLO = 'L' or 'l' Only the lower triangular part of the |
|---|
| 47 | * symmetric matrix is to be referenced. |
|---|
| 48 | * |
|---|
| 49 | * Unchanged on exit. |
|---|
| 50 | * |
|---|
| 51 | * M - INTEGER. |
|---|
| 52 | * On entry, M specifies the number of rows of the matrix C. |
|---|
| 53 | * M must be at least zero. |
|---|
| 54 | * Unchanged on exit. |
|---|
| 55 | * |
|---|
| 56 | * N - INTEGER. |
|---|
| 57 | * On entry, N specifies the number of columns of the matrix C. |
|---|
| 58 | * N must be at least zero. |
|---|
| 59 | * Unchanged on exit. |
|---|
| 60 | * |
|---|
| 61 | * ALPHA - DOUBLE PRECISION. |
|---|
| 62 | * On entry, ALPHA specifies the scalar alpha. |
|---|
| 63 | * Unchanged on exit. |
|---|
| 64 | * |
|---|
| 65 | * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is |
|---|
| 66 | * m when SIDE = 'L' or 'l' and is n otherwise. |
|---|
| 67 | * Before entry with SIDE = 'L' or 'l', the m by m part of |
|---|
| 68 | * the array A must contain the symmetric matrix, such that |
|---|
| 69 | * when UPLO = 'U' or 'u', the leading m by m upper triangular |
|---|
| 70 | * part of the array A must contain the upper triangular part |
|---|
| 71 | * of the symmetric matrix and the strictly lower triangular |
|---|
| 72 | * part of A is not referenced, and when UPLO = 'L' or 'l', |
|---|
| 73 | * the leading m by m lower triangular part of the array A |
|---|
| 74 | * must contain the lower triangular part of the symmetric |
|---|
| 75 | * matrix and the strictly upper triangular part of A is not |
|---|
| 76 | * referenced. |
|---|
| 77 | * Before entry with SIDE = 'R' or 'r', the n by n part of |
|---|
| 78 | * the array A must contain the symmetric matrix, such that |
|---|
| 79 | * when UPLO = 'U' or 'u', the leading n by n upper triangular |
|---|
| 80 | * part of the array A must contain the upper triangular part |
|---|
| 81 | * of the symmetric matrix and the strictly lower triangular |
|---|
| 82 | * part of A is not referenced, and when UPLO = 'L' or 'l', |
|---|
| 83 | * the leading n by n lower triangular part of the array A |
|---|
| 84 | * must contain the lower triangular part of the symmetric |
|---|
| 85 | * matrix and the strictly upper triangular part of A is not |
|---|
| 86 | * referenced. |
|---|
| 87 | * Unchanged on exit. |
|---|
| 88 | * |
|---|
| 89 | * LDA - INTEGER. |
|---|
| 90 | * On entry, LDA specifies the first dimension of A as declared |
|---|
| 91 | * in the calling (sub) program. When SIDE = 'L' or 'l' then |
|---|
| 92 | * LDA must be at least max( 1, m ), otherwise LDA must be at |
|---|
| 93 | * least max( 1, n ). |
|---|
| 94 | * Unchanged on exit. |
|---|
| 95 | * |
|---|
| 96 | * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). |
|---|
| 97 | * Before entry, the leading m by n part of the array B must |
|---|
| 98 | * contain the matrix B. |
|---|
| 99 | * Unchanged on exit. |
|---|
| 100 | * |
|---|
| 101 | * LDB - INTEGER. |
|---|
| 102 | * On entry, LDB specifies the first dimension of B as declared |
|---|
| 103 | * in the calling (sub) program. LDB must be at least |
|---|
| 104 | * max( 1, m ). |
|---|
| 105 | * Unchanged on exit. |
|---|
| 106 | * |
|---|
| 107 | * BETA - DOUBLE PRECISION. |
|---|
| 108 | * On entry, BETA specifies the scalar beta. When BETA is |
|---|
| 109 | * supplied as zero then C need not be set on input. |
|---|
| 110 | * Unchanged on exit. |
|---|
| 111 | * |
|---|
| 112 | * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). |
|---|
| 113 | * Before entry, the leading m by n part of the array C must |
|---|
| 114 | * contain the matrix C, except when beta is zero, in which |
|---|
| 115 | * case C need not be set on entry. |
|---|
| 116 | * On exit, the array C is overwritten by the m by n updated |
|---|
| 117 | * matrix. |
|---|
| 118 | * |
|---|
| 119 | * LDC - INTEGER. |
|---|
| 120 | * On entry, LDC specifies the first dimension of C as declared |
|---|
| 121 | * in the calling (sub) program. LDC must be at least |
|---|
| 122 | * max( 1, m ). |
|---|
| 123 | * Unchanged on exit. |
|---|
| 124 | * |
|---|
| 125 | * |
|---|
| 126 | * Level 3 Blas routine. |
|---|
| 127 | * |
|---|
| 128 | * -- Written on 8-February-1989. |
|---|
| 129 | * Jack Dongarra, Argonne National Laboratory. |
|---|
| 130 | * Iain Duff, AERE Harwell. |
|---|
| 131 | * Jeremy Du Croz, Numerical Algorithms Group Ltd. |
|---|
| 132 | * Sven Hammarling, Numerical Algorithms Group Ltd. |
|---|
| 133 | * |
|---|
| 134 | * |
|---|
| 135 | * .. External Functions .. |
|---|
| 136 | LOGICAL LSAME |
|---|
| 137 | EXTERNAL LSAME |
|---|
| 138 | * .. |
|---|
| 139 | * .. External Subroutines .. |
|---|
| 140 | EXTERNAL XERBLA |
|---|
| 141 | * .. |
|---|
| 142 | * .. Intrinsic Functions .. |
|---|
| 143 | INTRINSIC MAX |
|---|
| 144 | * .. |
|---|
| 145 | * .. Local Scalars .. |
|---|
| 146 | DOUBLE PRECISION TEMP1,TEMP2 |
|---|
| 147 | INTEGER I,INFO,J,K,NROWA |
|---|
| 148 | LOGICAL UPPER |
|---|
| 149 | * .. |
|---|
| 150 | * .. Parameters .. |
|---|
| 151 | DOUBLE PRECISION ONE,ZERO |
|---|
| 152 | PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) |
|---|
| 153 | * .. |
|---|
| 154 | * |
|---|
| 155 | * Set NROWA as the number of rows of A. |
|---|
| 156 | * |
|---|
| 157 | IF (LSAME(SIDE,'L')) THEN |
|---|
| 158 | NROWA = M |
|---|
| 159 | ELSE |
|---|
| 160 | NROWA = N |
|---|
| 161 | END IF |
|---|
| 162 | UPPER = LSAME(UPLO,'U') |
|---|
| 163 | * |
|---|
| 164 | * Test the input parameters. |
|---|
| 165 | * |
|---|
| 166 | INFO = 0 |
|---|
| 167 | IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN |
|---|
| 168 | INFO = 1 |
|---|
| 169 | ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN |
|---|
| 170 | INFO = 2 |
|---|
| 171 | ELSE IF (M.LT.0) THEN |
|---|
| 172 | INFO = 3 |
|---|
| 173 | ELSE IF (N.LT.0) THEN |
|---|
| 174 | INFO = 4 |
|---|
| 175 | ELSE IF (LDA.LT.MAX(1,NROWA)) THEN |
|---|
| 176 | INFO = 7 |
|---|
| 177 | ELSE IF (LDB.LT.MAX(1,M)) THEN |
|---|
| 178 | INFO = 9 |
|---|
| 179 | ELSE IF (LDC.LT.MAX(1,M)) THEN |
|---|
| 180 | INFO = 12 |
|---|
| 181 | END IF |
|---|
| 182 | IF (INFO.NE.0) THEN |
|---|
| 183 | CALL XERBLA('DSYMM ',INFO) |
|---|
| 184 | RETURN |
|---|
| 185 | END IF |
|---|
| 186 | * |
|---|
| 187 | * Quick return if possible. |
|---|
| 188 | * |
|---|
| 189 | IF ((M.EQ.0) .OR. (N.EQ.0) .OR. |
|---|
| 190 | + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN |
|---|
| 191 | * |
|---|
| 192 | * And when alpha.eq.zero. |
|---|
| 193 | * |
|---|
| 194 | IF (ALPHA.EQ.ZERO) THEN |
|---|
| 195 | IF (BETA.EQ.ZERO) THEN |
|---|
| 196 | DO 20 J = 1,N |
|---|
| 197 | DO 10 I = 1,M |
|---|
| 198 | C(I,J) = ZERO |
|---|
| 199 | 10 CONTINUE |
|---|
| 200 | 20 CONTINUE |
|---|
| 201 | ELSE |
|---|
| 202 | DO 40 J = 1,N |
|---|
| 203 | DO 30 I = 1,M |
|---|
| 204 | C(I,J) = BETA*C(I,J) |
|---|
| 205 | 30 CONTINUE |
|---|
| 206 | 40 CONTINUE |
|---|
| 207 | END IF |
|---|
| 208 | RETURN |
|---|
| 209 | END IF |
|---|
| 210 | * |
|---|
| 211 | * Start the operations. |
|---|
| 212 | * |
|---|
| 213 | IF (LSAME(SIDE,'L')) THEN |
|---|
| 214 | * |
|---|
| 215 | * Form C := alpha*A*B + beta*C. |
|---|
| 216 | * |
|---|
| 217 | IF (UPPER) THEN |
|---|
| 218 | DO 70 J = 1,N |
|---|
| 219 | DO 60 I = 1,M |
|---|
| 220 | TEMP1 = ALPHA*B(I,J) |
|---|
| 221 | TEMP2 = ZERO |
|---|
| 222 | DO 50 K = 1,I - 1 |
|---|
| 223 | C(K,J) = C(K,J) + TEMP1*A(K,I) |
|---|
| 224 | TEMP2 = TEMP2 + B(K,J)*A(K,I) |
|---|
| 225 | 50 CONTINUE |
|---|
| 226 | IF (BETA.EQ.ZERO) THEN |
|---|
| 227 | C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2 |
|---|
| 228 | ELSE |
|---|
| 229 | C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) + |
|---|
| 230 | + ALPHA*TEMP2 |
|---|
| 231 | END IF |
|---|
| 232 | 60 CONTINUE |
|---|
| 233 | 70 CONTINUE |
|---|
| 234 | ELSE |
|---|
| 235 | DO 100 J = 1,N |
|---|
| 236 | DO 90 I = M,1,-1 |
|---|
| 237 | TEMP1 = ALPHA*B(I,J) |
|---|
| 238 | TEMP2 = ZERO |
|---|
| 239 | DO 80 K = I + 1,M |
|---|
| 240 | C(K,J) = C(K,J) + TEMP1*A(K,I) |
|---|
| 241 | TEMP2 = TEMP2 + B(K,J)*A(K,I) |
|---|
| 242 | 80 CONTINUE |
|---|
| 243 | IF (BETA.EQ.ZERO) THEN |
|---|
| 244 | C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2 |
|---|
| 245 | ELSE |
|---|
| 246 | C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) + |
|---|
| 247 | + ALPHA*TEMP2 |
|---|
| 248 | END IF |
|---|
| 249 | 90 CONTINUE |
|---|
| 250 | 100 CONTINUE |
|---|
| 251 | END IF |
|---|
| 252 | ELSE |
|---|
| 253 | * |
|---|
| 254 | * Form C := alpha*B*A + beta*C. |
|---|
| 255 | * |
|---|
| 256 | DO 170 J = 1,N |
|---|
| 257 | TEMP1 = ALPHA*A(J,J) |
|---|
| 258 | IF (BETA.EQ.ZERO) THEN |
|---|
| 259 | DO 110 I = 1,M |
|---|
| 260 | C(I,J) = TEMP1*B(I,J) |
|---|
| 261 | 110 CONTINUE |
|---|
| 262 | ELSE |
|---|
| 263 | DO 120 I = 1,M |
|---|
| 264 | C(I,J) = BETA*C(I,J) + TEMP1*B(I,J) |
|---|
| 265 | 120 CONTINUE |
|---|
| 266 | END IF |
|---|
| 267 | DO 140 K = 1,J - 1 |
|---|
| 268 | IF (UPPER) THEN |
|---|
| 269 | TEMP1 = ALPHA*A(K,J) |
|---|
| 270 | ELSE |
|---|
| 271 | TEMP1 = ALPHA*A(J,K) |
|---|
| 272 | END IF |
|---|
| 273 | DO 130 I = 1,M |
|---|
| 274 | C(I,J) = C(I,J) + TEMP1*B(I,K) |
|---|
| 275 | 130 CONTINUE |
|---|
| 276 | 140 CONTINUE |
|---|
| 277 | DO 160 K = J + 1,N |
|---|
| 278 | IF (UPPER) THEN |
|---|
| 279 | TEMP1 = ALPHA*A(J,K) |
|---|
| 280 | ELSE |
|---|
| 281 | TEMP1 = ALPHA*A(K,J) |
|---|
| 282 | END IF |
|---|
| 283 | DO 150 I = 1,M |
|---|
| 284 | C(I,J) = C(I,J) + TEMP1*B(I,K) |
|---|
| 285 | 150 CONTINUE |
|---|
| 286 | 160 CONTINUE |
|---|
| 287 | 170 CONTINUE |
|---|
| 288 | END IF |
|---|
| 289 | * |
|---|
| 290 | RETURN |
|---|
| 291 | * |
|---|
| 292 | * End of DSYMM . |
|---|
| 293 | * |
|---|
| 294 | END |
|---|
