swiftnav: zhfrk.c Source File

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