001:       SUBROUTINE SSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
002: *     .. Scalar Arguments ..
003:       REAL ALPHA,BETA
004:       INTEGER K,LDA,LDB,LDC,N
005:       CHARACTER TRANS,UPLO
006: *     ..
007: *     .. Array Arguments ..
008:       REAL A(LDA,*),B(LDB,*),C(LDC,*)
009: *     ..
010: *
011: *  Purpose
012: *  =======
013: *
014: *  SSYR2K  performs one of the symmetric rank 2k operations
015: *
016: *     C := alpha*A*B' + alpha*B*A' + beta*C,
017: *
018: *  or
019: *
020: *     C := alpha*A'*B + alpha*B'*A + beta*C,
021: *
022: *  where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
023: *  and  A and B  are  n by k  matrices  in the  first  case  and  k by n
024: *  matrices in the second case.
025: *
026: *  Arguments
027: *  ==========
028: *
029: *  UPLO   - CHARACTER*1.
030: *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
031: *           triangular  part  of the  array  C  is to be  referenced  as
032: *           follows:
033: *
034: *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
035: *                                  is to be referenced.
036: *
037: *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
038: *                                  is to be referenced.
039: *
040: *           Unchanged on exit.
041: *
042: *  TRANS  - CHARACTER*1.
043: *           On entry,  TRANS  specifies the operation to be performed as
044: *           follows:
045: *
046: *              TRANS = 'N' or 'n'   C := alpha*A*B' + alpha*B*A' +
047: *                                        beta*C.
048: *
049: *              TRANS = 'T' or 't'   C := alpha*A'*B + alpha*B'*A +
050: *                                        beta*C.
051: *
052: *              TRANS = 'C' or 'c'   C := alpha*A'*B + alpha*B'*A +
053: *                                        beta*C.
054: *
055: *           Unchanged on exit.
056: *
057: *  N      - INTEGER.
058: *           On entry,  N specifies the order of the matrix C.  N must be
059: *           at least zero.
060: *           Unchanged on exit.
061: *
062: *  K      - INTEGER.
063: *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
064: *           of  columns  of the  matrices  A and B,  and on  entry  with
065: *           TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
066: *           of rows of the matrices  A and B.  K must be at least  zero.
067: *           Unchanged on exit.
068: *
069: *  ALPHA  - REAL            .
070: *           On entry, ALPHA specifies the scalar alpha.
071: *           Unchanged on exit.
072: *
073: *  A      - REAL             array of DIMENSION ( LDA, ka ), where ka is
074: *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
075: *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
076: *           part of the array  A  must contain the matrix  A,  otherwise
077: *           the leading  k by n  part of the array  A  must contain  the
078: *           matrix A.
079: *           Unchanged on exit.
080: *
081: *  LDA    - INTEGER.
082: *           On entry, LDA specifies the first dimension of A as declared
083: *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
084: *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
085: *           be at least  max( 1, k ).
086: *           Unchanged on exit.
087: *
088: *  B      - REAL             array of DIMENSION ( LDB, kb ), where kb is
089: *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
090: *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
091: *           part of the array  B  must contain the matrix  B,  otherwise
092: *           the leading  k by n  part of the array  B  must contain  the
093: *           matrix B.
094: *           Unchanged on exit.
095: *
096: *  LDB    - INTEGER.
097: *           On entry, LDB specifies the first dimension of B as declared
098: *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
099: *           then  LDB must be at least  max( 1, n ), otherwise  LDB must
100: *           be at least  max( 1, k ).
101: *           Unchanged on exit.
102: *
103: *  BETA   - REAL            .
104: *           On entry, BETA specifies the scalar beta.
105: *           Unchanged on exit.
106: *
107: *  C      - REAL             array of DIMENSION ( LDC, n ).
108: *           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
109: *           upper triangular part of the array C must contain the upper
110: *           triangular part  of the  symmetric matrix  and the strictly
111: *           lower triangular part of C is not referenced.  On exit, the
112: *           upper triangular part of the array  C is overwritten by the
113: *           upper triangular part of the updated matrix.
114: *           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
115: *           lower triangular part of the array C must contain the lower
116: *           triangular part  of the  symmetric matrix  and the strictly
117: *           upper triangular part of C is not referenced.  On exit, the
118: *           lower triangular part of the array  C is overwritten by the
119: *           lower triangular part of the updated matrix.
120: *
121: *  LDC    - INTEGER.
122: *           On entry, LDC specifies the first dimension of C as declared
123: *           in  the  calling  (sub)  program.   LDC  must  be  at  least
124: *           max( 1, n ).
125: *           Unchanged on exit.
126: *
127: *  Further Details
128: *  ===============
129: *
130: *  Level 3 Blas routine.
131: *
132: *
133: *  -- Written on 8-February-1989.
134: *     Jack Dongarra, Argonne National Laboratory.
135: *     Iain Duff, AERE Harwell.
136: *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
137: *     Sven Hammarling, Numerical Algorithms Group Ltd.
138: *
139: *  =====================================================================
140: *
141: *     .. External Functions ..
142:       LOGICAL LSAME
143:       EXTERNAL LSAME
144: *     ..
145: *     .. External Subroutines ..
146:       EXTERNAL XERBLA
147: *     ..
148: *     .. Intrinsic Functions ..
149:       INTRINSIC MAX
150: *     ..
151: *     .. Local Scalars ..
152:       REAL TEMP1,TEMP2
153:       INTEGER I,INFO,J,L,NROWA
154:       LOGICAL UPPER
155: *     ..
156: *     .. Parameters ..
157:       REAL ONE,ZERO
158:       PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
159: *     ..
160: *
161: *     Test the input parameters.
162: *
163:       IF (LSAME(TRANS,'N')) THEN
164:           NROWA = N
165:       ELSE
166:           NROWA = K
167:       END IF
168:       UPPER = LSAME(UPLO,'U')
169: *
170:       INFO = 0
171:       IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
172:           INFO = 1
173:       ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
174:      +         (.NOT.LSAME(TRANS,'T')) .AND.
175:      +         (.NOT.LSAME(TRANS,'C'))) THEN
176:           INFO = 2
177:       ELSE IF (N.LT.0) THEN
178:           INFO = 3
179:       ELSE IF (K.LT.0) THEN
180:           INFO = 4
181:       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
182:           INFO = 7
183:       ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
184:           INFO = 9
185:       ELSE IF (LDC.LT.MAX(1,N)) THEN
186:           INFO = 12
187:       END IF
188:       IF (INFO.NE.0) THEN
189:           CALL XERBLA('SSYR2K',INFO)
190:           RETURN
191:       END IF
192: *
193: *     Quick return if possible.
194: *
195:       IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
196:      +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
197: *
198: *     And when  alpha.eq.zero.
199: *
200:       IF (ALPHA.EQ.ZERO) THEN
201:           IF (UPPER) THEN
202:               IF (BETA.EQ.ZERO) THEN
203:                   DO 20 J = 1,N
204:                       DO 10 I = 1,J
205:                           C(I,J) = ZERO
206:    10                 CONTINUE
207:    20             CONTINUE
208:               ELSE
209:                   DO 40 J = 1,N
210:                       DO 30 I = 1,J
211:                           C(I,J) = BETA*C(I,J)
212:    30                 CONTINUE
213:    40             CONTINUE
214:               END IF
215:           ELSE
216:               IF (BETA.EQ.ZERO) THEN
217:                   DO 60 J = 1,N
218:                       DO 50 I = J,N
219:                           C(I,J) = ZERO
220:    50                 CONTINUE
221:    60             CONTINUE
222:               ELSE
223:                   DO 80 J = 1,N
224:                       DO 70 I = J,N
225:                           C(I,J) = BETA*C(I,J)
226:    70                 CONTINUE
227:    80             CONTINUE
228:               END IF
229:           END IF
230:           RETURN
231:       END IF
232: *
233: *     Start the operations.
234: *
235:       IF (LSAME(TRANS,'N')) THEN
236: *
237: *        Form  C := alpha*A*B' + alpha*B*A' + C.
238: *
239:           IF (UPPER) THEN
240:               DO 130 J = 1,N
241:                   IF (BETA.EQ.ZERO) THEN
242:                       DO 90 I = 1,J
243:                           C(I,J) = ZERO
244:    90                 CONTINUE
245:                   ELSE IF (BETA.NE.ONE) THEN
246:                       DO 100 I = 1,J
247:                           C(I,J) = BETA*C(I,J)
248:   100                 CONTINUE
249:                   END IF
250:                   DO 120 L = 1,K
251:                       IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
252:                           TEMP1 = ALPHA*B(J,L)
253:                           TEMP2 = ALPHA*A(J,L)
254:                           DO 110 I = 1,J
255:                               C(I,J) = C(I,J) + A(I,L)*TEMP1 +
256:      +                                 B(I,L)*TEMP2
257:   110                     CONTINUE
258:                       END IF
259:   120             CONTINUE
260:   130         CONTINUE
261:           ELSE
262:               DO 180 J = 1,N
263:                   IF (BETA.EQ.ZERO) THEN
264:                       DO 140 I = J,N
265:                           C(I,J) = ZERO
266:   140                 CONTINUE
267:                   ELSE IF (BETA.NE.ONE) THEN
268:                       DO 150 I = J,N
269:                           C(I,J) = BETA*C(I,J)
270:   150                 CONTINUE
271:                   END IF
272:                   DO 170 L = 1,K
273:                       IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
274:                           TEMP1 = ALPHA*B(J,L)
275:                           TEMP2 = ALPHA*A(J,L)
276:                           DO 160 I = J,N
277:                               C(I,J) = C(I,J) + A(I,L)*TEMP1 +
278:      +                                 B(I,L)*TEMP2
279:   160                     CONTINUE
280:                       END IF
281:   170             CONTINUE
282:   180         CONTINUE
283:           END IF
284:       ELSE
285: *
286: *        Form  C := alpha*A'*B + alpha*B'*A + C.
287: *
288:           IF (UPPER) THEN
289:               DO 210 J = 1,N
290:                   DO 200 I = 1,J
291:                       TEMP1 = ZERO
292:                       TEMP2 = ZERO
293:                       DO 190 L = 1,K
294:                           TEMP1 = TEMP1 + A(L,I)*B(L,J)
295:                           TEMP2 = TEMP2 + B(L,I)*A(L,J)
296:   190                 CONTINUE
297:                       IF (BETA.EQ.ZERO) THEN
298:                           C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
299:                       ELSE
300:                           C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
301:      +                             ALPHA*TEMP2
302:                       END IF
303:   200             CONTINUE
304:   210         CONTINUE
305:           ELSE
306:               DO 240 J = 1,N
307:                   DO 230 I = J,N
308:                       TEMP1 = ZERO
309:                       TEMP2 = ZERO
310:                       DO 220 L = 1,K
311:                           TEMP1 = TEMP1 + A(L,I)*B(L,J)
312:                           TEMP2 = TEMP2 + B(L,I)*A(L,J)
313:   220                 CONTINUE
314:                       IF (BETA.EQ.ZERO) THEN
315:                           C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
316:                       ELSE
317:                           C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
318:      +                             ALPHA*TEMP2
319:                       END IF
320:   230             CONTINUE
321:   240         CONTINUE
322:           END IF
323:       END IF
324: *
325:       RETURN
326: *
327: *     End of SSYR2K.
328: *
329:       END
330: