chpr2.c
Go to the documentation of this file.
00001 /* chpr2.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 chpr2_(char *uplo, integer *n, complex *alpha, complex *
00017         x, integer *incx, complex *y, integer *incy, complex *ap)
00018 {
00019     /* System generated locals */
00020     integer i__1, i__2, i__3, i__4, i__5, i__6;
00021     real r__1;
00022     complex q__1, q__2, q__3, q__4;
00023 
00024     /* Builtin functions */
00025     void r_cnjg(complex *, complex *);
00026 
00027     /* Local variables */
00028     integer i__, j, k, kk, ix, iy, jx, jy, kx, ky, info;
00029     complex temp1, temp2;
00030     extern logical lsame_(char *, char *);
00031     extern /* Subroutine */ int xerbla_(char *, integer *);
00032 
00033 /*     .. Scalar Arguments .. */
00034 /*     .. */
00035 /*     .. Array Arguments .. */
00036 /*     .. */
00037 
00038 /*  Purpose */
00039 /*  ======= */
00040 
00041 /*  CHPR2  performs the hermitian rank 2 operation */
00042 
00043 /*     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, */
00044 
00045 /*  where alpha is a scalar, x and y are n element vectors and A is an */
00046 /*  n by n hermitian matrix, supplied in packed form. */
00047 
00048 /*  Arguments */
00049 /*  ========== */
00050 
00051 /*  UPLO   - CHARACTER*1. */
00052 /*           On entry, UPLO specifies whether the upper or lower */
00053 /*           triangular part of the matrix A is supplied in the packed */
00054 /*           array AP as follows: */
00055 
00056 /*              UPLO = 'U' or 'u'   The upper triangular part of A is */
00057 /*                                  supplied in AP. */
00058 
00059 /*              UPLO = 'L' or 'l'   The lower triangular part of A is */
00060 /*                                  supplied in AP. */
00061 
00062 /*           Unchanged on exit. */
00063 
00064 /*  N      - INTEGER. */
00065 /*           On entry, N specifies the order of the matrix A. */
00066 /*           N must be at least zero. */
00067 /*           Unchanged on exit. */
00068 
00069 /*  ALPHA  - COMPLEX         . */
00070 /*           On entry, ALPHA specifies the scalar alpha. */
00071 /*           Unchanged on exit. */
00072 
00073 /*  X      - COMPLEX          array of dimension at least */
00074 /*           ( 1 + ( n - 1 )*abs( INCX ) ). */
00075 /*           Before entry, the incremented array X must contain the n */
00076 /*           element vector x. */
00077 /*           Unchanged on exit. */
00078 
00079 /*  INCX   - INTEGER. */
00080 /*           On entry, INCX specifies the increment for the elements of */
00081 /*           X. INCX must not be zero. */
00082 /*           Unchanged on exit. */
00083 
00084 /*  Y      - COMPLEX          array of dimension at least */
00085 /*           ( 1 + ( n - 1 )*abs( INCY ) ). */
00086 /*           Before entry, the incremented array Y must contain the n */
00087 /*           element vector y. */
00088 /*           Unchanged on exit. */
00089 
00090 /*  INCY   - INTEGER. */
00091 /*           On entry, INCY specifies the increment for the elements of */
00092 /*           Y. INCY must not be zero. */
00093 /*           Unchanged on exit. */
00094 
00095 /*  AP     - COMPLEX          array of DIMENSION at least */
00096 /*           ( ( n*( n + 1 ) )/2 ). */
00097 /*           Before entry with  UPLO = 'U' or 'u', the array AP must */
00098 /*           contain the upper triangular part of the hermitian matrix */
00099 /*           packed sequentially, column by column, so that AP( 1 ) */
00100 /*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
00101 /*           and a( 2, 2 ) respectively, and so on. On exit, the array */
00102 /*           AP is overwritten by the upper triangular part of the */
00103 /*           updated matrix. */
00104 /*           Before entry with UPLO = 'L' or 'l', the array AP must */
00105 /*           contain the lower triangular part of the hermitian matrix */
00106 /*           packed sequentially, column by column, so that AP( 1 ) */
00107 /*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
00108 /*           and a( 3, 1 ) respectively, and so on. On exit, the array */
00109 /*           AP is overwritten by the lower triangular part of the */
00110 /*           updated matrix. */
00111 /*           Note that the imaginary parts of the diagonal elements need */
00112 /*           not be set, they are assumed to be zero, and on exit they */
00113 /*           are set to zero. */
00114 
00115 
00116 /*  Level 2 Blas routine. */
00117 
00118 /*  -- Written on 22-October-1986. */
00119 /*     Jack Dongarra, Argonne National Lab. */
00120 /*     Jeremy Du Croz, Nag Central Office. */
00121 /*     Sven Hammarling, Nag Central Office. */
00122 /*     Richard Hanson, Sandia National Labs. */
00123 
00124 
00125 /*     .. Parameters .. */
00126 /*     .. */
00127 /*     .. Local Scalars .. */
00128 /*     .. */
00129 /*     .. External Functions .. */
00130 /*     .. */
00131 /*     .. External Subroutines .. */
00132 /*     .. */
00133 /*     .. Intrinsic Functions .. */
00134 /*     .. */
00135 
00136 /*     Test the input parameters. */
00137 
00138     /* Parameter adjustments */
00139     --ap;
00140     --y;
00141     --x;
00142 
00143     /* Function Body */
00144     info = 0;
00145     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00146         info = 1;
00147     } else if (*n < 0) {
00148         info = 2;
00149     } else if (*incx == 0) {
00150         info = 5;
00151     } else if (*incy == 0) {
00152         info = 7;
00153     }
00154     if (info != 0) {
00155         xerbla_("CHPR2 ", &info);
00156         return 0;
00157     }
00158 
00159 /*     Quick return if possible. */
00160 
00161     if (*n == 0 || alpha->r == 0.f && alpha->i == 0.f) {
00162         return 0;
00163     }
00164 
00165 /*     Set up the start points in X and Y if the increments are not both */
00166 /*     unity. */
00167 
00168     if (*incx != 1 || *incy != 1) {
00169         if (*incx > 0) {
00170             kx = 1;
00171         } else {
00172             kx = 1 - (*n - 1) * *incx;
00173         }
00174         if (*incy > 0) {
00175             ky = 1;
00176         } else {
00177             ky = 1 - (*n - 1) * *incy;
00178         }
00179         jx = kx;
00180         jy = ky;
00181     }
00182 
00183 /*     Start the operations. In this version the elements of the array AP */
00184 /*     are accessed sequentially with one pass through AP. */
00185 
00186     kk = 1;
00187     if (lsame_(uplo, "U")) {
00188 
00189 /*        Form  A  when upper triangle is stored in AP. */
00190 
00191         if (*incx == 1 && *incy == 1) {
00192             i__1 = *n;
00193             for (j = 1; j <= i__1; ++j) {
00194                 i__2 = j;
00195                 i__3 = j;
00196                 if (x[i__2].r != 0.f || x[i__2].i != 0.f || (y[i__3].r != 0.f 
00197                         || y[i__3].i != 0.f)) {
00198                     r_cnjg(&q__2, &y[j]);
00199                     q__1.r = alpha->r * q__2.r - alpha->i * q__2.i, q__1.i = 
00200                             alpha->r * q__2.i + alpha->i * q__2.r;
00201                     temp1.r = q__1.r, temp1.i = q__1.i;
00202                     i__2 = j;
00203                     q__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
00204                             q__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
00205                             .r;
00206                     r_cnjg(&q__1, &q__2);
00207                     temp2.r = q__1.r, temp2.i = q__1.i;
00208                     k = kk;
00209                     i__2 = j - 1;
00210                     for (i__ = 1; i__ <= i__2; ++i__) {
00211                         i__3 = k;
00212                         i__4 = k;
00213                         i__5 = i__;
00214                         q__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
00215                                 q__3.i = x[i__5].r * temp1.i + x[i__5].i * 
00216                                 temp1.r;
00217                         q__2.r = ap[i__4].r + q__3.r, q__2.i = ap[i__4].i + 
00218                                 q__3.i;
00219                         i__6 = i__;
00220                         q__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
00221                                 q__4.i = y[i__6].r * temp2.i + y[i__6].i * 
00222                                 temp2.r;
00223                         q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
00224                         ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
00225                         ++k;
00226 /* L10: */
00227                     }
00228                     i__2 = kk + j - 1;
00229                     i__3 = kk + j - 1;
00230                     i__4 = j;
00231                     q__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
00232                             q__2.i = x[i__4].r * temp1.i + x[i__4].i * 
00233                             temp1.r;
00234                     i__5 = j;
00235                     q__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
00236                             q__3.i = y[i__5].r * temp2.i + y[i__5].i * 
00237                             temp2.r;
00238                     q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00239                     r__1 = ap[i__3].r + q__1.r;
00240                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00241                 } else {
00242                     i__2 = kk + j - 1;
00243                     i__3 = kk + j - 1;
00244                     r__1 = ap[i__3].r;
00245                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00246                 }
00247                 kk += j;
00248 /* L20: */
00249             }
00250         } else {
00251             i__1 = *n;
00252             for (j = 1; j <= i__1; ++j) {
00253                 i__2 = jx;
00254                 i__3 = jy;
00255                 if (x[i__2].r != 0.f || x[i__2].i != 0.f || (y[i__3].r != 0.f 
00256                         || y[i__3].i != 0.f)) {
00257                     r_cnjg(&q__2, &y[jy]);
00258                     q__1.r = alpha->r * q__2.r - alpha->i * q__2.i, q__1.i = 
00259                             alpha->r * q__2.i + alpha->i * q__2.r;
00260                     temp1.r = q__1.r, temp1.i = q__1.i;
00261                     i__2 = jx;
00262                     q__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
00263                             q__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
00264                             .r;
00265                     r_cnjg(&q__1, &q__2);
00266                     temp2.r = q__1.r, temp2.i = q__1.i;
00267                     ix = kx;
00268                     iy = ky;
00269                     i__2 = kk + j - 2;
00270                     for (k = kk; k <= i__2; ++k) {
00271                         i__3 = k;
00272                         i__4 = k;
00273                         i__5 = ix;
00274                         q__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
00275                                 q__3.i = x[i__5].r * temp1.i + x[i__5].i * 
00276                                 temp1.r;
00277                         q__2.r = ap[i__4].r + q__3.r, q__2.i = ap[i__4].i + 
00278                                 q__3.i;
00279                         i__6 = iy;
00280                         q__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
00281                                 q__4.i = y[i__6].r * temp2.i + y[i__6].i * 
00282                                 temp2.r;
00283                         q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
00284                         ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
00285                         ix += *incx;
00286                         iy += *incy;
00287 /* L30: */
00288                     }
00289                     i__2 = kk + j - 1;
00290                     i__3 = kk + j - 1;
00291                     i__4 = jx;
00292                     q__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
00293                             q__2.i = x[i__4].r * temp1.i + x[i__4].i * 
00294                             temp1.r;
00295                     i__5 = jy;
00296                     q__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
00297                             q__3.i = y[i__5].r * temp2.i + y[i__5].i * 
00298                             temp2.r;
00299                     q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00300                     r__1 = ap[i__3].r + q__1.r;
00301                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00302                 } else {
00303                     i__2 = kk + j - 1;
00304                     i__3 = kk + j - 1;
00305                     r__1 = ap[i__3].r;
00306                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00307                 }
00308                 jx += *incx;
00309                 jy += *incy;
00310                 kk += j;
00311 /* L40: */
00312             }
00313         }
00314     } else {
00315 
00316 /*        Form  A  when lower triangle is stored in AP. */
00317 
00318         if (*incx == 1 && *incy == 1) {
00319             i__1 = *n;
00320             for (j = 1; j <= i__1; ++j) {
00321                 i__2 = j;
00322                 i__3 = j;
00323                 if (x[i__2].r != 0.f || x[i__2].i != 0.f || (y[i__3].r != 0.f 
00324                         || y[i__3].i != 0.f)) {
00325                     r_cnjg(&q__2, &y[j]);
00326                     q__1.r = alpha->r * q__2.r - alpha->i * q__2.i, q__1.i = 
00327                             alpha->r * q__2.i + alpha->i * q__2.r;
00328                     temp1.r = q__1.r, temp1.i = q__1.i;
00329                     i__2 = j;
00330                     q__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
00331                             q__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
00332                             .r;
00333                     r_cnjg(&q__1, &q__2);
00334                     temp2.r = q__1.r, temp2.i = q__1.i;
00335                     i__2 = kk;
00336                     i__3 = kk;
00337                     i__4 = j;
00338                     q__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
00339                             q__2.i = x[i__4].r * temp1.i + x[i__4].i * 
00340                             temp1.r;
00341                     i__5 = j;
00342                     q__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
00343                             q__3.i = y[i__5].r * temp2.i + y[i__5].i * 
00344                             temp2.r;
00345                     q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00346                     r__1 = ap[i__3].r + q__1.r;
00347                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00348                     k = kk + 1;
00349                     i__2 = *n;
00350                     for (i__ = j + 1; i__ <= i__2; ++i__) {
00351                         i__3 = k;
00352                         i__4 = k;
00353                         i__5 = i__;
00354                         q__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
00355                                 q__3.i = x[i__5].r * temp1.i + x[i__5].i * 
00356                                 temp1.r;
00357                         q__2.r = ap[i__4].r + q__3.r, q__2.i = ap[i__4].i + 
00358                                 q__3.i;
00359                         i__6 = i__;
00360                         q__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
00361                                 q__4.i = y[i__6].r * temp2.i + y[i__6].i * 
00362                                 temp2.r;
00363                         q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
00364                         ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
00365                         ++k;
00366 /* L50: */
00367                     }
00368                 } else {
00369                     i__2 = kk;
00370                     i__3 = kk;
00371                     r__1 = ap[i__3].r;
00372                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00373                 }
00374                 kk = kk + *n - j + 1;
00375 /* L60: */
00376             }
00377         } else {
00378             i__1 = *n;
00379             for (j = 1; j <= i__1; ++j) {
00380                 i__2 = jx;
00381                 i__3 = jy;
00382                 if (x[i__2].r != 0.f || x[i__2].i != 0.f || (y[i__3].r != 0.f 
00383                         || y[i__3].i != 0.f)) {
00384                     r_cnjg(&q__2, &y[jy]);
00385                     q__1.r = alpha->r * q__2.r - alpha->i * q__2.i, q__1.i = 
00386                             alpha->r * q__2.i + alpha->i * q__2.r;
00387                     temp1.r = q__1.r, temp1.i = q__1.i;
00388                     i__2 = jx;
00389                     q__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
00390                             q__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
00391                             .r;
00392                     r_cnjg(&q__1, &q__2);
00393                     temp2.r = q__1.r, temp2.i = q__1.i;
00394                     i__2 = kk;
00395                     i__3 = kk;
00396                     i__4 = jx;
00397                     q__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i, 
00398                             q__2.i = x[i__4].r * temp1.i + x[i__4].i * 
00399                             temp1.r;
00400                     i__5 = jy;
00401                     q__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i, 
00402                             q__3.i = y[i__5].r * temp2.i + y[i__5].i * 
00403                             temp2.r;
00404                     q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00405                     r__1 = ap[i__3].r + q__1.r;
00406                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00407                     ix = jx;
00408                     iy = jy;
00409                     i__2 = kk + *n - j;
00410                     for (k = kk + 1; k <= i__2; ++k) {
00411                         ix += *incx;
00412                         iy += *incy;
00413                         i__3 = k;
00414                         i__4 = k;
00415                         i__5 = ix;
00416                         q__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i, 
00417                                 q__3.i = x[i__5].r * temp1.i + x[i__5].i * 
00418                                 temp1.r;
00419                         q__2.r = ap[i__4].r + q__3.r, q__2.i = ap[i__4].i + 
00420                                 q__3.i;
00421                         i__6 = iy;
00422                         q__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i, 
00423                                 q__4.i = y[i__6].r * temp2.i + y[i__6].i * 
00424                                 temp2.r;
00425                         q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
00426                         ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
00427 /* L70: */
00428                     }
00429                 } else {
00430                     i__2 = kk;
00431                     i__3 = kk;
00432                     r__1 = ap[i__3].r;
00433                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00434                 }
00435                 jx += *incx;
00436                 jy += *incy;
00437                 kk = kk + *n - j + 1;
00438 /* L80: */
00439             }
00440         }
00441     }
00442 
00443     return 0;
00444 
00445 /*     End of CHPR2 . */
00446 
00447 } /* chpr2_ */

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