swiftnav: dsfrk.c Source File

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