swiftnav: zspr.c Source File

Go to the documentation of this file.
00001 /* zspr.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 zspr_(char *uplo, integer *n, doublecomplex *alpha, 
00017         doublecomplex *x, integer *incx, doublecomplex *ap)
00018 {
00019     /* System generated locals */
00020     integer i__1, i__2, i__3, i__4, i__5;
00021     doublecomplex z__1, z__2;
00022 
00023     /* Local variables */
00024     integer i__, j, k, kk, ix, jx, kx, info;
00025     doublecomplex temp;
00026     extern logical lsame_(char *, char *);
00027     extern /* Subroutine */ int xerbla_(char *, integer *);
00028 
00029 
00030 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00031 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00032 /*     November 2006 */
00033 
00034 /*     .. Scalar Arguments .. */
00035 /*     .. */
00036 /*     .. Array Arguments .. */
00037 /*     .. */
00038 
00039 /*  Purpose */
00040 /*  ======= */
00041 
00042 /*  ZSPR    performs the symmetric rank 1 operation */
00043 
00044 /*     A := alpha*x*conjg( x' ) + A, */
00045 
00046 /*  where alpha is a complex scalar, x is an n element vector and A is an */
00047 /*  n by n symmetric matrix, supplied in packed form. */
00048 
00049 /*  Arguments */
00050 /*  ========== */
00051 
00052 /*  UPLO     (input) CHARACTER*1 */
00053 /*           On entry, UPLO specifies whether the upper or lower */
00054 /*           triangular part of the matrix A is supplied in the packed */
00055 /*           array AP as follows: */
00056 
00057 /*              UPLO = 'U' or 'u'   The upper triangular part of A is */
00058 /*                                  supplied in AP. */
00059 
00060 /*              UPLO = 'L' or 'l'   The lower triangular part of A is */
00061 /*                                  supplied in AP. */
00062 
00063 /*           Unchanged on exit. */
00064 
00065 /*  N        (input) INTEGER */
00066 /*           On entry, N specifies the order of the matrix A. */
00067 /*           N must be at least zero. */
00068 /*           Unchanged on exit. */
00069 
00070 /*  ALPHA    (input) COMPLEX*16 */
00071 /*           On entry, ALPHA specifies the scalar alpha. */
00072 /*           Unchanged on exit. */
00073 
00074 /*  X        (input) COMPLEX*16 array, dimension at least */
00075 /*           ( 1 + ( N - 1 )*abs( INCX ) ). */
00076 /*           Before entry, the incremented array X must contain the N- */
00077 /*           element vector x. */
00078 /*           Unchanged on exit. */
00079 
00080 /*  INCX     (input) INTEGER */
00081 /*           On entry, INCX specifies the increment for the elements of */
00082 /*           X. INCX must not be zero. */
00083 /*           Unchanged on exit. */
00084 
00085 /*  AP       (input/output) COMPLEX*16 array, dimension at least */
00086 /*           ( ( N*( N + 1 ) )/2 ). */
00087 /*           Before entry, with  UPLO = 'U' or 'u', the array AP must */
00088 /*           contain the upper triangular part of the symmetric matrix */
00089 /*           packed sequentially, column by column, so that AP( 1 ) */
00090 /*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
00091 /*           and a( 2, 2 ) respectively, and so on. On exit, the array */
00092 /*           AP is overwritten by the upper triangular part of the */
00093 /*           updated matrix. */
00094 /*           Before entry, with UPLO = 'L' or 'l', the array AP must */
00095 /*           contain the lower triangular part of the symmetric matrix */
00096 /*           packed sequentially, column by column, so that AP( 1 ) */
00097 /*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
00098 /*           and a( 3, 1 ) respectively, and so on. On exit, the array */
00099 /*           AP is overwritten by the lower triangular part of the */
00100 /*           updated matrix. */
00101 /*           Note that the imaginary parts of the diagonal elements need */
00102 /*           not be set, they are assumed to be zero, and on exit they */
00103 /*           are set to zero. */
00104 
00105 /* ===================================================================== */
00106 
00107 /*     .. Parameters .. */
00108 /*     .. */
00109 /*     .. Local Scalars .. */
00110 /*     .. */
00111 /*     .. External Functions .. */
00112 /*     .. */
00113 /*     .. External Subroutines .. */
00114 /*     .. */
00115 /*     .. Executable Statements .. */
00116 
00117 /*     Test the input parameters. */
00118 
00119     /* Parameter adjustments */
00120     --ap;
00121     --x;
00122 
00123     /* Function Body */
00124     info = 0;
00125     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00126         info = 1;
00127     } else if (*n < 0) {
00128         info = 2;
00129     } else if (*incx == 0) {
00130         info = 5;
00131     }
00132     if (info != 0) {
00133         xerbla_("ZSPR  ", &info);
00134         return 0;
00135     }
00136 
00137 /*     Quick return if possible. */
00138 
00139     if (*n == 0 || alpha->r == 0. && alpha->i == 0.) {
00140         return 0;
00141     }
00142 
00143 /*     Set the start point in X if the increment is not unity. */
00144 
00145     if (*incx <= 0) {
00146         kx = 1 - (*n - 1) * *incx;
00147     } else if (*incx != 1) {
00148         kx = 1;
00149     }
00150 
00151 /*     Start the operations. In this version the elements of the array AP */
00152 /*     are accessed sequentially with one pass through AP. */
00153 
00154     kk = 1;
00155     if (lsame_(uplo, "U")) {
00156 
00157 /*        Form  A  when upper triangle is stored in AP. */
00158 
00159         if (*incx == 1) {
00160             i__1 = *n;
00161             for (j = 1; j <= i__1; ++j) {
00162                 i__2 = j;
00163                 if (x[i__2].r != 0. || x[i__2].i != 0.) {
00164                     i__2 = j;
00165                     z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
00166                             z__1.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
00167                             .r;
00168                     temp.r = z__1.r, temp.i = z__1.i;
00169                     k = kk;
00170                     i__2 = j - 1;
00171                     for (i__ = 1; i__ <= i__2; ++i__) {
00172                         i__3 = k;
00173                         i__4 = k;
00174                         i__5 = i__;
00175                         z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
00176                                 z__2.i = x[i__5].r * temp.i + x[i__5].i * 
00177                                 temp.r;
00178                         z__1.r = ap[i__4].r + z__2.r, z__1.i = ap[i__4].i + 
00179                                 z__2.i;
00180                         ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
00181                         ++k;
00182 /* L10: */
00183                     }
00184                     i__2 = kk + j - 1;
00185                     i__3 = kk + j - 1;
00186                     i__4 = j;
00187                     z__2.r = x[i__4].r * temp.r - x[i__4].i * temp.i, z__2.i =
00188                              x[i__4].r * temp.i + x[i__4].i * temp.r;
00189                     z__1.r = ap[i__3].r + z__2.r, z__1.i = ap[i__3].i + 
00190                             z__2.i;
00191                     ap[i__2].r = z__1.r, ap[i__2].i = z__1.i;
00192                 } else {
00193                     i__2 = kk + j - 1;
00194                     i__3 = kk + j - 1;
00195                     ap[i__2].r = ap[i__3].r, ap[i__2].i = ap[i__3].i;
00196                 }
00197                 kk += j;
00198 /* L20: */
00199             }
00200         } else {
00201             jx = kx;
00202             i__1 = *n;
00203             for (j = 1; j <= i__1; ++j) {
00204                 i__2 = jx;
00205                 if (x[i__2].r != 0. || x[i__2].i != 0.) {
00206                     i__2 = jx;
00207                     z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
00208                             z__1.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
00209                             .r;
00210                     temp.r = z__1.r, temp.i = z__1.i;
00211                     ix = kx;
00212                     i__2 = kk + j - 2;
00213                     for (k = kk; k <= i__2; ++k) {
00214                         i__3 = k;
00215                         i__4 = k;
00216                         i__5 = ix;
00217                         z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
00218                                 z__2.i = x[i__5].r * temp.i + x[i__5].i * 
00219                                 temp.r;
00220                         z__1.r = ap[i__4].r + z__2.r, z__1.i = ap[i__4].i + 
00221                                 z__2.i;
00222                         ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
00223                         ix += *incx;
00224 /* L30: */
00225                     }
00226                     i__2 = kk + j - 1;
00227                     i__3 = kk + j - 1;
00228                     i__4 = jx;
00229                     z__2.r = x[i__4].r * temp.r - x[i__4].i * temp.i, z__2.i =
00230                              x[i__4].r * temp.i + x[i__4].i * temp.r;
00231                     z__1.r = ap[i__3].r + z__2.r, z__1.i = ap[i__3].i + 
00232                             z__2.i;
00233                     ap[i__2].r = z__1.r, ap[i__2].i = z__1.i;
00234                 } else {
00235                     i__2 = kk + j - 1;
00236                     i__3 = kk + j - 1;
00237                     ap[i__2].r = ap[i__3].r, ap[i__2].i = ap[i__3].i;
00238                 }
00239                 jx += *incx;
00240                 kk += j;
00241 /* L40: */
00242             }
00243         }
00244     } else {
00245 
00246 /*        Form  A  when lower triangle is stored in AP. */
00247 
00248         if (*incx == 1) {
00249             i__1 = *n;
00250             for (j = 1; j <= i__1; ++j) {
00251                 i__2 = j;
00252                 if (x[i__2].r != 0. || x[i__2].i != 0.) {
00253                     i__2 = j;
00254                     z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
00255                             z__1.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
00256                             .r;
00257                     temp.r = z__1.r, temp.i = z__1.i;
00258                     i__2 = kk;
00259                     i__3 = kk;
00260                     i__4 = j;
00261                     z__2.r = temp.r * x[i__4].r - temp.i * x[i__4].i, z__2.i =
00262                              temp.r * x[i__4].i + temp.i * x[i__4].r;
00263                     z__1.r = ap[i__3].r + z__2.r, z__1.i = ap[i__3].i + 
00264                             z__2.i;
00265                     ap[i__2].r = z__1.r, ap[i__2].i = z__1.i;
00266                     k = kk + 1;
00267                     i__2 = *n;
00268                     for (i__ = j + 1; i__ <= i__2; ++i__) {
00269                         i__3 = k;
00270                         i__4 = k;
00271                         i__5 = i__;
00272                         z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
00273                                 z__2.i = x[i__5].r * temp.i + x[i__5].i * 
00274                                 temp.r;
00275                         z__1.r = ap[i__4].r + z__2.r, z__1.i = ap[i__4].i + 
00276                                 z__2.i;
00277                         ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
00278                         ++k;
00279 /* L50: */
00280                     }
00281                 } else {
00282                     i__2 = kk;
00283                     i__3 = kk;
00284                     ap[i__2].r = ap[i__3].r, ap[i__2].i = ap[i__3].i;
00285                 }
00286                 kk = kk + *n - j + 1;
00287 /* L60: */
00288             }
00289         } else {
00290             jx = kx;
00291             i__1 = *n;
00292             for (j = 1; j <= i__1; ++j) {
00293                 i__2 = jx;
00294                 if (x[i__2].r != 0. || x[i__2].i != 0.) {
00295                     i__2 = jx;
00296                     z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, 
00297                             z__1.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
00298                             .r;
00299                     temp.r = z__1.r, temp.i = z__1.i;
00300                     i__2 = kk;
00301                     i__3 = kk;
00302                     i__4 = jx;
00303                     z__2.r = temp.r * x[i__4].r - temp.i * x[i__4].i, z__2.i =
00304                              temp.r * x[i__4].i + temp.i * x[i__4].r;
00305                     z__1.r = ap[i__3].r + z__2.r, z__1.i = ap[i__3].i + 
00306                             z__2.i;
00307                     ap[i__2].r = z__1.r, ap[i__2].i = z__1.i;
00308                     ix = jx;
00309                     i__2 = kk + *n - j;
00310                     for (k = kk + 1; k <= i__2; ++k) {
00311                         ix += *incx;
00312                         i__3 = k;
00313                         i__4 = k;
00314                         i__5 = ix;
00315                         z__2.r = x[i__5].r * temp.r - x[i__5].i * temp.i, 
00316                                 z__2.i = x[i__5].r * temp.i + x[i__5].i * 
00317                                 temp.r;
00318                         z__1.r = ap[i__4].r + z__2.r, z__1.i = ap[i__4].i + 
00319                                 z__2.i;
00320                         ap[i__3].r = z__1.r, ap[i__3].i = z__1.i;
00321 /* L70: */
00322                     }
00323                 } else {
00324                     i__2 = kk;
00325                     i__3 = kk;
00326                     ap[i__2].r = ap[i__3].r, ap[i__2].i = ap[i__3].i;
00327                 }
00328                 jx += *incx;
00329                 kk = kk + *n - j + 1;
00330 /* L80: */
00331             }
00332         }
00333     }
00334 
00335     return 0;
00336 
00337 /*     End of ZSPR */
00338 
00339 } /* zspr_ */