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


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