1 /* gmvm.c */
2
3 #include "../InpMtx.h"
4
5 /*--------------------------------------------------------------------*/
6 static int checkInput ( InpMtx *A, double beta[], int ny, double y[],
7 double alpha[], int nx, double x[], char *methodname ) ;
8 /*--------------------------------------------------------------------*/
9 /*
10 --------------------------------------------------
11 purpose -- to compute y := beta*y + alpha*A*x
12 where x and y are double vectors.
13
14 return values ---
15 1 -- normal return
16 -1 -- A is NULL
17 -2 -- type of A is invalid
18 -3 -- indices or entries of A are NULL
19 -4 -- beta is NULL
20 -5 -- ny <= 0
21 -6 -- y is NULL
22 -7 -- alpha is NULL
23 -8 -- nx <= 0
24 -9 -- x is NULL
25
26 created -- 98nov14, cca
27 --------------------------------------------------
28 */
29 int
InpMtx_nonsym_gmvm(InpMtx * A,double beta[],int ny,double y[],double alpha[],int nx,double x[])30 InpMtx_nonsym_gmvm (
31 InpMtx *A,
32 double beta[],
33 int ny,
34 double y[],
35 double alpha[],
36 int nx,
37 double x[]
38 ) {
39 int nent, rc ;
40 int *ivec1, *ivec2 ;
41 double *dvec ;
42 /*
43 ---------------
44 check the input
45 ---------------
46 */
47 rc = checkInput(A, beta, ny, y, alpha, nx, x, "InpMtx_nonsym_gmvm") ;
48 if ( rc != 1 ) {
49 return(rc) ;
50 }
51 ivec1 = InpMtx_ivec1(A) ;
52 ivec2 = InpMtx_ivec2(A) ;
53 dvec = InpMtx_dvec(A) ;
54 /*
55 ----------------
56 scale y by beta
57 ----------------
58 */
59 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
60 if ( beta[0] == 0.0 ) {
61 DVzero(ny, y) ;
62 } else if ( beta[0] != 0.0 ) {
63 DVscale(ny, y, beta[0]) ;
64 }
65 } else {
66 if ( beta[0] == 0.0 && beta[1] == 0.0 ) {
67 ZVzero(ny, y) ;
68 } else if ( beta[0] != 1.0 || beta[1] != 0.0 ) {
69 ZVscale(ny, y, beta[0], beta[1]) ;
70 }
71 }
72 /*
73 --------------------------------
74 data is stored as triples
75 (deal with vector storage later)
76 --------------------------------
77 */
78 nent = A->nent ;
79 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
80 double rfac ;
81 int chev, col, ii, off, row ;
82
83 rfac = alpha[0] ;
84 if ( INPMTX_IS_BY_ROWS(A) ) {
85 if ( rfac == 1.0 ) {
86 for ( ii = 0 ; ii < nent ; ii++ ) {
87 row = ivec1[ii] ; col = ivec2[ii] ;
88 y[row] += dvec[ii]*x[col] ;
89 }
90 } else {
91 for ( ii = 0 ; ii < nent ; ii++ ) {
92 row = ivec1[ii] ; col = ivec2[ii] ;
93 y[row] += rfac * dvec[ii]*x[col] ;
94 }
95 }
96 } else if ( INPMTX_IS_BY_COLUMNS(A) ) {
97 if ( rfac == 1.0 ) {
98 for ( ii = 0 ; ii < nent ; ii++ ) {
99 col = ivec1[ii] ; row = ivec2[ii] ;
100 y[row] += dvec[ii]*x[col] ;
101 }
102 } else {
103 for ( ii = 0 ; ii < nent ; ii++ ) {
104 col = ivec1[ii] ; row = ivec2[ii] ;
105 y[row] += rfac * dvec[ii]*x[col] ;
106 }
107 }
108 } else if ( INPMTX_IS_BY_CHEVRONS(A) ) {
109 if ( rfac == 1.0 ) {
110 for ( ii = 0 ; ii < nent ; ii++ ) {
111 chev = ivec1[ii] ; off = ivec2[ii] ;
112 if ( off >= 0 ) {
113 row = chev ; col = chev + off ;
114 } else {
115 col = chev ; row = chev - off ;
116 }
117 y[row] += dvec[ii]*x[col] ;
118 }
119 } else {
120 for ( ii = 0 ; ii < nent ; ii++ ) {
121 chev = ivec1[ii] ; off = ivec2[ii] ;
122 if ( off >= 0 ) {
123 row = chev ; col = chev + off ;
124 } else {
125 col = chev ; row = chev - off ;
126 }
127 y[row] += rfac * dvec[ii]*x[col] ;
128 }
129 }
130 }
131 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
132 double aimag, areal, ifac, rfac, t1, t2, ximag, xreal ;
133 int chev, col, ii, jj, off, row ;
134
135 rfac = alpha[0] ; ifac = alpha[1] ;
136 if ( INPMTX_IS_BY_ROWS(A) ) {
137 if ( rfac == 1.0 && ifac == 0.0 ) {
138 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
139 row = ivec1[ii] ; col = ivec2[ii] ;
140 areal = dvec[jj] ; aimag = dvec[jj+1] ;
141 xreal = x[2*col] ; ximag = x[2*col+1] ;
142 y[2*row] += areal*xreal - aimag*ximag ;
143 y[2*row+1] += areal*ximag + aimag*xreal ;
144 }
145 } else if ( ifac == 0.0 ) {
146 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
147 row = ivec1[ii] ; col = ivec2[ii] ;
148 areal = dvec[jj] ; aimag = dvec[jj+1] ;
149 xreal = x[2*col] ; ximag = x[2*col+1] ;
150 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
151 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
152 }
153 } else {
154 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
155 row = ivec1[ii] ; col = ivec2[ii] ;
156 areal = dvec[jj] ; aimag = dvec[jj+1] ;
157 xreal = x[2*col] ; ximag = x[2*col+1] ;
158 t1 = areal*xreal - aimag*ximag ;
159 t2 = areal*ximag + aimag*xreal ;
160 y[2*row] += rfac*t1 - ifac*t2 ;
161 y[2*row+1] += rfac*t2 + ifac*t1 ;
162 }
163 }
164 } else if ( INPMTX_IS_BY_COLUMNS(A) ) {
165 if ( rfac == 1.0 && ifac == 0.0 ) {
166 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
167 col = ivec1[ii] ; row = ivec2[ii] ;
168 areal = dvec[jj] ; aimag = dvec[jj+1] ;
169 xreal = x[2*col] ; ximag = x[2*col+1] ;
170 y[2*row] += areal*xreal - aimag*ximag ;
171 y[2*row+1] += areal*ximag + aimag*xreal ;
172 }
173 } else if ( ifac == 0.0 ) {
174 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
175 col = ivec1[ii] ; row = ivec2[ii] ;
176 areal = dvec[jj] ; aimag = dvec[jj+1] ;
177 xreal = x[2*col] ; ximag = x[2*col+1] ;
178 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
179 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
180 }
181 } else {
182 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
183 col = ivec1[ii] ; row = ivec2[ii] ;
184 areal = dvec[jj] ; aimag = dvec[jj+1] ;
185 xreal = x[2*col] ; ximag = x[2*col+1] ;
186 t1 = areal*xreal - aimag*ximag ;
187 t2 = areal*ximag + aimag*xreal ;
188 y[2*row] += rfac*t1 - ifac*t2 ;
189 y[2*row+1] += rfac*t2 + ifac*t1 ;
190 }
191 }
192 } else if ( INPMTX_IS_BY_CHEVRONS(A) ) {
193 if ( rfac == 1.0 && ifac == 0.0 ) {
194 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
195 chev = ivec1[ii] ; off = ivec2[ii] ;
196 if ( off >= 0 ) {
197 row = chev ; col = chev + off ;
198 } else {
199 col = chev ; row = chev - off ;
200 }
201 areal = dvec[jj] ; aimag = dvec[jj+1] ;
202 xreal = x[2*col] ; ximag = x[2*col+1] ;
203 y[2*row] += areal*xreal - aimag*ximag ;
204 y[2*row+1] += areal*ximag + aimag*xreal ;
205 }
206 } else if ( ifac == 0.0 ) {
207 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
208 chev = ivec1[ii] ; off = ivec2[ii] ;
209 if ( off >= 0 ) {
210 row = chev ; col = chev + off ;
211 } else {
212 col = chev ; row = chev - off ;
213 }
214 areal = dvec[jj] ; aimag = dvec[jj+1] ;
215 xreal = x[2*col] ; ximag = x[2*col+1] ;
216 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
217 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
218 }
219 } else {
220 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
221 chev = ivec1[ii] ; off = ivec2[ii] ;
222 if ( off >= 0 ) {
223 row = chev ; col = chev + off ;
224 } else {
225 col = chev ; row = chev - off ;
226 }
227 areal = dvec[jj] ; aimag = dvec[jj+1] ;
228 xreal = x[2*col] ; ximag = x[2*col+1] ;
229 t1 = areal*xreal - aimag*ximag ;
230 t2 = areal*ximag + aimag*xreal ;
231 y[2*row] += rfac*t1 - ifac*t2 ;
232 y[2*row+1] += rfac*t2 + ifac*t1 ;
233 }
234 }
235 }
236 }
237 return(1) ; }
238
239 /*--------------------------------------------------------------------*/
240 /*
241 --------------------------------------------------
242 purpose -- to compute Y := beta*Y + alpha*A^T*X
243
244 return values ---
245 1 -- normal return
246 -1 -- A is NULL
247 -2 -- type of A is invalid
248 -3 -- indices or entries of A are NULL
249 -4 -- beta is NULL
250 -5 -- ny <= 0
251 -6 -- y is NULL
252 -7 -- alpha is NULL
253 -8 -- nx <= 0
254 -9 -- x is NULL
255
256 created -- 98may02, cca
257 --------------------------------------------------
258 */
259 int
InpMtx_nonsym_gmvm_T(InpMtx * A,double beta[],int ny,double y[],double alpha[],int nx,double x[])260 InpMtx_nonsym_gmvm_T (
261 InpMtx *A,
262 double beta[],
263 int ny,
264 double y[],
265 double alpha[],
266 int nx,
267 double x[]
268 ) {
269 int nent, rc ;
270 int *ivec1, *ivec2 ;
271 double *dvec ;
272 /*
273 ---------------
274 check the input
275 ---------------
276 */
277 rc = checkInput(A, beta, ny, y, alpha, nx, x, "InpMtx_nonsym_gmvm_T") ;
278 if ( rc != 1 ) {
279 return(rc) ;
280 }
281 ivec1 = InpMtx_ivec1(A) ;
282 ivec2 = InpMtx_ivec2(A) ;
283 dvec = InpMtx_dvec(A) ;
284 /*
285 ----------------
286 scale y by beta
287 ----------------
288 */
289 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
290 if ( beta[0] == 0.0 ) {
291 DVzero(ny, y) ;
292 } else if ( beta[0] != 0.0 ) {
293 DVscale(ny, y, beta[0]) ;
294 }
295 } else {
296 if ( beta[0] == 0.0 && beta[1] == 0.0 ) {
297 ZVzero(ny, y) ;
298 } else if ( beta[0] != 1.0 || beta[1] != 0.0 ) {
299 ZVscale(ny, y, beta[0], beta[1]) ;
300 }
301 }
302 /*
303 --------------------------------
304 data is stored as triples
305 (deal with vector storage later)
306 --------------------------------
307 */
308 nent = A->nent ;
309 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
310 double rfac ;
311 int chev, col, ii, jrhs, off, row ;
312
313 rfac = alpha[0] ;
314 if ( INPMTX_IS_BY_ROWS(A) ) {
315 if ( rfac == 1.0 ) {
316 for ( ii = 0 ; ii < nent ; ii++ ) {
317 row = ivec2[ii] ; col = ivec1[ii] ;
318 y[row] += dvec[ii]*x[col] ;
319 }
320 } else {
321 for ( ii = 0 ; ii < nent ; ii++ ) {
322 row = ivec2[ii] ; col = ivec1[ii] ;
323 y[row] += rfac * dvec[ii]*x[col] ;
324 }
325 }
326 } else if ( INPMTX_IS_BY_COLUMNS(A) ) {
327 if ( rfac == 1.0 ) {
328 for ( ii = 0 ; ii < nent ; ii++ ) {
329 col = ivec2[ii] ; row = ivec1[ii] ;
330 y[row] += dvec[ii]*x[col] ;
331 }
332 } else {
333 for ( ii = 0 ; ii < nent ; ii++ ) {
334 col = ivec2[ii] ; row = ivec1[ii] ;
335 y[row] += rfac * dvec[ii]*x[col] ;
336 }
337 }
338 } else if ( INPMTX_IS_BY_CHEVRONS(A) ) {
339 if ( rfac == 1.0 ) {
340 for ( ii = 0 ; ii < nent ; ii++ ) {
341 chev = ivec1[ii] ; off = ivec2[ii] ;
342 if ( off >= 0 ) {
343 col = chev ; row = chev + off ;
344 } else {
345 row = chev ; col = chev - off ;
346 }
347 y[row] += dvec[ii]*x[col] ;
348 }
349 } else {
350 for ( ii = 0 ; ii < nent ; ii++ ) {
351 chev = ivec1[ii] ; off = ivec2[ii] ;
352 if ( off >= 0 ) {
353 col = chev ; row = chev + off ;
354 } else {
355 row = chev ; col = chev - off ;
356 }
357 y[row] += rfac * dvec[ii]*x[col] ;
358 }
359 }
360 }
361 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
362 double aimag, areal, ifac, rfac, t1, t2, ximag, xreal ;
363 int chev, col, ii, jj, jrhs, off, row ;
364
365 rfac = alpha[0] ; ifac = alpha[1] ;
366 if ( INPMTX_IS_BY_ROWS(A) ) {
367 if ( rfac == 1.0 && ifac == 0.0 ) {
368 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
369 row = ivec2[ii] ; col = ivec1[ii] ;
370 areal = dvec[jj] ; aimag = dvec[jj+1] ;
371 xreal = x[2*col] ; ximag = x[2*col+1] ;
372 y[2*row] += areal*xreal - aimag*ximag ;
373 y[2*row+1] += areal*ximag + aimag*xreal ;
374 }
375 } else if ( ifac == 0.0 ) {
376 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
377 row = ivec2[ii] ; col = ivec1[ii] ;
378 areal = dvec[jj] ; aimag = dvec[jj+1] ;
379 xreal = x[2*col] ; ximag = x[2*col+1] ;
380 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
381 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
382 }
383 } else {
384 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
385 row = ivec2[ii] ; col = ivec1[ii] ;
386 areal = dvec[jj] ; aimag = dvec[jj+1] ;
387 xreal = x[2*col] ; ximag = x[2*col+1] ;
388 t1 = areal*xreal - aimag*ximag ;
389 t2 = areal*ximag + aimag*xreal ;
390 y[2*row] += rfac*t1 - ifac*t2 ;
391 y[2*row+1] += rfac*t2 + ifac*t1 ;
392 }
393 }
394 } else if ( INPMTX_IS_BY_COLUMNS(A) ) {
395 if ( rfac == 1.0 && ifac == 0.0 ) {
396 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
397 col = ivec2[ii] ; row = ivec1[ii] ;
398 areal = dvec[jj] ; aimag = dvec[jj+1] ;
399 xreal = x[2*col] ; ximag = x[2*col+1] ;
400 y[2*row] += areal*xreal - aimag*ximag ;
401 y[2*row+1] += areal*ximag + aimag*xreal ;
402 }
403 } else if ( ifac == 0.0 ) {
404 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
405 col = ivec2[ii] ; row = ivec1[ii] ;
406 areal = dvec[jj] ; aimag = dvec[jj+1] ;
407 xreal = x[2*col] ; ximag = x[2*col+1] ;
408 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
409 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
410 }
411 } else {
412 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
413 col = ivec2[ii] ; row = ivec1[ii] ;
414 areal = dvec[jj] ; aimag = dvec[jj+1] ;
415 xreal = x[2*col] ; ximag = x[2*col+1] ;
416 t1 = areal*xreal - aimag*ximag ;
417 t2 = areal*ximag + aimag*xreal ;
418 y[2*row] += rfac*t1 - ifac*t2 ;
419 y[2*row+1] += rfac*t2 + ifac*t1 ;
420 }
421 }
422 } else if ( INPMTX_IS_BY_CHEVRONS(A) ) {
423 if ( rfac == 1.0 && ifac == 0.0 ) {
424 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
425 chev = ivec1[ii] ; off = ivec2[ii] ;
426 if ( off >= 0 ) {
427 col = chev ; row = chev + off ;
428 } else {
429 row = chev ; col = chev - off ;
430 }
431 areal = dvec[jj] ; aimag = dvec[jj+1] ;
432 xreal = x[2*col] ; ximag = x[2*col+1] ;
433 y[2*row] += areal*xreal - aimag*ximag ;
434 y[2*row+1] += areal*ximag + aimag*xreal ;
435 }
436 } else if ( ifac == 0.0 ) {
437 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
438 chev = ivec1[ii] ; off = ivec2[ii] ;
439 if ( off >= 0 ) {
440 col = chev ; row = chev + off ;
441 } else {
442 row = chev ; col = chev - off ;
443 }
444 areal = dvec[jj] ; aimag = dvec[jj+1] ;
445 xreal = x[2*col] ; ximag = x[2*col+1] ;
446 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
447 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
448 }
449 } else {
450 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
451 chev = ivec1[ii] ; off = ivec2[ii] ;
452 if ( off >= 0 ) {
453 col = chev ; row = chev + off ;
454 } else {
455 row = chev ; col = chev - off ;
456 }
457 areal = dvec[jj] ; aimag = dvec[jj+1] ;
458 xreal = x[2*col] ; ximag = x[2*col+1] ;
459 t1 = areal*xreal - aimag*ximag ;
460 t2 = areal*ximag + aimag*xreal ;
461 y[2*row] += rfac*t1 - ifac*t2 ;
462 y[2*row+1] += rfac*t2 + ifac*t1 ;
463 }
464 }
465 }
466 }
467 return(1) ; }
468
469 /*--------------------------------------------------------------------*/
470 /*
471 ---------------------------------------------------
472 purpose -- to compute Y := beta*Y + alpha*A^H*X
473
474 return values ---
475 1 -- normal return
476 -1 -- A is NULL
477 -2 -- type of A is invalid
478 -3 -- indices or entries of A are NULL
479 -4 -- beta is NULL
480 -5 -- ny <= 0
481 -6 -- y is NULL
482 -7 -- alpha is NULL
483 -8 -- nx <= 0
484 -9 -- x is NULL
485 -10 -- A is real
486
487 created -- 98may02, cca
488 ---------------------------------------------------
489 */
490 int
InpMtx_nonsym_gmvm_H(InpMtx * A,double beta[],int ny,double y[],double alpha[],int nx,double x[])491 InpMtx_nonsym_gmvm_H (
492 InpMtx *A,
493 double beta[],
494 int ny,
495 double y[],
496 double alpha[],
497 int nx,
498 double x[]
499 ) {
500 int nent, rc ;
501 int *ivec1, *ivec2 ;
502 double *dvec ;
503 /*
504 ---------------
505 check the input
506 ---------------
507 */
508 rc = checkInput(A, beta, ny, y, alpha, nx, x, "InpMtx_nonsym_gmvm_H") ;
509 if ( rc != 1 ) {
510 return(rc) ;
511 }
512 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
513 fprintf(stderr, "\n fatal error in InpMtx_nonsym_gmvm_H()"
514 "\n A, X and Y are real\n") ;
515 return(-10) ;
516 }
517 ivec1 = InpMtx_ivec1(A) ;
518 ivec2 = InpMtx_ivec2(A) ;
519 dvec = InpMtx_dvec(A) ;
520 /*
521 ----------------
522 scale y by beta
523 ----------------
524 */
525 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
526 if ( beta[0] == 0.0 ) {
527 DVzero(ny, y) ;
528 } else if ( beta[0] != 0.0 ) {
529 DVscale(ny, y, beta[0]) ;
530 }
531 } else {
532 if ( beta[0] == 0.0 && beta[1] == 0.0 ) {
533 ZVzero(ny, y) ;
534 } else if ( beta[0] != 1.0 || beta[1] != 0.0 ) {
535 ZVscale(ny, y, beta[0], beta[1]) ;
536 }
537 }
538 /*
539 --------------------------------
540 data is stored as triples
541 (deal with vector storage later)
542 --------------------------------
543 */
544 nent = A->nent ;
545 if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
546 double aimag, areal, ifac, rfac, t1, t2, ximag, xreal ;
547 int chev, col, ii, jj, jrhs, off, row ;
548
549 rfac = alpha[0] ; ifac = alpha[1] ;
550 if ( INPMTX_IS_BY_ROWS(A) ) {
551 if ( rfac == 1.0 && ifac == 0.0 ) {
552 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
553 row = ivec2[ii] ; col = ivec1[ii] ;
554 areal = dvec[jj] ; aimag = dvec[jj+1] ;
555 xreal = x[2*col] ; ximag = x[2*col+1] ;
556 y[2*row] += areal*xreal + aimag*ximag ;
557 y[2*row+1] += areal*ximag - aimag*xreal ;
558 }
559 } else if ( ifac == 0.0 ) {
560 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
561 row = ivec2[ii] ; col = ivec1[ii] ;
562 areal = dvec[jj] ; aimag = dvec[jj+1] ;
563 xreal = x[2*col] ; ximag = x[2*col+1] ;
564 y[2*row] += rfac*(areal*xreal + aimag*ximag) ;
565 y[2*row+1] += rfac*(areal*ximag - aimag*xreal) ;
566 }
567 } else {
568 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
569 row = ivec2[ii] ; col = ivec1[ii] ;
570 areal = dvec[jj] ; aimag = dvec[jj+1] ;
571 xreal = x[2*col] ; ximag = x[2*col+1] ;
572 t1 = areal*xreal + aimag*ximag ;
573 t2 = areal*ximag - aimag*xreal ;
574 y[2*row] += rfac*t1 - ifac*t2 ;
575 y[2*row+1] += rfac*t2 + ifac*t1 ;
576 }
577 }
578 } else if ( INPMTX_IS_BY_COLUMNS(A) ) {
579 if ( rfac == 1.0 && ifac == 0.0 ) {
580 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
581 col = ivec2[ii] ; row = ivec1[ii] ;
582 areal = dvec[jj] ; aimag = dvec[jj+1] ;
583 xreal = x[2*col] ; ximag = x[2*col+1] ;
584 y[2*row] += areal*xreal + aimag*ximag ;
585 y[2*row+1] += areal*ximag - aimag*xreal ;
586 }
587 } else if ( ifac == 0.0 ) {
588 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
589 col = ivec2[ii] ; row = ivec1[ii] ;
590 areal = dvec[jj] ; aimag = dvec[jj+1] ;
591 xreal = x[2*col] ; ximag = x[2*col+1] ;
592 y[2*row] += rfac*(areal*xreal + aimag*ximag) ;
593 y[2*row+1] += rfac*(areal*ximag - aimag*xreal) ;
594 }
595 } else {
596 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
597 col = ivec2[ii] ; row = ivec1[ii] ;
598 areal = dvec[jj] ; aimag = dvec[jj+1] ;
599 xreal = x[2*col] ; ximag = x[2*col+1] ;
600 t1 = areal*xreal + aimag*ximag ;
601 t2 = areal*ximag - aimag*xreal ;
602 y[2*row] += rfac*t1 - ifac*t2 ;
603 y[2*row+1] += rfac*t2 + ifac*t1 ;
604 }
605 }
606 } else if ( INPMTX_IS_BY_CHEVRONS(A) ) {
607 if ( rfac == 1.0 && ifac == 0.0 ) {
608 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
609 chev = ivec1[ii] ; off = ivec2[ii] ;
610 if ( off >= 0 ) {
611 col = chev ; row = chev + off ;
612 } else {
613 row = chev ; col = chev - off ;
614 }
615 areal = dvec[jj] ; aimag = dvec[jj+1] ;
616 xreal = x[2*col] ; ximag = x[2*col+1] ;
617 y[2*row] += areal*xreal + aimag*ximag ;
618 y[2*row+1] += areal*ximag - aimag*xreal ;
619 }
620 } else if ( ifac == 0.0 ) {
621 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
622 chev = ivec1[ii] ; off = ivec2[ii] ;
623 if ( off >= 0 ) {
624 col = chev ; row = chev + off ;
625 } else {
626 row = chev ; col = chev - off ;
627 }
628 areal = dvec[jj] ; aimag = dvec[jj+1] ;
629 xreal = x[2*col] ; ximag = x[2*col+1] ;
630 y[2*row] += rfac*(areal*xreal + aimag*ximag) ;
631 y[2*row+1] += rfac*(areal*ximag - aimag*xreal) ;
632 }
633 } else {
634 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
635 chev = ivec1[ii] ; off = ivec2[ii] ;
636 if ( off >= 0 ) {
637 col = chev ; row = chev + off ;
638 } else {
639 row = chev ; col = chev - off ;
640 }
641 areal = dvec[jj] ; aimag = dvec[jj+1] ;
642 xreal = x[2*col] ; ximag = x[2*col+1] ;
643 t1 = areal*xreal + aimag*ximag ;
644 t2 = areal*ximag - aimag*xreal ;
645 y[2*row] += rfac*t1 - ifac*t2 ;
646 y[2*row+1] += rfac*t2 + ifac*t1 ;
647 }
648 }
649 }
650 }
651 return(1) ; }
652
653 /*--------------------------------------------------------------------*/
654 /*
655 ----------------------------------------------------
656 purpose -- to compute Y := beta*Y + alpha*A*X
657 where A is symmetric
658
659 return values ---
660 1 -- normal return
661 -1 -- A is NULL
662 -2 -- type of A is invalid
663 -3 -- indices or entries of A are NULL
664 -4 -- beta is NULL
665 -5 -- ny <= 0
666 -6 -- y is NULL
667 -7 -- alpha is NULL
668 -8 -- nx <= 0
669 -9 -- x is NULL
670
671 created -- 98nov06, cca
672 ----------------------------------------------------
673 */
674 int
InpMtx_sym_gmvm(InpMtx * A,double beta[],int ny,double y[],double alpha[],int nx,double x[])675 InpMtx_sym_gmvm (
676 InpMtx *A,
677 double beta[],
678 int ny,
679 double y[],
680 double alpha[],
681 int nx,
682 double x[]
683 ) {
684 int nent, rc ;
685 int *ivec1, *ivec2 ;
686 double *dvec ;
687 /*
688 ---------------
689 check the input
690 ---------------
691 */
692 rc = checkInput(A, beta, ny, y, alpha, nx, x, "InpMtx_sym_gmvm") ;
693 if ( rc != 1 ) {
694 return(rc) ;
695 }
696 ivec1 = InpMtx_ivec1(A) ;
697 ivec2 = InpMtx_ivec2(A) ;
698 dvec = InpMtx_dvec(A) ;
699 /*
700 ----------------
701 scale y by beta
702 ----------------
703 */
704 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
705 if ( beta[0] == 0.0 ) {
706 DVzero(ny, y) ;
707 } else if ( beta[0] != 0.0 ) {
708 DVscale(ny, y, beta[0]) ;
709 }
710 } else {
711 if ( beta[0] == 0.0 && beta[1] == 0.0 ) {
712 ZVzero(ny, y) ;
713 } else if ( beta[0] != 1.0 || beta[1] != 0.0 ) {
714 ZVscale(ny, y, beta[0], beta[1]) ;
715 }
716 }
717 /*
718 --------------------------------
719 data is stored as triples
720 (deal with vector storage later)
721 --------------------------------
722 */
723 nent = A->nent ;
724 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
725 double rfac ;
726 int chev, col, ii, jrhs, off, row ;
727
728 rfac = alpha[0] ;
729 if ( INPMTX_IS_BY_ROWS(A) ) {
730 if ( rfac == 1.0 ) {
731 for ( ii = 0 ; ii < nent ; ii++ ) {
732 row = ivec1[ii] ; col = ivec2[ii] ;
733 y[row] += dvec[ii]*x[col] ;
734 if ( row != col ) {
735 y[col] += dvec[ii]*x[row] ;
736 }
737 }
738 } else {
739 for ( ii = 0 ; ii < nent ; ii++ ) {
740 row = ivec1[ii] ; col = ivec2[ii] ;
741 y[row] += rfac * dvec[ii]*x[col] ;
742 if ( row != col ) {
743 y[col] += rfac * dvec[ii]*x[row] ;
744 }
745 }
746 }
747 } else if ( INPMTX_IS_BY_COLUMNS(A) ) {
748 if ( rfac == 1.0 ) {
749 for ( ii = 0 ; ii < nent ; ii++ ) {
750 col = ivec1[ii] ; row = ivec2[ii] ;
751 y[row] += dvec[ii]*x[col] ;
752 if ( row != col ) {
753 y[col] += dvec[ii]*x[row] ;
754 }
755 }
756 } else {
757 for ( ii = 0 ; ii < nent ; ii++ ) {
758 col = ivec1[ii] ; row = ivec2[ii] ;
759 y[row] += rfac * dvec[ii]*x[col] ;
760 if ( row != col ) {
761 y[col] += rfac * dvec[ii]*x[row] ;
762 }
763 }
764 }
765 } else if ( INPMTX_IS_BY_CHEVRONS(A) ) {
766 if ( rfac == 1.0 ) {
767 for ( ii = 0 ; ii < nent ; ii++ ) {
768 chev = ivec1[ii] ; off = ivec2[ii] ;
769 if ( off >= 0 ) {
770 row = chev ; col = chev + off ;
771 } else {
772 col = chev ; row = chev - off ;
773 }
774 y[row] += dvec[ii]*x[col] ;
775 if ( row != col ) {
776 y[col] += dvec[ii]*x[row] ;
777 }
778 }
779 } else {
780 for ( ii = 0 ; ii < nent ; ii++ ) {
781 chev = ivec1[ii] ; off = ivec2[ii] ;
782 if ( off >= 0 ) {
783 row = chev ; col = chev + off ;
784 } else {
785 col = chev ; row = chev - off ;
786 }
787 y[row] += rfac * dvec[ii]*x[col] ;
788 if ( row != col ) {
789 y[col] += rfac * dvec[ii]*x[row] ;
790 }
791 }
792 }
793 }
794 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
795 double aimag, areal, ifac, rfac, t1, t2, ximag, xreal ;
796 int chev, col, ii, jj, jrhs, off, row ;
797
798 rfac = alpha[0] ; ifac = alpha[1] ;
799 if ( INPMTX_IS_BY_ROWS(A) ) {
800 if ( rfac == 1.0 && ifac == 0.0 ) {
801 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
802 row = ivec1[ii] ; col = ivec2[ii] ;
803 areal = dvec[jj] ; aimag = dvec[jj+1] ;
804 xreal = x[2*col] ; ximag = x[2*col+1] ;
805 y[2*row] += areal*xreal - aimag*ximag ;
806 y[2*row+1] += areal*ximag + aimag*xreal ;
807 if ( row != col ) {
808 xreal = x[2*row] ; ximag = x[2*row+1] ;
809 y[2*col] += areal*xreal - aimag*ximag ;
810 y[2*col+1] += areal*ximag + aimag*xreal ;
811 }
812 }
813 } else if ( ifac == 0.0 ) {
814 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
815 row = ivec1[ii] ; col = ivec2[ii] ;
816 areal = dvec[jj] ; aimag = dvec[jj+1] ;
817 xreal = x[2*col] ; ximag = x[2*col+1] ;
818 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
819 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
820 if ( row != col ) {
821 xreal = x[2*row] ; ximag = x[2*row+1] ;
822 y[2*col] += rfac*(areal*xreal - aimag*ximag) ;
823 y[2*col+1] += rfac*(areal*ximag + aimag*xreal) ;
824 }
825 }
826 } else {
827 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
828 row = ivec1[ii] ; col = ivec2[ii] ;
829 areal = dvec[jj] ; aimag = dvec[jj+1] ;
830 xreal = x[2*col] ; ximag = x[2*col+1] ;
831 t1 = areal*xreal - aimag*ximag ;
832 t2 = areal*ximag + aimag*xreal ;
833 y[2*row] += rfac*t1 - ifac*t2 ;
834 y[2*row+1] += rfac*t2 + ifac*t1 ;
835 if ( row != col ) {
836 xreal = x[2*row] ; ximag = x[2*row+1] ;
837 t1 = areal*xreal - aimag*ximag ;
838 t2 = areal*ximag + aimag*xreal ;
839 y[2*col] += rfac*t1 - ifac*t2 ;
840 y[2*col+1] += rfac*t2 + ifac*t1 ;
841 }
842 }
843 }
844 } else if ( INPMTX_IS_BY_COLUMNS(A) ) {
845 if ( rfac == 1.0 && ifac == 0.0 ) {
846 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
847 col = ivec1[ii] ; row = ivec2[ii] ;
848 areal = dvec[jj] ; aimag = dvec[jj+1] ;
849 xreal = x[2*col] ; ximag = x[2*col+1] ;
850 y[2*row] += areal*xreal - aimag*ximag ;
851 y[2*row+1] += areal*ximag + aimag*xreal ;
852 if ( row != col ) {
853 xreal = x[2*row] ; ximag = x[2*row+1] ;
854 y[2*col] += areal*xreal - aimag*ximag ;
855 y[2*col+1] += areal*ximag + aimag*xreal ;
856 }
857 }
858 } else if ( ifac == 0.0 ) {
859 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
860 col = ivec1[ii] ; row = ivec2[ii] ;
861 areal = dvec[jj] ; aimag = dvec[jj+1] ;
862 xreal = x[2*col] ; ximag = x[2*col+1] ;
863 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
864 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
865 if ( row != col ) {
866 xreal = x[2*row] ; ximag = x[2*row+1] ;
867 y[2*col] += rfac*(areal*xreal - aimag*ximag) ;
868 y[2*col+1] += rfac*(areal*ximag + aimag*xreal) ;
869 }
870 }
871 } else {
872 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
873 col = ivec1[ii] ; row = ivec2[ii] ;
874 areal = dvec[jj] ; aimag = dvec[jj+1] ;
875 xreal = x[2*col] ; ximag = x[2*col+1] ;
876 t1 = areal*xreal - aimag*ximag ;
877 t2 = areal*ximag + aimag*xreal ;
878 y[2*row] += rfac*t1 - ifac*t2 ;
879 y[2*row+1] += rfac*t2 + ifac*t1 ;
880 if ( row != col ) {
881 xreal = x[2*row] ; ximag = x[2*row+1] ;
882 t1 = areal*xreal - aimag*ximag ;
883 t2 = areal*ximag + aimag*xreal ;
884 y[2*col] += rfac*t1 - ifac*t2 ;
885 y[2*col+1] += rfac*t2 + ifac*t1 ;
886 }
887 }
888 }
889 } else if ( INPMTX_IS_BY_CHEVRONS(A) ) {
890 if ( rfac == 1.0 && ifac == 0.0 ) {
891 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
892 chev = ivec1[ii] ; off = ivec2[ii] ;
893 if ( off >= 0 ) {
894 row = chev ; col = chev + off ;
895 } else {
896 col = chev ; row = chev - off ;
897 }
898 areal = dvec[jj] ; aimag = dvec[jj+1] ;
899 xreal = x[2*col] ; ximag = x[2*col+1] ;
900 y[2*row] += areal*xreal - aimag*ximag ;
901 y[2*row+1] += areal*ximag + aimag*xreal ;
902 if ( row != col ) {
903 xreal = x[2*row] ; ximag = x[2*row+1] ;
904 y[2*col] += areal*xreal - aimag*ximag ;
905 y[2*col+1] += areal*ximag + aimag*xreal ;
906 }
907 }
908 } else if ( ifac == 0.0 ) {
909 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
910 chev = ivec1[ii] ; off = ivec2[ii] ;
911 if ( off >= 0 ) {
912 row = chev ; col = chev + off ;
913 } else {
914 col = chev ; row = chev - off ;
915 }
916 areal = dvec[jj] ; aimag = dvec[jj+1] ;
917 xreal = x[2*col] ; ximag = x[2*col+1] ;
918 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
919 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
920 if ( row != col ) {
921 xreal = x[2*row] ; ximag = x[2*row+1] ;
922 y[2*col] += rfac*(areal*xreal - aimag*ximag) ;
923 y[2*col+1] += rfac*(areal*ximag + aimag*xreal) ;
924 }
925 }
926 } else {
927 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
928 chev = ivec1[ii] ; off = ivec2[ii] ;
929 if ( off >= 0 ) {
930 row = chev ; col = chev + off ;
931 } else {
932 col = chev ; row = chev - off ;
933 }
934 areal = dvec[jj] ; aimag = dvec[jj+1] ;
935 xreal = x[2*col] ; ximag = x[2*col+1] ;
936 t1 = areal*xreal - aimag*ximag ;
937 t2 = areal*ximag + aimag*xreal ;
938 y[2*row] += rfac*t1 - ifac*t2 ;
939 y[2*row+1] += rfac*t2 + ifac*t1 ;
940 if ( row != col ) {
941 xreal = x[2*row] ; ximag = x[2*row+1] ;
942 t1 = areal*xreal - aimag*ximag ;
943 t2 = areal*ximag + aimag*xreal ;
944 y[2*col] += rfac*t1 - ifac*t2 ;
945 y[2*col+1] += rfac*t2 + ifac*t1 ;
946 }
947 }
948 }
949 }
950 }
951 return(1) ; }
952
953 /*--------------------------------------------------------------------*/
954 /*
955 --------------------------------------------------
956 purpose -- to compute Y := beta*Y + alpha*A*X
957 where A is hermitian
958
959 return values ---
960 1 -- normal return
961 -1 -- A is NULL
962 -2 -- type of A is invalid
963 -3 -- indices or entries of A are NULL
964 -4 -- beta is NULL
965 -5 -- ny <= 0
966 -6 -- y is NULL
967 -7 -- alpha is NULL
968 -8 -- nx <= 0
969 -9 -- x is NULL
970 -10 -- A, X and Y are real
971
972 created -- 98nov06, cca
973 --------------------------------------------------
974 */
975 int
InpMtx_herm_gmvm(InpMtx * A,double beta[],int ny,double y[],double alpha[],int nx,double x[])976 InpMtx_herm_gmvm (
977 InpMtx *A,
978 double beta[],
979 int ny,
980 double y[],
981 double alpha[],
982 int nx,
983 double x[]
984 ) {
985 int nent, rc ;
986 int *ivec1, *ivec2 ;
987 double *dvec ;
988 /*
989 ---------------
990 check the input
991 ---------------
992 */
993 rc = checkInput(A, beta, ny, y, alpha, nx, x, "InpMtx_herm_gmvm") ;
994 if ( rc != 1 ) {
995 return(rc) ;
996 }
997 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
998 fprintf(stderr, "\n fatal error in InpMtx_herm_gmvm()"
999 "\n A, X and Y are real\n") ;
1000 return(-10) ;
1001 }
1002 ivec1 = InpMtx_ivec1(A) ;
1003 ivec2 = InpMtx_ivec2(A) ;
1004 dvec = InpMtx_dvec(A) ;
1005 /*
1006 ----------------
1007 scale y by beta
1008 ----------------
1009 */
1010 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
1011 if ( beta[0] == 0.0 ) {
1012 DVzero(ny, y) ;
1013 } else if ( beta[0] != 0.0 ) {
1014 DVscale(ny, y, beta[0]) ;
1015 }
1016 } else {
1017 if ( beta[0] == 0.0 && beta[1] == 0.0 ) {
1018 ZVzero(ny, y) ;
1019 } else if ( beta[0] != 1.0 || beta[1] != 0.0 ) {
1020 ZVscale(ny, y, beta[0], beta[1]) ;
1021 }
1022 }
1023 /*
1024 --------------------------------
1025 data is stored as triples
1026 (deal with vector storage later)
1027 --------------------------------
1028 */
1029 nent = A->nent ;
1030 if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
1031 double aimag, areal, ifac, rfac, t1, t2, ximag, xreal ;
1032 int chev, col, ii, jj, jrhs, off, row ;
1033
1034 rfac = alpha[0] ; ifac = alpha[1] ;
1035 if ( INPMTX_IS_BY_ROWS(A) ) {
1036 if ( rfac == 1.0 && ifac == 0.0 ) {
1037 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
1038 row = ivec1[ii] ; col = ivec2[ii] ;
1039 areal = dvec[jj] ; aimag = dvec[jj+1] ;
1040 xreal = x[2*col] ; ximag = x[2*col+1] ;
1041 y[2*row] += areal*xreal - aimag*ximag ;
1042 y[2*row+1] += areal*ximag + aimag*xreal ;
1043 if ( row != col ) {
1044 xreal = x[2*row] ; ximag = x[2*row+1] ;
1045 y[2*col] += areal*xreal + aimag*ximag ;
1046 y[2*col+1] += areal*ximag - aimag*xreal ;
1047 }
1048 }
1049 } else if ( ifac == 0.0 ) {
1050 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
1051 row = ivec1[ii] ; col = ivec2[ii] ;
1052 areal = dvec[jj] ; aimag = dvec[jj+1] ;
1053 xreal = x[2*col] ; ximag = x[2*col+1] ;
1054 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
1055 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
1056 if ( row != col ) {
1057 xreal = x[2*row] ; ximag = x[2*row+1] ;
1058 y[2*col] += rfac*(areal*xreal + aimag*ximag) ;
1059 y[2*col+1] += rfac*(areal*ximag - aimag*xreal) ;
1060 }
1061 }
1062 } else {
1063 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
1064 row = ivec1[ii] ; col = ivec2[ii] ;
1065 areal = dvec[jj] ; aimag = dvec[jj+1] ;
1066 xreal = x[2*col] ; ximag = x[2*col+1] ;
1067 t1 = areal*xreal - aimag*ximag ;
1068 t2 = areal*ximag + aimag*xreal ;
1069 y[2*row] += rfac*t1 - ifac*t2 ;
1070 y[2*row+1] += rfac*t2 + ifac*t1 ;
1071 if ( row != col ) {
1072 xreal = x[2*row] ; ximag = x[2*row+1] ;
1073 t1 = areal*xreal + aimag*ximag ;
1074 t2 = areal*ximag - aimag*xreal ;
1075 y[2*col] += rfac*t1 - ifac*t2 ;
1076 y[2*col+1] += rfac*t2 + ifac*t1 ;
1077 }
1078 }
1079 }
1080 } else if ( INPMTX_IS_BY_COLUMNS(A) ) {
1081 if ( rfac == 1.0 && ifac == 0.0 ) {
1082 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
1083 col = ivec1[ii] ; row = ivec2[ii] ;
1084 areal = dvec[jj] ; aimag = dvec[jj+1] ;
1085 xreal = x[2*col] ; ximag = x[2*col+1] ;
1086 y[2*row] += areal*xreal - aimag*ximag ;
1087 y[2*row+1] += areal*ximag + aimag*xreal ;
1088 if ( row != col ) {
1089 xreal = x[2*row] ; ximag = x[2*row+1] ;
1090 y[2*col] += areal*xreal + aimag*ximag ;
1091 y[2*col+1] += areal*ximag - aimag*xreal ;
1092 }
1093 }
1094 } else if ( ifac == 0.0 ) {
1095 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
1096 col = ivec1[ii] ; row = ivec2[ii] ;
1097 areal = dvec[jj] ; aimag = dvec[jj+1] ;
1098 xreal = x[2*col] ; ximag = x[2*col+1] ;
1099 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
1100 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
1101 if ( row != col ) {
1102 xreal = x[2*row] ; ximag = x[2*row+1] ;
1103 y[2*col] += rfac*(areal*xreal + aimag*ximag) ;
1104 y[2*col+1] += rfac*(areal*ximag - aimag*xreal) ;
1105 }
1106 }
1107 } else {
1108 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
1109 col = ivec1[ii] ; row = ivec2[ii] ;
1110 areal = dvec[jj] ; aimag = dvec[jj+1] ;
1111 xreal = x[2*col] ; ximag = x[2*col+1] ;
1112 t1 = areal*xreal - aimag*ximag ;
1113 t2 = areal*ximag + aimag*xreal ;
1114 y[2*row] += rfac*t1 - ifac*t2 ;
1115 y[2*row+1] += rfac*t2 + ifac*t1 ;
1116 if ( row != col ) {
1117 xreal = x[2*row] ; ximag = x[2*row+1] ;
1118 t1 = areal*xreal + aimag*ximag ;
1119 t2 = areal*ximag - aimag*xreal ;
1120 y[2*col] += rfac*t1 - ifac*t2 ;
1121 y[2*col+1] += rfac*t2 + ifac*t1 ;
1122 }
1123 }
1124 }
1125 } else if ( INPMTX_IS_BY_CHEVRONS(A) ) {
1126 if ( rfac == 1.0 && ifac == 0.0 ) {
1127 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
1128 chev = ivec1[ii] ; off = ivec2[ii] ;
1129 if ( off >= 0 ) {
1130 row = chev ; col = chev + off ;
1131 } else {
1132 col = chev ; row = chev - off ;
1133 }
1134 areal = dvec[jj] ; aimag = dvec[jj+1] ;
1135 xreal = x[2*col] ; ximag = x[2*col+1] ;
1136 y[2*row] += areal*xreal - aimag*ximag ;
1137 y[2*row+1] += areal*ximag + aimag*xreal ;
1138 if ( row != col ) {
1139 xreal = x[2*row] ; ximag = x[2*row+1] ;
1140 y[2*col] += areal*xreal + aimag*ximag ;
1141 y[2*col+1] += areal*ximag - aimag*xreal ;
1142 }
1143 }
1144 } else if ( ifac == 0.0 ) {
1145 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
1146 chev = ivec1[ii] ; off = ivec2[ii] ;
1147 if ( off >= 0 ) {
1148 row = chev ; col = chev + off ;
1149 } else {
1150 col = chev ; row = chev - off ;
1151 }
1152 areal = dvec[jj] ; aimag = dvec[jj+1] ;
1153 xreal = x[2*col] ; ximag = x[2*col+1] ;
1154 y[2*row] += rfac*(areal*xreal - aimag*ximag) ;
1155 y[2*row+1] += rfac*(areal*ximag + aimag*xreal) ;
1156 if ( row != col ) {
1157 xreal = x[2*row] ; ximag = x[2*row+1] ;
1158 y[2*col] += rfac*(areal*xreal + aimag*ximag) ;
1159 y[2*col+1] += rfac*(areal*ximag - aimag*xreal) ;
1160 }
1161 }
1162 } else {
1163 for ( ii = jj = 0 ; ii < nent ; ii++, jj += 2 ) {
1164 chev = ivec1[ii] ; off = ivec2[ii] ;
1165 if ( off >= 0 ) {
1166 row = chev ; col = chev + off ;
1167 } else {
1168 col = chev ; row = chev - off ;
1169 }
1170 areal = dvec[jj] ; aimag = dvec[jj+1] ;
1171 xreal = x[2*col] ; ximag = x[2*col+1] ;
1172 t1 = areal*xreal - aimag*ximag ;
1173 t2 = areal*ximag + aimag*xreal ;
1174 y[2*row] += rfac*t1 - ifac*t2 ;
1175 y[2*row+1] += rfac*t2 + ifac*t1 ;
1176 if ( row != col ) {
1177 xreal = x[2*row] ; ximag = x[2*row+1] ;
1178 t1 = areal*xreal + aimag*ximag ;
1179 t2 = areal*ximag - aimag*xreal ;
1180 y[2*col] += rfac*t1 - ifac*t2 ;
1181 y[2*col+1] += rfac*t2 + ifac*t1 ;
1182 }
1183 }
1184 }
1185 }
1186 }
1187 return(1) ; }
1188
1189 /*--------------------------------------------------------------------*/
1190 /*
1191 --------------------------------------------------
1192 purpose -- to check the input
1193
1194 return values ---
1195 1 -- normal return
1196 -1 -- A is NULL
1197 -2 -- type of A is invalid
1198 -3 -- indices or entries of A are NULL
1199 -4 -- beta is NULL
1200 -5 -- ny <= 0
1201 -6 -- y is NULL
1202 -7 -- alpha is NULL
1203 -8 -- nx <= 0
1204 -9 -- x is NULL
1205
1206 created -- 98may02, cca
1207 --------------------------------------------------
1208 */
1209 static int
checkInput(InpMtx * A,double beta[],int ny,double y[],double alpha[],int nx,double x[],char * methodname)1210 checkInput (
1211 InpMtx *A,
1212 double beta[],
1213 int ny,
1214 double y[],
1215 double alpha[],
1216 int nx,
1217 double x[],
1218 char *methodname
1219 ) {
1220 double *dvec ;
1221 int typeX, typeY ;
1222 int *ivec1, *ivec2 ;
1223
1224 if ( A == NULL ) {
1225 fprintf(stderr, "\n fatal error in %s()"
1226 "\n A is NULL\n", methodname) ;
1227 return(-1) ;
1228 }
1229 if ( ! INPMTX_IS_REAL_ENTRIES(A) && ! INPMTX_IS_COMPLEX_ENTRIES(A) ) {
1230 fprintf(stderr, "\n fatal error in %s()"
1231 "\n type of A is %d, invalid\n", methodname, A->inputMode) ;
1232 return(-2) ;
1233 }
1234 ivec1 = InpMtx_ivec1(A) ;
1235 ivec2 = InpMtx_ivec2(A) ;
1236 dvec = InpMtx_dvec(A) ;
1237 if ( ivec1 == NULL || ivec2 == NULL || dvec == NULL ) {
1238 fprintf(stderr, "\n fatal error in %s()"
1239 "\n ivec1 %p, ivec2 %p, dvec %p\n",
1240 methodname, ivec1, ivec2, dvec) ;
1241 return(-3) ;
1242 }
1243 if ( beta == NULL ) {
1244 fprintf(stderr, "\n fatal error in %s()"
1245 "\n beta is NULL\n", methodname) ;
1246 return(-4) ;
1247 }
1248 if ( ny <= 0 ) {
1249 fprintf(stderr, "\n fatal error in %s()"
1250 "\n ny = %d\n", methodname, ny) ;
1251 return(-5) ;
1252 }
1253 if ( y == NULL ) {
1254 fprintf(stderr, "\n fatal error in %s()"
1255 "\n y is NULL\n", methodname) ;
1256 return(-6) ;
1257 }
1258 if ( alpha == NULL ) {
1259 fprintf(stderr, "\n fatal error in %s()"
1260 "\n alpha is NULL\n", methodname) ;
1261 return(-7) ;
1262 }
1263 if ( nx <= 0 ) {
1264 fprintf(stderr, "\n fatal error in %s()"
1265 "\n nx = %d\n", methodname, nx) ;
1266 return(-8) ;
1267 }
1268 if ( x == NULL ) {
1269 fprintf(stderr, "\n fatal error in %s()"
1270 "\n x is NULL\n", methodname) ;
1271 return(-9) ;
1272 }
1273 return(1) ; }
1274
1275 /*--------------------------------------------------------------------*/
1276