swiftnav: ssfrk.c Source File

Go to the documentation of this file.
00001 /* ssfrk.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Subroutine */ int ssfrk_(char *transr, char *uplo, char *trans, integer *n, 
00017          integer *k, real *alpha, real *a, integer *lda, real *beta, real *
00018         c__)
00019 {
00020     /* System generated locals */
00021     integer a_dim1, a_offset, i__1;
00022 
00023     /* Local variables */
00024     integer j, n1, n2, nk, info;
00025     logical normaltransr;
00026     extern logical lsame_(char *, char *);
00027     extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 
00028             integer *, real *, real *, integer *, real *, integer *, real *, 
00029             real *, integer *);
00030     integer nrowa;
00031     logical lower;
00032     extern /* Subroutine */ int ssyrk_(char *, char *, integer *, integer *, 
00033             real *, real *, integer *, real *, real *, integer *), xerbla_(char *, integer *);
00034     logical nisodd, notrans;
00035 
00036 
00037 /*  -- LAPACK routine (version 3.2)                                    -- */
00038 
00039 /*  -- Contributed by Julien Langou of the Univ. of Colorado Denver    -- */
00040 /*  -- November 2008                                                   -- */
00041 
00042 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
00043 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
00044 
00045 /*     .. */
00046 /*     .. Scalar Arguments .. */
00047 /*     .. */
00048 /*     .. Array Arguments .. */
00049 /*     .. */
00050 
00051 /*  Purpose */
00052 /*  ======= */
00053 
00054 /*  Level 3 BLAS like routine for C in RFP Format. */
00055 
00056 /*  SSFRK performs one of the symmetric rank--k operations */
00057 
00058 /*     C := alpha*A*A' + beta*C, */
00059 
00060 /*  or */
00061 
00062 /*     C := alpha*A'*A + beta*C, */
00063 
00064 /*  where alpha and beta are real scalars, C is an n--by--n symmetric */
00065 /*  matrix and A is an n--by--k matrix in the first case and a k--by--n */
00066 /*  matrix in the second case. */
00067 
00068 /*  Arguments */
00069 /*  ========== */
00070 
00071 /*  TRANSR    (input) CHARACTER */
00072 /*          = 'N':  The Normal Form of RFP A is stored; */
00073 /*          = 'T':  The Transpose Form of RFP A is stored. */
00074 
00075 /*  UPLO   - (input) CHARACTER */
00076 /*           On  entry, UPLO specifies whether the upper or lower */
00077 /*           triangular part of the array C is to be referenced as */
00078 /*           follows: */
00079 
00080 /*              UPLO = 'U' or 'u'   Only the upper triangular part of C */
00081 /*                                  is to be referenced. */
00082 
00083 /*              UPLO = 'L' or 'l'   Only the lower triangular part of C */
00084 /*                                  is to be referenced. */
00085 
00086 /*           Unchanged on exit. */
00087 
00088 /*  TRANS  - (input) CHARACTER */
00089 /*           On entry, TRANS specifies the operation to be performed as */
00090 /*           follows: */
00091 
00092 /*              TRANS = 'N' or 'n'   C := alpha*A*A' + beta*C. */
00093 
00094 /*              TRANS = 'T' or 't'   C := alpha*A'*A + beta*C. */
00095 
00096 /*           Unchanged on exit. */
00097 
00098 /*  N      - (input) INTEGER. */
00099 /*           On entry, N specifies the order of the matrix C. N must be */
00100 /*           at least zero. */
00101 /*           Unchanged on exit. */
00102 
00103 /*  K      - (input) INTEGER. */
00104 /*           On entry with TRANS = 'N' or 'n', K specifies the number */
00105 /*           of  columns of the matrix A, and on entry with TRANS = 'T' */
00106 /*           or 't', K specifies the number of rows of the matrix A. K */
00107 /*           must be at least zero. */
00108 /*           Unchanged on exit. */
00109 
00110 /*  ALPHA  - (input) REAL. */
00111 /*           On entry, ALPHA specifies the scalar alpha. */
00112 /*           Unchanged on exit. */
00113 
00114 /*  A      - (input) REAL array of DIMENSION ( LDA, ka ), where KA */
00115 /*           is K  when TRANS = 'N' or 'n', and is N otherwise. Before */
00116 /*           entry with TRANS = 'N' or 'n', the leading N--by--K part of */
00117 /*           the array A must contain the matrix A, otherwise the leading */
00118 /*           K--by--N part of the array A must contain the matrix A. */
00119 /*           Unchanged on exit. */
00120 
00121 /*  LDA    - (input) INTEGER. */
00122 /*           On entry, LDA specifies the first dimension of A as declared */
00123 /*           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n' */
00124 /*           then  LDA must be at least  max( 1, n ), otherwise  LDA must */
00125 /*           be at least  max( 1, k ). */
00126 /*           Unchanged on exit. */
00127 
00128 /*  BETA   - (input) REAL. */
00129 /*           On entry, BETA specifies the scalar beta. */
00130 /*           Unchanged on exit. */
00131 
00132 
00133 /*  C      - (input/output) REAL array, dimension ( NT ); */
00134 /*           NT = N*(N+1)/2. On entry, the symmetric matrix C in RFP */
00135 /*           Format. RFP Format is described by TRANSR, UPLO and N. */
00136 
00137 /*  Arguments */
00138 /*  ========== */
00139 
00140 /*     .. */
00141 /*     .. Parameters .. */
00142 /*     .. */
00143 /*     .. Local Scalars .. */
00144 /*     .. */
00145 /*     .. External Functions .. */
00146 /*     .. */
00147 /*     .. External Subroutines .. */
00148 /*     .. */
00149 /*     .. Intrinsic Functions .. */
00150 /*     .. */
00151 /*     .. Executable Statements .. */
00152 
00153 /*     Test the input parameters. */
00154 
00155     /* Parameter adjustments */
00156     a_dim1 = *lda;
00157     a_offset = 1 + a_dim1;
00158     a -= a_offset;
00159     --c__;
00160 
00161     /* Function Body */
00162     info = 0;
00163     normaltransr = lsame_(transr, "N");
00164     lower = lsame_(uplo, "L");
00165     notrans = lsame_(trans, "N");
00166 
00167     if (notrans) {
00168         nrowa = *n;
00169     } else {
00170         nrowa = *k;
00171     }
00172 
00173     if (! normaltransr && ! lsame_(transr, "T")) {
00174         info = -1;
00175     } else if (! lower && ! lsame_(uplo, "U")) {
00176         info = -2;
00177     } else if (! notrans && ! lsame_(trans, "T")) {
00178         info = -3;
00179     } else if (*n < 0) {
00180         info = -4;
00181     } else if (*k < 0) {
00182         info = -5;
00183     } else if (*lda < max(1,nrowa)) {
00184         info = -8;
00185     }
00186     if (info != 0) {
00187         i__1 = -info;
00188         xerbla_("SSFRK ", &i__1);
00189         return 0;
00190     }
00191 
00192 /*     Quick return if possible. */
00193 
00194 /*     The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not */
00195 /*     done (it is in SSYRK for example) and left in the general case. */
00196 
00197     if (*n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
00198         return 0;
00199     }
00200 
00201     if (*alpha == 0.f && *beta == 0.f) {
00202         i__1 = *n * (*n + 1) / 2;
00203         for (j = 1; j <= i__1; ++j) {
00204             c__[j] = 0.f;
00205         }
00206         return 0;
00207     }
00208 
00209 /*     C is N-by-N. */
00210 /*     If N is odd, set NISODD = .TRUE., and N1 and N2. */
00211 /*     If N is even, NISODD = .FALSE., and NK. */
00212 
00213     if (*n % 2 == 0) {
00214         nisodd = FALSE_;
00215         nk = *n / 2;
00216     } else {
00217         nisodd = TRUE_;
00218         if (lower) {
00219             n2 = *n / 2;
00220             n1 = *n - n2;
00221         } else {
00222             n1 = *n / 2;
00223             n2 = *n - n1;
00224         }
00225     }
00226 
00227     if (nisodd) {
00228 
00229 /*        N is odd */
00230 
00231         if (normaltransr) {
00232 
00233 /*           N is odd and TRANSR = 'N' */
00234 
00235             if (lower) {
00236 
00237 /*              N is odd, TRANSR = 'N', and UPLO = 'L' */
00238 
00239                 if (notrans) {
00240 
00241 /*                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' */
00242 
00243                     ssyrk_("L", "N", &n1, k, alpha, &a[a_dim1 + 1], lda, beta, 
00244                              &c__[1], n);
00245                     ssyrk_("U", "N", &n2, k, alpha, &a[n1 + 1 + a_dim1], lda, 
00246                             beta, &c__[*n + 1], n);
00247                     sgemm_("N", "T", &n2, &n1, k, alpha, &a[n1 + 1 + a_dim1], 
00248                             lda, &a[a_dim1 + 1], lda, beta, &c__[n1 + 1], n);
00249 
00250                 } else {
00251 
00252 /*                 N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T' */
00253 
00254                     ssyrk_("L", "T", &n1, k, alpha, &a[a_dim1 + 1], lda, beta, 
00255                              &c__[1], n);
00256                     ssyrk_("U", "T", &n2, k, alpha, &a[(n1 + 1) * a_dim1 + 1], 
00257                              lda, beta, &c__[*n + 1], n)
00258                             ;
00259                     sgemm_("T", "N", &n2, &n1, k, alpha, &a[(n1 + 1) * a_dim1 
00260                             + 1], lda, &a[a_dim1 + 1], lda, beta, &c__[n1 + 1]
00261 , n);
00262 
00263                 }
00264 
00265             } else {
00266 
00267 /*              N is odd, TRANSR = 'N', and UPLO = 'U' */
00268 
00269                 if (notrans) {
00270 
00271 /*                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' */
00272 
00273                     ssyrk_("L", "N", &n1, k, alpha, &a[a_dim1 + 1], lda, beta, 
00274                              &c__[n2 + 1], n);
00275                     ssyrk_("U", "N", &n2, k, alpha, &a[n2 + a_dim1], lda, 
00276                             beta, &c__[n1 + 1], n);
00277                     sgemm_("N", "T", &n1, &n2, k, alpha, &a[a_dim1 + 1], lda, 
00278                             &a[n2 + a_dim1], lda, beta, &c__[1], n);
00279 
00280                 } else {
00281 
00282 /*                 N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T' */
00283 
00284                     ssyrk_("L", "T", &n1, k, alpha, &a[a_dim1 + 1], lda, beta, 
00285                              &c__[n2 + 1], n);
00286                     ssyrk_("U", "T", &n2, k, alpha, &a[n2 * a_dim1 + 1], lda, 
00287                             beta, &c__[n1 + 1], n);
00288                     sgemm_("T", "N", &n1, &n2, k, alpha, &a[a_dim1 + 1], lda, 
00289                             &a[n2 * a_dim1 + 1], lda, beta, &c__[1], n);
00290 
00291                 }
00292 
00293             }
00294 
00295         } else {
00296 
00297 /*           N is odd, and TRANSR = 'T' */
00298 
00299             if (lower) {
00300 
00301 /*              N is odd, TRANSR = 'T', and UPLO = 'L' */
00302 
00303                 if (notrans) {
00304 
00305 /*                 N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N' */
00306 
00307                     ssyrk_("U", "N", &n1, k, alpha, &a[a_dim1 + 1], lda, beta, 
00308                              &c__[1], &n1);
00309                     ssyrk_("L", "N", &n2, k, alpha, &a[n1 + 1 + a_dim1], lda, 
00310                             beta, &c__[2], &n1);
00311                     sgemm_("N", "T", &n1, &n2, k, alpha, &a[a_dim1 + 1], lda, 
00312                             &a[n1 + 1 + a_dim1], lda, beta, &c__[n1 * n1 + 1], 
00313                              &n1);
00314 
00315                 } else {
00316 
00317 /*                 N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T' */
00318 
00319                     ssyrk_("U", "T", &n1, k, alpha, &a[a_dim1 + 1], lda, beta, 
00320                              &c__[1], &n1);
00321                     ssyrk_("L", "T", &n2, k, alpha, &a[(n1 + 1) * a_dim1 + 1], 
00322                              lda, beta, &c__[2], &n1);
00323                     sgemm_("T", "N", &n1, &n2, k, alpha, &a[a_dim1 + 1], lda, 
00324                             &a[(n1 + 1) * a_dim1 + 1], lda, beta, &c__[n1 * 
00325                             n1 + 1], &n1);
00326 
00327                 }
00328 
00329             } else {
00330 
00331 /*              N is odd, TRANSR = 'T', and UPLO = 'U' */
00332 
00333                 if (notrans) {
00334 
00335 /*                 N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N' */
00336 
00337                     ssyrk_("U", "N", &n1, k, alpha, &a[a_dim1 + 1], lda, beta, 
00338                              &c__[n2 * n2 + 1], &n2);
00339                     ssyrk_("L", "N", &n2, k, alpha, &a[n1 + 1 + a_dim1], lda, 
00340                             beta, &c__[n1 * n2 + 1], &n2);
00341                     sgemm_("N", "T", &n2, &n1, k, alpha, &a[n1 + 1 + a_dim1], 
00342                             lda, &a[a_dim1 + 1], lda, beta, &c__[1], &n2);
00343 
00344                 } else {
00345 
00346 /*                 N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T' */
00347 
00348                     ssyrk_("U", "T", &n1, k, alpha, &a[a_dim1 + 1], lda, beta, 
00349                              &c__[n2 * n2 + 1], &n2);
00350                     ssyrk_("L", "T", &n2, k, alpha, &a[(n1 + 1) * a_dim1 + 1], 
00351                              lda, beta, &c__[n1 * n2 + 1], &n2);
00352                     sgemm_("T", "N", &n2, &n1, k, alpha, &a[(n1 + 1) * a_dim1 
00353                             + 1], lda, &a[a_dim1 + 1], lda, beta, &c__[1], &
00354                             n2);
00355 
00356                 }
00357 
00358             }
00359 
00360         }
00361 
00362     } else {
00363 
00364 /*        N is even */
00365 
00366         if (normaltransr) {
00367 
00368 /*           N is even and TRANSR = 'N' */
00369 
00370             if (lower) {
00371 
00372 /*              N is even, TRANSR = 'N', and UPLO = 'L' */
00373 
00374                 if (notrans) {
00375 
00376 /*                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' */
00377 
00378                     i__1 = *n + 1;
00379                     ssyrk_("L", "N", &nk, k, alpha, &a[a_dim1 + 1], lda, beta, 
00380                              &c__[2], &i__1);
00381                     i__1 = *n + 1;
00382                     ssyrk_("U", "N", &nk, k, alpha, &a[nk + 1 + a_dim1], lda, 
00383                             beta, &c__[1], &i__1);
00384                     i__1 = *n + 1;
00385                     sgemm_("N", "T", &nk, &nk, k, alpha, &a[nk + 1 + a_dim1], 
00386                             lda, &a[a_dim1 + 1], lda, beta, &c__[nk + 2], &
00387                             i__1);
00388 
00389                 } else {
00390 
00391 /*                 N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T' */
00392 
00393                     i__1 = *n + 1;
00394                     ssyrk_("L", "T", &nk, k, alpha, &a[a_dim1 + 1], lda, beta, 
00395                              &c__[2], &i__1);
00396                     i__1 = *n + 1;
00397                     ssyrk_("U", "T", &nk, k, alpha, &a[(nk + 1) * a_dim1 + 1], 
00398                              lda, beta, &c__[1], &i__1);
00399                     i__1 = *n + 1;
00400                     sgemm_("T", "N", &nk, &nk, k, alpha, &a[(nk + 1) * a_dim1 
00401                             + 1], lda, &a[a_dim1 + 1], lda, beta, &c__[nk + 2]
00402 , &i__1);
00403 
00404                 }
00405 
00406             } else {
00407 
00408 /*              N is even, TRANSR = 'N', and UPLO = 'U' */
00409 
00410                 if (notrans) {
00411 
00412 /*                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' */
00413 
00414                     i__1 = *n + 1;
00415                     ssyrk_("L", "N", &nk, k, alpha, &a[a_dim1 + 1], lda, beta, 
00416                              &c__[nk + 2], &i__1);
00417                     i__1 = *n + 1;
00418                     ssyrk_("U", "N", &nk, k, alpha, &a[nk + 1 + a_dim1], lda, 
00419                             beta, &c__[nk + 1], &i__1);
00420                     i__1 = *n + 1;
00421                     sgemm_("N", "T", &nk, &nk, k, alpha, &a[a_dim1 + 1], lda, 
00422                             &a[nk + 1 + a_dim1], lda, beta, &c__[1], &i__1);
00423 
00424                 } else {
00425 
00426 /*                 N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T' */
00427 
00428                     i__1 = *n + 1;
00429                     ssyrk_("L", "T", &nk, k, alpha, &a[a_dim1 + 1], lda, beta, 
00430                              &c__[nk + 2], &i__1);
00431                     i__1 = *n + 1;
00432                     ssyrk_("U", "T", &nk, k, alpha, &a[(nk + 1) * a_dim1 + 1], 
00433                              lda, beta, &c__[nk + 1], &i__1);
00434                     i__1 = *n + 1;
00435                     sgemm_("T", "N", &nk, &nk, k, alpha, &a[a_dim1 + 1], lda, 
00436                             &a[(nk + 1) * a_dim1 + 1], lda, beta, &c__[1], &
00437                             i__1);
00438 
00439                 }
00440 
00441             }
00442 
00443         } else {
00444 
00445 /*           N is even, and TRANSR = 'T' */
00446 
00447             if (lower) {
00448 
00449 /*              N is even, TRANSR = 'T', and UPLO = 'L' */
00450 
00451                 if (notrans) {
00452 
00453 /*                 N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N' */
00454 
00455                     ssyrk_("U", "N", &nk, k, alpha, &a[a_dim1 + 1], lda, beta, 
00456                              &c__[nk + 1], &nk);
00457                     ssyrk_("L", "N", &nk, k, alpha, &a[nk + 1 + a_dim1], lda, 
00458                             beta, &c__[1], &nk);
00459                     sgemm_("N", "T", &nk, &nk, k, alpha, &a[a_dim1 + 1], lda, 
00460                             &a[nk + 1 + a_dim1], lda, beta, &c__[(nk + 1) * 
00461                             nk + 1], &nk);
00462 
00463                 } else {
00464 
00465 /*                 N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T' */
00466 
00467                     ssyrk_("U", "T", &nk, k, alpha, &a[a_dim1 + 1], lda, beta, 
00468                              &c__[nk + 1], &nk);
00469                     ssyrk_("L", "T", &nk, k, alpha, &a[(nk + 1) * a_dim1 + 1], 
00470                              lda, beta, &c__[1], &nk);
00471                     sgemm_("T", "N", &nk, &nk, k, alpha, &a[a_dim1 + 1], lda, 
00472                             &a[(nk + 1) * a_dim1 + 1], lda, beta, &c__[(nk + 
00473                             1) * nk + 1], &nk);
00474 
00475                 }
00476 
00477             } else {
00478 
00479 /*              N is even, TRANSR = 'T', and UPLO = 'U' */
00480 
00481                 if (notrans) {
00482 
00483 /*                 N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N' */
00484 
00485                     ssyrk_("U", "N", &nk, k, alpha, &a[a_dim1 + 1], lda, beta, 
00486                              &c__[nk * (nk + 1) + 1], &nk);
00487                     ssyrk_("L", "N", &nk, k, alpha, &a[nk + 1 + a_dim1], lda, 
00488                             beta, &c__[nk * nk + 1], &nk);
00489                     sgemm_("N", "T", &nk, &nk, k, alpha, &a[nk + 1 + a_dim1], 
00490                             lda, &a[a_dim1 + 1], lda, beta, &c__[1], &nk);
00491 
00492                 } else {
00493 
00494 /*                 N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T' */
00495 
00496                     ssyrk_("U", "T", &nk, k, alpha, &a[a_dim1 + 1], lda, beta, 
00497                              &c__[nk * (nk + 1) + 1], &nk);
00498                     ssyrk_("L", "T", &nk, k, alpha, &a[(nk + 1) * a_dim1 + 1], 
00499                              lda, beta, &c__[nk * nk + 1], &nk);
00500                     sgemm_("T", "N", &nk, &nk, k, alpha, &a[(nk + 1) * a_dim1 
00501                             + 1], lda, &a[a_dim1 + 1], lda, beta, &c__[1], &
00502                             nk);
00503 
00504                 }
00505 
00506             }
00507 
00508         }
00509 
00510     }
00511 
00512     return 0;
00513 
00514 /*     End of SSFRK */
00515 
00516 } /* ssfrk_ */