1 /* chpmv.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 chpmv_(char *uplo, integer *n, complex *alpha, complex * 16 ap, complex *x, integer *incx, complex *beta, complex *y, integer * 17 incy, ftnlen uplo_len) 18 { 19 /* System generated locals */ 20 integer i__1, i__2, i__3, i__4, i__5; 21 real r__1; 22 complex q__1, q__2, q__3, q__4; 23 24 /* Builtin functions */ 25 void r_cnjg(complex *, complex *); 26 27 /* Local variables */ 28 integer i__, j, k, kk, ix, iy, jx, jy, kx, ky, info; 29 complex 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 /* CHPMV 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 . */ 70 /* On entry, ALPHA specifies the scalar alpha. */ 71 /* Unchanged on exit. */ 72 73 /* AP - COMPLEX 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 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 . */ 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 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_("CHPMV ", &info, (ftnlen)6); 161 return 0; 162 } 163 164 /* Quick return if possible. */ 165 166 if (*n == 0 || (alpha->r == 0.f && alpha->i == 0.f && (beta->r == 1.f && 167 beta->i == 0.f))) { 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.f || beta->i != 0.f) { 190 if (*incy == 1) { 191 if (beta->r == 0.f && beta->i == 0.f) { 192 i__1 = *n; 193 for (i__ = 1; i__ <= i__1; ++i__) { 194 i__2 = i__; 195 y[i__2].r = 0.f, y[i__2].i = 0.f; 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 q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 204 q__1.i = beta->r * y[i__3].i + beta->i * y[i__3] 205 .r; 206 y[i__2].r = q__1.r, y[i__2].i = q__1.i; 207 /* L20: */ 208 } 209 } 210 } else { 211 iy = ky; 212 if (beta->r == 0.f && beta->i == 0.f) { 213 i__1 = *n; 214 for (i__ = 1; i__ <= i__1; ++i__) { 215 i__2 = iy; 216 y[i__2].r = 0.f, y[i__2].i = 0.f; 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 q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i, 226 q__1.i = beta->r * y[i__3].i + beta->i * y[i__3] 227 .r; 228 y[i__2].r = q__1.r, y[i__2].i = q__1.i; 229 iy += *incy; 230 /* L40: */ 231 } 232 } 233 } 234 } 235 if (alpha->r == 0.f && alpha->i == 0.f) { 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 q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i = 248 alpha->r * x[i__2].i + alpha->i * x[i__2].r; 249 temp1.r = q__1.r, temp1.i = q__1.i; 250 temp2.r = 0.f, temp2.i = 0.f; 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 q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 258 q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5] 259 .r; 260 q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i; 261 y[i__3].r = q__1.r, y[i__3].i = q__1.i; 262 r_cnjg(&q__3, &ap[k]); 263 i__3 = i__; 264 q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i = 265 q__3.r * x[i__3].i + q__3.i * x[i__3].r; 266 q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i; 267 temp2.r = q__1.r, temp2.i = q__1.i; 268 ++k; 269 /* L50: */ 270 } 271 i__2 = j; 272 i__3 = j; 273 i__4 = kk + j - 1; 274 r__1 = ap[i__4].r; 275 q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i; 276 q__2.r = y[i__3].r + q__3.r, q__2.i = y[i__3].i + q__3.i; 277 q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i = 278 alpha->r * temp2.i + alpha->i * temp2.r; 279 q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i; 280 y[i__2].r = q__1.r, y[i__2].i = q__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 q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i = 291 alpha->r * x[i__2].i + alpha->i * x[i__2].r; 292 temp1.r = q__1.r, temp1.i = q__1.i; 293 temp2.r = 0.f, temp2.i = 0.f; 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 q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 302 q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5] 303 .r; 304 q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i; 305 y[i__3].r = q__1.r, y[i__3].i = q__1.i; 306 r_cnjg(&q__3, &ap[k]); 307 i__3 = ix; 308 q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i = 309 q__3.r * x[i__3].i + q__3.i * x[i__3].r; 310 q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i; 311 temp2.r = q__1.r, temp2.i = q__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 r__1 = ap[i__4].r; 320 q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i; 321 q__2.r = y[i__3].r + q__3.r, q__2.i = y[i__3].i + q__3.i; 322 q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i = 323 alpha->r * temp2.i + alpha->i * temp2.r; 324 q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i; 325 y[i__2].r = q__1.r, y[i__2].i = q__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 q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i = 341 alpha->r * x[i__2].i + alpha->i * x[i__2].r; 342 temp1.r = q__1.r, temp1.i = q__1.i; 343 temp2.r = 0.f, temp2.i = 0.f; 344 i__2 = j; 345 i__3 = j; 346 i__4 = kk; 347 r__1 = ap[i__4].r; 348 q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i; 349 q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i; 350 y[i__2].r = q__1.r, y[i__2].i = q__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 q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 358 q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5] 359 .r; 360 q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i; 361 y[i__3].r = q__1.r, y[i__3].i = q__1.i; 362 r_cnjg(&q__3, &ap[k]); 363 i__3 = i__; 364 q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i = 365 q__3.r * x[i__3].i + q__3.i * x[i__3].r; 366 q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i; 367 temp2.r = q__1.r, temp2.i = q__1.i; 368 ++k; 369 /* L90: */ 370 } 371 i__2 = j; 372 i__3 = j; 373 q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i = 374 alpha->r * temp2.i + alpha->i * temp2.r; 375 q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i; 376 y[i__2].r = q__1.r, y[i__2].i = q__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 q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i = 387 alpha->r * x[i__2].i + alpha->i * x[i__2].r; 388 temp1.r = q__1.r, temp1.i = q__1.i; 389 temp2.r = 0.f, temp2.i = 0.f; 390 i__2 = jy; 391 i__3 = jy; 392 i__4 = kk; 393 r__1 = ap[i__4].r; 394 q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i; 395 q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i; 396 y[i__2].r = q__1.r, y[i__2].i = q__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 q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i, 407 q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5] 408 .r; 409 q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i; 410 y[i__3].r = q__1.r, y[i__3].i = q__1.i; 411 r_cnjg(&q__3, &ap[k]); 412 i__3 = ix; 413 q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i = 414 q__3.r * x[i__3].i + q__3.i * x[i__3].r; 415 q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i; 416 temp2.r = q__1.r, temp2.i = q__1.i; 417 /* L110: */ 418 } 419 i__2 = jy; 420 i__3 = jy; 421 q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i = 422 alpha->r * temp2.i + alpha->i * temp2.r; 423 q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i; 424 y[i__2].r = q__1.r, y[i__2].i = q__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 CHPMV . */ 436 437 } /* chpmv_ */ 438 439
