chpr.c
Go to the documentation of this file.
00001 /* chpr.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 chpr_(char *uplo, integer *n, real *alpha, complex *x, 
00017         integer *incx, complex *ap)
00018 {
00019     /* System generated locals */
00020     integer i__1, i__2, i__3, i__4, i__5;
00021     real r__1;
00022     complex q__1, q__2;
00023 
00024     /* Builtin functions */
00025     void r_cnjg(complex *, complex *);
00026 
00027     /* Local variables */
00028     integer i__, j, k, kk, ix, jx, kx, info;
00029     complex temp;
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 /*  CHPR    performs the hermitian rank 1 operation */
00042 
00043 /*     A := alpha*x*conjg( x' ) + A, */
00044 
00045 /*  where alpha is a real scalar, x is an n element vector 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  - REAL            . */
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 /*  AP     - COMPLEX          array of DIMENSION at least */
00085 /*           ( ( n*( n + 1 ) )/2 ). */
00086 /*           Before entry with  UPLO = 'U' or 'u', the array AP must */
00087 /*           contain the upper triangular part of the hermitian matrix */
00088 /*           packed sequentially, column by column, so that AP( 1 ) */
00089 /*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
00090 /*           and a( 2, 2 ) respectively, and so on. On exit, the array */
00091 /*           AP is overwritten by the upper triangular part of the */
00092 /*           updated matrix. */
00093 /*           Before entry with UPLO = 'L' or 'l', the array AP must */
00094 /*           contain the lower triangular part of the hermitian matrix */
00095 /*           packed sequentially, column by column, so that AP( 1 ) */
00096 /*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
00097 /*           and a( 3, 1 ) respectively, and so on. On exit, the array */
00098 /*           AP is overwritten by the lower triangular part of the */
00099 /*           updated matrix. */
00100 /*           Note that the imaginary parts of the diagonal elements need */
00101 /*           not be set, they are assumed to be zero, and on exit they */
00102 /*           are set to zero. */
00103 
00104 
00105 /*  Level 2 Blas routine. */
00106 
00107 /*  -- Written on 22-October-1986. */
00108 /*     Jack Dongarra, Argonne National Lab. */
00109 /*     Jeremy Du Croz, Nag Central Office. */
00110 /*     Sven Hammarling, Nag Central Office. */
00111 /*     Richard Hanson, Sandia National Labs. */
00112 
00113 
00114 /*     .. Parameters .. */
00115 /*     .. */
00116 /*     .. Local Scalars .. */
00117 /*     .. */
00118 /*     .. External Functions .. */
00119 /*     .. */
00120 /*     .. External Subroutines .. */
00121 /*     .. */
00122 /*     .. Intrinsic Functions .. */
00123 /*     .. */
00124 
00125 /*     Test the input parameters. */
00126 
00127     /* Parameter adjustments */
00128     --ap;
00129     --x;
00130 
00131     /* Function Body */
00132     info = 0;
00133     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00134         info = 1;
00135     } else if (*n < 0) {
00136         info = 2;
00137     } else if (*incx == 0) {
00138         info = 5;
00139     }
00140     if (info != 0) {
00141         xerbla_("CHPR  ", &info);
00142         return 0;
00143     }
00144 
00145 /*     Quick return if possible. */
00146 
00147     if (*n == 0 || *alpha == 0.f) {
00148         return 0;
00149     }
00150 
00151 /*     Set the start point in X if the increment is not unity. */
00152 
00153     if (*incx <= 0) {
00154         kx = 1 - (*n - 1) * *incx;
00155     } else if (*incx != 1) {
00156         kx = 1;
00157     }
00158 
00159 /*     Start the operations. In this version the elements of the array AP */
00160 /*     are accessed sequentially with one pass through AP. */
00161 
00162     kk = 1;
00163     if (lsame_(uplo, "U")) {
00164 
00165 /*        Form  A  when upper triangle is stored in AP. */
00166 
00167         if (*incx == 1) {
00168             i__1 = *n;
00169             for (j = 1; j <= i__1; ++j) {
00170                 i__2 = j;
00171                 if (x[i__2].r != 0.f || x[i__2].i != 0.f) {
00172                     r_cnjg(&q__2, &x[j]);
00173                     q__1.r = *alpha * q__2.r, q__1.i = *alpha * q__2.i;
00174                     temp.r = q__1.r, temp.i = q__1.i;
00175                     k = kk;
00176                     i__2 = j - 1;
00177                     for (i__ = 1; i__ <= i__2; ++i__) {
00178                         i__3 = k;
00179                         i__4 = k;
00180                         i__5 = i__;
00181                         q__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
00182                                 q__2.i = x[i__5].r * temp.i + x[i__5].i * 
00183                                 temp.r;
00184                         q__1.r = ap[i__4].r + q__2.r, q__1.i = ap[i__4].i + 
00185                                 q__2.i;
00186                         ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
00187                         ++k;
00188 /* L10: */
00189                     }
00190                     i__2 = kk + j - 1;
00191                     i__3 = kk + j - 1;
00192                     i__4 = j;
00193                     q__1.r = x[i__4].r * temp.r - x[i__4].i * temp.i, q__1.i =
00194                              x[i__4].r * temp.i + x[i__4].i * temp.r;
00195                     r__1 = ap[i__3].r + q__1.r;
00196                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00197                 } else {
00198                     i__2 = kk + j - 1;
00199                     i__3 = kk + j - 1;
00200                     r__1 = ap[i__3].r;
00201                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00202                 }
00203                 kk += j;
00204 /* L20: */
00205             }
00206         } else {
00207             jx = kx;
00208             i__1 = *n;
00209             for (j = 1; j <= i__1; ++j) {
00210                 i__2 = jx;
00211                 if (x[i__2].r != 0.f || x[i__2].i != 0.f) {
00212                     r_cnjg(&q__2, &x[jx]);
00213                     q__1.r = *alpha * q__2.r, q__1.i = *alpha * q__2.i;
00214                     temp.r = q__1.r, temp.i = q__1.i;
00215                     ix = kx;
00216                     i__2 = kk + j - 2;
00217                     for (k = kk; k <= i__2; ++k) {
00218                         i__3 = k;
00219                         i__4 = k;
00220                         i__5 = ix;
00221                         q__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
00222                                 q__2.i = x[i__5].r * temp.i + x[i__5].i * 
00223                                 temp.r;
00224                         q__1.r = ap[i__4].r + q__2.r, q__1.i = ap[i__4].i + 
00225                                 q__2.i;
00226                         ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
00227                         ix += *incx;
00228 /* L30: */
00229                     }
00230                     i__2 = kk + j - 1;
00231                     i__3 = kk + j - 1;
00232                     i__4 = jx;
00233                     q__1.r = x[i__4].r * temp.r - x[i__4].i * temp.i, q__1.i =
00234                              x[i__4].r * temp.i + x[i__4].i * temp.r;
00235                     r__1 = ap[i__3].r + q__1.r;
00236                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00237                 } else {
00238                     i__2 = kk + j - 1;
00239                     i__3 = kk + j - 1;
00240                     r__1 = ap[i__3].r;
00241                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00242                 }
00243                 jx += *incx;
00244                 kk += j;
00245 /* L40: */
00246             }
00247         }
00248     } else {
00249 
00250 /*        Form  A  when lower triangle is stored in AP. */
00251 
00252         if (*incx == 1) {
00253             i__1 = *n;
00254             for (j = 1; j <= i__1; ++j) {
00255                 i__2 = j;
00256                 if (x[i__2].r != 0.f || x[i__2].i != 0.f) {
00257                     r_cnjg(&q__2, &x[j]);
00258                     q__1.r = *alpha * q__2.r, q__1.i = *alpha * q__2.i;
00259                     temp.r = q__1.r, temp.i = q__1.i;
00260                     i__2 = kk;
00261                     i__3 = kk;
00262                     i__4 = j;
00263                     q__1.r = temp.r * x[i__4].r - temp.i * x[i__4].i, q__1.i =
00264                              temp.r * x[i__4].i + temp.i * x[i__4].r;
00265                     r__1 = ap[i__3].r + q__1.r;
00266                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00267                     k = kk + 1;
00268                     i__2 = *n;
00269                     for (i__ = j + 1; i__ <= i__2; ++i__) {
00270                         i__3 = k;
00271                         i__4 = k;
00272                         i__5 = i__;
00273                         q__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
00274                                 q__2.i = x[i__5].r * temp.i + x[i__5].i * 
00275                                 temp.r;
00276                         q__1.r = ap[i__4].r + q__2.r, q__1.i = ap[i__4].i + 
00277                                 q__2.i;
00278                         ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
00279                         ++k;
00280 /* L50: */
00281                     }
00282                 } else {
00283                     i__2 = kk;
00284                     i__3 = kk;
00285                     r__1 = ap[i__3].r;
00286                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00287                 }
00288                 kk = kk + *n - j + 1;
00289 /* L60: */
00290             }
00291         } else {
00292             jx = kx;
00293             i__1 = *n;
00294             for (j = 1; j <= i__1; ++j) {
00295                 i__2 = jx;
00296                 if (x[i__2].r != 0.f || x[i__2].i != 0.f) {
00297                     r_cnjg(&q__2, &x[jx]);
00298                     q__1.r = *alpha * q__2.r, q__1.i = *alpha * q__2.i;
00299                     temp.r = q__1.r, temp.i = q__1.i;
00300                     i__2 = kk;
00301                     i__3 = kk;
00302                     i__4 = jx;
00303                     q__1.r = temp.r * x[i__4].r - temp.i * x[i__4].i, q__1.i =
00304                              temp.r * x[i__4].i + temp.i * x[i__4].r;
00305                     r__1 = ap[i__3].r + q__1.r;
00306                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00307                     ix = jx;
00308                     i__2 = kk + *n - j;
00309                     for (k = kk + 1; k <= i__2; ++k) {
00310                         ix += *incx;
00311                         i__3 = k;
00312                         i__4 = k;
00313                         i__5 = ix;
00314                         q__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
00315                                 q__2.i = x[i__5].r * temp.i + x[i__5].i * 
00316                                 temp.r;
00317                         q__1.r = ap[i__4].r + q__2.r, q__1.i = ap[i__4].i + 
00318                                 q__2.i;
00319                         ap[i__3].r = q__1.r, ap[i__3].i = q__1.i;
00320 /* L70: */
00321                     }
00322                 } else {
00323                     i__2 = kk;
00324                     i__3 = kk;
00325                     r__1 = ap[i__3].r;
00326                     ap[i__2].r = r__1, ap[i__2].i = 0.f;
00327                 }
00328                 jx += *incx;
00329                 kk = kk + *n - j + 1;
00330 /* L80: */
00331             }
00332         }
00333     }
00334 
00335     return 0;
00336 
00337 /*     End of CHPR  . */
00338 
00339 } /* chpr_ */


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