001:       SUBROUTINE SSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA,
002:      +                  C )
003: *
004: *  -- LAPACK routine (version 3.2)                                    --
005: *
006: *  -- Contributed by Julien Langou of the Univ. of Colorado Denver    --
007: *  -- November 2008                                                   --
008: *
009: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
010: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
011: *
012: *     ..
013: *     .. Scalar Arguments ..
014:       REAL               ALPHA, BETA
015:       INTEGER            K, LDA, N
016:       CHARACTER          TRANS, TRANSR, UPLO
017: *     ..
018: *     .. Array Arguments ..
019:       REAL               A( LDA, * ), C( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  Level 3 BLAS like routine for C in RFP Format.
026: *
027: *  SSFRK performs one of the symmetric rank--k operations
028: *
029: *     C := alpha*A*A' + beta*C,
030: *
031: *  or
032: *
033: *     C := alpha*A'*A + beta*C,
034: *
035: *  where alpha and beta are real scalars, C is an n--by--n symmetric
036: *  matrix and A is an n--by--k matrix in the first case and a k--by--n
037: *  matrix in the second case.
038: *
039: *  Arguments
040: *  ==========
041: *
042: *  TRANSR    (input) CHARACTER
043: *          = 'N':  The Normal Form of RFP A is stored;
044: *          = 'T':  The Transpose Form of RFP A is stored.
045: *
046: *  UPLO   - (input) CHARACTER
047: *           On  entry, UPLO specifies whether the upper or lower
048: *           triangular part of the array C is to be referenced as
049: *           follows:
050: *
051: *              UPLO = 'U' or 'u'   Only the upper triangular part of C
052: *                                  is to be referenced.
053: *
054: *              UPLO = 'L' or 'l'   Only the lower triangular part of C
055: *                                  is to be referenced.
056: *
057: *           Unchanged on exit.
058: *
059: *  TRANS  - (input) CHARACTER
060: *           On entry, TRANS specifies the operation to be performed as
061: *           follows:
062: *
063: *              TRANS = 'N' or 'n'   C := alpha*A*A' + beta*C.
064: *
065: *              TRANS = 'T' or 't'   C := alpha*A'*A + beta*C.
066: *
067: *           Unchanged on exit.
068: *
069: *  N      - (input) INTEGER.
070: *           On entry, N specifies the order of the matrix C. N must be
071: *           at least zero.
072: *           Unchanged on exit.
073: *
074: *  K      - (input) INTEGER.
075: *           On entry with TRANS = 'N' or 'n', K specifies the number
076: *           of  columns of the matrix A, and on entry with TRANS = 'T'
077: *           or 't', K specifies the number of rows of the matrix A. K
078: *           must be at least zero.
079: *           Unchanged on exit.
080: *
081: *  ALPHA  - (input) REAL.
082: *           On entry, ALPHA specifies the scalar alpha.
083: *           Unchanged on exit.
084: *
085: *  A      - (input) REAL array of DIMENSION ( LDA, ka ), where KA
086: *           is K  when TRANS = 'N' or 'n', and is N otherwise. Before
087: *           entry with TRANS = 'N' or 'n', the leading N--by--K part of
088: *           the array A must contain the matrix A, otherwise the leading
089: *           K--by--N part of the array A must contain the matrix A.
090: *           Unchanged on exit.
091: *
092: *  LDA    - (input) INTEGER.
093: *           On entry, LDA specifies the first dimension of A as declared
094: *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
095: *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
096: *           be at least  max( 1, k ).
097: *           Unchanged on exit.
098: *
099: *  BETA   - (input) REAL.
100: *           On entry, BETA specifies the scalar beta.
101: *           Unchanged on exit.
102: *
103: *
104: *  C      - (input/output) REAL array, dimension ( NT );
105: *           NT = N*(N+1)/2. On entry, the symmetric matrix C in RFP
106: *           Format. RFP Format is described by TRANSR, UPLO and N.
107: *
108: *  Arguments
109: *  ==========
110: *
111: *     ..
112: *     .. Parameters ..
113:       REAL               ONE, ZERO
114:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
115: *     ..
116: *     .. Local Scalars ..
117:       LOGICAL            LOWER, NORMALTRANSR, NISODD, NOTRANS
118:       INTEGER            INFO, NROWA, J, NK, N1, N2
119: *     ..
120: *     .. External Functions ..
121:       LOGICAL            LSAME
122:       EXTERNAL           LSAME
123: *     ..
124: *     .. External Subroutines ..
125:       EXTERNAL           SGEMM, SSYRK, XERBLA
126: *     ..
127: *     .. Intrinsic Functions ..
128:       INTRINSIC          MAX
129: *     ..
130: *     .. Executable Statements ..
131: *
132: *     Test the input parameters.
133: *
134:       INFO = 0
135:       NORMALTRANSR = LSAME( TRANSR, 'N' )
136:       LOWER = LSAME( UPLO, 'L' )
137:       NOTRANS = LSAME( TRANS, 'N' )
138: *
139:       IF( NOTRANS ) THEN
140:          NROWA = N
141:       ELSE
142:          NROWA = K
143:       END IF
144: *
145:       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
146:          INFO = -1
147:       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
148:          INFO = -2
149:       ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
150:          INFO = -3
151:       ELSE IF( N.LT.0 ) THEN
152:          INFO = -4
153:       ELSE IF( K.LT.0 ) THEN
154:          INFO = -5
155:       ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
156:          INFO = -8
157:       END IF
158:       IF( INFO.NE.0 ) THEN
159:          CALL XERBLA( 'SSFRK ', -INFO )
160:          RETURN
161:       END IF
162: *
163: *     Quick return if possible.
164: *
165: *     The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not
166: *     done (it is in SSYRK for example) and left in the general case.
167: *
168:       IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
169:      +    ( BETA.EQ.ONE ) ) )RETURN
170: *
171:       IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN
172:          DO J = 1, ( ( N*( N+1 ) ) / 2 )
173:             C( J ) = ZERO
174:          END DO
175:          RETURN
176:       END IF
177: *
178: *     C is N-by-N.
179: *     If N is odd, set NISODD = .TRUE., and N1 and N2.
180: *     If N is even, NISODD = .FALSE., and NK.
181: *
182:       IF( MOD( N, 2 ).EQ.0 ) THEN
183:          NISODD = .FALSE.
184:          NK = N / 2
185:       ELSE
186:          NISODD = .TRUE.
187:          IF( LOWER ) THEN
188:             N2 = N / 2
189:             N1 = N - N2
190:          ELSE
191:             N1 = N / 2
192:             N2 = N - N1
193:          END IF
194:       END IF
195: *
196:       IF( NISODD ) THEN
197: *
198: *        N is odd
199: *
200:          IF( NORMALTRANSR ) THEN
201: *
202: *           N is odd and TRANSR = 'N'
203: *
204:             IF( LOWER ) THEN
205: *
206: *              N is odd, TRANSR = 'N', and UPLO = 'L'
207: *
208:                IF( NOTRANS ) THEN
209: *
210: *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
211: *
212:                   CALL SSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
213:      +                        BETA, C( 1 ), N )
214:                   CALL SSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
215:      +                        BETA, C( N+1 ), N )
216:                   CALL SGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
217:      +                        LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
218: *
219:                ELSE
220: *
221: *                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
222: *
223:                   CALL SSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
224:      +                        BETA, C( 1 ), N )
225:                   CALL SSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
226:      +                        BETA, C( N+1 ), N )
227:                   CALL SGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
228:      +                        LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N )
229: *
230:                END IF
231: *
232:             ELSE
233: *
234: *              N is odd, TRANSR = 'N', and UPLO = 'U'
235: *
236:                IF( NOTRANS ) THEN
237: *
238: *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
239: *
240:                   CALL SSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
241:      +                        BETA, C( N2+1 ), N )
242:                   CALL SSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA,
243:      +                        BETA, C( N1+1 ), N )
244:                   CALL SGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
245:      +                        LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N )
246: *
247:                ELSE
248: *
249: *                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
250: *
251:                   CALL SSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
252:      +                        BETA, C( N2+1 ), N )
253:                   CALL SSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA,
254:      +                        BETA, C( N1+1 ), N )
255:                   CALL SGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
256:      +                        LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N )
257: *
258:                END IF
259: *
260:             END IF
261: *
262:          ELSE
263: *
264: *           N is odd, and TRANSR = 'T'
265: *
266:             IF( LOWER ) THEN
267: *
268: *              N is odd, TRANSR = 'T', and UPLO = 'L'
269: *
270:                IF( NOTRANS ) THEN
271: *
272: *                 N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
273: *
274:                   CALL SSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
275:      +                        BETA, C( 1 ), N1 )
276:                   CALL SSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
277:      +                        BETA, C( 2 ), N1 )
278:                   CALL SGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ),
279:      +                        LDA, A( N1+1, 1 ), LDA, BETA,
280:      +                        C( N1*N1+1 ), N1 )
281: *
282:                ELSE
283: *
284: *                 N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
285: *
286:                   CALL SSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
287:      +                        BETA, C( 1 ), N1 )
288:                   CALL SSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
289:      +                        BETA, C( 2 ), N1 )
290:                   CALL SGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ),
291:      +                        LDA, A( 1, N1+1 ), LDA, BETA,
292:      +                        C( N1*N1+1 ), N1 )
293: *
294:                END IF
295: *
296:             ELSE
297: *
298: *              N is odd, TRANSR = 'T', and UPLO = 'U'
299: *
300:                IF( NOTRANS ) THEN
301: *
302: *                 N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
303: *
304:                   CALL SSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA,
305:      +                        BETA, C( N2*N2+1 ), N2 )
306:                   CALL SSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA,
307:      +                        BETA, C( N1*N2+1 ), N2 )
308:                   CALL SGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ),
309:      +                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
310: *
311:                ELSE
312: *
313: *                 N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
314: *
315:                   CALL SSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA,
316:      +                        BETA, C( N2*N2+1 ), N2 )
317:                   CALL SSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA,
318:      +                        BETA, C( N1*N2+1 ), N2 )
319:                   CALL SGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ),
320:      +                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 )
321: *
322:                END IF
323: *
324:             END IF
325: *
326:          END IF
327: *
328:       ELSE
329: *
330: *        N is even
331: *
332:          IF( NORMALTRANSR ) THEN
333: *
334: *           N is even and TRANSR = 'N'
335: *
336:             IF( LOWER ) THEN
337: *
338: *              N is even, TRANSR = 'N', and UPLO = 'L'
339: *
340:                IF( NOTRANS ) THEN
341: *
342: *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N'
343: *
344:                   CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
345:      +                        BETA, C( 2 ), N+1 )
346:                   CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
347:      +                        BETA, C( 1 ), N+1 )
348:                   CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
349:      +                        LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
350:      +                        N+1 )
351: *
352:                ELSE
353: *
354: *                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T'
355: *
356:                   CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
357:      +                        BETA, C( 2 ), N+1 )
358:                   CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
359:      +                        BETA, C( 1 ), N+1 )
360:                   CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
361:      +                        LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ),
362:      +                        N+1 )
363: *
364:                END IF
365: *
366:             ELSE
367: *
368: *              N is even, TRANSR = 'N', and UPLO = 'U'
369: *
370:                IF( NOTRANS ) THEN
371: *
372: *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N'
373: *
374:                   CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
375:      +                        BETA, C( NK+2 ), N+1 )
376:                   CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
377:      +                        BETA, C( NK+1 ), N+1 )
378:                   CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
379:      +                        LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ),
380:      +                        N+1 )
381: *
382:                ELSE
383: *
384: *                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T'
385: *
386:                   CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
387:      +                        BETA, C( NK+2 ), N+1 )
388:                   CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
389:      +                        BETA, C( NK+1 ), N+1 )
390:                   CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
391:      +                        LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ),
392:      +                        N+1 )
393: *
394:                END IF
395: *
396:             END IF
397: *
398:          ELSE
399: *
400: *           N is even, and TRANSR = 'T'
401: *
402:             IF( LOWER ) THEN
403: *
404: *              N is even, TRANSR = 'T', and UPLO = 'L'
405: *
406:                IF( NOTRANS ) THEN
407: *
408: *                 N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N'
409: *
410:                   CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
411:      +                        BETA, C( NK+1 ), NK )
412:                   CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
413:      +                        BETA, C( 1 ), NK )
414:                   CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ),
415:      +                        LDA, A( NK+1, 1 ), LDA, BETA,
416:      +                        C( ( ( NK+1 )*NK )+1 ), NK )
417: *
418:                ELSE
419: *
420: *                 N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T'
421: *
422:                   CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
423:      +                        BETA, C( NK+1 ), NK )
424:                   CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
425:      +                        BETA, C( 1 ), NK )
426:                   CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ),
427:      +                        LDA, A( 1, NK+1 ), LDA, BETA,
428:      +                        C( ( ( NK+1 )*NK )+1 ), NK )
429: *
430:                END IF
431: *
432:             ELSE
433: *
434: *              N is even, TRANSR = 'T', and UPLO = 'U'
435: *
436:                IF( NOTRANS ) THEN
437: *
438: *                 N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N'
439: *
440:                   CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA,
441:      +                        BETA, C( NK*( NK+1 )+1 ), NK )
442:                   CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA,
443:      +                        BETA, C( NK*NK+1 ), NK )
444:                   CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ),
445:      +                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
446: *
447:                ELSE
448: *
449: *                 N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T'
450: *
451:                   CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA,
452:      +                        BETA, C( NK*( NK+1 )+1 ), NK )
453:                   CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA,
454:      +                        BETA, C( NK*NK+1 ), NK )
455:                   CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ),
456:      +                        LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK )
457: *
458:                END IF
459: *
460:             END IF
461: *
462:          END IF
463: *
464:       END IF
465: *
466:       RETURN
467: *
468: *     End of SSFRK
469: *
470:       END
471: