Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
BlaILU.c
Go to the documentation of this file.
1
16#include <math.h>
17#include <time.h>
18
19#include "fasp.h"
20#include "fasp_functs.h"
21
22/*---------------------------------*/
23/*-- Declare Private Functions --*/
24/*---------------------------------*/
25
26static void fasp_qsplit (REAL *a, INT *ind, INT n, INT ncut);
27static void fasp_sortrow (INT num,INT *q);
28static void fasp_check_col_index (INT row, INT num, INT *q);
29
30/*---------------------------------*/
31/*-- Public Functions --*/
32/*---------------------------------*/
33
72void fasp_iluk (INT n,
73 REAL *a,
74 INT *ja,
75 INT *ia,
76 INT lfil,
77 REAL *alu,
78 INT *jlu,
79 INT iwk,
80 INT *ierr,
81 INT *nzlu)
82{
83#if DEBUG_MODE > 0
84 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
85#endif
86
87 /*----------------------------------------------------------------------
88 SPARSKIT ROUTINE ILUK -- ILU WITH LEVEL OF FILL-IN OF K (ILU(k))
89 ----------------------------------------------------------------------
90
91 on entry:
92 ==========
93 n = integer. The row dimension of the matrix A. The matrix
94
95 a,ja,ia = matrix stored in Compressed Sparse Row format.
96
97 lfil = integer. The fill-in parameter. Each element whose
98 leve-of-fill exceeds lfil during the ILU process is dropped.
99 lfil must be .ge. 0
100
101 iwk = integer. The minimum length of arrays alu, jlu, and levs.
102
103 On return:
104 ===========
105
106 alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
107 the L and U factors together. The diagonal (stored in
108 alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix
109 contains the i-th row of L (excluding the diagonal entry=1)
110 followed by the i-th row of U.
111
112 jlu = integer array of length n containing the pointers to
113 the beginning of each row of U in the matrix alu,jlu.
114
115 levs = integer (work) array of size iwk -- which contains the
116 levels of each element in alu, jlu.
117
118 ierr = integer. Error message with the following meaning.
119 ierr = 0 --> successful return.
120 ierr .gt. 0 --> zero pivot encountered at step number ierr.
121 ierr = -1 --> Error. input matrix may be wrong.
122 (The elimination process has generated a
123 row in L or U whose length is .gt. n.)
124 ierr = -2 --> The matrix L overflows the array al.
125 ierr = -3 --> The matrix U overflows the array alu.
126 ierr = -4 --> Illegal value for lfil.
127 ierr = -5 --> zero row encountered in A or U.
128
129 work arrays:
130 =============
131 jw = integer work array of length 3*n.
132 w = real work array of length n
133
134 ----------------------------------------------------------------------
135 w, ju (1:n) store the working array [1:ii-1 = L-part, ii:n = U-part]
136 jw(n+1:2n) stores the nonzero indicator.
137
138 Notes:
139 ------
140 All the diagonal elements of the input matrix must be nonzero.
141 ---------------------------------------------------------------------- */
142
143 // locals
144 INT ju0, k, j1, j2, j, ii, i, lenl, lenu, jj, jrow, jpos, n2, jlev, NE;
145 REAL t, s, fact;
146 SHORT cinindex=0;
147 REAL *w;
148 INT *ju, *jw, *levs;
149
150 if (lfil < 0) goto F998;
151
152 w = (REAL *)fasp_mem_calloc(n, sizeof(REAL));
153 ju = (INT *)fasp_mem_calloc(n, sizeof(INT));
154 jw = (INT *)fasp_mem_calloc(3*n, sizeof(INT));
155 levs = (INT *)fasp_mem_calloc(iwk, sizeof(INT));
156
157 --jw;
158 --w;
159 --ju;
160 --jlu;
161 --alu;
162 --ia;
163 --ja;
164 --a;
165 --levs;
166
167 /*-----------------------------------------------------------------------
168 shift index for C routines
169 -----------------------------------------------------------------------*/
170 if (ia[1] == 0) cinindex=1 ;
171 if (cinindex)
172 {
173 NE = n + 1; //modify by chunsheng 2012, Sep, 1;
174 for (i=1; i<=NE; ++i) ++ia[i];
175 NE = ia[n+1] - 1;
176 for (i=1; i<=NE; ++i) ++ja[i];
177 }
178
179 /*-----------------------------------------------------------------------
180 initialize ju0 (points to next element to be added to alu,jlu)
181 and pointer array.
182 ------------------------------------------------------------------------*/
183 n2 = n + n;
184 ju0 = n + 2;
185 jlu[1] = ju0;
186
187 // initialize nonzero indicator array + levs array --
188 for(j = 1; j<=2*n; ++j) jw[j] = 0;
189
190 /*-----------------------------------------------------------------------
191 beginning of main loop.
192 ------------------------------------------------------------------------*/
193 for(ii = 1; ii <= n; ++ii) { //500
194 j1 = ia[ii];
195 j2 = ia[ii + 1] - 1;
196
197 // unpack L-part and U-part of row of A in arrays w
198 lenu = 1;
199 lenl = 0;
200 jw[ii] = ii;
201 w[ii] = 0.0;
202 jw[n + ii] = ii;
203
204 //
205 for(j = j1; j <= j2; ++j) { //170
206 k = ja[j];
207 t = a[j];
208 if (t == 0.0) continue; //goto g170;
209 if (k < ii) {
210 ++lenl;
211 jw[lenl] = k;
212 w[lenl] = t;
213 jw[n2 + lenl] = 0;
214 jw[n + k] = lenl;
215 } else if (k == ii) {
216 w[ii] = t;
217 jw[n2 + ii] = 0;
218 } else {
219 ++lenu;
220 jpos = ii + lenu - 1;
221 jw[jpos] = k;
222 w[jpos] = t;
223 jw[n2 + jpos] = 0;
224 jw[n + k] = jpos;
225 }
226
227 } //170
228
229 jj = 0;
230 // eliminate previous rows
231
232 F150:
233 ++jj;
234 if (jj > lenl) goto F160;
235
236 /*----------------------------------------------------------------------
237 in order to do the elimination in the correct order we must select
238 the smallest column index among jw(k), k=jj+1, ..., lenl.
239 -----------------------------------------------------------------------*/
240
241 jrow = jw[jj];
242 k = jj;
243
244 // determine smallest column index
245 for(j = jj + 1; j <= lenl; ++j) { //151
246 if (jw[j] < jrow) {
247 jrow = jw[j];
248 k = j;
249 }
250 } //151
251
252 if (k != jj) {
253 // exchange in jw
254 j = jw[jj];
255 jw[jj] = jw[k];
256 jw[k] = j;
257 // exchange in jw(n+ (pointers/ nonzero indicator).
258 jw[n + jrow] = jj;
259 jw[n + j] = k;
260 // exchange in jw(n2+ (levels)
261 j = jw[n2 + jj];
262 jw[n2 + jj] = jw[n2 + k];
263 jw[n2 + k] = j;
264 // exchange in w
265 s = w[jj];
266 w[jj] = w[k];
267 w[k] = s;
268 }
269
270 // zero out element in row by resetting jw(n+jrow) to zero.
271 jw[n + jrow] = 0;
272
273 // get the multiplier for row to be eliminated (jrow) + its level
274 fact = w[jj]*alu[jrow];
275 jlev = jw[n2 + jj];
276 if (jlev > lfil) goto F150;
277
278 // combine current row and row jrow
279 for(k = ju[jrow]; k <= jlu[jrow + 1] - 1; ++k ) { // 203
280 s = fact*alu[k];
281 j = jlu[k];
282 jpos = jw[n + j];
283 if (j >= ii) {
284 // dealing with upper part.
285 if (jpos == 0) {
286 // this is a fill-in element
287 ++lenu;
288 if (lenu > n) goto F995;
289 i = ii + lenu - 1;
290 jw[i] = j;
291 jw[n + j] = i;
292 w[i] = -s;
293 jw[n2 + i] = jlev + levs[k] + 1;
294 } else {
295 // this is not a fill-in element
296 w[jpos] = w[jpos] - s;
297 jw[n2 + jpos] = MIN(jw[n2 + jpos], jlev + levs[k] + 1);
298 }
299 } else {
300 // dealing with lower part.
301 if (jpos == 0) {
302 // this is a fill-in element
303 ++lenl;
304 if (lenl > n) goto F995;
305 jw[lenl] = j;
306 jw[n + j] = lenl;
307 w[lenl] = -s;
308 jw[n2 + lenl] = jlev + levs[k] + 1;
309 } else {
310 // this is not a fill-in element
311 w[jpos] = w[jpos] - s;
312 jw[n2 + jpos] = MIN(jw[n2 + jpos], jlev + levs[k] + 1);
313 }
314 }
315
316 } //203
317 w[jj] = fact;
318 jw[jj] = jrow;
319 goto F150;
320
321 F160:
322 // reset double-pointer to zero (U-part)
323 for(k = 1; k <= lenu; ++k) jw[n + jw[ii + k - 1]] = 0;
324
325 // update l-matrix
326 for(k = 1; k <= lenl; ++k ) { //204
327 if (ju0 > iwk) goto F996;
328 if (jw[n2 + k] <= lfil) {
329 alu[ju0] = w[k];
330 jlu[ju0] = jw[k];
331 ++ju0;
332 }
333 } //204
334
335 // save pointer to beginning of row ii of U
336 ju[ii] = ju0;
337
338 // update u-matrix
339 for(k = ii + 1; k <= ii + lenu - 1; ++k ) { //302
340 if (ju0 > iwk) goto F997;
341
342 if (jw[n2 + k] <= lfil) {
343 jlu[ju0] = jw[k];
344 alu[ju0] = w[k];
345 levs[ju0] = jw[n2 + k];
346 ++ju0;
347 }
348
349 } //302
350
351 if (w[ii] == 0.0) goto F999;
352 //
353 alu[ii] = 1.0/w[ii];
354
355 // update pointer to beginning of next row of U.
356 jlu[ii + 1] = ju0;
357 /*----------------------------------------------------------------------
358 end main loop
359 -----------------------------------------------------------------------*/
360 } //500
361
362 *nzlu = ju[n] - 1;
363
364 if (cinindex) {
365 for ( i = 1; i <= *nzlu; ++i ) --jlu[i];
366 }
367
368 *ierr = 0;
369
370F100:
371 ++jw;
372 ++w;
373 ++ju;
374 ++jlu;
375 ++alu;
376 ++ia;
377 ++ja;
378 ++a;
379 ++levs;
380
381 fasp_mem_free(w); w = NULL;
382 fasp_mem_free(ju); ju = NULL;
383 fasp_mem_free(jw); jw = NULL;
384 fasp_mem_free(levs); levs = NULL;
385
386#if DEBUG_MODE > 0
387 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
388#endif
389
390 return;
391
392 // incomprehensible error. Matrix must be wrong.
393F995:
394 printf("### ERROR: Incomprehensible error. [%s]\n", __FUNCTION__);
395 *ierr = -1;
396 goto F100;
397
398 // insufficient storage in L.
399F996:
400 printf("### ERROR: Insufficient storage in L. [%s]\n", __FUNCTION__);
401 *ierr = -2;
402 goto F100;
403
404 // insufficient storage in U.
405F997:
406 printf("### ERROR: Insufficient storage in U. [%s]\n", __FUNCTION__);
407 *ierr = -3;
408 goto F100;
409
410 // illegal lfil entered.
411F998:
412 printf("### ERROR: Illegal lfil entered. [%s]\n", __FUNCTION__);
413 *ierr = -4;
414 return;
415
416 // zero row encountered in A or U.
417F999:
418 printf("### ERROR: Zero row encountered in A or U. [%s]\n", __FUNCTION__);
419 *ierr = -5;
420 goto F100;
421 /*----------------end-of-iluk--------------------------------------------
422 ---------------------------------------------------------------------- */
423}
424
467void fasp_ilut (INT n,
468 REAL *a,
469 INT *ja,
470 INT *ia,
471 INT lfil,
472 REAL droptol,
473 REAL *alu,
474 INT *jlu,
475 INT iwk,
476 INT *ierr,
477 INT *nz)
478{
479#if DEBUG_MODE > 0
480 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
481#endif
482
483 /*--------------------------------------------------------------------*
484 *** ILUT preconditioner *** *
485 incomplete LU factorization with dual truncation mechanism *
486 ----------------------------------------------------------------------*
487 Author: Yousef Saad *May, 5, 1990, Latest revision, August 1996 *
488 ----------------------------------------------------------------------*
489 PARAMETERS
490 -----------
491
492 on entry:
493 ==========
494 n = integer. The row dimension of the matrix A. The matrix
495
496 a,ja,ia = matrix stored in Compressed Sparse Row format.
497
498 lfil = integer. The fill-in parameter. Each row of L and each row
499 of U will have a maximum of lfil elements (excluding the diagonal
500 element). lfil must be .ge. 0.
501
502 droptol = real*8. Sets the threshold for dropping small terms in the
503 factorization. See below for details on dropping strategy.
504
505 iwk = integer. The lengths of arrays alu and jlu. If the arrays
506 are not big enough to store the ILU factorizations, ilut
507 will stop with an error message.
508
509 On return:
510 ===========
511
512 alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
513 the L and U factors together. The diagonal (stored in
514 alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix
515 contains the i-th row of L (excluding the diagonal entry=1)
516 followed by the i-th row of U.
517
518 ju = integer array of length n containing the pointers to
519 the beginning of each row of U in the matrix alu,jlu.
520
521 ierr = integer. Error message with the following meaning.
522 ierr = 0 --> successful return.
523 ierr .gt. 0 --> zero pivot encountered at step number ierr.
524 ierr = -1 --> Error. input matrix may be wrong.
525 (The elimination process has generated a
526 row in L or U whose length is .gt. n.)
527 ierr = -2 --> The matrix L overflows the array al.
528 ierr = -3 --> The matrix U overflows the array alu.
529 ierr = -4 --> Illegal value for lfil.
530 ierr = -5 --> zero row encountered.
531
532 work arrays:
533 =============
534 jw = integer work array of length 2*n.
535 w = real work array of length n+1.
536
537 ----------------------------------------------------------------------
538 w, ju (1:n) store the working array [1:ii-1 = L-part, ii:n = u]
539 jw(n+1:2n) stores nonzero indicators
540
541 Notes:
542 ------
543 The diagonal elements of the input matrix must be nonzero (at least
544 'structurally').
545
546 ----------------------------------------------------------------------*
547 ---- Dual drop strategy works as follows. *
548 *
549 1) Theresholding in L and U as set by droptol. Any element whose *
550 magnitude is less than some tolerance (relative to the abs *
551 value of diagonal element in u) is dropped. *
552 *
553 2) Keeping only the largest lfil elements in the i-th row of L *
554 and the largest lfil elements in the i-th row of U (excluding *
555 diagonal elements). *
556 *
557 Flexibility: one can use droptol=0 to get a strategy based on *
558 keeping the largest elements in each row of L and U. Taking *
559 droptol .ne. 0 but lfil=n will give the usual threshold strategy *
560 (however, fill-in is then unpredictible). *
561 ---------------------------------------------------------------------- */
562
563 // locals
564 INT ju0, k, j1, j2, j, ii, i, lenl, lenu, jj, jrow, jpos, NE, len;
565 REAL t, s, fact, tmp;
566 SHORT cinindex=0;
567 REAL *w, *tnorm;
568 INT *ju, *jw;
569
570 if (lfil < 0) goto F998;
571
572 ju = (INT *)fasp_mem_calloc(n, sizeof(INT));
573 jw = (INT *)fasp_mem_calloc(2*n, sizeof(INT));
574 w = (REAL *)fasp_mem_calloc(n+1, sizeof(REAL));
575 tnorm = (REAL *)fasp_mem_calloc(n, sizeof(REAL));
576
577 --jw;
578 --ju;
579 --w;
580 --tnorm;
581 --jlu;
582 --alu;
583 --ia;
584 --ja;
585 --a;
586
587 if (ia[1] == 0) cinindex=1 ;
588
589 if (cinindex)
590 {
591 NE = n + 1; //modify by chunsheng 2012, Sep, 1;
592 for (i=1; i<=NE; ++i) ++ia[i];
593 NE = ia[n+1] - 1;
594 for (i=1; i<=NE; ++i) ++ja[i];
595 }
596
597 /*-----------------------------------------------------------------------
598 initialize ju0 (points to next element to be added to alu,jlu)
599 and pointer array.
600 -----------------------------------------------------------------------*/
601 ju0 = n + 2;
602 jlu[1] = ju0;
603
604 // initialize nonzero indicator array.
605 for (j = 1; j<=n; ++j) jw[n + j] = 0;
606
607 /*-----------------------------------------------------------------------
608 beginning of main loop.
609 -----------------------------------------------------------------------*/
610 for (ii = 1; ii <= n; ++ii ) {
611 j1 = ia[ii];
612 j2 = ia[ii + 1] - 1;
613 tmp = 0.0;
614 for ( k = j1; k<= j2; ++k) tmp = tmp + ABS(a[k]);
615 tmp = tmp/(REAL)(j2 - j1 + 1);
616 tnorm[ii] = tmp*droptol;;
617 }
618
619 for (ii = 1; ii<=n; ++ii) {
620 j1 = ia[ii];
621 j2 = ia[ii + 1] - 1;
622
623 // unpack L-part and U-part of row of A in arrays w
624 lenu = 1;
625 lenl = 0;
626 jw[ii] = ii;
627 w[ii] = 0.0;
628 jw[n + ii] = ii;
629
630 for(j = j1; j<=j2; ++j) {
631 k = ja[j];
632 t = a[j];
633 if (k < ii) {
634 ++lenl;
635 jw[lenl] = k;
636 w[lenl] = t;
637 jw[n + k] = lenl;
638 } else if (k == ii) {
639 w[ii] = t;
640 } else {
641 ++lenu ;
642 jpos = ii + lenu - 1;
643 jw[jpos] = k;
644 w[jpos] = t;
645 jw[n + k] = jpos;
646 }
647 }
648 jj = 0;
649 len = 0;
650
651 // eliminate previous rows
652 F150:
653 ++jj;
654 if (jj > lenl) goto F160;
655
656 /*-----------------------------------------------------------------------
657 in order to do the elimination in the correct order we must select
658 the smallest column index among jw(k), k=jj+1, ..., lenl.
659 -----------------------------------------------------------------------*/
660 jrow = jw[jj];
661 k = jj;
662
663 /*
664 determine smallest column index
665 */
666 for(j = jj + 1; j<=lenl; ++j) { //151
667 if (jw[j] < jrow) {
668 jrow = jw[j];
669 k = j;
670 }
671 } //151
672
673 if (k != jj) {
674 // exchange in jw
675 j = jw[jj];
676 jw[jj] = jw[k];
677 jw[k] = j;
678 // exchange in jr
679 jw[n + jrow] = jj;
680 jw[n + j] = k;
681 // exchange in w
682 s = w[jj];
683 w[jj] = w[k];
684 w[k] = s;
685 }
686
687 // zero out element in row by setting jw(n+jrow) to zero.
688 jw[n + jrow] = 0;
689
690 // get the multiplier for row to be eliminated (jrow).
691 fact = w[jj]*alu[jrow];
692
693 if (ABS(fact) <= droptol) goto F150;
694
695 // combine current row and row jrow
696 for ( k = ju[jrow]; k <= jlu[jrow + 1] - 1; ++k) { //203
697 s = fact*alu[k];
698 j = jlu[k];
699 jpos = jw[n + j];
700 if (j >= ii) {
701 // dealing with upper part.
702 if (jpos == 0)
703 {
704 // this is a fill-in element
705 ++lenu;
706 if (lenu > n) goto F995;
707 i = ii + lenu - 1;
708 jw[i] = j;
709 jw[n + j] = i;
710 w[i] = -s;
711 } else {
712 // this is not a fill-in element
713 w[jpos] = w[jpos] - s;
714 }
715 } else {
716 // dealing with lower part.
717 if (jpos == 0) {
718 // this is a fill-in element
719 ++lenl;
720 if (lenl > n) goto F995;
721 jw[lenl] = j;
722 jw[n + j] = lenl;
723 w[lenl] = -s;
724 } else {
725 // this is not a fill-in element
726 w[jpos] = w[jpos] - s;
727 }
728 }
729 } //203
730
731 /*
732 store this pivot element -- (from left to right -- no danger of
733 overlap with the working elements in L (pivots).
734 */
735 ++len;
736 w[len] = fact;
737 jw[len] = jrow;
738 goto F150;
739
740 F160:
741 // reset double-pointer to zero (U-part)
742 for (k = 1; k <= lenu; ++k ) jw[n + jw[ii + k - 1]] = 0; //308
743
744 // update L-matrix
745 lenl = len;
746 len = MIN(lenl, lfil);
747
748 // sort by quick-split
749 fasp_qsplit(&w[1], &jw[1], lenl, len);
750
751 // store L-part
752 for (k = 1; k <= len; ++k ) { //204
753 if (ju0 > iwk) goto F996;
754 alu[ju0] = w[k];
755 jlu[ju0] = jw[k];
756 ++ju0;
757 }
758
759 // save pointer to beginning of row ii of U
760 ju[ii] = ju0;
761
762 // update U-matrix -- first apply dropping strategy
763 len = 0;
764 for (k = 1; k <= lenu - 1; ++k) {
765 // if ( ABS(w[ii + k]) > droptol*tnorm )
766 if ( ABS(w[ii + k]) > tnorm[ii] ) {
767 ++len;
768 w[ii + len] = w[ii + k];
769 jw[ii + len] = jw[ii + k];
770 }
771 }
772
773 lenu = len + 1;
774 len = MIN(lenu, lfil);
775
776 fasp_qsplit(&w[ii + 1], &jw[ii + 1], lenu - 1, len);
777
778 // copy
779 t = ABS(w[ii]);
780 if (len + ju0 > iwk) goto F997;
781 for (k = ii + 1; k<=ii + len - 1; ++k) { //302
782 jlu[ju0] = jw[k];
783 alu[ju0] = w[k];
784 t = t + ABS(w[k]);
785 ++ju0;
786 }
787
788 // store inverse of diagonal element of u
789 // if (w(ii) .eq. 0.0) w(ii) = (0.0001 + droptol)*tnorm
790 if (w[ii] == 0.0) w[ii] = tnorm[ii];
791
792 alu[ii] = 1.0/w[ii];
793
794 // update pointer to beginning of next row of U.
795 jlu[ii + 1] = ju0;
796 /*-----------------------------------------------------------------------
797 end main loop
798 ----------------------------------------------------------------------- */
799 }
800
801 *nz = ju[n] - 1;
802
803 if (cinindex) {
804 for(i = 1; i <= *nz; ++i) --jlu[i];
805 }
806
807 *ierr = 0;
808
809F100:
810 ++jw;
811 ++ju;
812 ++w;
813 ++tnorm;
814 ++jlu;
815 ++alu;
816 ++ia;
817 ++ja;
818 ++a;
819
820 fasp_mem_free(ju); ju = NULL;
821 fasp_mem_free(jw); jw = NULL;
822 fasp_mem_free(w); w = NULL;
823 fasp_mem_free(tnorm); tnorm = NULL;
824
825#if DEBUG_MODE > 0
826 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
827#endif
828
829 return;
830
831F995: // incomprehensible error. Matrix must be wrong.
832 printf("### ERROR: Input matrix may be wrong. [%s]\n", __FUNCTION__);
833 *ierr = -1;
834 goto F100;
835
836F996: // insufficient storage in L.
837 printf("### ERROR: Insufficient storage in L. [%s]\n", __FUNCTION__);
838 *ierr = -2;
839 goto F100;
840
841F997: // insufficient storage in U.
842 printf("### ERROR: Insufficient storage in U. [%s]\n", __FUNCTION__);
843 *ierr = -3;
844 goto F100;
845
846F998: // illegal lfil entered.
847 *ierr = -4;
848 printf("### ERROR: Illegal lfil entered. [%s]\n", __FUNCTION__);
849 return;
850 /*----------------end-of-ilut--------------------------------------------
851 -----------------------------------------------------------------------*/
852}
853
906void fasp_ilutp (INT n,
907 REAL *a,
908 INT *ja,
909 INT *ia,
910 INT lfil,
911 REAL droptol,
912 REAL permtol,
913 INT mbloc,
914 REAL *alu,
915 INT *jlu,
916 INT *iperm,
917 INT iwk,
918 INT *ierr,
919 INT *nz)
920{
921#if DEBUG_MODE > 0
922 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
923#endif
924
925 /*----------------------------------------------------------------------*
926 *** ILUTP preconditioner -- ILUT with pivoting *** *
927 incomplete LU factorization with dual truncation mechanism *
928 ----------------------------------------------------------------------*
929 author Yousef Saad *Sep 8, 1993 -- Latest revision, August 1996. *
930 ----------------------------------------------------------------------*
931 on entry:
932 ==========
933 n = integer. The dimension of the matrix A.
934
935 a,ja,ia = matrix stored in Compressed Sparse Row format.
936 ON RETURN THE COLUMNS OF A ARE PERMUTED. SEE BELOW FOR
937 DETAILS.
938
939 lfil = integer. The fill-in parameter. Each row of L and each row
940 of U will have a maximum of lfil elements (excluding the
941 diagonal element). lfil must be .ge. 0.
942 ** WARNING: THE MEANING OF LFIL HAS CHANGED WITH RESPECT TO
943 EARLIER VERSIONS.
944
945 droptol = real*8. Sets the threshold for dropping small terms in the
946 factorization. See below for details on dropping strategy.
947
948 lfil = integer. The fill-in parameter. Each row of L and
949 each row of U will have a maximum of lfil elements.
950
951 permtol = tolerance ratio used to determne whether or not to permute
952 two columns. At step i columns i and j are permuted when
953
954 abs(a(i,j))*permtol .gt. abs(a(i,i))
955
956 [0 --> never permute; good values 0.1 to 0.01]
957
958 mbloc = if desired, permuting can be done only within the diagonal
959 blocks of size mbloc. Useful for PDE problems with several
960 degrees of freedom.. If feature not wanted take mbloc=n.
961
962 iwk = integer. The lengths of arrays alu and jlu. If the arrays
963 are not big enough to store the ILU factorizations, ilut
964 will stop with an error message.
965
966 On return:
967 ===========
968
969 alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
970 the L and U factors together. The diagonal (stored in
971 alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix
972 contains the i-th row of L (excluding the diagonal entry=1)
973 followed by the i-th row of U.
974
975 ju = integer array of length n containing the pointers to
976 the beginning of each row of U in the matrix alu,jlu.
977
978 iperm = contains the permutation arrays.
979 iperm(1:n) = old numbers of unknowns
980 iperm(n+1:2*n) = reverse permutation = new unknowns.
981
982 ierr = integer. Error message with the following meaning.
983 ierr = 0 --> successful return.
984 ierr .gt. 0 --> zero pivot encountered at step number ierr.
985 ierr = -1 --> Error. input matrix may be wrong.
986 (The elimination process has generated a
987 row in L or U whose length is .gt. n.)
988 ierr = -2 --> The matrix L overflows the array al.
989 ierr = -3 --> The matrix U overflows the array alu.
990 ierr = -4 --> Illegal value for lfil.
991 ierr = -5 --> zero row encountered.
992
993 work arrays:
994 =============
995 jw = integer work array of length 2*n.
996 w = real work array of length n
997
998 IMPORTANR NOTE:
999 --------------
1000 TO AVOID PERMUTING THE SOLUTION VECTORS ARRAYS FOR EACH LU-SOLVE,
1001 THE MATRIX A IS PERMUTED ON RETURN. [all column indices are
1002 changed]. SIMILARLY FOR THE U MATRIX.
1003 To permute the matrix back to its original state use the loop:
1004
1005 do k=ia(1), ia(n+1)-1
1006 ja(k) = iperm(ja(k))
1007 enddo
1008
1009 -----------------------------------------------------------------------*/
1010
1011 // local variables
1012 INT k, i, j, jrow, ju0, ii, j1, j2, jpos, len, imax, lenu, lenl, jj, icut,NE;
1013 REAL s, tmp, tnorm, xmax, xmax0, fact, t;
1014 SHORT cinindex=0;
1015 REAL *w;
1016 INT *ju, *jw;
1017
1018 if (lfil < 0) goto F998;
1019
1020 ju = (INT *) fasp_mem_calloc(n, sizeof(INT));
1021 jw = (INT *) fasp_mem_calloc(2*n, sizeof(INT));
1022 w = (REAL *)fasp_mem_calloc(n+1, sizeof(REAL));
1023
1024 --ju;
1025 --jw;
1026 --iperm;
1027 --w;
1028 --jlu;
1029 --alu;
1030 --ia;
1031 --ja;
1032 --a;
1033
1034 /*-----------------------------------------------------------------------
1035 shift index for C routines
1036 -----------------------------------------------------------------------*/
1037 if (ia[1] == 0) cinindex=1 ;
1038
1039 if (cinindex)
1040 {
1041 NE = n + 1; //modify by chunsheng 2012, Sep, 1;
1042 for (i=1; i<=NE; ++i) ++ia[i];
1043 NE = ia[n+1] - 1;
1044 for (i=1; i<=NE; ++i) ++ja[i];
1045 }
1046
1047 /*-----------------------------------------------------------------------
1048 initialize ju0 (points to next element to be added to alu,jlu)
1049 and pointer array.
1050 -----------------------------------------------------------------------*/
1051 ju0 = n + 2;
1052 jlu[1] = ju0;
1053
1054
1055 // integer double pointer array.
1056 for ( j = 1; j <= n; ++j ) { //1
1057 jw[n + j] = 0;
1058 iperm[j] = j;
1059 iperm[n + j] = j;
1060 } //1
1061
1062 /*-----------------------------------------------------------------------
1063 beginning of main loop.
1064 -----------------------------------------------------------------------*/
1065 for (ii = 1; ii <= n; ++ii ) { //500
1066 j1 = ia[ii];
1067 j2 = ia[ii + 1] - 1;
1068
1069 tnorm = 0.0;
1070 for (k = j1; k <= j2; ++k ) tnorm = tnorm + ABS( a[k] ); //501
1071 if (tnorm == 0.0) goto F999;
1072 tnorm = tnorm/(REAL)(j2 - j1 + 1);
1073
1074 // unpack L-part and U-part of row of A in arrays w --
1075 lenu = 1;
1076 lenl = 0;
1077 jw[ii] = ii;
1078 w[ii] = 0.0;
1079 jw[n + ii] = ii;
1080 //
1081 for (j = j1; j <= j2; ++j ) { // 170
1082 k = iperm[n + ja[j]];
1083 t = a[j];
1084 if (k < ii) {
1085 ++lenl;
1086 jw[lenl] = k;
1087 w[lenl] = t;
1088 jw[n + k] = lenl;
1089 } else if (k == ii) {
1090 w[ii] = t;
1091 } else {
1092 ++lenu;
1093 jpos = ii + lenu - 1;
1094 jw[jpos] = k;
1095 w[jpos] = t;
1096 jw[n + k] = jpos;
1097 }
1098 } //170
1099
1100 jj = 0;
1101 len = 0;
1102
1103
1104 // eliminate previous rows
1105 F150:
1106 ++jj;
1107 if (jj > lenl) goto F160;
1108
1109 /*-----------------------------------------------------------------------
1110 in order to do the elimination in the correct order we must select
1111 the smallest column index among jw(k), k=jj+1, ..., lenl.
1112 -----------------------------------------------------------------------*/
1113 jrow = jw[jj];
1114 k = jj;
1115
1116 // determine smallest column index
1117 for (j = jj + 1; j <= lenl; ++j) { //151
1118 if (jw[j] < jrow) {
1119 jrow = jw[j];
1120 k = j;
1121 }
1122 }
1123
1124 if (k != jj) {
1125 // exchange in jw
1126 j = jw[jj];
1127 jw[jj] = jw[k];
1128 jw[k] = j;
1129 // exchange in jr
1130 jw[n + jrow] = jj;
1131 jw[n + j] = k;
1132 // exchange in w
1133 s = w[jj];
1134 w[jj] = w[k];
1135 w[k] = s;
1136 }
1137
1138 // zero out element in row by resetting jw(n+jrow) to zero.
1139 jw[n + jrow] = 0;
1140
1141 // get the multiplier for row to be eliminated: jrow
1142 fact = w[jj]*alu[jrow];
1143
1144 // drop term if small
1145 if (ABS(fact) <= droptol) goto F150;
1146
1147 // combine current row and row jrow
1148
1149 for ( k = ju[jrow]; k <= jlu[jrow + 1] - 1; ++k ) { //203
1150 s = fact*alu[k];
1151 // new column number
1152 j = iperm[n + jlu[k]];
1153 jpos = jw[n + j];
1154 if (j >= ii) {
1155 // dealing with upper part.
1156 if (jpos == 0) {
1157 // this is a fill-in element
1158 ++lenu;
1159 i = ii + lenu - 1;
1160 if (lenu > n) goto F995;
1161 jw[i] = j;
1162 jw[n + j] = i;
1163 w[i] = -s;
1164 } else {
1165 // no fill-in element --
1166 w[jpos] = w[jpos] - s;
1167 }
1168
1169 } else {
1170 // dealing with lower part.
1171 if (jpos == 0) {
1172 // this is a fill-in element
1173 ++lenl;
1174 if (lenl > n) goto F995;
1175 jw[lenl] = j;
1176 jw[n + j] = lenl;
1177 w[lenl] = -s;
1178 } else {
1179 // this is not a fill-in element
1180 w[jpos] = w[jpos] - s;
1181 }
1182 }
1183 } //203
1184
1185 /*
1186 store this pivot element -- (from left to right -- no danger of
1187 overlap with the working elements in L (pivots).
1188 */
1189
1190 ++len;
1191 w[len] = fact;
1192 jw[len] = jrow;
1193 goto F150;
1194
1195 F160:
1196 // reset double-pointer to zero (U-part)
1197 for ( k = 1; k <= lenu; ++k ) jw[n + jw[ii + k - 1]] = 0; //308
1198
1199 // update L-matrix
1200 lenl = len;
1201 len = MIN(lenl, lfil);
1202
1203 // sort by quick-split
1204 fasp_qsplit(&w[1], &jw[1], lenl, len);
1205
1206 // store L-part -- in original coordinates ..
1207 for ( k = 1; k <= len; ++k ) { // 204
1208 if (ju0 > iwk) goto F996;
1209 alu[ju0] = w[k];
1210 jlu[ju0] = iperm[jw[k]];
1211 ++ju0;
1212 } //204
1213
1214 // save pointer to beginning of row ii of U
1215 ju[ii] = ju0;
1216
1217 // update U-matrix -- first apply dropping strategy
1218 len = 0;
1219 for(k = 1; k <= lenu - 1; ++k ) {
1220 if ( ABS(w[ii + k]) > droptol*tnorm) {
1221 ++len;
1222 w[ii + len] = w[ii + k];
1223 jw[ii + len] = jw[ii + k];
1224 }
1225 }
1226
1227 lenu = len + 1;
1228 len = MIN(lenu, lfil);
1229 fasp_qsplit(&w[ii + 1], &jw[ii + 1], lenu-1, len);
1230
1231 // determine next pivot --
1232 imax = ii;
1233 xmax = ABS(w[imax]);
1234 xmax0 = xmax;
1235 icut = ii - 1 + mbloc - (ii - 1)%mbloc;
1236
1237 for ( k = ii + 1; k <= ii + len - 1; ++k ) {
1238 t = ABS(w[k]);
1239 if ((t > xmax) && (t*permtol > xmax0) && (jw[k] <= icut)) {
1240 imax = k;
1241 xmax = t;
1242 }
1243 }
1244
1245 // exchange w's
1246 tmp = w[ii];
1247 w[ii] = w[imax];
1248 w[imax] = tmp;
1249
1250 // update iperm and reverse iperm
1251 j = jw[imax];
1252 i = iperm[ii];
1253 iperm[ii] = iperm[j];
1254 iperm[j] = i;
1255
1256 // reverse iperm
1257 iperm[n + iperm[ii]] = ii;
1258 iperm[n + iperm[j]] = j;
1259
1260 //-----------------------------------------------------------------------
1261 if (len + ju0 > iwk) goto F997;
1262
1263
1264 // copy U-part in original coordinates
1265 for ( k = ii + 1; k <= ii + len - 1; ++k ) { //302
1266 jlu[ju0] = iperm[jw[k]];
1267 alu[ju0] = w[k];
1268 ++ju0;
1269 }
1270
1271 // store inverse of diagonal element of u
1272 if (w[ii] == 0.0) w[ii] = (1.0e-4 + droptol)*tnorm;
1273 alu[ii] = 1.0/w[ii];
1274
1275 // update pointer to beginning of next row of U.
1276 jlu[ii + 1] = ju0;
1277
1278 /*-----------------------------------------------------------------------
1279 end main loop
1280 -----------------------------------------------------------------------*/
1281 } //500
1282
1283 // permute all column indices of LU ...
1284 for ( k = jlu[1]; k <= jlu[n + 1] - 1; ++k ) jlu[k] = iperm[n + jlu[k]];
1285
1286 // ...and of A
1287 for ( k = ia[1]; k <= ia[n + 1] - 1; ++k ) ja[k] = iperm[n + ja[k]];
1288
1289 *nz = ju[n]- 1;
1290
1291 if (cinindex) {
1292 for (i = 1; i <= *nz; ++i ) --jlu[i];
1293 }
1294
1295 *ierr = 0;
1296
1297F100:
1298 ++jw;
1299 ++ju;
1300 ++iperm;
1301 ++w;
1302 ++jlu;
1303 ++alu;
1304 ++ia;
1305 ++ja;
1306 ++a;
1307
1308 fasp_mem_free(ju); ju = NULL;
1309 fasp_mem_free(jw); jw = NULL;
1310 fasp_mem_free(w); w = NULL;
1311
1312#if DEBUG_MODE > 0
1313 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1314#endif
1315
1316 return;
1317
1318F995: // incomprehensible error. Matrix must be wrong.
1319 printf("### ERROR: Input matrix may be wrong. [%s]\n", __FUNCTION__);
1320 *ierr = -1;
1321 goto F100;
1322
1323F996: // insufficient storage in L.
1324 printf("### ERROR: Insufficient storage in L. [%s]\n", __FUNCTION__);
1325 *ierr = -2;
1326 goto F100;
1327
1328F997: // insufficient storage in U.
1329 printf("### ERROR: Insufficient storage in U. [%s]\n", __FUNCTION__);
1330 *ierr = -3;
1331 goto F100;
1332
1333F998: // illegal lfil entered.
1334 printf("### ERROR: Illegal lfil entered. [%s]\n", __FUNCTION__);
1335 *ierr = -4;
1336 // goto F100;
1337 return;
1338
1339F999: // zero row encountered
1340 printf("### ERROR: Zero row encountered. [%s]\n", __FUNCTION__);
1341 *ierr = -5;
1342 goto F100;
1343 //----------------end-of-ilutp-------------------------------------------
1344}
1345
1372void fasp_symbfactor (INT n,
1373 INT *colind,
1374 INT *rwptr,
1375 INT levfill,
1376 INT nzmax,
1377 INT *nzlu,
1378 INT *ijlu,
1379 INT *uptr,
1380 INT *ierr)
1381{
1382#if DEBUG_MODE > 0
1383 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1384#endif
1385
1579 INT icolindj, ijlum, i, j, k, m, ibegin, iend, Ujbeg, Ujend,NE;
1580 INT head, prev, lm, actlev, lowct, k1, k2, levp1, lmk, nzi, rowct;
1581 SHORT cinindex=0;
1582 INT *rowll, *lastcol, *levels;
1583
1584 rowll = (INT *)fasp_mem_calloc(n, sizeof(INT));
1585 lastcol = (INT *)fasp_mem_calloc(n, sizeof(INT));
1586 levels = (INT *)fasp_mem_calloc(nzmax, sizeof(INT));
1587
1588 //========================================================================
1589 // Beginning of Executable Statements
1590 //========================================================================
1591
1592 /*-----------------------------------------------------------------------
1593 shift index for C routines
1594 -----------------------------------------------------------------------*/
1595 --rowll;
1596 --lastcol;
1597 --levels;
1598 --colind;
1599 --rwptr;
1600 --ijlu;
1601 --uptr;
1602
1603 if (rwptr[1] == 0) cinindex=1 ;
1604 if (cinindex) {
1605 NE = n + 1;
1606 for (i=1; i<=NE; ++i) ++rwptr[i];
1607 NE = rwptr[n+1] - 1;
1608 for (i=1; i<=NE; ++i) ++colind[i];
1609 }
1610
1611 // --------------------------------------------------------------
1612 // Because the first row of the factor contains no strictly lower
1613 // triangular parts (parts of L), uptr(1) = ijlu(1) = n+2:
1614 // --------------------------------------------------------------
1615 ijlu[1] = n + 2;
1616 uptr[1] = n + 2;
1617
1618 // --------------------------------------------------------
1619 // The storage for the nonzeros of LU must be at least n+1,
1620 // for a diagonal matrix:
1621 // --------------------------------------------------------
1622 *nzlu = n + 1;
1623
1624 // --------------------------------------------------------------------
1625 // Number of allowed levels plus 1; used for the test of accept/reject.
1626 // See the notes about the methodology above.
1627 // --------------------------------------------------------------------
1628 levp1 = levfill + 1;
1629
1630 // -------------------------------------------------------------
1631 // Initially, for all columns there were no nonzeros in the rows
1632 // above, because there are no rows above the first one.
1633 // -------------------------------------------------------------
1634 for (i = 1; i<=n; ++i) lastcol[i] = 0;
1635
1636 // -------------------
1637 // Proceed row by row:
1638 // -------------------
1639 for (i = 1; i <= n; ++i) { // 100
1640
1641 // ----------------------------------------------------------
1642 // Because the matrix diagonal entry is nonzero, the level of
1643 // fill for that diagonal entry is zero:
1644 // ----------------------------------------------------------
1645 levels[i] = 0;
1646
1647 // ----------------------------------------------------------
1648 // ibegin and iend are the beginning of rows i and i+1, resp.
1649 // ----------------------------------------------------------
1650 ibegin = rwptr[i];
1651 iend = rwptr[i + 1];
1652
1653 // -------------------------------------------------------------
1654 // Number of offdiagonal nonzeros in the original matrix's row i
1655 // -------------------------------------------------------------
1656 nzi = iend - ibegin;
1657
1658 // --------------------------------------------------------
1659 // If only the diagonal entry in row i is nonzero, skip the
1660 // fancy stuff; nothing need be done:
1661 // --------------------------------------------------------
1662 if (nzi > 1) {
1663 // ----------------------------------------------------------
1664 // Decrement iend, so that it can be used as the ending index
1665 // in icolind of row i:
1666 // ----------------------------------------------------------
1667 iend = iend - 1;
1668
1669 // ---------------------------------------------------------
1670 // rowct keeps count of the number of nondiagonal entries in
1671 // the current row:
1672 // ---------------------------------------------------------
1673 rowct = 0;
1674
1675 // ------------------------------------------------------------
1676 // For nonzeros in the current row from the original matrix A,
1677 // set lastcol to be the current row number, and the levels of
1678 // the entry to be 1. Note that this is really the true level
1679 // of the element, plus 1. At the same time, load up the work
1680 // array rowll with the column numbers for the original entries
1681 // from row i:
1682 // ------------------------------------------------------------
1683#if DEBUG_MODE > 0
1684 printf("### DEBUG: %s %d row\n", __FUNCTION__, i);
1685#endif
1686
1687 for ( j = ibegin; j <= iend; ++j) {
1688 icolindj = colind[j];
1689 lastcol[icolindj] = i;
1690 if (icolindj != i) {
1691 levels[icolindj] = 1;
1692 rowct = rowct + 1;
1693 rowll[rowct] = icolindj;
1694 }
1695#if DEBUG_MODE > 0
1696 printf("### DEBUG: %d\n", icolindj);
1697#endif
1698 }
1699
1700 // ---------------------------------------------------------
1701 // Sort the entries in rowll, so that the row has its column
1702 // entries in increasing order.
1703 // ---------------------------------------------------------
1704 fasp_sortrow(nzi - 1, &rowll[1]);
1705
1706 //check col index
1707 fasp_check_col_index(i, nzi-1, &rowll[1]);
1708 // ---------------------------------------------------------
1709 // Now set up rowll as a linked list containing the original
1710 // nonzero column numbers, as described in the methods section:
1711 // ---------------------------------------------------------
1712 head = rowll[1];
1713 k1 = n + 1;
1714 for (j = nzi - 1; j >= 1; --j) {
1715 k2 = rowll[j];
1716 rowll[k2] = k1;
1717 k1 = k2;
1718 }
1719
1720 // ------------------------------------------------------------
1721 // Increment count of nonzeros in the LU factors by the number
1722 // of nonzeros in the original matrix's row i. Further
1723 // incrementing will be necessary if any fill-in actually occurs
1724 // ------------------------------------------------------------
1725 *nzlu = *nzlu + nzi - 1;
1726
1727 // ------------------------------------------------------------
1728 // The integer j will be used as a pointer to track through the
1729 // linked list rowll:
1730 // ------------------------------------------------------------
1731 j = head;
1732
1733 // ------------------------------------------------------------
1734 // The integer lowct is used to keep count of the number of
1735 // nonzeros in the current row's strictly lower triangular part,
1736 // for setting uptr pointers to indicate where in ijlu the upperc
1737 // triangular part starts.
1738 // ------------------------------------------------------------
1739 lowct = 0;
1740
1741 // ------------------------------------------------------------
1742 // Fill-in could only have resulted from rows preceding row i,
1743 // so we only need check those rows with index j < i.
1744 // Furthermore, if the current row has a zero in column j,
1745 // there is no need to check the preceding rows; there clearly
1746 // could not be any fill-in from those rows to this entry.
1747 // ------------------------------------------------------------
1748 while (j < i) { //80
1749 // ------------------------------------------------------------
1750 // Increment lower triangular part count, since in this case
1751 // (j<i) we got another entry in L:
1752 // ------------------------------------------------------------
1753 lowct = lowct + 1;
1754
1755 // ---------------------------------------------------------
1756 // If the fill level is zero, there is no way to get fill in
1757 // occuring.
1758 // ---------------------------------------------------------
1759 if (levfill != 0) {
1760
1761 // -----------------------------------------------------
1762 // Ujbeg is beginning index of strictly upper triangular
1763 // part of U's j-th row, and Ujend is the ending index
1764 // of it, in ijlu().
1765 // -----------------------------------------------------
1766 Ujbeg = uptr[j];
1767 Ujend = ijlu[j + 1] - 1;
1768
1769 // -----------------------------------------------------
1770 // Need to set pointer to previous entry before working
1771 // segment of rowll, because if fill occurs that will be
1772 // a moving segment.
1773 // -----------------------------------------------------
1774 prev = j;
1775
1776 // -----------------------------------------------------
1777 // lm is the next nonzero pointer in linked list rowll:
1778 // -----------------------------------------------------
1779 lm = rowll[j];
1780
1781 // -------------------------------------------------------
1782 // lmk is the fill level in this row, caused by
1783 // eliminating column entry j. That is, level s1 from the
1784 // methodology explanation above.
1785 // -------------------------------------------------------
1786 lmk = levels[j];
1787
1788 // -------------------------------------------------------
1789 // Now proceed through the j-th row of U, because in the
1790 // elimination we add a multiple of it to row i to zero
1791 // out entry (i,j). If a column entry in row j of U is
1792 // zero, there is no need to worry about fill, because it
1793 // cannot cause a fill in the corresponding entry of row i
1794 // -------------------------------------------------------
1795 for (m = Ujbeg; m <= Ujend; ++m) { //60
1796 // ----------------------------------------------------
1797 // ijlum is the column number of the current nonzero in
1798 // row j of U:
1799 // ----------------------------------------------------
1800 ijlum = ijlu[m];
1801
1802 // ---------------------------------------------------
1803 // actlev is the actual level (plus 1) of column entry
1804 // j in row i, from summing the level contributions
1805 // s1 and s2 as explained in the methods section.
1806 // Note that the next line could reasonably be
1807 // replaced by, e.g., actlev = max(lmk, levels(m)),
1808 // but this would cause greater fill-in:
1809 // ---------------------------------------------------
1810 actlev = lmk + levels[m];
1811
1812 // ---------------------------------------------------
1813 // If lastcol of the current column entry in U is not
1814 // equal to the current row number i, then the current
1815 // row has a zero in column j, and the earlier row j
1816 // in U has a nonzero, so possible fill can occur.
1817 // ---------------------------------------------------
1818 if (lastcol[ijlum] != i) {
1819
1820 // --------------------------------------------------
1821 // If actlev < levfill + 1, then the new entry has an
1822 // acceptable fill level and needs to be added to the
1823 // data structure.
1824 // --------------------------------------------------
1825 if (actlev <= levp1) {
1826
1827 // -------------------------------------------
1828 // Since the column entry ijlum in the current
1829 // row i is to be filled, we need to update
1830 // lastcol for that column number. Also, the
1831 // level number of the current entry needs to be
1832 // set to actlev. Note that when we finish
1833 // processing this row, the n-vector levels(1:n)
1834 // will be copied over to the corresponding
1835 // trailing part of levels, so that it can be
1836 // used in subsequent rows:
1837 // -------------------------------------------
1838 lastcol[ijlum] = i;
1839 levels[ijlum] = actlev;
1840
1841 // -------------------------------------------
1842 // Now find location in the linked list rowll
1843 // where the fillin entry should be placed.
1844 // Chase through the linked list until the next
1845 // nonzero column is to the right of the fill
1846 // column number.
1847 // -------------------------------------------
1848 while (lm <= ijlum) { //50
1849 prev = lm;
1850 lm = rowll[lm];
1851 } //50
1852
1853 // -------------------------------------------
1854 // Insert new entry into the linked list for
1855 // row i, and increase the nonzero count for LU
1856 // -------------------------------------------
1857 rowll[prev] = ijlum;
1858 rowll[ijlum] = lm;
1859 prev = ijlum;
1860 *nzlu = *nzlu + 1;
1861 }
1862
1863 // -------------------------------------------------
1864 // Else clause is for when lastcol(ijlum) = i. In
1865 // this case, the current column has a nonzero, but
1866 // it resulted from an earlier fill-in or from an
1867 // original matrix entry. In this case, need to
1868 // update the level number for this column to be the
1869 // smaller of the two possible fill contributors,
1870 // the current fill number or the computed one from
1871 // updating this entry from a previous row.
1872 // -------------------------------------------------
1873 } else {
1874 levels[ijlum] = MIN(levels[ijlum], actlev);
1875 }
1876
1877 // -------------------------------------------------
1878 // Now go and pick up the next column entry from row
1879 // j of U:
1880 // -------------------------------------------------
1881
1882 } //60
1883 // -------------------------------------------
1884 // End if clause for levfill not equal to zero
1885 // -------------------------------------------
1886 }
1887
1888 // ------------------------------------------------------
1889 // Pick up next nonzero column index from the linked
1890 // list, and continue processing the i-th row's nonzeros.
1891 // This ends the first while loop (j < i).
1892 // ------------------------------------------------------
1893 j = rowll[j];
1894 } //80
1895
1896 // ---------------------------------------------------------
1897 // Check to see if we have exceeded the allowed memory
1898 // storage before storing the results of computing row i's
1899 // sparsity pattern into the ijlu and uptr data structures.
1900 // ---------------------------------------------------------
1901 if (*nzlu > nzmax) {
1902 printf("### ERROR: More storage needed! [%s]\n", __FUNCTION__);
1903 *ierr = 1;
1904 goto F100;
1905 }
1906
1907 // ---------------------------------------------------------
1908 // Storage is adequate, so update ijlu data structure.
1909 // Row i ends at nzlu + 1:
1910 // ---------------------------------------------------------
1911 ijlu[i + 1] = *nzlu + 1;
1912
1913 // ---------------------------------------------------------
1914 // ... and the upper triangular part of LU begins at
1915 // lowct entries to right of where row i begins.
1916 // ---------------------------------------------------------
1917 uptr[i] = ijlu[i] + lowct;
1918
1919 // -----------------------------------------------------
1920 // Now chase through linked list for row i, recording
1921 // information into ijlu. At same time, put level data
1922 // into the levels array for use on later rows:
1923 // -----------------------------------------------------
1924 j = head;
1925 k1 = ijlu[i];
1926 for (k = k1; k <= *nzlu; ++k) {
1927 ijlu[k] = j;
1928 levels[k] = levels[j];
1929 j = rowll[j];
1930 }
1931
1932 } else {
1933
1934 // ---------------------------------------------------------
1935 // This else clause ends the (nzi > 1) if. If nzi = 1, then
1936 // the update of ijlu and uptr is trivial:
1937 // ---------------------------------------------------------
1938 ijlu[i + 1] = *nzlu + 1;
1939 uptr[i] = ijlu[i];
1940 }
1941
1942 // ----------------------------------------------
1943 // And you thought we would never get through....
1944 // ----------------------------------------------
1945 } //100
1946
1947 if (cinindex) {
1948 for ( i = 1; i <= *nzlu; ++i ) --ijlu[i];
1949 for ( i = 1; i <= n; ++i ) --uptr[i];
1950 NE = rwptr[n + 1] - 1;
1951 for ( i = 1; i <= NE; ++i ) --colind[i];
1952 NE = n + 1;
1953 for ( i = 1; i <= NE; ++i ) --rwptr[i];
1954 }
1955
1956 *ierr = 0;
1957
1958F100:
1959 ++rowll;
1960 ++lastcol;
1961 ++levels;
1962 ++colind;
1963 ++rwptr;
1964 ++ijlu;
1965 ++uptr;
1966
1967 fasp_mem_free(rowll); rowll = NULL;
1968 fasp_mem_free(lastcol); lastcol = NULL;
1969 fasp_mem_free(levels); levels = NULL;
1970
1971#if DEBUG_MODE > 0
1972 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1973#endif
1974
1975 return;
1976 //======================== End of symbfac ==============================
1977}
1978
1979/*---------------------------------*/
1980/*-- Private Functions --*/
1981/*---------------------------------*/
1982
1998static void fasp_qsplit (REAL *a,
1999 INT *ind,
2000 INT n,
2001 INT ncut)
2002{
2003 /*-----------------------------------------------------------------------
2004 does a quick-sort split of a real array.
2005 on input a(1:n). is a real array
2006 on output a(1:n) is permuted such that its elements satisfy:
2007 abs(a(i)) .ge. abs(a(ncut)) for i .lt. ncut and
2008 abs(a(i)) .le. abs(a(ncut)) for i .gt. ncut
2009 ind(1:n) is an integer array which permuted in the same way as a(*).
2010 -----------------------------------------------------------------------*/
2011 REAL tmp, abskey;
2012 INT itmp, first, last, mid, j;
2013
2014 /* Parameter adjustments */
2015 --ind;
2016 --a;
2017
2018 first = 1;
2019 last = n;
2020 if ((ncut < first) || (ncut > last)) return;
2021
2022 // outer loop -- while mid .ne. ncut do
2023F161:
2024 mid = first;
2025 abskey = ABS(a[mid]);
2026 for (j = first + 1; j <= last; ++j ) {
2027 if (ABS(a[j]) > abskey) {
2028 ++mid;
2029 // interchange
2030 tmp = a[mid];
2031 itmp = ind[mid];
2032 a[mid] = a[j];
2033 ind[mid] = ind[j];
2034 a[j] = tmp;
2035 ind[j] = itmp;
2036 }
2037 }
2038
2039 // interchange
2040 tmp = a[mid];
2041 a[mid] = a[first];
2042 a[first] = tmp;
2043 //
2044 itmp = ind[mid];
2045 ind[mid] = ind[first];
2046 ind[first] = itmp;
2047
2048 // test for while loop
2049 if (mid == ncut) {
2050 ++ind;
2051 ++a;
2052 return;
2053 }
2054
2055 if (mid > ncut) {
2056 last = mid - 1;
2057 } else {
2058 first = mid + 1;
2059 }
2060
2061 goto F161;
2062 /*----------------end-of-qsplit------------------------------------------*/
2063}
2064
2077static void fasp_sortrow (INT num,
2078 INT *q)
2079{
2080#if DEBUG_MODE > 0
2081 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
2082 #endif
2119 INT key, icn, ih, ii, i, j, jj;
2120 INT iinc[6] = {0,1, 4, 13, 40, 121};
2121 //data iinc/1, 4, 13, 40, 121/;
2122
2123 --q;
2124 if (num == 0)
2125 icn = 0;
2126 else if (num < 14)
2127 icn = 1;
2128 else if (num < 41)
2129 icn = 2;
2130 else if (num < 122)
2131 icn = 3;
2132 else if (num < 365)
2133 icn = 4;
2134 else
2135 icn = 5;
2136
2137 for(ii = 1; ii <= icn; ++ii) { // 40
2138 ih = iinc[icn + 1 - ii];
2139 for(j = ih + 1; j <= num; ++j) { // 30
2140 i = j - ih;
2141 key = q[j];
2142 for(jj = 1; jj <= j - ih; jj += ih) { // 10
2143 if (key >= q[i]) {
2144 goto F20;
2145 } else {
2146 q[i + ih] = q[i];
2147 i = i - ih;
2148 }
2149 } // 10
2150 F20:
2151 q[i + ih] = key;
2152 } // 30
2153 } // 40
2154
2155 ++q;
2156
2157#if DEBUG_MODE > 0
2158 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
2159#endif
2160 return;
2161}
2162
2175static void fasp_check_col_index (INT row,
2176 INT num,
2177 INT *q)
2178{
2179#if DEBUG_MODE > 0
2180 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
2181#endif
2182
2183 INT ii;
2184 INT num_1 = num - 1;
2185
2186 for ( ii = 0; ii < num_1; ++ii ) {
2187 if ( q[ii] == q[ii+1] ) {
2188 printf("### ERROR: Multiple entries with same col indices!\n");
2189 printf("### ERROR: row = %d, col = %d, %d!\n", row, q[ii], q[ii+1]);
2190 fasp_chkerr(ERROR_SOLVER_ILUSETUP, __FUNCTION__);
2191 }
2192 }
2193
2194#if DEBUG_MODE > 0
2195 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
2196#endif
2197
2198 return;
2199}
2200
2201/*---------------------------------*/
2202/*-- End of File --*/
2203/*---------------------------------*/
fasp_mem_free
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:152
fasp_mem_calloc
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
fasp_chkerr
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
fasp_iluk
void fasp_iluk(INT n, REAL *a, INT *ja, INT *ia, INT lfil, REAL *alu, INT *jlu, INT iwk, INT *ierr, INT *nzlu)
Get ILU factorization with level of fill-in k (ilu(k)) for a CSR matrix A.
Definition: BlaILU.c:72
fasp_ilutp
void fasp_ilutp(INT n, REAL *a, INT *ja, INT *ia, INT lfil, REAL droptol, REAL permtol, INT mbloc, REAL *alu, INT *jlu, INT *iperm, INT iwk, INT *ierr, INT *nz)
Get incomplete LU factorization with pivoting dual truncations of a CSR matrix A.
Definition: BlaILU.c:906
fasp_ilut
void fasp_ilut(INT n, REAL *a, INT *ja, INT *ia, INT lfil, REAL droptol, REAL *alu, INT *jlu, INT iwk, INT *ierr, INT *nz)
Get incomplete LU factorization with dual truncations of a CSR matrix A.
Definition: BlaILU.c:467
fasp_symbfactor
void fasp_symbfactor(INT n, INT *colind, INT *rwptr, INT levfill, INT nzmax, INT *nzlu, INT *ijlu, INT *uptr, INT *ierr)
Symbolic factorization of a CSR matrix A in compressed sparse row format, with resulting factors stor...
Definition: BlaILU.c:1372
fasp.h
Main header file for the FASP project.
MIN
#define MIN(a, b)
Definition: fasp.h:83
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
ABS
#define ABS(a)
Definition: fasp.h:84
INT
#define INT
Definition: fasp.h:72
ERROR_SOLVER_ILUSETUP
#define ERROR_SOLVER_ILUSETUP
Definition: fasp_const.h:46