zla_heamv.c
Go to the documentation of this file.
00001 /* zla_heamv.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 zla_heamv__(integer *uplo, integer *n, doublereal *alpha,
00017          doublecomplex *a, integer *lda, doublecomplex *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, i__3;
00022     doublereal d__1, d__2;
00023 
00024     /* Builtin functions */
00025     double d_imag(doublecomplex *), 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 /*  ZLA_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      - COMPLEX*16         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      - COMPLEX*16         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 /*     .. Statement Functions .. */
00156 /*     .. */
00157 /*     .. Statement Function Definitions .. */
00158 /*     .. */
00159 /*     .. Executable Statements .. */
00160 
00161 /*     Test the input parameters. */
00162 
00163     /* Parameter adjustments */
00164     a_dim1 = *lda;
00165     a_offset = 1 + a_dim1;
00166     a -= a_offset;
00167     --x;
00168     --y;
00169 
00170     /* Function Body */
00171     info = 0;
00172     if (*uplo != ilauplo_("U") && *uplo != ilauplo_("L")
00173             ) {
00174         info = 1;
00175     } else if (*n < 0) {
00176         info = 2;
00177     } else if (*lda < max(1,*n)) {
00178         info = 5;
00179     } else if (*incx == 0) {
00180         info = 7;
00181     } else if (*incy == 0) {
00182         info = 10;
00183     }
00184     if (info != 0) {
00185         xerbla_("ZHEMV ", &info);
00186         return 0;
00187     }
00188 
00189 /*     Quick return if possible. */
00190 
00191     if (*n == 0 || *alpha == 0. && *beta == 1.) {
00192         return 0;
00193     }
00194 
00195 /*     Set up the start points in  X  and  Y. */
00196 
00197     if (*incx > 0) {
00198         kx = 1;
00199     } else {
00200         kx = 1 - (*n - 1) * *incx;
00201     }
00202     if (*incy > 0) {
00203         ky = 1;
00204     } else {
00205         ky = 1 - (*n - 1) * *incy;
00206     }
00207 
00208 /*     Set SAFE1 essentially to be the underflow threshold times the */
00209 /*     number of additions in each row. */
00210 
00211     safe1 = dlamch_("Safe minimum");
00212     safe1 = (*n + 1) * safe1;
00213 
00214 /*     Form  y := alpha*abs(A)*abs(x) + beta*abs(y). */
00215 
00216 /*     The O(N^2) SYMB_ZERO tests could be replaced by O(N) queries to */
00217 /*     the inexact flag.  Still doesn't help change the iteration order */
00218 /*     to per-column. */
00219 
00220     iy = ky;
00221     if (*incx == 1) {
00222         i__1 = *n;
00223         for (i__ = 1; i__ <= i__1; ++i__) {
00224             if (*beta == 0.) {
00225                 symb_zero__ = TRUE_;
00226                 y[iy] = 0.;
00227             } else if (y[iy] == 0.) {
00228                 symb_zero__ = TRUE_;
00229             } else {
00230                 symb_zero__ = FALSE_;
00231                 y[iy] = *beta * (d__1 = y[iy], abs(d__1));
00232             }
00233             if (*alpha != 0.) {
00234                 i__2 = *n;
00235                 for (j = 1; j <= i__2; ++j) {
00236                     if (*uplo == ilauplo_("U")) {
00237                         if (i__ <= j) {
00238                             i__3 = i__ + j * a_dim1;
00239                             temp = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = 
00240                                     d_imag(&a[i__ + j * a_dim1]), abs(d__2));
00241                         } else {
00242                             i__3 = j + i__ * a_dim1;
00243                             temp = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = 
00244                                     d_imag(&a[j + i__ * a_dim1]), abs(d__2));
00245                         }
00246                     } else {
00247                         if (i__ >= j) {
00248                             i__3 = i__ + j * a_dim1;
00249                             temp = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = 
00250                                     d_imag(&a[i__ + j * a_dim1]), abs(d__2));
00251                         } else {
00252                             i__3 = j + i__ * a_dim1;
00253                             temp = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = 
00254                                     d_imag(&a[j + i__ * a_dim1]), abs(d__2));
00255                         }
00256                     }
00257                     i__3 = j;
00258                     symb_zero__ = symb_zero__ && (x[i__3].r == 0. && x[i__3]
00259                             .i == 0. || temp == 0.);
00260                     i__3 = j;
00261                     y[iy] += *alpha * ((d__1 = x[i__3].r, abs(d__1)) + (d__2 =
00262                              d_imag(&x[j]), abs(d__2))) * temp;
00263                 }
00264             }
00265             if (! symb_zero__) {
00266                 y[iy] += d_sign(&safe1, &y[iy]);
00267             }
00268             iy += *incy;
00269         }
00270     } else {
00271         i__1 = *n;
00272         for (i__ = 1; i__ <= i__1; ++i__) {
00273             if (*beta == 0.) {
00274                 symb_zero__ = TRUE_;
00275                 y[iy] = 0.;
00276             } else if (y[iy] == 0.) {
00277                 symb_zero__ = TRUE_;
00278             } else {
00279                 symb_zero__ = FALSE_;
00280                 y[iy] = *beta * (d__1 = y[iy], abs(d__1));
00281             }
00282             jx = kx;
00283             if (*alpha != 0.) {
00284                 i__2 = *n;
00285                 for (j = 1; j <= i__2; ++j) {
00286                     if (*uplo == ilauplo_("U")) {
00287                         if (i__ <= j) {
00288                             i__3 = i__ + j * a_dim1;
00289                             temp = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = 
00290                                     d_imag(&a[i__ + j * a_dim1]), abs(d__2));
00291                         } else {
00292                             i__3 = j + i__ * a_dim1;
00293                             temp = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = 
00294                                     d_imag(&a[j + i__ * a_dim1]), abs(d__2));
00295                         }
00296                     } else {
00297                         if (i__ >= j) {
00298                             i__3 = i__ + j * a_dim1;
00299                             temp = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = 
00300                                     d_imag(&a[i__ + j * a_dim1]), abs(d__2));
00301                         } else {
00302                             i__3 = j + i__ * a_dim1;
00303                             temp = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = 
00304                                     d_imag(&a[j + i__ * a_dim1]), abs(d__2));
00305                         }
00306                     }
00307                     i__3 = j;
00308                     symb_zero__ = symb_zero__ && (x[i__3].r == 0. && x[i__3]
00309                             .i == 0. || temp == 0.);
00310                     i__3 = jx;
00311                     y[iy] += *alpha * ((d__1 = x[i__3].r, abs(d__1)) + (d__2 =
00312                              d_imag(&x[jx]), abs(d__2))) * temp;
00313                     jx += *incx;
00314                 }
00315             }
00316             if (! symb_zero__) {
00317                 y[iy] += d_sign(&safe1, &y[iy]);
00318             }
00319             iy += *incy;
00320         }
00321     }
00322 
00323     return 0;
00324 
00325 /*     End of ZLA_HEAMV */
00326 
00327 } /* zla_heamv__ */


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