1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 #include <math.h>
9 #include "_hypre_utilities.h"
10
11 /*--------------------------------------------------------------------------
12 *--------------------------------------------------------------------------*/
13
hypre_swap(HYPRE_Int * v,HYPRE_Int i,HYPRE_Int j)14 void hypre_swap( HYPRE_Int *v,
15 HYPRE_Int i,
16 HYPRE_Int j )
17 {
18 HYPRE_Int temp;
19
20 temp = v[i];
21 v[i] = v[j];
22 v[j] = temp;
23 }
24
25 /*--------------------------------------------------------------------------
26 *--------------------------------------------------------------------------*/
27
hypre_swap_c(HYPRE_Complex * v,HYPRE_Int i,HYPRE_Int j)28 void hypre_swap_c( HYPRE_Complex *v,
29 HYPRE_Int i,
30 HYPRE_Int j )
31 {
32 HYPRE_Complex temp;
33
34 temp = v[i];
35 v[i] = v[j];
36 v[j] = temp;
37 }
38
39 /*--------------------------------------------------------------------------
40 *--------------------------------------------------------------------------*/
41
hypre_swap2(HYPRE_Int * v,HYPRE_Real * w,HYPRE_Int i,HYPRE_Int j)42 void hypre_swap2( HYPRE_Int *v,
43 HYPRE_Real *w,
44 HYPRE_Int i,
45 HYPRE_Int j )
46 {
47 HYPRE_Int temp;
48 HYPRE_Real temp2;
49
50 temp = v[i];
51 v[i] = v[j];
52 v[j] = temp;
53 temp2 = w[i];
54 w[i] = w[j];
55 w[j] = temp2;
56 }
57
58 /*--------------------------------------------------------------------------
59 *--------------------------------------------------------------------------*/
60
hypre_BigSwap2(HYPRE_BigInt * v,HYPRE_Real * w,HYPRE_Int i,HYPRE_Int j)61 void hypre_BigSwap2( HYPRE_BigInt *v,
62 HYPRE_Real *w,
63 HYPRE_Int i,
64 HYPRE_Int j )
65 {
66 HYPRE_BigInt temp;
67 HYPRE_Real temp2;
68
69 temp = v[i];
70 v[i] = v[j];
71 v[j] = temp;
72 temp2 = w[i];
73 w[i] = w[j];
74 w[j] = temp2;
75 }
76
77 /*--------------------------------------------------------------------------
78 *--------------------------------------------------------------------------*/
79
hypre_swap2i(HYPRE_Int * v,HYPRE_Int * w,HYPRE_Int i,HYPRE_Int j)80 void hypre_swap2i( HYPRE_Int *v,
81 HYPRE_Int *w,
82 HYPRE_Int i,
83 HYPRE_Int j )
84 {
85 HYPRE_Int temp;
86
87 temp = v[i];
88 v[i] = v[j];
89 v[j] = temp;
90 temp = w[i];
91 w[i] = w[j];
92 w[j] = temp;
93 }
94
hypre_BigSwap2i(HYPRE_BigInt * v,HYPRE_Int * w,HYPRE_Int i,HYPRE_Int j)95 void hypre_BigSwap2i( HYPRE_BigInt *v,
96 HYPRE_Int *w,
97 HYPRE_Int i,
98 HYPRE_Int j )
99 {
100 HYPRE_BigInt big_temp;
101 HYPRE_Int temp;
102
103 big_temp = v[i];
104 v[i] = v[j];
105 v[j] = big_temp;
106 temp = w[i];
107 w[i] = w[j];
108 w[j] = temp;
109 }
110
111 /*--------------------------------------------------------------------------
112 *--------------------------------------------------------------------------*/
113
114
115 /* AB 11/04 */
116
hypre_swap3i(HYPRE_Int * v,HYPRE_Int * w,HYPRE_Int * z,HYPRE_Int i,HYPRE_Int j)117 void hypre_swap3i( HYPRE_Int *v,
118 HYPRE_Int *w,
119 HYPRE_Int *z,
120 HYPRE_Int i,
121 HYPRE_Int j )
122 {
123 HYPRE_Int temp;
124
125 temp = v[i];
126 v[i] = v[j];
127 v[j] = temp;
128 temp = w[i];
129 w[i] = w[j];
130 w[j] = temp;
131 temp = z[i];
132 z[i] = z[j];
133 z[j] = temp;
134 }
135
136 /*--------------------------------------------------------------------------
137 *--------------------------------------------------------------------------*/
138
hypre_swap3_d(HYPRE_Real * v,HYPRE_Int * w,HYPRE_Int * z,HYPRE_Int i,HYPRE_Int j)139 void hypre_swap3_d( HYPRE_Real *v,
140 HYPRE_Int *w,
141 HYPRE_Int *z,
142 HYPRE_Int i,
143 HYPRE_Int j )
144 {
145 HYPRE_Int temp;
146 HYPRE_Real temp_d;
147
148 temp_d = v[i];
149 v[i] = v[j];
150 v[j] = temp_d;
151 temp = w[i];
152 w[i] = w[j];
153 w[j] = temp;
154 temp = z[i];
155 z[i] = z[j];
156 z[j] = temp;
157 }
158
159 /* swap (v[i], v[j]), (w[i], w[j]), and (z[v[i]], z[v[j]]) - DOK */
hypre_swap3_d_perm(HYPRE_Int * v,HYPRE_Real * w,HYPRE_Int * z,HYPRE_Int i,HYPRE_Int j)160 void hypre_swap3_d_perm( HYPRE_Int *v,
161 HYPRE_Real *w,
162 HYPRE_Int *z,
163 HYPRE_Int i,
164 HYPRE_Int j )
165 {
166 HYPRE_Int temp;
167 HYPRE_Real temp_d;
168
169 temp = v[i];
170 v[i] = v[j];
171 v[j] = temp;
172 temp_d = w[i];
173 w[i] = w[j];
174 w[j] = temp_d;
175 temp = z[v[i]];
176 z[v[i]] = z[v[j]];
177 z[v[j]] = temp;
178 }
179 /*--------------------------------------------------------------------------
180 *--------------------------------------------------------------------------*/
181
hypre_BigSwap4_d(HYPRE_Real * v,HYPRE_BigInt * w,HYPRE_Int * z,HYPRE_Int * y,HYPRE_Int i,HYPRE_Int j)182 void hypre_BigSwap4_d( HYPRE_Real *v,
183 HYPRE_BigInt *w,
184 HYPRE_Int *z,
185 HYPRE_Int *y,
186 HYPRE_Int i,
187 HYPRE_Int j )
188 {
189 HYPRE_Int temp;
190 HYPRE_BigInt big_temp;
191 HYPRE_Real temp_d;
192
193 temp_d = v[i];
194 v[i] = v[j];
195 v[j] = temp_d;
196 big_temp = w[i];
197 w[i] = w[j];
198 w[j] = big_temp;
199 temp = z[i];
200 z[i] = z[j];
201 z[j] = temp;
202 temp = y[i];
203 y[i] = y[j];
204 y[j] = temp;
205
206 }
207
208 /*--------------------------------------------------------------------------
209 *--------------------------------------------------------------------------*/
210
hypre_swap_d(HYPRE_Real * v,HYPRE_Int i,HYPRE_Int j)211 void hypre_swap_d( HYPRE_Real *v,
212 HYPRE_Int i,
213 HYPRE_Int j )
214 {
215 HYPRE_Real temp;
216
217 temp = v[i];
218 v[i] = v[j];
219 v[j] = temp;
220 }
221
222 /*--------------------------------------------------------------------------
223 *--------------------------------------------------------------------------*/
224
hypre_qsort0(HYPRE_Int * v,HYPRE_Int left,HYPRE_Int right)225 void hypre_qsort0( HYPRE_Int *v,
226 HYPRE_Int left,
227 HYPRE_Int right )
228 {
229 HYPRE_Int i, last;
230
231 if (left >= right)
232 {
233 return;
234 }
235 hypre_swap(v, left, (left+right)/2);
236 last = left;
237 for (i = left+1; i <= right; i++)
238 {
239 if (v[i] < v[left])
240 {
241 hypre_swap(v, ++last, i);
242 }
243 }
244 hypre_swap(v, left, last);
245 hypre_qsort0(v, left, last-1);
246 hypre_qsort0(v, last+1, right);
247 }
248
249 /*--------------------------------------------------------------------------
250 *--------------------------------------------------------------------------*/
251
hypre_qsort1(HYPRE_Int * v,HYPRE_Real * w,HYPRE_Int left,HYPRE_Int right)252 void hypre_qsort1( HYPRE_Int *v,
253 HYPRE_Real *w,
254 HYPRE_Int left,
255 HYPRE_Int right )
256 {
257 HYPRE_Int i, last;
258
259 if (left >= right)
260 {
261 return;
262 }
263 hypre_swap2( v, w, left, (left+right)/2);
264 last = left;
265 for (i = left+1; i <= right; i++)
266 {
267 if (v[i] < v[left])
268 {
269 hypre_swap2(v, w, ++last, i);
270 }
271 }
272 hypre_swap2(v, w, left, last);
273 hypre_qsort1(v, w, left, last-1);
274 hypre_qsort1(v, w, last+1, right);
275 }
276
277 /*--------------------------------------------------------------------------
278 *--------------------------------------------------------------------------*/
279
hypre_BigQsort1(HYPRE_BigInt * v,HYPRE_Real * w,HYPRE_Int left,HYPRE_Int right)280 void hypre_BigQsort1( HYPRE_BigInt *v,
281 HYPRE_Real *w,
282 HYPRE_Int left,
283 HYPRE_Int right )
284 {
285 HYPRE_Int i, last;
286
287 if (left >= right)
288 {
289 return;
290 }
291 hypre_BigSwap2(v, w, left, (left+right)/2);
292 last = left;
293 for (i = left+1; i <= right; i++)
294 {
295 if (v[i] < v[left])
296 {
297 hypre_BigSwap2(v, w, ++last, i);
298 }
299 }
300 hypre_BigSwap2(v, w, left, last);
301 hypre_BigQsort1(v, w, left, last-1);
302 hypre_BigQsort1(v, w, last+1, right);
303 }
304
305 /*--------------------------------------------------------------------------
306 *--------------------------------------------------------------------------*/
307
hypre_qsort2i(HYPRE_Int * v,HYPRE_Int * w,HYPRE_Int left,HYPRE_Int right)308 void hypre_qsort2i( HYPRE_Int *v,
309 HYPRE_Int *w,
310 HYPRE_Int left,
311 HYPRE_Int right )
312 {
313 HYPRE_Int i, last;
314
315 if (left >= right)
316 {
317 return;
318 }
319 hypre_swap2i( v, w, left, (left+right)/2);
320 last = left;
321 for (i = left+1; i <= right; i++)
322 {
323 if (v[i] < v[left])
324 {
325 hypre_swap2i(v, w, ++last, i);
326 }
327 }
328 hypre_swap2i(v, w, left, last);
329 hypre_qsort2i(v, w, left, last-1);
330 hypre_qsort2i(v, w, last+1, right);
331 }
332
333 /*--------------------------------------------------------------------------
334 *--------------------------------------------------------------------------*/
335
hypre_BigQsort2i(HYPRE_BigInt * v,HYPRE_Int * w,HYPRE_Int left,HYPRE_Int right)336 void hypre_BigQsort2i( HYPRE_BigInt *v,
337 HYPRE_Int *w,
338 HYPRE_Int left,
339 HYPRE_Int right )
340 {
341 HYPRE_Int i, last;
342
343 if (left >= right)
344 {
345 return;
346 }
347 hypre_BigSwap2i( v, w, left, (left+right)/2);
348 last = left;
349 for (i = left+1; i <= right; i++)
350 {
351 if (v[i] < v[left])
352 {
353 hypre_BigSwap2i(v, w, ++last, i);
354 }
355 }
356 hypre_BigSwap2i(v, w, left, last);
357 hypre_BigQsort2i(v, w, left, last-1);
358 hypre_BigQsort2i(v, w, last+1, right);
359 }
360
361 /*--------------------------------------------------------------------------
362 *--------------------------------------------------------------------------*/
363
364 /* sort on w (HYPRE_Real), move v (AB 11/04) */
365
hypre_qsort2(HYPRE_Int * v,HYPRE_Real * w,HYPRE_Int left,HYPRE_Int right)366 void hypre_qsort2( HYPRE_Int *v,
367 HYPRE_Real *w,
368 HYPRE_Int left,
369 HYPRE_Int right )
370 {
371 HYPRE_Int i, last;
372
373 if (left >= right)
374 {
375 return;
376 }
377 hypre_swap2( v, w, left, (left+right)/2);
378 last = left;
379 for (i = left+1; i <= right; i++)
380 {
381 if (w[i] < w[left])
382 {
383 hypre_swap2(v, w, ++last, i);
384 }
385 }
386 hypre_swap2(v, w, left, last);
387 hypre_qsort2(v, w, left, last-1);
388 hypre_qsort2(v, w, last+1, right);
389 }
390
391 /*--------------------------------------------------------------------------
392 *--------------------------------------------------------------------------*/
393
394 /* qsort2 based on absolute value of entries in w. */
hypre_qsort2_abs(HYPRE_Int * v,HYPRE_Real * w,HYPRE_Int left,HYPRE_Int right)395 void hypre_qsort2_abs( HYPRE_Int *v,
396 HYPRE_Real *w,
397 HYPRE_Int left,
398 HYPRE_Int right )
399 {
400 HYPRE_Int i, last;
401 if (left >= right)
402 {
403 return;
404 }
405 hypre_swap2( v, w, left, (left+right)/2);
406 last = left;
407 for (i = left+1; i <= right; i++)
408 {
409 if (fabs(w[i]) > fabs(w[left]))
410 {
411 hypre_swap2(v, w, ++last, i);
412 }
413 }
414 hypre_swap2(v, w, left, last);
415 hypre_qsort2_abs(v, w, left, last-1);
416 hypre_qsort2_abs(v, w, last+1, right);
417 }
418 /*--------------------------------------------------------------------------
419 *--------------------------------------------------------------------------*/
420
421 /* sort on v, move w and z (AB 11/04) */
422
hypre_qsort3i(HYPRE_Int * v,HYPRE_Int * w,HYPRE_Int * z,HYPRE_Int left,HYPRE_Int right)423 void hypre_qsort3i( HYPRE_Int *v,
424 HYPRE_Int *w,
425 HYPRE_Int *z,
426 HYPRE_Int left,
427 HYPRE_Int right )
428 {
429 HYPRE_Int i, last;
430
431 if (left >= right)
432 {
433 return;
434 }
435 hypre_swap3i( v, w, z, left, (left+right)/2);
436 last = left;
437 for (i = left+1; i <= right; i++)
438 {
439 if (v[i] < v[left])
440 {
441 hypre_swap3i(v, w, z, ++last, i);
442 }
443 }
444 hypre_swap3i(v, w, z, left, last);
445 hypre_qsort3i(v, w, z, left, last-1);
446 hypre_qsort3i(v, w, z, last+1, right);
447 }
448
449 /* sort on v, move w and z DOK */
hypre_qsort3ir(HYPRE_Int * v,HYPRE_Real * w,HYPRE_Int * z,HYPRE_Int left,HYPRE_Int right)450 void hypre_qsort3ir( HYPRE_Int *v,
451 HYPRE_Real *w,
452 HYPRE_Int *z,
453 HYPRE_Int left,
454 HYPRE_Int right )
455 {
456 HYPRE_Int i, last;
457
458 if (left >= right)
459 {
460 return;
461 }
462 hypre_swap3_d_perm( v, w, z, left, (left+right)/2);
463 last = left;
464 for (i = left+1; i <= right; i++)
465 {
466 if (v[i] < v[left])
467 {
468 hypre_swap3_d_perm(v, w, z, ++last, i);
469 }
470 }
471 hypre_swap3_d_perm(v, w, z, left, last);
472 hypre_qsort3ir(v, w, z, left, last-1);
473 hypre_qsort3ir(v, w, z, last+1, right);
474 }
475
476 /*--------------------------------------------------------------------------
477 *--------------------------------------------------------------------------*/
478
479 /* sort min to max based on real array v */
hypre_qsort3(HYPRE_Real * v,HYPRE_Int * w,HYPRE_Int * z,HYPRE_Int left,HYPRE_Int right)480 void hypre_qsort3( HYPRE_Real *v,
481 HYPRE_Int *w,
482 HYPRE_Int *z,
483 HYPRE_Int left,
484 HYPRE_Int right )
485 {
486 HYPRE_Int i, last;
487
488 if (left >= right)
489 {
490 return;
491 }
492 hypre_swap3_d( v, w, z, left, (left+right)/2);
493 last = left;
494 for (i = left+1; i <= right; i++)
495 {
496 if (v[i] < v[left])
497 {
498 hypre_swap3_d(v,w, z, ++last, i);
499 }
500 }
501 hypre_swap3_d(v, w, z, left, last);
502 hypre_qsort3(v, w, z, left, last-1);
503 hypre_qsort3(v, w, z, last+1, right);
504 }
505
506 /*--------------------------------------------------------------------------
507 *--------------------------------------------------------------------------*/
508
509 /* sort min to max based on absolute value */
510
hypre_qsort3_abs(HYPRE_Real * v,HYPRE_Int * w,HYPRE_Int * z,HYPRE_Int left,HYPRE_Int right)511 void hypre_qsort3_abs(HYPRE_Real *v,
512 HYPRE_Int *w,
513 HYPRE_Int *z,
514 HYPRE_Int left,
515 HYPRE_Int right )
516 {
517 HYPRE_Int i, last;
518
519 if (left >= right)
520 {
521 return;
522 }
523 hypre_swap3_d( v, w, z, left, (left+right)/2);
524 last = left;
525 for (i = left+1; i <= right; i++)
526 {
527 if (fabs(v[i]) < fabs(v[left]))
528 {
529 hypre_swap3_d(v,w, z, ++last, i);
530 }
531 }
532 hypre_swap3_d(v, w, z, left, last);
533 hypre_qsort3_abs(v, w, z, left, last-1);
534 hypre_qsort3_abs(v, w, z, last+1, right);
535 }
536
537 /*--------------------------------------------------------------------------
538 *--------------------------------------------------------------------------*/
539
540 /* sort min to max based on absolute value */
541
hypre_BigQsort4_abs(HYPRE_Real * v,HYPRE_BigInt * w,HYPRE_Int * z,HYPRE_Int * y,HYPRE_Int left,HYPRE_Int right)542 void hypre_BigQsort4_abs( HYPRE_Real *v,
543 HYPRE_BigInt *w,
544 HYPRE_Int *z,
545 HYPRE_Int *y,
546 HYPRE_Int left,
547 HYPRE_Int right )
548 {
549 HYPRE_Int i, last;
550
551 if (left >= right)
552 {
553 return;
554 }
555 hypre_BigSwap4_d( v, w, z, y, left, (left+right)/2);
556 last = left;
557 for (i = left+1; i <= right; i++)
558 {
559 if (fabs(v[i]) < fabs(v[left]))
560 {
561 hypre_BigSwap4_d(v,w, z, y, ++last, i);
562 }
563 }
564 hypre_BigSwap4_d(v, w, z, y, left, last);
565 hypre_BigQsort4_abs(v, w, z, y, left, last-1);
566 hypre_BigQsort4_abs(v, w, z, y, last+1, right);
567 }
568
569 /*--------------------------------------------------------------------------
570 *--------------------------------------------------------------------------*/
571 /* sort min to max based on absolute value */
572
hypre_qsort_abs(HYPRE_Real * w,HYPRE_Int left,HYPRE_Int right)573 void hypre_qsort_abs( HYPRE_Real *w,
574 HYPRE_Int left,
575 HYPRE_Int right )
576 {
577 HYPRE_Int i, last;
578
579 if (left >= right)
580 {
581 return;
582 }
583 hypre_swap_d( w, left, (left+right)/2);
584 last = left;
585 for (i = left+1; i <= right; i++)
586 {
587 if (fabs(w[i]) < fabs(w[left]))
588 {
589 hypre_swap_d(w, ++last, i);
590 }
591 }
592 hypre_swap_d(w, left, last);
593 hypre_qsort_abs(w, left, last-1);
594 hypre_qsort_abs(w, last+1, right);
595 }
596
597
598 /*--------------------------------------------------------------------------
599 *--------------------------------------------------------------------------*/
600
hypre_BigSwapbi(HYPRE_BigInt * v,HYPRE_Int * w,HYPRE_Int i,HYPRE_Int j)601 void hypre_BigSwapbi( HYPRE_BigInt *v,
602 HYPRE_Int *w,
603 HYPRE_Int i,
604 HYPRE_Int j )
605 {
606 HYPRE_BigInt big_temp;
607 HYPRE_Int temp;
608
609 big_temp = v[i];
610 v[i] = v[j];
611 v[j] = big_temp;
612 temp = w[i];
613 w[i] = w[j];
614 w[j] = temp;
615 }
616
617 /*--------------------------------------------------------------------------
618 *--------------------------------------------------------------------------*/
619
hypre_BigQsortbi(HYPRE_BigInt * v,HYPRE_Int * w,HYPRE_Int left,HYPRE_Int right)620 void hypre_BigQsortbi( HYPRE_BigInt *v,
621 HYPRE_Int *w,
622 HYPRE_Int left,
623 HYPRE_Int right )
624 {
625 HYPRE_Int i, last;
626
627 if (left >= right)
628 {
629 return;
630 }
631 hypre_BigSwapbi( v, w, left, (left+right)/2);
632 last = left;
633 for (i = left+1; i <= right; i++)
634 {
635 if (v[i] < v[left])
636 {
637 hypre_BigSwapbi(v, w, ++last, i);
638 }
639 }
640 hypre_BigSwapbi(v, w, left, last);
641 hypre_BigQsortbi(v, w, left, last-1);
642 hypre_BigQsortbi(v, w, last+1, right);
643 }
644
645 /*--------------------------------------------------------------------------
646 *--------------------------------------------------------------------------*/
647
hypre_BigSwapLoc(HYPRE_BigInt * v,HYPRE_Int * w,HYPRE_Int i,HYPRE_Int j)648 void hypre_BigSwapLoc( HYPRE_BigInt *v,
649 HYPRE_Int *w,
650 HYPRE_Int i,
651 HYPRE_Int j )
652 {
653 HYPRE_BigInt big_temp;
654
655 big_temp = v[i];
656 v[i] = v[j];
657 v[j] = big_temp;
658 w[i] = j;
659 w[j] = i;
660 }
661
662 /*--------------------------------------------------------------------------
663 *--------------------------------------------------------------------------*/
664
hypre_BigQsortbLoc(HYPRE_BigInt * v,HYPRE_Int * w,HYPRE_Int left,HYPRE_Int right)665 void hypre_BigQsortbLoc( HYPRE_BigInt *v,
666 HYPRE_Int *w,
667 HYPRE_Int left,
668 HYPRE_Int right )
669 {
670 HYPRE_Int i, last;
671
672 if (left >= right)
673 {
674 return;
675 }
676 hypre_BigSwapLoc( v, w, left, (left+right)/2);
677 last = left;
678 for (i = left+1; i <= right; i++)
679 {
680 if (v[i] < v[left])
681 {
682 hypre_BigSwapLoc(v, w, ++last, i);
683 }
684 }
685 hypre_BigSwapLoc(v, w, left, last);
686 hypre_BigQsortbLoc(v, w, left, last-1);
687 hypre_BigQsortbLoc(v, w, last+1, right);
688 }
689
690 /*--------------------------------------------------------------------------
691 *--------------------------------------------------------------------------*/
692
693
hypre_BigSwapb2i(HYPRE_BigInt * v,HYPRE_Int * w,HYPRE_Int * z,HYPRE_Int i,HYPRE_Int j)694 void hypre_BigSwapb2i( HYPRE_BigInt *v,
695 HYPRE_Int *w,
696 HYPRE_Int *z,
697 HYPRE_Int i,
698 HYPRE_Int j )
699 {
700 HYPRE_BigInt big_temp;
701 HYPRE_Int temp;
702
703 big_temp = v[i];
704 v[i] = v[j];
705 v[j] = big_temp;
706 temp = w[i];
707 w[i] = w[j];
708 w[j] = temp;
709 temp = z[i];
710 z[i] = z[j];
711 z[j] = temp;
712 }
713
714 /*--------------------------------------------------------------------------
715 *--------------------------------------------------------------------------*/
716
hypre_BigQsortb2i(HYPRE_BigInt * v,HYPRE_Int * w,HYPRE_Int * z,HYPRE_Int left,HYPRE_Int right)717 void hypre_BigQsortb2i( HYPRE_BigInt *v,
718 HYPRE_Int *w,
719 HYPRE_Int *z,
720 HYPRE_Int left,
721 HYPRE_Int right )
722 {
723 HYPRE_Int i, last;
724
725 if (left >= right)
726 {
727 return;
728 }
729 hypre_BigSwapb2i( v, w, z, left, (left+right)/2);
730 last = left;
731 for (i = left+1; i <= right; i++)
732 {
733 if (v[i] < v[left])
734 {
735 hypre_BigSwapb2i(v, w, z, ++last, i);
736 }
737 }
738 hypre_BigSwapb2i(v, w, z, left, last);
739 hypre_BigQsortb2i(v, w, z, left, last-1);
740 hypre_BigQsortb2i(v, w, z, last+1, right);
741 }
742
743 /*--------------------------------------------------------------------------
744 *--------------------------------------------------------------------------*/
745
hypre_BigSwap(HYPRE_BigInt * v,HYPRE_Int i,HYPRE_Int j)746 void hypre_BigSwap( HYPRE_BigInt *v,
747 HYPRE_Int i,
748 HYPRE_Int j )
749 {
750 HYPRE_BigInt temp;
751
752 temp = v[i];
753 v[i] = v[j];
754 v[j] = temp;
755 }
756
757 /*--------------------------------------------------------------------------
758 *--------------------------------------------------------------------------*/
759
hypre_BigQsort0(HYPRE_BigInt * v,HYPRE_Int left,HYPRE_Int right)760 void hypre_BigQsort0( HYPRE_BigInt *v,
761 HYPRE_Int left,
762 HYPRE_Int right )
763 {
764 HYPRE_Int i, last;
765
766 if (left >= right)
767 {
768 return;
769 }
770 hypre_BigSwap( v, left, (left+right)/2);
771 last = left;
772 for (i = left+1; i <= right; i++)
773 {
774 if (v[i] < v[left])
775 {
776 hypre_BigSwap(v, ++last, i);
777 }
778 }
779 hypre_BigSwap(v, left, last);
780 hypre_BigQsort0(v, left, last-1);
781 hypre_BigQsort0(v, last+1, right);
782 }
783
784 // Recursive DFS search.
hypre_search_row(HYPRE_Int row,const HYPRE_Int * row_ptr,const HYPRE_Int * col_inds,const HYPRE_Complex * data,HYPRE_Int * visited,HYPRE_Int * ordering,HYPRE_Int * order_ind)785 static void hypre_search_row(HYPRE_Int row,
786 const HYPRE_Int *row_ptr,
787 const HYPRE_Int *col_inds,
788 const HYPRE_Complex *data,
789 HYPRE_Int *visited,
790 HYPRE_Int *ordering,
791 HYPRE_Int *order_ind)
792 {
793 // If this row has not been visited, call recursive DFS on nonzero
794 // column entries
795 if (!visited[row])
796 {
797 HYPRE_Int j;
798 visited[row] = 1;
799 for (j=row_ptr[row]; j<row_ptr[row+1]; j++)
800 {
801 HYPRE_Int col = col_inds[j];
802 hypre_search_row(col, row_ptr, col_inds, data,
803 visited, ordering, order_ind);
804 }
805 // Add node to ordering *after* it has been searched
806 ordering[*order_ind] = row;
807 *order_ind += 1;
808 }
809 }
810
811
812 // Find topological ordering on acyclic CSR matrix. That is, find ordering
813 // of matrix to be triangular.
814 //
815 // INPUT
816 // -----
817 // - rowptr[], colinds[], data[] form a CSR structure for nxn matrix
818 // - ordering[] should be empty array of length n
hypre_topo_sort(const HYPRE_Int * row_ptr,const HYPRE_Int * col_inds,const HYPRE_Complex * data,HYPRE_Int * ordering,HYPRE_Int n)819 void hypre_topo_sort( const HYPRE_Int *row_ptr,
820 const HYPRE_Int *col_inds,
821 const HYPRE_Complex *data,
822 HYPRE_Int *ordering,
823 HYPRE_Int n)
824 {
825 HYPRE_Int *visited = hypre_CTAlloc(HYPRE_Int, n, HYPRE_MEMORY_HOST);
826 HYPRE_Int order_ind = 0;
827 HYPRE_Int temp_row = 0;
828 while (order_ind < n)
829 {
830 hypre_search_row(temp_row, row_ptr, col_inds, data,
831 visited, ordering, &order_ind);
832 temp_row += 1;
833 if (temp_row == n)
834 {
835 temp_row = 0;
836 }
837 }
838 hypre_TFree(visited, HYPRE_MEMORY_HOST);
839 }
840
841
842 // Recursive DFS search.
hypre_dense_search_row(HYPRE_Int row,const HYPRE_Complex * L,HYPRE_Int * visited,HYPRE_Int * ordering,HYPRE_Int * order_ind,HYPRE_Int n,HYPRE_Int is_col_major)843 static void hypre_dense_search_row(HYPRE_Int row,
844 const HYPRE_Complex *L,
845 HYPRE_Int *visited,
846 HYPRE_Int *ordering,
847 HYPRE_Int *order_ind,
848 HYPRE_Int n,
849 HYPRE_Int is_col_major)
850 {
851 // If this row has not been visited, call recursive DFS on nonzero
852 // column entries
853 if (!visited[row])
854 {
855 HYPRE_Int col;
856 visited[row] = 1;
857 for (col=0; col<n; col++)
858 {
859 HYPRE_Complex val;
860 if (is_col_major)
861 {
862 val = L[col*n + row];
863 }
864 else
865 {
866 val = L[row*n + col];
867 }
868 if (fabs(val) > 1e-14)
869 {
870 hypre_dense_search_row(col, L, visited, ordering, order_ind, n, is_col_major);
871 }
872 }
873 // Add node to ordering *after* it has been searched
874 ordering[*order_ind] = row;
875 *order_ind += 1;
876 }
877 }
878
879
880 // Find topological ordering of acyclic dense matrix in column major
881 // format. That is, find ordering of matrix to be triangular.
882 //
883 // INPUT
884 // -----
885 // - L[] : dense nxn matrix in column major format
886 // - ordering[] should be empty array of length n
887 // - row is the row to start the search from
hypre_dense_topo_sort(const HYPRE_Complex * L,HYPRE_Int * ordering,HYPRE_Int n,HYPRE_Int is_col_major)888 void hypre_dense_topo_sort(const HYPRE_Complex *L,
889 HYPRE_Int *ordering,
890 HYPRE_Int n,
891 HYPRE_Int is_col_major)
892 {
893 HYPRE_Int *visited = hypre_CTAlloc(HYPRE_Int, n, HYPRE_MEMORY_HOST);
894 HYPRE_Int order_ind = 0;
895 HYPRE_Int temp_row = 0;
896 while (order_ind < n)
897 {
898 hypre_dense_search_row(temp_row, L, visited, ordering, &order_ind, n, is_col_major);
899 temp_row += 1;
900 if (temp_row == n)
901 {
902 temp_row = 0;
903 }
904 }
905 hypre_TFree(visited, HYPRE_MEMORY_HOST);
906 }
907