1 /* search.c */
2
3 #include "../Chv.h"
4
5 #define MYDEBUG 0
6
7 /*--------------------------------------------------------------------*/
8 /*
9 -----------------------------------
10 find the first unmarked entry in
11 the diagonal with largest magnitude
12 if ( mark[jj] == tag ) then
13 we can compare this entry
14 endif
15
16 created -- 98apr30, cca
17 -----------------------------------
18 */
19 int
Chv_maxabsInDiagonal11(Chv * chv,int mark[],int tag,double * pmaxval)20 Chv_maxabsInDiagonal11 (
21 Chv *chv,
22 int mark[],
23 int tag,
24 double *pmaxval
25 ) {
26 double maxval, val ;
27 double *entries ;
28 int jcol, jj, nD, nL, nU, off, stride ;
29 /*
30 --------------
31 check the data
32 --------------
33 */
34 if ( chv == NULL || mark == NULL || pmaxval == NULL ) {
35 fprintf(stderr,
36 "\n fatal error in Chv_maxabsInDiagonal11(%p,%p,%d,%p)"
37 "\n bad input\n", chv, mark, tag, pmaxval) ;
38 exit(-1) ;
39 }
40 Chv_dimensions(chv, &nD, &nL, &nU) ;
41 entries = Chv_entries(chv) ;
42 if ( CHV_IS_REAL(chv) ) {
43 if ( CHV_IS_NONSYMMETRIC(chv) ) {
44 /*
45 --------------------
46 nonsymmetric chevron
47 --------------------
48 */
49 jcol = -1 ;
50 maxval = 0.0 ;
51 off = nD + nL - 1 ;
52 stride = 2*nD + nL + nU - 2 ;
53 for ( jj = 0 ; jj < nD ; jj++ ) {
54 if ( mark[jj] == tag ) {
55 val = fabs(entries[off]) ;
56 if ( jcol == -1 || maxval < val ) {
57 jcol = jj ; maxval = val ;
58 }
59 }
60 off += stride ;
61 stride -= 2 ;
62 }
63 } else if ( CHV_IS_SYMMETRIC(chv) ) {
64 /*
65 -----------------
66 symmetric chevron
67 -----------------
68 */
69 jcol = -1 ;
70 maxval = 0.0 ;
71 off = 0 ;
72 stride = nD + nU ;
73 for ( jj = 0 ; jj < nD ; jj++ ) {
74 if ( mark[jj] == tag ) {
75 val = fabs(entries[off]) ;
76 if ( jcol == -1 || maxval < val ) {
77 jcol = jj ; maxval = val ;
78 }
79 }
80 off += stride ;
81 stride-- ;
82 }
83 } else {
84 fprintf(stderr,
85 "\n fatal error in Chv_maxabsInDiagonal11(%p,%p,%d,%p)"
86 "\n type = SPOOLES_REAL, bad symflag %d \n",
87 chv, mark, tag, pmaxval, chv->symflag) ;
88 exit(-1) ;
89 }
90 } else if ( CHV_IS_COMPLEX(chv) ) {
91 if ( CHV_IS_NONSYMMETRIC(chv) ) {
92 /*
93 --------------------
94 nonsymmetric chevron
95 --------------------
96 */
97 jcol = -1 ;
98 maxval = 0.0 ;
99 off = nD + nL - 1 ;
100 stride = 2*nD + nL + nU - 2 ;
101 for ( jj = 0 ; jj < nD ; jj++ ) {
102 if ( mark[jj] == tag ) {
103 val = Zabs(entries[2*off], entries[2*off+1]) ;
104 if ( jcol == -1 || maxval < val ) {
105 jcol = jj ; maxval = val ;
106 }
107 }
108 off += stride ;
109 stride -= 2 ;
110 }
111 } else if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
112 /*
113 ------------------------------
114 hermitian or symmetric chevron
115 ------------------------------
116 */
117 jcol = -1 ;
118 maxval = 0.0 ;
119 off = 0 ;
120 stride = nD + nU ;
121 for ( jj = 0 ; jj < nD ; jj++ ) {
122 if ( mark[jj] == tag ) {
123 val = Zabs(entries[2*off], entries[2*off+1]) ;
124 if ( jcol == -1 || maxval < val ) {
125 jcol = jj ; maxval = val ;
126 }
127 }
128 off += stride ;
129 stride-- ;
130 }
131 } else {
132 fprintf(stderr,
133 "\n fatal error in Chv_maxabsInDiagonal11(%p,%p,%d,%p)"
134 "\n type = SPOOLES_COMPLEX, bad symflag %d \n",
135 chv, mark, tag, pmaxval, chv->symflag) ;
136 exit(-1) ;
137 }
138 } else {
139 fprintf(stderr,
140 "\n fatal error in Chv_maxabsInDiagonal11(%p,%p,%d,%p)"
141 "\n bad type, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
142 chv, mark, tag, pmaxval) ;
143 exit(-1) ;
144 }
145 *pmaxval = maxval ;
146
147 return(jcol) ; }
148
149 /*--------------------------------------------------------------------*/
150 /*
151 --------------------------------------------
152 find the first unmarked entry in
153 row irow with largest magnitude
154 if ( colmark[jj] == tag ) then
155 we can examined this entry
156 endif
157 only entries in the (1,1) block are examined
158
159 created -- 98apr30, cca
160 --------------------------------------------
161 */
162 int
Chv_maxabsInRow11(Chv * chv,int irow,int colmark[],int tag,double * pmaxval)163 Chv_maxabsInRow11 (
164 Chv *chv,
165 int irow,
166 int colmark[],
167 int tag,
168 double *pmaxval
169 ) {
170 double maxval, val ;
171 double *entries ;
172 int jcol, jj, nD, nL, nU, off, stride ;
173 /*
174 --------------
175 check the data
176 --------------
177 */
178 if ( chv == NULL || irow < 0 || colmark == NULL || pmaxval == NULL ) {
179 fprintf(stderr,
180 "\n fatal error in Chv_maxabsInRow11(%p,%d,%p,%d,%p)"
181 "\n bad input\n", chv, irow, colmark, tag, pmaxval) ;
182 exit(-1) ;
183 }
184 Chv_dimensions(chv, &nD, &nL, &nU) ;
185 entries = Chv_entries(chv) ;
186 if ( CHV_IS_REAL(chv) ) {
187 if ( CHV_IS_NONSYMMETRIC(chv) ) {
188 /*
189 --------------------
190 nonsymmetric chevron
191 --------------------
192 */
193 jcol = -1 ;
194 maxval = 0.0 ;
195 off = nD + nL - 1 - irow ;
196 stride = 2*nD + nL + nU - 1 ;
197 for ( jj = 0 ; jj < irow ; jj++ ) {
198 if ( colmark[jj] == tag ) {
199 val = fabs(entries[off]) ;
200 if ( jcol == -1 || maxval < val ) {
201 jcol = jj ; maxval = val ;
202 }
203 }
204 off += stride ;
205 stride -= 2 ;
206 }
207 for ( jj = irow ; jj < nD ; jj++, off++ ) {
208 if ( colmark[jj] == tag ) {
209 val = fabs(entries[off]) ;
210 if ( jcol == -1 || maxval < val ) {
211 jcol = jj ; maxval = val ;
212 }
213 }
214 }
215 } else if ( CHV_IS_SYMMETRIC(chv) ) {
216 /*
217 -----------------
218 symmetric chevron
219 -----------------
220 */
221 jcol = -1 ;
222 maxval = 0.0 ;
223 off = irow ;
224 stride = nD + nU - 1 ;
225 for ( jj = 0 ; jj < irow ; jj++ ) {
226 if ( colmark[jj] == tag ) {
227 val = fabs(entries[off]) ;
228 if ( jcol == -1 || maxval < val ) {
229 jcol = jj ; maxval = val ;
230 }
231 }
232 off += stride ;
233 stride-- ;
234 }
235 for ( jj = irow ; jj < nD ; jj++, off++ ) {
236 if ( colmark[jj] == tag ) {
237 val = fabs(entries[off]) ;
238 if ( jcol == -1 || maxval < val ) {
239 jcol = jj ; maxval = val ;
240 }
241 }
242 }
243 } else {
244 fprintf(stderr,
245 "\n fatal error in Chv_maxabsInRow11(%p,%d,%p,%d,%p)"
246 "\n type is SPOOLES_REAL, bad symflag %d \n",
247 chv, irow, colmark, tag, pmaxval, chv->symflag) ;
248 exit(-1) ;
249 }
250 } else if ( CHV_IS_COMPLEX(chv) ) {
251 if ( CHV_IS_NONSYMMETRIC(chv) ) {
252 /*
253 --------------------
254 nonsymmetric chevron
255 --------------------
256 */
257 jcol = -1 ;
258 maxval = 0.0 ;
259 off = nD + nL - 1 - irow ;
260 stride = 2*nD + nL + nU - 1 ;
261 for ( jj = 0 ; jj < irow ; jj++ ) {
262 if ( colmark[jj] == tag ) {
263 val = Zabs(entries[2*off], entries[2*off+1]) ;
264 if ( jcol == -1 || maxval < val ) {
265 jcol = jj ; maxval = val ;
266 }
267 }
268 off += stride ;
269 stride -= 2 ;
270 }
271 for ( jj = irow ; jj < nD ; jj++, off++ ) {
272 if ( colmark[jj] == tag ) {
273 val = Zabs(entries[2*off], entries[2*off+1]) ;
274 if ( jcol == -1 || maxval < val ) {
275 jcol = jj ; maxval = val ;
276 }
277 }
278 }
279 } else if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
280 /*
281 ------------------------------
282 hermitian or symmetric chevron
283 ------------------------------
284 */
285 jcol = -1 ;
286 maxval = 0.0 ;
287 off = irow ;
288 stride = nD + nU - 1 ;
289 for ( jj = 0 ; jj < irow ; jj++ ) {
290 if ( colmark[jj] == tag ) {
291 val = Zabs(entries[2*off], entries[2*off+1]) ;
292 if ( jcol == -1 || maxval < val ) {
293 jcol = jj ; maxval = val ;
294 }
295 }
296 off += stride ;
297 stride-- ;
298 }
299 for ( jj = irow ; jj < nD ; jj++, off++ ) {
300 if ( colmark[jj] == tag ) {
301 val = Zabs(entries[2*off], entries[2*off+1]) ;
302 if ( jcol == -1 || maxval < val ) {
303 jcol = jj ; maxval = val ;
304 }
305 }
306 }
307 } else {
308 fprintf(stderr,
309 "\n fatal error in Chv_maxabsInRow11(%p,%d,%p,%d,%p)"
310 "\n type is SPOOLES_COMPLEX, bad symflag %d \n",
311 chv, irow, colmark, tag, pmaxval, chv->symflag) ;
312 exit(-1) ;
313 }
314 } else {
315 fprintf(stderr,
316 "\n fatal error in Chv_maxabsInRow11(%p,%d,%p,%d,%p)"
317 "\n bad type, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
318 chv, irow, colmark, tag, pmaxval) ;
319 exit(-1) ;
320 }
321 *pmaxval = maxval ;
322
323 return(jcol) ; }
324
325 /*--------------------------------------------------------------------*/
326 /*
327 --------------------------------------------
328 find the first unmarked entry in
329 column jcol with largest magnitude
330 if ( rowmark[ii] == tag ) then
331 we can examined this entry
332 endif
333 only entries in the (1,1) block are examined
334
335 created -- 98apr30, cca
336 --------------------------------------------
337 */
338 int
Chv_maxabsInColumn11(Chv * chv,int jcol,int rowmark[],int tag,double * pmaxval)339 Chv_maxabsInColumn11 (
340 Chv *chv,
341 int jcol,
342 int rowmark[],
343 int tag,
344 double *pmaxval
345 ) {
346 double maxval, val ;
347 double *entries ;
348 int irow, ii, nD, nL, nU, off, stride ;
349 /*
350 --------------
351 check the data
352 --------------
353 */
354 if ( chv == NULL || jcol < 0 || rowmark == NULL || pmaxval == NULL ) {
355 fprintf(stderr,
356 "\n fatal error in Chv_maxabsInColumn11(%p,%d,%p,%d,%p)"
357 "\n bad input\n", chv, jcol, rowmark, tag, pmaxval) ;
358 exit(-1) ;
359 }
360 Chv_dimensions(chv, &nD, &nL, &nU) ;
361 entries = Chv_entries(chv) ;
362 irow = -1 ;
363 maxval = 0.0 ;
364 if ( CHV_IS_REAL(chv) ) {
365 if ( CHV_IS_NONSYMMETRIC(chv) ) {
366 /*
367 --------------------
368 nonsymmetric chevron
369 --------------------
370 */
371 maxval = 0.0 ;
372 off = nD + nL + jcol - 1 ;
373 stride = 2*nD + nL + nU - 3 ;
374 for ( ii = 0 ; ii < jcol ; ii++ ) {
375 if ( rowmark[ii] == tag ) {
376 val = fabs(entries[off]) ;
377 if ( irow == -1 || maxval < val ) {
378 irow = ii ; maxval = val ;
379 }
380 }
381 off += stride ;
382 stride -= 2 ;
383 }
384 for ( ii = jcol ; ii < nD ; ii++, off-- ) {
385 if ( rowmark[ii] == tag ) {
386 val = fabs(entries[off]) ;
387 if ( irow == -1 || maxval < val ) {
388 irow = ii ; maxval = val ;
389 }
390 }
391 }
392 } else if ( CHV_IS_SYMMETRIC(chv) ) {
393 /*
394 -----------------
395 symmetric chevron
396 -----------------
397 */
398 maxval = 0.0 ;
399 off = jcol ;
400 stride = nD + nU - 1 ;
401 for ( ii = 0 ; ii < jcol ; ii++ ) {
402 if ( rowmark[ii] == tag ) {
403 val = fabs(entries[off]) ;
404 if ( irow == -1 || maxval < val ) {
405 irow = ii ; maxval = val ;
406 }
407 }
408 off += stride ;
409 stride-- ;
410 }
411 for ( ii = jcol ; ii < nD ; ii++, off++ ) {
412 if ( rowmark[ii] == tag ) {
413 val = fabs(entries[off]) ;
414 if ( irow == -1 || maxval < val ) {
415 irow = ii ; maxval = val ;
416 }
417 }
418 }
419 }
420 } else if ( CHV_IS_COMPLEX(chv) ) {
421 if ( CHV_IS_NONSYMMETRIC(chv) ) {
422 /*
423 --------------------
424 nonsymmetric chevron
425 --------------------
426 */
427 maxval = 0.0 ;
428 off = nD + nL + jcol - 1 ;
429 stride = 2*nD + nL + nU - 3 ;
430 for ( ii = 0 ; ii < jcol ; ii++ ) {
431 if ( rowmark[ii] == tag ) {
432 val = Zabs(entries[2*off], entries[2*off+1]) ;
433 if ( irow == -1 || maxval < val ) {
434 irow = ii ; maxval = val ;
435 }
436 }
437 off += stride ;
438 stride -= 2 ;
439 }
440 for ( ii = jcol ; ii < nD ; ii++, off-- ) {
441 if ( rowmark[ii] == tag ) {
442 val = Zabs(entries[2*off], entries[2*off+1]) ;
443 if ( irow == -1 || maxval < val ) {
444 irow = ii ; maxval = val ;
445 }
446 }
447 }
448 } else if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
449 /*
450 ------------------------------
451 hermitian or symmetric chevron
452 ------------------------------
453 */
454 maxval = 0.0 ;
455 off = jcol ;
456 stride = nD + nU - 1 ;
457 for ( ii = 0 ; ii < jcol ; ii++ ) {
458 if ( rowmark[ii] == tag ) {
459 val = Zabs(entries[2*off], entries[2*off+1]) ;
460 if ( irow == -1 || maxval < val ) {
461 irow = ii ; maxval = val ;
462 }
463 }
464 off += stride ;
465 stride-- ;
466 }
467 for ( ii = jcol ; ii < nD ; ii++, off++ ) {
468 if ( rowmark[ii] == tag ) {
469 val = Zabs(entries[2*off], entries[2*off+1]) ;
470 if ( irow == -1 || maxval < val ) {
471 irow = ii ; maxval = val ;
472 }
473 }
474 }
475 }
476 } else {
477 fprintf(stderr,
478 "\n fatal error in Chv_maxabsInColumn11(%p,%d,%p,%d,%p)"
479 "\n bad type, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
480 chv, jcol, rowmark, tag, pmaxval) ;
481 exit(-1) ;
482 }
483 *pmaxval = maxval ;
484
485 return(irow) ; }
486
487 /*--------------------------------------------------------------------*/
488 /*
489 --------------------------------------
490 return the location of the first entry
491 with largest magnitude in row irow.
492 *pmaxval is filled with its magnitude.
493
494 created -- 98apr30, cca
495 --------------------------------------
496 */
497 int
Chv_maxabsInRow(Chv * chv,int irow,double * pmaxval)498 Chv_maxabsInRow (
499 Chv *chv,
500 int irow,
501 double *pmaxval
502 ) {
503 double maxval, val ;
504 double *entries ;
505 int jcol, jj, ncol, nD, nL, nU, off, stride ;
506 /*
507 --------------
508 check the data
509 --------------
510 */
511 if ( chv == NULL || irow < 0 || pmaxval == NULL ) {
512 fprintf(stderr, "\n fatal error in Chv_maxabsInRow(%p,%d,%p)"
513 "\n bad input\n", chv, irow, pmaxval) ;
514 exit(-1) ;
515 }
516 Chv_dimensions(chv, &nD, &nL, &nU) ;
517 entries = Chv_entries(chv) ;
518 ncol = nD + nU ;
519 jcol = -1 ;
520 maxval = 0.0 ;
521 if ( CHV_IS_REAL(chv) ) {
522 if ( CHV_IS_NONSYMMETRIC(chv) ) {
523 /*
524 --------------------
525 nonsymmetric chevron
526 --------------------
527 */
528 jcol = -1 ;
529 maxval = 0.0 ;
530 off = nD + nL - 1 - irow ;
531 stride = 2*nD + nL + nU - 1 ;
532 for ( jj = 0 ; jj < irow ; jj++ ) {
533 val = fabs(entries[off]) ;
534 if ( jcol == -1 || maxval < val ) {
535 jcol = jj ; maxval = val ;
536 }
537 off += stride ;
538 stride -= 2 ;
539 }
540 for ( jj = irow ; jj < ncol ; jj++, off++ ) {
541 val = fabs(entries[off]) ;
542 if ( jcol == -1 || maxval < val ) {
543 jcol = jj ; maxval = val ;
544 }
545 }
546 } else if ( CHV_IS_SYMMETRIC(chv) ) {
547 /*
548 -----------------
549 symmetric chevron
550 -----------------
551 */
552 jcol = -1 ;
553 maxval = 0.0 ;
554 off = irow ;
555 stride = nD + nU - 1 ;
556 for ( jj = 0 ; jj < irow ; jj++ ) {
557 val = fabs(entries[off]) ;
558 if ( jcol == -1 || maxval < val ) {
559 jcol = jj ; maxval = val ;
560 }
561 off += stride ;
562 stride-- ;
563 }
564 for ( jj = irow ; jj < ncol ; jj++, off++ ) {
565 val = fabs(entries[off]) ;
566 if ( jcol == -1 || maxval < val ) {
567 jcol = jj ; maxval = val ;
568 }
569 }
570 }
571 } else if ( CHV_IS_COMPLEX(chv) ) {
572 if ( CHV_IS_NONSYMMETRIC(chv) ) {
573 /*
574 --------------------
575 nonsymmetric chevron
576 --------------------
577 */
578 jcol = -1 ;
579 maxval = 0.0 ;
580 off = nD + nL - 1 - irow ;
581 stride = 2*nD + nL + nU - 1 ;
582 for ( jj = 0 ; jj < irow ; jj++ ) {
583 val = Zabs(entries[2*off], entries[2*off+1]) ;
584 if ( jcol == -1 || maxval < val ) {
585 jcol = jj ; maxval = val ;
586 }
587 off += stride ;
588 stride -= 2 ;
589 }
590 for ( jj = irow ; jj < ncol ; jj++, off++ ) {
591 val = Zabs(entries[2*off], entries[2*off+1]) ;
592 if ( jcol == -1 || maxval < val ) {
593 jcol = jj ; maxval = val ;
594 }
595 }
596 } else if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
597 /*
598 ------------------------------
599 hermitian or symmetric chevron
600 ------------------------------
601 */
602 jcol = -1 ;
603 maxval = 0.0 ;
604 off = irow ;
605 stride = nD + nU - 1 ;
606 for ( jj = 0 ; jj < irow ; jj++ ) {
607 val = Zabs(entries[2*off], entries[2*off+1]) ;
608 if ( jcol == -1 || maxval < val ) {
609 jcol = jj ; maxval = val ;
610 }
611 off += stride ;
612 stride-- ;
613 }
614 for ( jj = irow ; jj < ncol ; jj++, off++ ) {
615 val = Zabs(entries[2*off], entries[2*off+1]) ;
616 if ( jcol == -1 || maxval < val ) {
617 jcol = jj ; maxval = val ;
618 }
619 }
620 }
621 } else {
622 fprintf(stderr,
623 "\n fatal error in Chv_maxabsInRow(%p,%d,%p)"
624 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX \n",
625 chv, irow, pmaxval, chv->symflag) ;
626 exit(-1) ;
627 }
628 *pmaxval = maxval ;
629
630 return(jcol) ; }
631
632 /*--------------------------------------------------------------------*/
633 /*
634 --------------------------------------
635 return the location of the first entry
636 with largest magnitude in column jcol.
637 *pmaxval is filled with its magnitude.
638
639 created -- 98apr30, cca
640 --------------------------------------
641 */
642 int
Chv_maxabsInColumn(Chv * chv,int jcol,double * pmaxval)643 Chv_maxabsInColumn (
644 Chv *chv,
645 int jcol,
646 double *pmaxval
647 ) {
648 double maxval, val ;
649 double *entries ;
650 int irow, ii, nD, nL, nrow, nU, off, stride ;
651 /*
652 --------------
653 check the data
654 --------------
655 */
656 if ( chv == NULL || jcol < 0 || pmaxval == NULL ) {
657 fprintf(stderr, "\n fatal error in Chv_maxabsInColumn(%p,%d,%p)"
658 "\n bad input\n", chv, jcol, pmaxval) ;
659 exit(-1) ;
660 }
661 Chv_dimensions(chv, &nD, &nL, &nU) ;
662 entries = Chv_entries(chv) ;
663 nrow = nD + nL ;
664 irow = -1 ;
665 maxval = 0.0 ;
666 if ( CHV_IS_REAL(chv) ) {
667 if ( CHV_IS_NONSYMMETRIC(chv) ) {
668 /*
669 --------------------
670 nonsymmetric chevron
671 --------------------
672 */
673 maxval = 0.0 ;
674 off = nD + nL + jcol - 1 ;
675 stride = 2*nD + nL + nU - 3 ;
676 for ( ii = 0 ; ii < jcol ; ii++ ) {
677 val = fabs(entries[off]) ;
678 if ( irow == -1 || maxval < val ) {
679 irow = ii ; maxval = val ;
680 }
681 off += stride ;
682 stride -= 2 ;
683 }
684 for ( ii = jcol ; ii < nrow ; ii++, off-- ) {
685 val = fabs(entries[off]) ;
686 if ( irow == -1 || maxval < val ) {
687 irow = ii ; maxval = val ;
688 }
689 }
690 } else if ( CHV_IS_SYMMETRIC(chv) ) {
691 /*
692 -----------------
693 symmetric chevron
694 -----------------
695 */
696 maxval = 0.0 ;
697 off = jcol ;
698 stride = nD + nU - 1 ;
699 for ( ii = 0 ; ii < jcol ; ii++ ) {
700 val = fabs(entries[off]) ;
701 if ( irow == -1 || maxval < val ) {
702 irow = ii ; maxval = val ;
703 }
704 off += stride ;
705 stride-- ;
706 }
707 for ( ii = jcol ; ii < nrow ; ii++, off++ ) {
708 val = fabs(entries[off]) ;
709 if ( irow == -1 || maxval < val ) {
710 irow = ii ; maxval = val ;
711 }
712 }
713 }
714 } else if ( CHV_IS_COMPLEX(chv) ) {
715 if ( CHV_IS_NONSYMMETRIC(chv) ) {
716 /*
717 --------------------
718 nonsymmetric chevron
719 --------------------
720 */
721 maxval = 0.0 ;
722 off = nD + nL + jcol - 1 ;
723 stride = 2*nD + nL + nU - 3 ;
724 for ( ii = 0 ; ii < jcol ; ii++ ) {
725 val = Zabs(entries[2*off], entries[2*off+1]) ;
726 if ( irow == -1 || maxval < val ) {
727 irow = ii ; maxval = val ;
728 }
729 off += stride ;
730 stride -= 2 ;
731 }
732 for ( ii = jcol ; ii < nrow ; ii++, off-- ) {
733 val = Zabs(entries[2*off], entries[2*off+1]) ;
734 if ( irow == -1 || maxval < val ) {
735 irow = ii ; maxval = val ;
736 }
737 }
738 } else if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
739 /*
740 ------------------------------
741 hermitian or symmetric chevron
742 ------------------------------
743 */
744 maxval = 0.0 ;
745 off = jcol ;
746 stride = nD + nU - 1 ;
747 for ( ii = 0 ; ii < jcol ; ii++ ) {
748 val = Zabs(entries[2*off], entries[2*off+1]) ;
749 if ( irow == -1 || maxval < val ) {
750 irow = ii ; maxval = val ;
751 }
752 off += stride ;
753 stride-- ;
754 }
755 for ( ii = jcol ; ii < nrow ; ii++, off++ ) {
756 val = Zabs(entries[2*off], entries[2*off+1]) ;
757 if ( irow == -1 || maxval < val ) {
758 irow = ii ; maxval = val ;
759 }
760 }
761 }
762 } else {
763 fprintf(stderr,
764 "\n fatal error in Chv_maxabsInColumn(%p,%d,%p)"
765 "\n bad symflag %d \n", chv, jcol, pmaxval, chv->symflag) ;
766 exit(-1) ;
767 }
768 *pmaxval = maxval ;
769
770 return(irow) ; }
771
772 /*--------------------------------------------------------------------*/
773 /*
774 -------------------------------------------------------------
775 return the magnitude of a quasimax entry from the unmarked
776 rows and columns and fill *pirow and *pjcol with its location
777
778 created -- 98apr30, cca
779 -------------------------------------------------------------
780 */
781 double
Chv_quasimax(Chv * chv,int rowmark[],int colmark[],int tag,int * pirow,int * pjcol)782 Chv_quasimax (
783 Chv *chv,
784 int rowmark[],
785 int colmark[],
786 int tag,
787 int *pirow,
788 int *pjcol
789 ) {
790 double maxval ;
791 int irow, jcol, nD, qcol, qrow ;
792 /*
793 ---------------
794 check the input
795 ---------------
796 */
797 if ( chv == NULL || rowmark == NULL || colmark == NULL
798 || pirow == NULL || pjcol == NULL ) {
799 fprintf(stderr,
800 "\n fatal error in Chv_quasimax(%p,%p,%p,%d,%p,%p)"
801 "\n bad input\n", chv, rowmark, colmark, tag, pirow, pjcol) ;
802 exit(-1) ;
803 }
804 if ( ! CHV_IS_NONSYMMETRIC(chv) ) {
805 fprintf(stderr,
806 "\n fatal error in Chv_quasimax(%p,%p,%p,%d,%p,%p)"
807 "\n chv->symflag = %d"
808 "\n chevron is not symmetric or hermitian"
809 "\n method cannot be used \n",
810 chv, rowmark, colmark, tag, pirow, pjcol, chv->symflag) ;
811 exit(-1) ;
812 }
813 nD = chv->nD ;
814 /*
815 ----------------------
816 set the default values
817 ----------------------
818 */
819 *pirow = *pjcol = -1 ;
820 maxval = 0.0 ;
821 /*
822 -----------------------
823 find an unmarked column
824 -----------------------
825 */
826 for ( jcol = 0 ; jcol < nD ; jcol++ ) {
827 if ( colmark[jcol] == tag ) {
828 break ;
829 }
830 }
831 if ( jcol == nD ) {
832 /*
833 ---------------------------
834 no unmarked columns, return
835 ---------------------------
836 */
837 return(maxval) ;
838 }
839 /*
840 ----------------------------------------
841 find a maxabs entry in the unmarked rows
842 ----------------------------------------
843 */
844 irow = Chv_maxabsInColumn11(chv, jcol, rowmark, tag, &maxval) ;
845 #if MYDEBUG > 0
846 fprintf(stdout, "\n first maxabs = %12.4e in (%d,%d)",
847 Chv_entry(chv, irow, jcol), irow, jcol) ;
848 fflush(stdout) ;
849 #endif
850 if ( irow == -1 ) {
851 /*
852 ------------------------
853 no unmarked rows, return
854 ------------------------
855 */
856 return(maxval) ;
857 }
858 while ( 1 ) {
859 /*
860 ----------------------------------
861 find a new maxabs entry in the row
862 ----------------------------------
863 */
864 qcol = Chv_maxabsInRow11(chv, irow, colmark, tag, &maxval) ;
865 #if MYDEBUG > 0
866 fprintf(stdout, "\n new maxabs = %12.4e in (%d,%d)",
867 Chv_entry(chv, irow, qcol), irow, qcol) ;
868 fflush(stdout) ;
869 #endif
870 if ( qcol == jcol ) {
871 /*
872 --------------------------------------------
873 same as before, break out of the search loop
874 --------------------------------------------
875 */
876 break ;
877 }
878 jcol = qcol ;
879 /*
880 -------------------------------------
881 find a new maxabs entry in the column
882 -------------------------------------
883 */
884 qrow = Chv_maxabsInColumn11(chv, jcol, rowmark, tag, &maxval) ;
885 #if MYDEBUG > 0
886 fprintf(stdout, "\n new maxabs = %12.4e in (%d,%d)",
887 Chv_entry(chv, qrow, jcol), qrow, jcol) ;
888 fflush(stdout) ;
889 #endif
890 if ( qrow == irow ) {
891 /*
892 --------------------------------------------
893 same as before, break out of the search loop
894 --------------------------------------------
895 */
896 break ;
897 }
898 irow = qrow ;
899 }
900 /*
901 --------------------------------------------------------
902 set the row and column where the quasimax entry is found
903 --------------------------------------------------------
904 */
905 *pjcol = jcol ;
906 *pirow = irow ;
907
908 return(maxval) ; }
909
910 /*--------------------------------------------------------------------*/
911 /*
912 ---------------------------------------------------------------
913 find a 1x1 or 2x2 pivot using the fast Bunch-Parlett algorithm.
914 used only with symmetric chevrons.
915
916 created -- 98apr30, cca
917 ---------------------------------------------------------------
918 */
919 void
Chv_fastBunchParlettPivot(Chv * chv,int mark[],int tag,int * pirow,int * pjcol)920 Chv_fastBunchParlettPivot (
921 Chv *chv,
922 int mark[],
923 int tag,
924 int *pirow,
925 int *pjcol
926 ) {
927 double maxdiag, gamma_r, gamma_s ;
928 double *entries ;
929 int nD, nL, nU, r, s, t ;
930 /*
931 --------------
932 check the data
933 --------------
934 */
935 if ( chv == NULL || mark == NULL || pirow == NULL || pjcol == NULL ) {
936 fprintf(stderr,
937 "\n fatal error in Chv_fastBunchParlettPivot(%p,%p,%d,%p,%p)"
938 "\n bad input\n",
939 chv, mark, tag, pirow, pjcol) ;
940 exit(-1) ;
941 }
942 Chv_dimensions(chv, &nD, &nL, &nU) ;
943 entries = Chv_entries(chv) ;
944 /*
945 ----------------------
946 set the default values
947 ----------------------
948 */
949 *pirow = *pjcol = -1 ;
950 /*
951 ------------------------------------------
952 find an unmarked entry of maximum magitude
953 ------------------------------------------
954 */
955 r = Chv_maxabsInDiagonal11(chv, mark, tag, &maxdiag) ;
956 if ( r == -1 ) {
957 /*
958 -------------------------------------------------------
959 all rows and columns are marked, return without success
960 -------------------------------------------------------
961 */
962 *pirow = *pjcol = -1 ;
963 return ;
964 }
965 /*
966 -------------------------------------------------------------------
967 find the offdiagonal entry of maximum magnitude in row and column r
968 -------------------------------------------------------------------
969 */
970 s = -1 ;
971 gamma_r = 0.0 ;
972 s = Chv_maxabsInRow11(chv, r, mark, tag, &gamma_r) ;
973 if ( s == -1 ) {
974 /*
975 ---------------------------------------------
976 r is the only unmarked row and column, return
977 ---------------------------------------------
978 */
979 *pirow = *pjcol = r ;
980 return ;
981 }
982 if ( maxdiag >= 0.6404 * gamma_r ) {
983 /*
984 -------------------------
985 1 x 1 pivot is acceptable
986 -------------------------
987 */
988 *pirow = *pjcol = r ;
989 return ;
990 } else {
991 /*
992 ---------------
993 loop until done
994 ---------------
995 */
996 while ( 1 ) {
997 /*
998 ----------------------------------------------
999 find
1000 t = index of max off diag entry in column s
1001 gamma_s is its magnitude
1002 ----------------------------------------------
1003 */
1004 t = Chv_maxabsInRow11(chv, s, mark, tag, &gamma_s) ;
1005 if ( t == r || gamma_s == gamma_r ) {
1006 /*
1007 -------------------------
1008 2 x 2 pivot is acceptable
1009 -------------------------
1010 */
1011 *pirow = r ;
1012 *pjcol = s ;
1013 return ;
1014 } else {
1015 /*
1016 --------------------------------
1017 keep looking for a local maximum
1018 --------------------------------
1019 */
1020 r = s ;
1021 gamma_r = gamma_s ;
1022 s = t ;
1023 }
1024 }
1025 }
1026 return ; }
1027
1028 /*--------------------------------------------------------------------*/
1029