ssbmv.c
Go to the documentation of this file.
00001 /* ssbmv.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 ssbmv_(char *uplo, integer *n, integer *k, real *alpha, 
00017         real *a, integer *lda, real *x, integer *incx, real *beta, real *y, 
00018         integer *incy)
00019 {
00020     /* System generated locals */
00021     integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
00022 
00023     /* Local variables */
00024     integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
00025     real temp1, temp2;
00026     extern logical lsame_(char *, char *);
00027     integer kplus1;
00028     extern /* Subroutine */ int xerbla_(char *, integer *);
00029 
00030 /*     .. Scalar Arguments .. */
00031 /*     .. */
00032 /*     .. Array Arguments .. */
00033 /*     .. */
00034 
00035 /*  Purpose */
00036 /*  ======= */
00037 
00038 /*  SSBMV  performs the matrix-vector  operation */
00039 
00040 /*     y := alpha*A*x + beta*y, */
00041 
00042 /*  where alpha and beta are scalars, x and y are n element vectors and */
00043 /*  A is an n by n symmetric band matrix, with k super-diagonals. */
00044 
00045 /*  Arguments */
00046 /*  ========== */
00047 
00048 /*  UPLO   - CHARACTER*1. */
00049 /*           On entry, UPLO specifies whether the upper or lower */
00050 /*           triangular part of the band matrix A is being supplied as */
00051 /*           follows: */
00052 
00053 /*              UPLO = 'U' or 'u'   The upper triangular part of A is */
00054 /*                                  being supplied. */
00055 
00056 /*              UPLO = 'L' or 'l'   The lower triangular part of A is */
00057 /*                                  being supplied. */
00058 
00059 /*           Unchanged on exit. */
00060 
00061 /*  N      - INTEGER. */
00062 /*           On entry, N specifies the order of the matrix A. */
00063 /*           N must be at least zero. */
00064 /*           Unchanged on exit. */
00065 
00066 /*  K      - INTEGER. */
00067 /*           On entry, K specifies the number of super-diagonals of the */
00068 /*           matrix A. K must satisfy  0 .le. K. */
00069 /*           Unchanged on exit. */
00070 
00071 /*  ALPHA  - REAL            . */
00072 /*           On entry, ALPHA specifies the scalar alpha. */
00073 /*           Unchanged on exit. */
00074 
00075 /*  A      - REAL             array of DIMENSION ( LDA, n ). */
00076 /*           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
00077 /*           by n part of the array A must contain the upper triangular */
00078 /*           band part of the symmetric matrix, supplied column by */
00079 /*           column, with the leading diagonal of the matrix in row */
00080 /*           ( k + 1 ) of the array, the first super-diagonal starting at */
00081 /*           position 2 in row k, and so on. The top left k by k triangle */
00082 /*           of the array A is not referenced. */
00083 /*           The following program segment will transfer the upper */
00084 /*           triangular part of a symmetric band matrix from conventional */
00085 /*           full matrix storage to band storage: */
00086 
00087 /*                 DO 20, J = 1, N */
00088 /*                    M = K + 1 - J */
00089 /*                    DO 10, I = MAX( 1, J - K ), J */
00090 /*                       A( M + I, J ) = matrix( I, J ) */
00091 /*              10    CONTINUE */
00092 /*              20 CONTINUE */
00093 
00094 /*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
00095 /*           by n part of the array A must contain the lower triangular */
00096 /*           band part of the symmetric matrix, supplied column by */
00097 /*           column, with the leading diagonal of the matrix in row 1 of */
00098 /*           the array, the first sub-diagonal starting at position 1 in */
00099 /*           row 2, and so on. The bottom right k by k triangle of the */
00100 /*           array A is not referenced. */
00101 /*           The following program segment will transfer the lower */
00102 /*           triangular part of a symmetric band matrix from conventional */
00103 /*           full matrix storage to band storage: */
00104 
00105 /*                 DO 20, J = 1, N */
00106 /*                    M = 1 - J */
00107 /*                    DO 10, I = J, MIN( N, J + K ) */
00108 /*                       A( M + I, J ) = matrix( I, J ) */
00109 /*              10    CONTINUE */
00110 /*              20 CONTINUE */
00111 
00112 /*           Unchanged on exit. */
00113 
00114 /*  LDA    - INTEGER. */
00115 /*           On entry, LDA specifies the first dimension of A as declared */
00116 /*           in the calling (sub) program. LDA must be at least */
00117 /*           ( k + 1 ). */
00118 /*           Unchanged on exit. */
00119 
00120 /*  X      - REAL             array of DIMENSION at least */
00121 /*           ( 1 + ( n - 1 )*abs( INCX ) ). */
00122 /*           Before entry, the incremented array X must contain the */
00123 /*           vector x. */
00124 /*           Unchanged on exit. */
00125 
00126 /*  INCX   - INTEGER. */
00127 /*           On entry, INCX specifies the increment for the elements of */
00128 /*           X. INCX must not be zero. */
00129 /*           Unchanged on exit. */
00130 
00131 /*  BETA   - REAL            . */
00132 /*           On entry, BETA specifies the scalar beta. */
00133 /*           Unchanged on exit. */
00134 
00135 /*  Y      - REAL             array of DIMENSION at least */
00136 /*           ( 1 + ( n - 1 )*abs( INCY ) ). */
00137 /*           Before entry, the incremented array Y must contain the */
00138 /*           vector y. On exit, Y is overwritten by the updated vector y. */
00139 
00140 /*  INCY   - INTEGER. */
00141 /*           On entry, INCY specifies the increment for the elements of */
00142 /*           Y. INCY must not be zero. */
00143 /*           Unchanged on exit. */
00144 
00145 
00146 /*  Level 2 Blas routine. */
00147 
00148 /*  -- Written on 22-October-1986. */
00149 /*     Jack Dongarra, Argonne National Lab. */
00150 /*     Jeremy Du Croz, Nag Central Office. */
00151 /*     Sven Hammarling, Nag Central Office. */
00152 /*     Richard Hanson, Sandia National Labs. */
00153 
00154 
00155 /*     .. Parameters .. */
00156 /*     .. */
00157 /*     .. Local Scalars .. */
00158 /*     .. */
00159 /*     .. External Functions .. */
00160 /*     .. */
00161 /*     .. External Subroutines .. */
00162 /*     .. */
00163 /*     .. Intrinsic Functions .. */
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     --x;
00173     --y;
00174 
00175     /* Function Body */
00176     info = 0;
00177     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00178         info = 1;
00179     } else if (*n < 0) {
00180         info = 2;
00181     } else if (*k < 0) {
00182         info = 3;
00183     } else if (*lda < *k + 1) {
00184         info = 6;
00185     } else if (*incx == 0) {
00186         info = 8;
00187     } else if (*incy == 0) {
00188         info = 11;
00189     }
00190     if (info != 0) {
00191         xerbla_("SSBMV ", &info);
00192         return 0;
00193     }
00194 
00195 /*     Quick return if possible. */
00196 
00197     if (*n == 0 || *alpha == 0.f && *beta == 1.f) {
00198         return 0;
00199     }
00200 
00201 /*     Set up the start points in  X  and  Y. */
00202 
00203     if (*incx > 0) {
00204         kx = 1;
00205     } else {
00206         kx = 1 - (*n - 1) * *incx;
00207     }
00208     if (*incy > 0) {
00209         ky = 1;
00210     } else {
00211         ky = 1 - (*n - 1) * *incy;
00212     }
00213 
00214 /*     Start the operations. In this version the elements of the array A */
00215 /*     are accessed sequentially with one pass through A. */
00216 
00217 /*     First form  y := beta*y. */
00218 
00219     if (*beta != 1.f) {
00220         if (*incy == 1) {
00221             if (*beta == 0.f) {
00222                 i__1 = *n;
00223                 for (i__ = 1; i__ <= i__1; ++i__) {
00224                     y[i__] = 0.f;
00225 /* L10: */
00226                 }
00227             } else {
00228                 i__1 = *n;
00229                 for (i__ = 1; i__ <= i__1; ++i__) {
00230                     y[i__] = *beta * y[i__];
00231 /* L20: */
00232                 }
00233             }
00234         } else {
00235             iy = ky;
00236             if (*beta == 0.f) {
00237                 i__1 = *n;
00238                 for (i__ = 1; i__ <= i__1; ++i__) {
00239                     y[iy] = 0.f;
00240                     iy += *incy;
00241 /* L30: */
00242                 }
00243             } else {
00244                 i__1 = *n;
00245                 for (i__ = 1; i__ <= i__1; ++i__) {
00246                     y[iy] = *beta * y[iy];
00247                     iy += *incy;
00248 /* L40: */
00249                 }
00250             }
00251         }
00252     }
00253     if (*alpha == 0.f) {
00254         return 0;
00255     }
00256     if (lsame_(uplo, "U")) {
00257 
00258 /*        Form  y  when upper triangle of A is stored. */
00259 
00260         kplus1 = *k + 1;
00261         if (*incx == 1 && *incy == 1) {
00262             i__1 = *n;
00263             for (j = 1; j <= i__1; ++j) {
00264                 temp1 = *alpha * x[j];
00265                 temp2 = 0.f;
00266                 l = kplus1 - j;
00267 /* Computing MAX */
00268                 i__2 = 1, i__3 = j - *k;
00269                 i__4 = j - 1;
00270                 for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
00271                     y[i__] += temp1 * a[l + i__ + j * a_dim1];
00272                     temp2 += a[l + i__ + j * a_dim1] * x[i__];
00273 /* L50: */
00274                 }
00275                 y[j] = y[j] + temp1 * a[kplus1 + j * a_dim1] + *alpha * temp2;
00276 /* L60: */
00277             }
00278         } else {
00279             jx = kx;
00280             jy = ky;
00281             i__1 = *n;
00282             for (j = 1; j <= i__1; ++j) {
00283                 temp1 = *alpha * x[jx];
00284                 temp2 = 0.f;
00285                 ix = kx;
00286                 iy = ky;
00287                 l = kplus1 - j;
00288 /* Computing MAX */
00289                 i__4 = 1, i__2 = j - *k;
00290                 i__3 = j - 1;
00291                 for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
00292                     y[iy] += temp1 * a[l + i__ + j * a_dim1];
00293                     temp2 += a[l + i__ + j * a_dim1] * x[ix];
00294                     ix += *incx;
00295                     iy += *incy;
00296 /* L70: */
00297                 }
00298                 y[jy] = y[jy] + temp1 * a[kplus1 + j * a_dim1] + *alpha * 
00299                         temp2;
00300                 jx += *incx;
00301                 jy += *incy;
00302                 if (j > *k) {
00303                     kx += *incx;
00304                     ky += *incy;
00305                 }
00306 /* L80: */
00307             }
00308         }
00309     } else {
00310 
00311 /*        Form  y  when lower triangle of A is stored. */
00312 
00313         if (*incx == 1 && *incy == 1) {
00314             i__1 = *n;
00315             for (j = 1; j <= i__1; ++j) {
00316                 temp1 = *alpha * x[j];
00317                 temp2 = 0.f;
00318                 y[j] += temp1 * a[j * a_dim1 + 1];
00319                 l = 1 - j;
00320 /* Computing MIN */
00321                 i__4 = *n, i__2 = j + *k;
00322                 i__3 = min(i__4,i__2);
00323                 for (i__ = j + 1; i__ <= i__3; ++i__) {
00324                     y[i__] += temp1 * a[l + i__ + j * a_dim1];
00325                     temp2 += a[l + i__ + j * a_dim1] * x[i__];
00326 /* L90: */
00327                 }
00328                 y[j] += *alpha * temp2;
00329 /* L100: */
00330             }
00331         } else {
00332             jx = kx;
00333             jy = ky;
00334             i__1 = *n;
00335             for (j = 1; j <= i__1; ++j) {
00336                 temp1 = *alpha * x[jx];
00337                 temp2 = 0.f;
00338                 y[jy] += temp1 * a[j * a_dim1 + 1];
00339                 l = 1 - j;
00340                 ix = jx;
00341                 iy = jy;
00342 /* Computing MIN */
00343                 i__4 = *n, i__2 = j + *k;
00344                 i__3 = min(i__4,i__2);
00345                 for (i__ = j + 1; i__ <= i__3; ++i__) {
00346                     ix += *incx;
00347                     iy += *incy;
00348                     y[iy] += temp1 * a[l + i__ + j * a_dim1];
00349                     temp2 += a[l + i__ + j * a_dim1] * x[ix];
00350 /* L110: */
00351                 }
00352                 y[jy] += *alpha * temp2;
00353                 jx += *incx;
00354                 jy += *incy;
00355 /* L120: */
00356             }
00357         }
00358     }
00359 
00360     return 0;
00361 
00362 /*     End of SSBMV . */
00363 
00364 } /* ssbmv_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:13