Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
XtrUmfpack.c
Go to the documentation of this file.
1
14#include <time.h>
15
16#include "fasp.h"
17#include "fasp_functs.h"
18
19#if WITH_UMFPACK
20#include "umfpack.h"
21#endif
22
23/*---------------------------------*/
24/*-- Public Functions --*/
25/*---------------------------------*/
26
44INT fasp_solver_umfpack(dCSRmat* ptrA, dvector* b, dvector* u, const SHORT prtlvl)
45{
46
47#if WITH_UMFPACK
48
49 const INT n = ptrA->col;
50 INT* Ap = ptrA->IA;
51 INT* Ai = ptrA->JA;
52 double* Ax = ptrA->val;
53
54 INT status = FASP_SUCCESS;
55 void* Symbolic;
56 void* Numeric;
57 double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO];
58
59 if (prtlvl > PRINT_SOME) {
60 Control[UMFPACK_PRL] = prtlvl;
61 } else {
62 Control[UMFPACK_PRL] = 0;
63 }
64
65#if DEBUG_MODE
66 const INT m = ptrA->row;
67 const INT nnz = ptrA->nnz;
68 printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
69 printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", m, n, nnz);
70#endif
71
72 REAL start_time, end_time;
73 fasp_gettime(&start_time);
74
75 /* Symbolic factorization */
76 status = umfpack_di_symbolic(n, n, Ap, Ai, Ax, &Symbolic, Control, Info);
77 if (prtlvl > PRINT_MORE) umfpack_di_report_status(Control, status);
78 if (prtlvl > PRINT_MOST) umfpack_di_report_control(Control);
79 if (prtlvl > PRINT_MORE) umfpack_di_report_info(Control, Info);
80 if (status < 0) {
81 printf("### ERROR: %d, %s %d\n", status, __FUNCTION__, __LINE__);
82 printf("### ERROR: UMFPACK symbolic factorization failed!\n");
83 exit(ERROR_SOLVER_MISC);
84 }
85
86 /* Numeric factorization */
87 status = umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &Numeric, Control, Info);
88 if (prtlvl > PRINT_MORE) umfpack_di_report_status(Control, status);
89 if (prtlvl > PRINT_MOST) umfpack_di_report_control(Control);
90 if (prtlvl > PRINT_MORE) umfpack_di_report_info(Control, Info);
91 if (status < 0) {
92 printf("### ERROR: %d, %s %d\n", status, __FUNCTION__, __LINE__);
93 printf("### ERROR: UMFPACK numeric factorization failed!\n");
94 exit(ERROR_SOLVER_MISC);
95 }
96
97 /* Clean up symbolic factorization */
98 umfpack_di_free_symbolic(&Symbolic);
99
100 /* Solve using after numeric factorization */
101 status =
102 umfpack_di_solve(UMFPACK_A, Ap, Ai, Ax, u->val, b->val, Numeric, Control, Info);
103 if (prtlvl > PRINT_MORE) umfpack_di_report_status(Control, status);
104 if (prtlvl > PRINT_MOST) umfpack_di_report_control(Control);
105 if (prtlvl > PRINT_MORE) umfpack_di_report_info(Control, Info);
106 if (status < 0) {
107 printf("### ERROR: %d, %s %d\n", status, __FUNCTION__, __LINE__);
108 printf("### ERROR: UMFPACK solve failed!\n");
109 exit(ERROR_SOLVER_MISC);
110 }
111
112 /* Clean up numeric factorization */
113 umfpack_di_free_numeric(&Numeric);
114
115 if (prtlvl > PRINT_MIN) {
116 fasp_gettime(&end_time);
117 fasp_cputime("UMFPACK costs", end_time - start_time);
118 }
119
120#if DEBUG_MODE
121 printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
122#endif
123
124 return status;
125
126#else
127
128 printf("### ERROR: UMFPACK is not available!\n");
129 return ERROR_SOLVER_EXIT;
130
131#endif
132}
133
134#if WITH_UMFPACK
145void* fasp_umfpack_factorize(dCSRmat* ptrA, const SHORT prtlvl)
146{
147 INT status = FASP_SUCCESS;
148 const INT n = ptrA->col;
149 INT* Ap = ptrA->IA;
150 INT* Ai = ptrA->JA;
151 REAL* Ax = ptrA->val;
152 void* Symbolic;
153 void* Numeric;
154 double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO];
155
156 if (prtlvl > PRINT_SOME) {
157 Control[UMFPACK_PRL] = prtlvl;
158 } else {
159 Control[UMFPACK_PRL] = 0;
160 }
161
162#if DEBUG_MODE
163 const INT m = ptrA->row;
164 const INT nnz = ptrA->nnz;
165 printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
166 printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", m, n, nnz);
167#endif
168
169 REAL start_time, end_time;
170 fasp_gettime(&start_time);
171
172 /* Symbolic factorization */
173 status = umfpack_di_symbolic(n, n, Ap, Ai, Ax, &Symbolic, Control, Info);
174 if (prtlvl > PRINT_MORE) umfpack_di_report_status(Control, status);
175 if (prtlvl > PRINT_MOST) umfpack_di_report_control(Control);
176 if (prtlvl > PRINT_MORE) umfpack_di_report_info(Control, Info);
177 if (status < 0) {
178 printf("### ERROR: %d, %s %d\n", status, __FUNCTION__, __LINE__);
179 printf("### ERROR: UMFPACK symbolic factorization failed!\n");
180 exit(ERROR_SOLVER_MISC);
181 }
182
183 /* Numeric factorization */
184 status = umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &Numeric, Control, Info);
185 if (prtlvl > PRINT_MORE) umfpack_di_report_status(Control, status);
186 if (prtlvl > PRINT_MOST) umfpack_di_report_control(Control);
187 if (prtlvl > PRINT_MORE) umfpack_di_report_info(Control, Info);
188 if (status < 0) {
189 printf("### ERROR: %d, %s %d\n", status, __FUNCTION__, __LINE__);
190 printf("### ERROR: UMFPACK numeric factorization failed!\n");
191 exit(ERROR_SOLVER_MISC);
192 }
193
194 umfpack_di_free_symbolic(&Symbolic);
195
196 if (prtlvl > PRINT_MIN) {
197 fasp_gettime(&end_time);
198 fasp_cputime("UMFPACK setup", end_time - start_time);
199 }
200
201#if DEBUG_MODE
202 printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
203#endif
204
205 return Numeric;
206}
207
222INT fasp_umfpack_solve(dCSRmat* ptrA, dvector* b, dvector* u, void* Numeric,
223 const SHORT prtlvl)
224{
225 INT status = FASP_SUCCESS;
226 INT* Ap = ptrA->IA;
227 INT* Ai = ptrA->JA;
228 double* Ax = ptrA->val;
229 double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO];
230
231 if (prtlvl > PRINT_SOME) {
232 Control[UMFPACK_PRL] = prtlvl;
233 } else {
234 Control[UMFPACK_PRL] = 0;
235 }
236
237#if DEBUG_MODE
238 const INT m = ptrA->row;
239 const INT n = ptrA->col;
240 const INT nnz = ptrA->nnz;
241 printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
242 printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", m, n, nnz);
243#endif
244
245 REAL start_time, end_time;
246 fasp_gettime(&start_time);
247
248 /* Solve using after numeric factorization */
249 status =
250 umfpack_di_solve(UMFPACK_A, Ap, Ai, Ax, u->val, b->val, Numeric, Control, Info);
251 if (prtlvl > PRINT_MORE) umfpack_di_report_status(Control, status);
252 if (prtlvl > PRINT_MOST) umfpack_di_report_control(Control);
253 if (prtlvl > PRINT_MORE) umfpack_di_report_info(Control, Info);
254 if (status < 0) {
255 printf("### ERROR: %d, %s %d\n", status, __FUNCTION__, __LINE__);
256 printf("### ERROR: UMFPACK solve failed!\n");
257 exit(ERROR_SOLVER_MISC);
258 }
259
260 if (prtlvl > PRINT_NONE) {
261 fasp_gettime(&end_time);
262 fasp_cputime("UMFPACK solve", end_time - start_time);
263 }
264
265#if DEBUG_MODE
266 printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
267#endif
268
269 return status;
270}
271
281INT fasp_umfpack_free_numeric(void* Numeric)
282{
283 INT status = FASP_SUCCESS;
284
285#if DEBUG_MODE
286 printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
287#endif
288
289 umfpack_di_free_numeric(&Numeric);
290
291#if DEBUG_MODE
292 printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
293#endif
294
295 return status;
296}
297
298#endif
299
300/*---------------------------------*/
301/*-- End of File --*/
302/*---------------------------------*/
fasp_cputime
void fasp_cputime(const char *message, const REAL cputime)
Print CPU walltime.
Definition: AuxMessage.c:179
fasp_gettime
void fasp_gettime(REAL *time)
Get system time.
Definition: AuxTiming.c:37
fasp_solver_umfpack
INT fasp_solver_umfpack(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Au=b by UMFpack.
Definition: XtrUmfpack.c:44
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
INT
#define INT
Definition: fasp.h:72
PRINT_MOST
#define PRINT_MOST
Definition: fasp_const.h:77
FASP_SUCCESS
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:19
ERROR_SOLVER_MISC
#define ERROR_SOLVER_MISC
Definition: fasp_const.h:47
PRINT_SOME
#define PRINT_SOME
Definition: fasp_const.h:75
PRINT_MORE
#define PRINT_MORE
Definition: fasp_const.h:76
PRINT_NONE
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
ERROR_SOLVER_EXIT
#define ERROR_SOLVER_EXIT
Definition: fasp_const.h:49
PRINT_MIN
#define PRINT_MIN
Definition: fasp_const.h:74
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::col
INT col
column of matrix A, n
Definition: fasp.h:157
dCSRmat::val
REAL * val
nonzero entries of A
Definition: fasp.h:169
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
dCSRmat::IA
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:163
dCSRmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:160
dCSRmat::JA
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:166
dvector
Vector with n entries of REAL type.
Definition: fasp.h:354
dvector::val
REAL * val
actual vector entries
Definition: fasp.h:360