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


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:28