3 /* Subroutine */ int ssyrk_(char *uplo, char *trans, integer *n, integer *k,
4 real *alpha, real *a, integer *lda, real *beta, real *c__, integer *
7 /* System generated locals */
8 integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3;
11 integer i__, j, l, info;
13 extern logical lsame_(char *, char *);
16 extern /* Subroutine */ int xerbla_(char *, integer *);
18 /* .. Scalar Arguments .. */
20 /* .. Array Arguments .. */
26 /* SSYRK performs one of the symmetric rank k operations */
28 /* C := alpha*A*A' + beta*C, */
32 /* C := alpha*A'*A + beta*C, */
34 /* where alpha and beta are scalars, C is an n by n symmetric matrix */
35 /* and A is an n by k matrix in the first case and a k by n matrix */
36 /* in the second case. */
41 /* UPLO - CHARACTER*1. */
42 /* On entry, UPLO specifies whether the upper or lower */
43 /* triangular part of the array C is to be referenced as */
46 /* UPLO = 'U' or 'u' Only the upper triangular part of C */
47 /* is to be referenced. */
49 /* UPLO = 'L' or 'l' Only the lower triangular part of C */
50 /* is to be referenced. */
52 /* Unchanged on exit. */
54 /* TRANS - CHARACTER*1. */
55 /* On entry, TRANS specifies the operation to be performed as */
58 /* TRANS = 'N' or 'n' C := alpha*A*A' + beta*C. */
60 /* TRANS = 'T' or 't' C := alpha*A'*A + beta*C. */
62 /* TRANS = 'C' or 'c' C := alpha*A'*A + beta*C. */
64 /* Unchanged on exit. */
67 /* On entry, N specifies the order of the matrix C. N must be */
69 /* Unchanged on exit. */
72 /* On entry with TRANS = 'N' or 'n', K specifies the number */
73 /* of columns of the matrix A, and on entry with */
74 /* TRANS = 'T' or 't' or 'C' or 'c', K specifies the number */
75 /* of rows of the matrix A. K must be at least zero. */
76 /* Unchanged on exit. */
79 /* On entry, ALPHA specifies the scalar alpha. */
80 /* Unchanged on exit. */
82 /* A - REAL array of DIMENSION ( LDA, ka ), where ka is */
83 /* k when TRANS = 'N' or 'n', and is n otherwise. */
84 /* Before entry with TRANS = 'N' or 'n', the leading n by k */
85 /* part of the array A must contain the matrix A, otherwise */
86 /* the leading k by n part of the array A must contain the */
88 /* Unchanged on exit. */
91 /* On entry, LDA specifies the first dimension of A as declared */
92 /* in the calling (sub) program. When TRANS = 'N' or 'n' */
93 /* then LDA must be at least max( 1, n ), otherwise LDA must */
94 /* be at least max( 1, k ). */
95 /* Unchanged on exit. */
98 /* On entry, BETA specifies the scalar beta. */
99 /* Unchanged on exit. */
101 /* C - REAL array of DIMENSION ( LDC, n ). */
102 /* Before entry with UPLO = 'U' or 'u', the leading n by n */
103 /* upper triangular part of the array C must contain the upper */
104 /* triangular part of the symmetric matrix and the strictly */
105 /* lower triangular part of C is not referenced. On exit, the */
106 /* upper triangular part of the array C is overwritten by the */
107 /* upper triangular part of the updated matrix. */
108 /* Before entry with UPLO = 'L' or 'l', the leading n by n */
109 /* lower triangular part of the array C must contain the lower */
110 /* triangular part of the symmetric matrix and the strictly */
111 /* upper triangular part of C is not referenced. On exit, the */
112 /* lower triangular part of the array C is overwritten by the */
113 /* lower triangular part of the updated matrix. */
116 /* On entry, LDC specifies the first dimension of C as declared */
117 /* in the calling (sub) program. LDC must be at least */
119 /* Unchanged on exit. */
122 /* Level 3 Blas routine. */
124 /* -- Written on 8-February-1989. */
125 /* Jack Dongarra, Argonne National Laboratory. */
126 /* Iain Duff, AERE Harwell. */
127 /* Jeremy Du Croz, Numerical Algorithms Group Ltd. */
128 /* Sven Hammarling, Numerical Algorithms Group Ltd. */
131 /* .. External Functions .. */
133 /* .. External Subroutines .. */
135 /* .. Intrinsic Functions .. */
137 /* .. Local Scalars .. */
139 /* .. Parameters .. */
142 /* Test the input parameters. */
144 /* Parameter adjustments */
146 a_offset = 1 + a_dim1;
149 c_offset = 1 + c_dim1;
153 if (lsame_(trans, "N")) {
158 upper = lsame_(uplo, "U");
161 if (! upper && ! lsame_(uplo, "L")) {
163 } else if (! lsame_(trans, "N") && ! lsame_(trans,
164 "T") && ! lsame_(trans, "C")) {
170 } else if (*lda < max(1,nrowa)) {
172 } else if (*ldc < max(1,*n)) {
176 xerbla_("SSYRK ", &info);
180 /* Quick return if possible. */
182 if (*n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
186 /* And when alpha.eq.zero. */
192 for (j = 1; j <= i__1; ++j) {
194 for (i__ = 1; i__ <= i__2; ++i__) {
195 c__[i__ + j * c_dim1] = 0.f;
202 for (j = 1; j <= i__1; ++j) {
204 for (i__ = 1; i__ <= i__2; ++i__) {
205 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
214 for (j = 1; j <= i__1; ++j) {
216 for (i__ = j; i__ <= i__2; ++i__) {
217 c__[i__ + j * c_dim1] = 0.f;
224 for (j = 1; j <= i__1; ++j) {
226 for (i__ = j; i__ <= i__2; ++i__) {
227 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
237 /* Start the operations. */
239 if (lsame_(trans, "N")) {
241 /* Form C := alpha*A*A' + beta*C. */
245 for (j = 1; j <= i__1; ++j) {
248 for (i__ = 1; i__ <= i__2; ++i__) {
249 c__[i__ + j * c_dim1] = 0.f;
252 } else if (*beta != 1.f) {
254 for (i__ = 1; i__ <= i__2; ++i__) {
255 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
260 for (l = 1; l <= i__2; ++l) {
261 if (a[j + l * a_dim1] != 0.f) {
262 temp = *alpha * a[j + l * a_dim1];
264 for (i__ = 1; i__ <= i__3; ++i__) {
265 c__[i__ + j * c_dim1] += temp * a[i__ + l *
276 for (j = 1; j <= i__1; ++j) {
279 for (i__ = j; i__ <= i__2; ++i__) {
280 c__[i__ + j * c_dim1] = 0.f;
283 } else if (*beta != 1.f) {
285 for (i__ = j; i__ <= i__2; ++i__) {
286 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
291 for (l = 1; l <= i__2; ++l) {
292 if (a[j + l * a_dim1] != 0.f) {
293 temp = *alpha * a[j + l * a_dim1];
295 for (i__ = j; i__ <= i__3; ++i__) {
296 c__[i__ + j * c_dim1] += temp * a[i__ + l *
308 /* Form C := alpha*A'*A + beta*C. */
312 for (j = 1; j <= i__1; ++j) {
314 for (i__ = 1; i__ <= i__2; ++i__) {
317 for (l = 1; l <= i__3; ++l) {
318 temp += a[l + i__ * a_dim1] * a[l + j * a_dim1];
322 c__[i__ + j * c_dim1] = *alpha * temp;
324 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
333 for (j = 1; j <= i__1; ++j) {
335 for (i__ = j; i__ <= i__2; ++i__) {
338 for (l = 1; l <= i__3; ++l) {
339 temp += a[l + i__ * a_dim1] * a[l + j * a_dim1];
343 c__[i__ + j * c_dim1] = *alpha * temp;
345 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[