dla_syamv.c
Go to the documentation of this file.
00001 /* dla_syamv.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 dla_syamv__(integer *uplo, integer *n, doublereal *alpha,
00017          doublereal *a, integer *lda, doublereal *x, integer *incx, 
00018         doublereal *beta, doublereal *y, integer *incy)
00019 {
00020     /* System generated locals */
00021     integer a_dim1, a_offset, i__1, i__2;
00022     doublereal d__1;
00023 
00024     /* Builtin functions */
00025     double d_sign(doublereal *, doublereal *);
00026 
00027     /* Local variables */
00028     integer i__, j;
00029     logical symb_zero__;
00030     integer iy, jx, kx, ky, info;
00031     doublereal temp, safe1;
00032     extern doublereal dlamch_(char *);
00033     extern /* Subroutine */ int xerbla_(char *, integer *);
00034     extern integer ilauplo_(char *);
00035 
00036 
00037 /*     -- LAPACK routine (version 3.2)                                 -- */
00038 /*     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- */
00039 /*     -- Jason Riedy of Univ. of California Berkeley.                 -- */
00040 /*     -- November 2008                                                -- */
00041 
00042 /*     -- LAPACK is a software package provided by Univ. of Tennessee, -- */
00043 /*     -- Univ. of California Berkeley and NAG Ltd.                    -- */
00044 
00045 /*     .. */
00046 /*     .. Scalar Arguments .. */
00047 /*     .. */
00048 /*     .. Array Arguments .. */
00049 /*     .. */
00050 
00051 /*  Purpose */
00052 /*  ======= */
00053 
00054 /*  DLA_SYAMV  performs the matrix-vector operation */
00055 
00056 /*          y := alpha*abs(A)*abs(x) + beta*abs(y), */
00057 
00058 /*  where alpha and beta are scalars, x and y are vectors and A is an */
00059 /*  n by n symmetric matrix. */
00060 
00061 /*  This function is primarily used in calculating error bounds. */
00062 /*  To protect against underflow during evaluation, components in */
00063 /*  the resulting vector are perturbed away from zero by (N+1) */
00064 /*  times the underflow threshold.  To prevent unnecessarily large */
00065 /*  errors for block-structure embedded in general matrices, */
00066 /*  "symbolically" zero components are not perturbed.  A zero */
00067 /*  entry is considered "symbolic" if all multiplications involved */
00068 /*  in computing that entry have at least one zero multiplicand. */
00069 
00070 /*  Parameters */
00071 /*  ========== */
00072 
00073 /*  UPLO   - INTEGER */
00074 /*           On entry, UPLO specifies whether the upper or lower */
00075 /*           triangular part of the array A is to be referenced as */
00076 /*           follows: */
00077 
00078 /*              UPLO = BLAS_UPPER   Only the upper triangular part of A */
00079 /*                                  is to be referenced. */
00080 
00081 /*              UPLO = BLAS_LOWER   Only the lower triangular part of A */
00082 /*                                  is to be referenced. */
00083 
00084 /*           Unchanged on exit. */
00085 
00086 /*  N      - INTEGER. */
00087 /*           On entry, N specifies the number of columns of the matrix A. */
00088 /*           N must be at least zero. */
00089 /*           Unchanged on exit. */
00090 
00091 /*  ALPHA  - DOUBLE PRECISION   . */
00092 /*           On entry, ALPHA specifies the scalar alpha. */
00093 /*           Unchanged on exit. */
00094 
00095 /*  A      - DOUBLE PRECISION   array of DIMENSION ( LDA, n ). */
00096 /*           Before entry, the leading m by n part of the array A must */
00097 /*           contain the matrix of coefficients. */
00098 /*           Unchanged on exit. */
00099 
00100 /*  LDA    - INTEGER. */
00101 /*           On entry, LDA specifies the first dimension of A as declared */
00102 /*           in the calling (sub) program. LDA must be at least */
00103 /*           max( 1, n ). */
00104 /*           Unchanged on exit. */
00105 
00106 /*  X      - DOUBLE PRECISION   array of DIMENSION at least */
00107 /*           ( 1 + ( n - 1 )*abs( INCX ) ) */
00108 /*           Before entry, the incremented array X must contain the */
00109 /*           vector x. */
00110 /*           Unchanged on exit. */
00111 
00112 /*  INCX   - INTEGER. */
00113 /*           On entry, INCX specifies the increment for the elements of */
00114 /*           X. INCX must not be zero. */
00115 /*           Unchanged on exit. */
00116 
00117 /*  BETA   - DOUBLE PRECISION   . */
00118 /*           On entry, BETA specifies the scalar beta. When BETA is */
00119 /*           supplied as zero then Y need not be set on input. */
00120 /*           Unchanged on exit. */
00121 
00122 /*  Y      - DOUBLE PRECISION   array of DIMENSION at least */
00123 /*           ( 1 + ( n - 1 )*abs( INCY ) ) */
00124 /*           Before entry with BETA non-zero, the incremented array Y */
00125 /*           must contain the vector y. On exit, Y is overwritten by the */
00126 /*           updated vector y. */
00127 
00128 /*  INCY   - INTEGER. */
00129 /*           On entry, INCY specifies the increment for the elements of */
00130 /*           Y. INCY must not be zero. */
00131 /*           Unchanged on exit. */
00132 
00133 
00134 /*  Level 2 Blas routine. */
00135 
00136 /*  -- Written on 22-October-1986. */
00137 /*     Jack Dongarra, Argonne National Lab. */
00138 /*     Jeremy Du Croz, Nag Central Office. */
00139 /*     Sven Hammarling, Nag Central Office. */
00140 /*     Richard Hanson, Sandia National Labs. */
00141 /*  -- Modified for the absolute-value product, April 2006 */
00142 /*     Jason Riedy, UC Berkeley */
00143 
00144 /*     .. */
00145 /*     .. Parameters .. */
00146 /*     .. */
00147 /*     .. Local Scalars .. */
00148 /*     .. */
00149 /*     .. External Subroutines .. */
00150 /*     .. */
00151 /*     .. External Functions .. */
00152 /*     .. */
00153 /*     .. Intrinsic Functions .. */
00154 /*     .. */
00155 /*     .. Executable Statements .. */
00156 
00157 /*     Test the input parameters. */
00158 
00159     /* Parameter adjustments */
00160     a_dim1 = *lda;
00161     a_offset = 1 + a_dim1;
00162     a -= a_offset;
00163     --x;
00164     --y;
00165 
00166     /* Function Body */
00167     info = 0;
00168     if (*uplo != ilauplo_("U") && *uplo != ilauplo_("L")
00169             ) {
00170         info = 1;
00171     } else if (*n < 0) {
00172         info = 2;
00173     } else if (*lda < max(1,*n)) {
00174         info = 5;
00175     } else if (*incx == 0) {
00176         info = 7;
00177     } else if (*incy == 0) {
00178         info = 10;
00179     }
00180     if (info != 0) {
00181         xerbla_("DSYMV ", &info);
00182         return 0;
00183     }
00184 
00185 /*     Quick return if possible. */
00186 
00187     if (*n == 0 || *alpha == 0. && *beta == 1.) {
00188         return 0;
00189     }
00190 
00191 /*     Set up the start points in  X  and  Y. */
00192 
00193     if (*incx > 0) {
00194         kx = 1;
00195     } else {
00196         kx = 1 - (*n - 1) * *incx;
00197     }
00198     if (*incy > 0) {
00199         ky = 1;
00200     } else {
00201         ky = 1 - (*n - 1) * *incy;
00202     }
00203 
00204 /*     Set SAFE1 essentially to be the underflow threshold times the */
00205 /*     number of additions in each row. */
00206 
00207     safe1 = dlamch_("Safe minimum");
00208     safe1 = (*n + 1) * safe1;
00209 
00210 /*     Form  y := alpha*abs(A)*abs(x) + beta*abs(y). */
00211 
00212 /*     The O(N^2) SYMB_ZERO tests could be replaced by O(N) queries to */
00213 /*     the inexact flag.  Still doesn't help change the iteration order */
00214 /*     to per-column. */
00215 
00216     iy = ky;
00217     if (*incx == 1) {
00218         i__1 = *n;
00219         for (i__ = 1; i__ <= i__1; ++i__) {
00220             if (*beta == 0.) {
00221                 symb_zero__ = TRUE_;
00222                 y[iy] = 0.;
00223             } else if (y[iy] == 0.) {
00224                 symb_zero__ = TRUE_;
00225             } else {
00226                 symb_zero__ = FALSE_;
00227                 y[iy] = *beta * (d__1 = y[iy], abs(d__1));
00228             }
00229             if (*alpha != 0.) {
00230                 i__2 = *n;
00231                 for (j = 1; j <= i__2; ++j) {
00232                     if (*uplo == ilauplo_("U")) {
00233                         if (i__ <= j) {
00234                             temp = (d__1 = a[i__ + j * a_dim1], abs(d__1));
00235                         } else {
00236                             temp = (d__1 = a[j + i__ * a_dim1], abs(d__1));
00237                         }
00238                     } else {
00239                         if (i__ >= j) {
00240                             temp = (d__1 = a[i__ + j * a_dim1], abs(d__1));
00241                         } else {
00242                             temp = (d__1 = a[j + i__ * a_dim1], abs(d__1));
00243                         }
00244                     }
00245                     symb_zero__ = symb_zero__ && (x[j] == 0. || temp == 0.);
00246                     y[iy] += *alpha * (d__1 = x[j], abs(d__1)) * temp;
00247                 }
00248             }
00249             if (! symb_zero__) {
00250                 y[iy] += d_sign(&safe1, &y[iy]);
00251             }
00252             iy += *incy;
00253         }
00254     } else {
00255         i__1 = *n;
00256         for (i__ = 1; i__ <= i__1; ++i__) {
00257             if (*beta == 0.) {
00258                 symb_zero__ = TRUE_;
00259                 y[iy] = 0.;
00260             } else if (y[iy] == 0.) {
00261                 symb_zero__ = TRUE_;
00262             } else {
00263                 symb_zero__ = FALSE_;
00264                 y[iy] = *beta * (d__1 = y[iy], abs(d__1));
00265             }
00266             jx = kx;
00267             if (*alpha != 0.) {
00268                 i__2 = *n;
00269                 for (j = 1; j <= i__2; ++j) {
00270                     if (*uplo == ilauplo_("U")) {
00271                         if (i__ <= j) {
00272                             temp = (d__1 = a[i__ + j * a_dim1], abs(d__1));
00273                         } else {
00274                             temp = (d__1 = a[j + i__ * a_dim1], abs(d__1));
00275                         }
00276                     } else {
00277                         if (i__ >= j) {
00278                             temp = (d__1 = a[i__ + j * a_dim1], abs(d__1));
00279                         } else {
00280                             temp = (d__1 = a[j + i__ * a_dim1], abs(d__1));
00281                         }
00282                     }
00283                     symb_zero__ = symb_zero__ && (x[j] == 0. || temp == 0.);
00284                     y[iy] += *alpha * (d__1 = x[jx], abs(d__1)) * temp;
00285                     jx += *incx;
00286                 }
00287             }
00288             if (! symb_zero__) {
00289                 y[iy] += d_sign(&safe1, &y[iy]);
00290             }
00291             iy += *incy;
00292         }
00293     }
00294 
00295     return 0;
00296 
00297 /*     End of DLA_SYAMV */
00298 
00299 } /* dla_syamv__ */


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