1/* zhpmv.f -- translated by f2c (version 20100827). 2 You must link the resulting object file with libf2c: 3 on Microsoft Windows system, link with libf2c.lib; 4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 5 or, if you install libf2c.a in a standard place, with -lf2c -lm 6 -- in that order, at the end of the command line, as in 7 cc *.o -lf2c -lm 8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 9 10 http://www.netlib.org/f2c/libf2c.zip 11*/ 12 13#include "datatypes.h" 14 15/* Subroutine */ int zhpmv_(char *uplo, integer *n, doublecomplex *alpha, 16 doublecomplex *ap, doublecomplex *x, integer *incx, doublecomplex * 17 beta, doublecomplex *y, integer *incy, ftnlen uplo_len) 18{ 19 /* System generated locals */ 20 integer i__1, i__2, i__3, i__4, i__5; 21 doublereal d__1; 22 doublecomplex z__1, z__2, z__3, z__4; 23 24 /* Builtin functions */ 25 void d_cnjg(doublecomplex *, doublecomplex *); 26 27 /* Local variables */ 28 integer i__, j, k, kk, ix, iy, jx, jy, kx, ky, info; 29 doublecomplex temp1, temp2; 30 extern logical lsame_(char *, char *, ftnlen, ftnlen); 31 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen); 32 33/* .. Scalar Arguments .. */ 34/* .. */ 35/* .. Array Arguments .. */ 36/* .. */ 37 38/* Purpose */ 39/* ======= */ 40 41/* ZHPMV performs the matrix-vector operation */ 42 43/* y := alpha*A*x + beta*y, */ 44 45/* where alpha and beta are scalars, x and y are n element vectors and */ 46/* A is an n by n hermitian matrix, supplied in packed form. */ 47 48/* Arguments */ 49/* ========== */ 50 51/* UPLO - CHARACTER*1. */ 52/* On entry, UPLO specifies whether the upper or lower */ 53/* triangular part of the matrix A is supplied in the packed */ 54/* array AP as follows: */ 55 56/* UPLO = 'U' or 'u' The upper triangular part of A is */ 57/* supplied in AP. */ 58 59/* UPLO = 'L' or 'l' The lower triangular part of A is */ 60/* supplied in AP. */ 61 62/* Unchanged on exit. */ 63 64/* N - INTEGER. */ 65/* On entry, N specifies the order of the matrix A. */ 66/* N must be at least zero. */ 67/* Unchanged on exit. */ 68 69/* ALPHA - COMPLEX*16 . */ 70/* On entry, ALPHA specifies the scalar alpha. */ 71/* Unchanged on exit. */ 72 73/* AP - COMPLEX*16 array of DIMENSION at least */ 74/* ( ( n*( n + 1 ) )/2 ). */ 75/* Before entry with UPLO = 'U' or 'u', the array AP must */ 76/* contain the upper triangular part of the hermitian matrix */ 77/* packed sequentially, column by column, so that AP( 1 ) */ 78/* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */ 79/* and a( 2, 2 ) respectively, and so on. */ 80/* Before entry with UPLO = 'L' or 'l', the array AP must */ 81/* contain the lower triangular part of the hermitian matrix */ 82/* packed sequentially, column by column, so that AP( 1 ) */ 83/* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */ 84/* and a( 3, 1 ) respectively, and so on. */ 85/* Note that the imaginary parts of the diagonal elements need */ 86/* not be set and are assumed to be zero. */ 87/* Unchanged on exit. */ 88 89/* X - COMPLEX*16 array of dimension at least */ 90/* ( 1 + ( n - 1 )*abs( INCX ) ). */ 91/* Before entry, the incremented array X must contain the n */ 92/* element vector x. */ 93/* Unchanged on exit. */ 94 95/* INCX - INTEGER. */ 96/* On entry, INCX specifies the increment for the elements of */ 97/* X. INCX must not be zero. */ 98/* Unchanged on exit. */ 99 100/* BETA - COMPLEX*16 . */ 101/* On entry, BETA specifies the scalar beta. When BETA is */ 102/* supplied as zero then Y need not be set on input. */ 103/* Unchanged on exit. */ 104 105/* Y - COMPLEX*16 array of dimension at least */ 106/* ( 1 + ( n - 1 )*abs( INCY ) ). */ 107/* Before entry, the incremented array Y must contain the n */ 108/* element vector y. On exit, Y is overwritten by the updated */ 109/* vector y. */ 110 111/* INCY - INTEGER. */ 112/* On entry, INCY specifies the increment for the elements of */ 113/* Y. INCY must not be zero. */ 114/* Unchanged on exit. */ 115 116/* Further Details */ 117/* =============== */ 118 119/* Level 2 Blas routine. */ 120 121/* -- Written on 22-October-1986. */ 122/* Jack Dongarra, Argonne National Lab. */ 123/* Jeremy Du Croz, Nag Central Office. */ 124/* Sven Hammarling, Nag Central Office. */ 125/* Richard Hanson, Sandia National Labs. */ 126 127/* ===================================================================== */ 128 129/* .. Parameters .. */ 130/* .. */ 131/* .. Local Scalars .. */ 132/* .. */ 133/* .. External Functions .. */ 134/* .. */ 135/* .. External Subroutines .. */ 136/* .. */ 137/* .. Intrinsic Functions .. */ 138/* .. */ 139 140/* Test the input parameters. */ 141 142 /* Parameter adjustments */ 143 --y; 144 --x; 145 --ap; 146 147 /* Function Body */ 148 info = 0; 149 if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", ( 150 ftnlen)1, (ftnlen)1)) { 151 info = 1; 152 } else if (*n < 0) { 153 info = 2; 154 } else if (*incx == 0) { 155 info = 6; 156 } else if (*incy == 0) { 157 info = 9; 158 } 159 if (info != 0) { 160 xerbla_("ZHPMV ", &info, (ftnlen)6); 161 return 0; 162 } 163 164/* Quick return if possible. */ 165 166 if (*n == 0 || (alpha->r == 0. && alpha->i == 0. && (beta->r == 1. && 167 beta->i == 0.))) { 168 return 0; 169 } 170 171/* Set up the start points in X and Y. */ 172 173 if (*incx > 0) { 174 kx = 1; 175 } else { 176 kx = 1 - (*n - 1) * *incx; 177 } 178 if (*incy > 0) { 179 ky = 1; 180 } else { 181 ky = 1 - (*n - 1) * *incy; 182 } 183 184/* Start the operations. In this version the elements of the array AP */ 185/* are accessed sequentially with one pass through AP. */ 186 187/* First form y := beta*y. */ 188 189 if (beta->r != 1. || beta->i != 0.) { 190 if (*incy == 1) { 191 if (beta->r == 0. && beta->i == 0.) { 192 i__1 = *n; 193 for (i__ = 1; i__ <= i__1; ++i__) { 194 i__2 = i__; 195 y[i__2].r = 0., y[i__2].i = 0.; 196/* L10: */ 197 } 198 } else { 199 i__1 = *n; 200 for (i__ = 1; i__ <= i__1; ++i__) { 201 i__2 = i__; 202 i__3 = i__; 203 z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 204 z__1.i = beta->r * y[i__3].i + beta->i * y[i__3] 205 .r; 206 y[i__2].r = z__1.r, y[i__2].i = z__1.i; 207/* L20: */ 208 } 209 } 210 } else { 211 iy = ky; 212 if (beta->r == 0. && beta->i == 0.) { 213 i__1 = *n; 214 for (i__ = 1; i__ <= i__1; ++i__) { 215 i__2 = iy; 216 y[i__2].r = 0., y[i__2].i = 0.; 217 iy += *incy; 218/* L30: */ 219 } 220 } else { 221 i__1 = *n; 222 for (i__ = 1; i__ <= i__1; ++i__) { 223 i__2 = iy; 224 i__3 = iy; 225 z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 226 z__1.i = beta->r * y[i__3].i + beta->i * y[i__3] 227 .r; 228 y[i__2].r = z__1.r, y[i__2].i = z__1.i; 229 iy += *incy; 230/* L40: */ 231 } 232 } 233 } 234 } 235 if (alpha->r == 0. && alpha->i == 0.) { 236 return 0; 237 } 238 kk = 1; 239 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) { 240 241/* Form y when AP contains the upper triangle. */ 242 243 if (*incx == 1 && *incy == 1) { 244 i__1 = *n; 245 for (j = 1; j <= i__1; ++j) { 246 i__2 = j; 247 z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i = 248 alpha->r * x[i__2].i + alpha->i * x[i__2].r; 249 temp1.r = z__1.r, temp1.i = z__1.i; 250 temp2.r = 0., temp2.i = 0.; 251 k = kk; 252 i__2 = j - 1; 253 for (i__ = 1; i__ <= i__2; ++i__) { 254 i__3 = i__; 255 i__4 = i__; 256 i__5 = k; 257 z__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 258 z__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5] 259 .r; 260 z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i; 261 y[i__3].r = z__1.r, y[i__3].i = z__1.i; 262 d_cnjg(&z__3, &ap[k]); 263 i__3 = i__; 264 z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i, z__2.i = 265 z__3.r * x[i__3].i + z__3.i * x[i__3].r; 266 z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i; 267 temp2.r = z__1.r, temp2.i = z__1.i; 268 ++k; 269/* L50: */ 270 } 271 i__2 = j; 272 i__3 = j; 273 i__4 = kk + j - 1; 274 d__1 = ap[i__4].r; 275 z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i; 276 z__2.r = y[i__3].r + z__3.r, z__2.i = y[i__3].i + z__3.i; 277 z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i = 278 alpha->r * temp2.i + alpha->i * temp2.r; 279 z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i; 280 y[i__2].r = z__1.r, y[i__2].i = z__1.i; 281 kk += j; 282/* L60: */ 283 } 284 } else { 285 jx = kx; 286 jy = ky; 287 i__1 = *n; 288 for (j = 1; j <= i__1; ++j) { 289 i__2 = jx; 290 z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i = 291 alpha->r * x[i__2].i + alpha->i * x[i__2].r; 292 temp1.r = z__1.r, temp1.i = z__1.i; 293 temp2.r = 0., temp2.i = 0.; 294 ix = kx; 295 iy = ky; 296 i__2 = kk + j - 2; 297 for (k = kk; k <= i__2; ++k) { 298 i__3 = iy; 299 i__4 = iy; 300 i__5 = k; 301 z__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 302 z__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5] 303 .r; 304 z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i; 305 y[i__3].r = z__1.r, y[i__3].i = z__1.i; 306 d_cnjg(&z__3, &ap[k]); 307 i__3 = ix; 308 z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i, z__2.i = 309 z__3.r * x[i__3].i + z__3.i * x[i__3].r; 310 z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i; 311 temp2.r = z__1.r, temp2.i = z__1.i; 312 ix += *incx; 313 iy += *incy; 314/* L70: */ 315 } 316 i__2 = jy; 317 i__3 = jy; 318 i__4 = kk + j - 1; 319 d__1 = ap[i__4].r; 320 z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i; 321 z__2.r = y[i__3].r + z__3.r, z__2.i = y[i__3].i + z__3.i; 322 z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i = 323 alpha->r * temp2.i + alpha->i * temp2.r; 324 z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i; 325 y[i__2].r = z__1.r, y[i__2].i = z__1.i; 326 jx += *incx; 327 jy += *incy; 328 kk += j; 329/* L80: */ 330 } 331 } 332 } else { 333 334/* Form y when AP contains the lower triangle. */ 335 336 if (*incx == 1 && *incy == 1) { 337 i__1 = *n; 338 for (j = 1; j <= i__1; ++j) { 339 i__2 = j; 340 z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i = 341 alpha->r * x[i__2].i + alpha->i * x[i__2].r; 342 temp1.r = z__1.r, temp1.i = z__1.i; 343 temp2.r = 0., temp2.i = 0.; 344 i__2 = j; 345 i__3 = j; 346 i__4 = kk; 347 d__1 = ap[i__4].r; 348 z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i; 349 z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i; 350 y[i__2].r = z__1.r, y[i__2].i = z__1.i; 351 k = kk + 1; 352 i__2 = *n; 353 for (i__ = j + 1; i__ <= i__2; ++i__) { 354 i__3 = i__; 355 i__4 = i__; 356 i__5 = k; 357 z__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 358 z__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5] 359 .r; 360 z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i; 361 y[i__3].r = z__1.r, y[i__3].i = z__1.i; 362 d_cnjg(&z__3, &ap[k]); 363 i__3 = i__; 364 z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i, z__2.i = 365 z__3.r * x[i__3].i + z__3.i * x[i__3].r; 366 z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i; 367 temp2.r = z__1.r, temp2.i = z__1.i; 368 ++k; 369/* L90: */ 370 } 371 i__2 = j; 372 i__3 = j; 373 z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i = 374 alpha->r * temp2.i + alpha->i * temp2.r; 375 z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i; 376 y[i__2].r = z__1.r, y[i__2].i = z__1.i; 377 kk += *n - j + 1; 378/* L100: */ 379 } 380 } else { 381 jx = kx; 382 jy = ky; 383 i__1 = *n; 384 for (j = 1; j <= i__1; ++j) { 385 i__2 = jx; 386 z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i = 387 alpha->r * x[i__2].i + alpha->i * x[i__2].r; 388 temp1.r = z__1.r, temp1.i = z__1.i; 389 temp2.r = 0., temp2.i = 0.; 390 i__2 = jy; 391 i__3 = jy; 392 i__4 = kk; 393 d__1 = ap[i__4].r; 394 z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i; 395 z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i; 396 y[i__2].r = z__1.r, y[i__2].i = z__1.i; 397 ix = jx; 398 iy = jy; 399 i__2 = kk + *n - j; 400 for (k = kk + 1; k <= i__2; ++k) { 401 ix += *incx; 402 iy += *incy; 403 i__3 = iy; 404 i__4 = iy; 405 i__5 = k; 406 z__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 407 z__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5] 408 .r; 409 z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i; 410 y[i__3].r = z__1.r, y[i__3].i = z__1.i; 411 d_cnjg(&z__3, &ap[k]); 412 i__3 = ix; 413 z__2.r = z__3.r * x[i__3].r - z__3.i * x[i__3].i, z__2.i = 414 z__3.r * x[i__3].i + z__3.i * x[i__3].r; 415 z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i; 416 temp2.r = z__1.r, temp2.i = z__1.i; 417/* L110: */ 418 } 419 i__2 = jy; 420 i__3 = jy; 421 z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i = 422 alpha->r * temp2.i + alpha->i * temp2.r; 423 z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i; 424 y[i__2].r = z__1.r, y[i__2].i = z__1.i; 425 jx += *incx; 426 jy += *incy; 427 kk += *n - j + 1; 428/* L120: */ 429 } 430 } 431 } 432 433 return 0; 434 435/* End of ZHPMV . */ 436 437} /* zhpmv_ */ 438 439