Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
ilu_setup_csr.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 #include <time.h>
8 
9 #include "fasp.h"
10 #include "fasp_functs.h"
11 
12 /* declarations for ilu.for */
13 #ifdef __cplusplus
14 extern "C" {void iluk_(const INT *n,REAL *a,INT *ja,INT *ia,INT *lfil,REAL *alu,
15  INT *jlu,INT *iwk,INT *ierr,INT *nzlu);}
16 extern "C" {void ilut_(const INT *n,REAL *a,INT *ja,INT *ia,INT *lfil,
17  const REAL *droptol,REAL *alu,INT *jlu,INT *iwk,
18  INT *ierr,INT *nz);}
19 extern "C" {void ilutp_(const INT *n,REAL *a,INT *ja,INT *ia,INT *lfil,
20  const REAL *droptol,const REAL *permtol,const INT *mbloc,
21  REAL *alu,INT *jlu,INT *iwk,INT *ierr,INT *nz);}
22 #else
23 extern void iluk_(const INT *n,REAL *a,INT *ja,INT *ia,INT *lfil,REAL *alu,
24  INT *jlu,INT *iwk,INT *ierr,INT *nzlu);
25 extern void ilut_(const INT *n,REAL *a,INT *ja,INT *ia,INT *lfil,const REAL *droptol,
26  REAL *alu,INT *jlu,INT *iwk,INT *ierr,INT *nz);
27 extern void ilutp_(const INT *n,REAL *a,INT *ja,INT *ia,INT *lfil,const REAL *droptol,
28  const REAL *permtol,const INT *mbloc,REAL *alu,INT *jlu,INT *iwk,
29  INT *ierr,INT *nz);
30 #endif
31 
32 /*---------------------------------*/
33 /*-- Public Functions --*/
34 /*---------------------------------*/
35 
51  ILU_data *iludata,
52  ILU_param *iluparam)
53 {
54 #if FASP_USE_ILU
55  const INT type=iluparam->ILU_type, print_level=iluparam->print_level;
56  const INT n=A->col, nnz=A->nnz, mbloc=n;
57  const REAL ILU_droptol=iluparam->ILU_droptol;
58  const REAL permtol=iluparam->ILU_permtol;
59 
60  // local variable
61  INT lfil=iluparam->ILU_lfil, lfilt=iluparam->ILU_lfil;
62  INT ierr, iwk, nzlu, nwork, *ijlu;
63  REAL *luval;
64 
65  REAL setup_start, setup_end, setup_duration;
66  SHORT status = FASP_SUCCESS;
67 
68 #if DEBUG_MODE > 0
69  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
70  printf("### DEBUG: m=%d, n=%d, nnz=%d\n", A->row, n, nnz);
71 #endif
72 
73  fasp_gettime(&setup_start);
74 
75  // Expected amount of memory for ILU needed and allocate memory
76  switch (type) {
77  case ILUt:
78  iwk=10*nnz; // iwk is the maxim possible nnz for ILU
79  lfilt=floor(n*0.5)+1;
80  break;
81  case ILUtp:
82  iwk=10*nnz; // iwk is the maxim possible nnz for ILU
83  lfilt=floor(n*0.5)+1;
84  break;
85  default: // ILUk
86  if (lfil == 0) iwk=nnz+500;
87  else iwk=(lfil+5)*nnz;
88  break;
89  }
90 
91  nwork = 4*n;
92 
93 #if DEBUG_MODE > 1
94  printf("### DEBUG: fill-in = %d, iwk = %d, nwork = %d\n", lfil, iwk, nwork);
95 #endif
96 
97  // setup ILU preconditioner
98  iludata->row=iludata->col=n;
99  fasp_ilu_data_alloc(iwk, nwork, iludata);
100 
101 #if CHMEM_MODE
102  printf("### DEBUG: memory usage after %s: \n", __FUNCTION__);
103  fasp_mem_usage();
104 #endif
105 
106  // ILU decomposition
107  ijlu=iludata->ijlu;
108  luval=iludata->luval;
109 
110  switch (type) {
111  case ILUt:
112 #if DEBUG_MODE > 0
113  printf("### DEBUG: %s (ILUt) ...... [Start]\n", __FUNCTION__);
114 #endif
115  ilut_(&n,A->val,A->JA,A->IA,&lfilt,&ILU_droptol,luval,ijlu,&iwk,&ierr,&nzlu);
116  break;
117  case ILUtp:
118 #if DEBUG_MODE > 0
119  printf("### DEBUG: %s (ILUp) ...... [Start]\n", __FUNCTION__);
120 #endif
121  ilutp_(&n,A->val,A->JA,A->IA,&lfilt,&ILU_droptol,&permtol,
122  &mbloc,luval,ijlu,&iwk,&ierr,&nzlu);
123  break;
124  default: // ILUk
125 #if DEBUG_MODE > 0
126  printf("### DEBUG: %s (ILUk) ...... [Start]\n", __FUNCTION__);
127 #endif
128  iluk_(&n,A->val,A->JA,A->IA,&lfil,luval,ijlu,&iwk,&ierr,&nzlu);
129  break;
130  }
131 
132  fasp_dcsr_shift(A, -1);
133 
134 #if CHMEM_MODE
135  printf("### DEBUG: memory usage after ILU setup: \n");
136  fasp_mem_usage();
137 #endif
138 
139  iludata->nzlu=nzlu;
140  iludata->nwork=nwork;
141 
142 #if DEBUG_MODE > 1
143  printf("### DEBUG: iwk = %d, nzlu = %d\n",iwk,nzlu);
144 #endif
145 
146  if (ierr!=0) {
147  printf("### ERROR: ILU setup failed (ierr=%d)!\n", ierr);
148  status = ERROR_SOLVER_ILUSETUP;
149  goto FINISHED;
150  }
151 
152  if (iwk<nzlu) {
153  printf("### ERROR: Need more memory for ILU %d!\n", iwk-nzlu);
154  status = ERROR_SOLVER_ILUSETUP;
155  goto FINISHED;
156  }
157 
158  if (print_level>PRINT_NONE) {
159  fasp_gettime(&setup_end);
160  setup_duration = setup_end - setup_start;
161 
162  switch (type) {
163  case ILUt:
164  printf("ILUt setup costs %f seconds.\n", setup_duration);
165  break;
166  case ILUtp:
167  printf("ILUtp setup costs %f seconds.\n", setup_duration);
168  break;
169  default: // ILUk
170  printf("ILUk setup costs %f seconds.\n", setup_duration);
171  break;
172  }
173  }
174 
175 FINISHED:
176 
177 #if DEBUG_MODE > 0
178  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
179 #endif
180 
181  return status;
182 
183 #else // WITH_ILU
184 
185  printf("### ERROR: ILU is not enabled!\n");
186  fasp_chkerr(ERROR_MISC, __FUNCTION__);
187 
188 #endif
189 }
190 
191 /*---------------------------------*/
192 /*-- End of File --*/
193 /*---------------------------------*/
INT col
column of matrix LU, n
Definition: fasp.h:406
INT * ijlu
integer array of row pointers and column indexes, the size is nzlu
Definition: fasp.h:412
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: message.c:199
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
#define REAL
Definition: fasp.h:67
SHORT ILU_type
ILU type for decomposition.
Definition: fasp.h:380
void fasp_mem_usage()
Show total allocated memory currently.
Definition: memory.c:175
REAL * val
nonzero entries of A
Definition: fasp.h:166
REAL ILU_permtol
permuted if permtol*|a(i,j)| > |a(i,i)|
Definition: fasp.h:392
INT ILU_lfil
level of fill-in for ILUk
Definition: fasp.h:383
#define INT
Definition: fasp.h:64
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
Data for ILU setup.
Definition: fasp.h:400
INT nwork
work space size
Definition: fasp.h:421
void fasp_dcsr_shift(dCSRmat *A, INT offset)
Re-index a REAL matrix in CSR format to make the index starting from 0 or 1.
Definition: sparse_csr.c:1085
INT nnz
number of nonzero entries
Definition: fasp.h:157
REAL * luval
nonzero entries of LU
Definition: fasp.h:415
#define ERROR_MISC
Definition: fasp_const.h:35
INT col
column of matrix A, n
Definition: fasp.h:154
#define ERROR_SOLVER_ILUSETUP
Definition: fasp_const.h:52
Main header file for FASP.
void fasp_gettime(REAL *time)
Get system time.
Definition: timing.c:28
INT row
row number of matrix LU, m
Definition: fasp.h:403
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
REAL ILU_droptol
drop tolerance for ILUt
Definition: fasp.h:386
INT row
row number of matrix A, m
Definition: fasp.h:151
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
void fasp_ilu_data_alloc(const INT iwk, const INT nwork, ILU_data *iludata)
Allocate workspace for ILU factorization.
Definition: init.c:118
SHORT fasp_ilu_dcsr_setup(dCSRmat *A, ILU_data *iludata, ILU_param *iluparam)
Get ILU decomposition of a CSR matrix A.
Definition: ilu_setup_csr.c:50
Parameters for ILU.
Definition: fasp.h:374
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
INT nzlu
number of nonzero entries
Definition: fasp.h:409
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160
#define ILUtp
Definition: fasp_const.h:150
#define ILUt
Definition: fasp_const.h:149
SHORT print_level
print level
Definition: fasp.h:377