Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
XtrPardiso.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_PARDISO
20#include "mkl_pardiso.h"
21#include "mkl_spblas.h"
22#include "mkl_types.h"
23#endif
24
25/*---------------------------------*/
26/*-- Public Functions --*/
27/*---------------------------------*/
28
45INT fasp_solver_pardiso(dCSRmat* ptrA, dvector* b, dvector* u, const SHORT prtlvl)
46{
47#if WITH_PARDISO
48
49 INT status = FASP_SUCCESS;
50
51 MKL_INT n = ptrA->col;
52 MKL_INT* ia = ptrA->IA;
53 MKL_INT* ja = ptrA->JA;
54 REAL* a = ptrA->val;
55
56 MKL_INT mtype = 11; /* Real unsymmetric matrix */
57 MKL_INT nrhs = 1; /* Number of right hand sides */
58 MKL_INT idum; /* Integer dummy */
59 MKL_INT iparm[64]; /* Pardiso control parameters */
60 MKL_INT maxfct, mnum, phase, error, msglvl; /* Auxiliary variables */
61
62 REAL* f = b->val; /* RHS vector */
63 REAL* x = u->val; /* Solution vector */
64 void* pt[64]; /* Internal solver memory pointer pt */
65 double ddum; /* Double dummy */
66
67#if DEBUG_MODE
68 printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
69 printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", ptrA->row, ptrA->col, ptrA->nnz);
70#endif
71
72 PARDISOINIT(pt, &mtype, iparm); /* Initialize */
73 iparm[34] = 1; /* Use 0-based indexing */
74 maxfct = 1; /* Maximum number of numerical factorizations */
75 mnum = 1; /* Which factorization to use */
76 error = 0; /* Initialize error flag */
77
78 if (prtlvl > PRINT_MORE) {
79 msglvl = 1;
80 iparm[26] = 1; /* Turn on matrix checker */
81 } else {
82 msglvl = 0; /* Do not print statistical information in file */
83 }
84
85 REAL start_time, end_time;
86 fasp_gettime(&start_time);
87
88 phase = 11; /* Reordering and symbolic factorization */
89 PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm,
90 &msglvl, &ddum, &ddum, &error);
91 if (error != 0) {
92 printf("### ERROR: %d, %s %d\n", error, __FUNCTION__, __LINE__);
93 printf("### ERROR: PARDISO symbolic factorization failed!\n");
94 exit(ERROR_SOLVER_MISC);
95 }
96
97 phase = 22; /* Numeric factorization */
98 PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm,
99 &msglvl, &ddum, &ddum, &error);
100 if (error != 0) {
101 printf("### ERROR: %d, %s %d\n", error, __FUNCTION__, __LINE__);
102 printf("### ERROR: PARDISO numeric factorization failed!\n");
103 exit(ERROR_SOLVER_MISC);
104 }
105
106 phase = 33; /* Back substitution and iterative refinement */
107 PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm,
108 &msglvl, f, x, &error);
109
110 if (error != 0) {
111 printf("### ERROR: %d, %s %d\n", error, __FUNCTION__, __LINE__);
112 printf("### ERROR: Solution failed!\n");
113 exit(ERROR_SOLVER_MISC);
114 }
115
116 if (prtlvl > PRINT_MIN) {
117 fasp_gettime(&end_time);
118 fasp_cputime("PARDISO solver", end_time - start_time);
119 }
120
121 phase = -1; /* Release internal memory */
122 PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, &ddum, ia, ja, &idum, &nrhs, iparm,
123 &msglvl, &ddum, &ddum, &error);
124
125#if DEBUG_MODE
126 printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
127#endif
128
129 return status;
130
131#else
132
133 printf("### ERROR: PARDISO is not available!\n");
134 return ERROR_SOLVER_EXIT;
135
136#endif
137}
138
139#if WITH_PARDISO
153INT fasp_pardiso_factorize(dCSRmat* ptrA, Pardiso_data* pdata, const SHORT prtlvl)
154{
155 INT status = FASP_SUCCESS;
156
157 MKL_INT n = ptrA->col;
158 MKL_INT* ia = ptrA->IA;
159 MKL_INT* ja = ptrA->JA;
160 REAL* a = ptrA->val;
161
162 double ddum; /* Double dummy */
163 MKL_INT nrhs = 1; /* Number of right hand sides */
164 MKL_INT idum; /* Integer dummy */
165 MKL_INT phase, error, msglvl; /* Auxiliary variables */
166
167#if DEBUG_MODE
168 printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
169 printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", ptrA->row, n, ptrA->nnz);
170#endif
171
172 pdata->mtype = 11; /* Real unsymmetric matrix */
173
174 PARDISOINIT(pdata->pt, &(pdata->mtype), pdata->iparm); /* Initialize */
175 pdata->iparm[34] = 1; /* Use 0-based indexing */
176
177 /* Numbers of processors, value of OMP_NUM_THREADS */
178#ifdef _OPENMP
179 pdata->iparm[2] = fasp_get_num_threads();
180#endif
181
182 pdata->maxfct = 1; /* Maximum number of numerical factorizations */
183 pdata->mnum = 1; /* Which factorization to use */
184 error = 0; /* Initialize error flag */
185
186 if (prtlvl > PRINT_MORE) {
187 msglvl = 1;
188 } else {
189 msglvl = 0; /* Do not print statistical information in file */
190 }
191
192 REAL start_time, end_time;
193 fasp_gettime(&start_time);
194
195 phase = 11; /* Reordering and symbolic factorization */
196 PARDISO(pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase, &n, a,
197 ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, &ddum, &ddum, &error);
198 if (error != 0) {
199 printf("### ERROR: %d, %s %d\n", error, __FUNCTION__, __LINE__);
200 printf("### ERROR: PARDISO symbolic factorization failed!\n");
201 exit(ERROR_SOLVER_MISC);
202 }
203
204 phase = 22; /* Numerical factorization */
205 PARDISO(pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase, &n, a,
206 ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, &ddum, &ddum, &error);
207
208 if (error != 0) {
209 printf("### ERROR: %d, %s %d\n", error, __FUNCTION__, __LINE__);
210 printf("### ERROR: PARDISO numeric factorization failed!\n");
211 exit(ERROR_SOLVER_MISC);
212 }
213
214 if (prtlvl > PRINT_MIN) {
215 fasp_gettime(&end_time);
216 fasp_cputime("PARDISO setup", end_time - start_time);
217 }
218
219#if DEBUG_MODE
220 printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
221#endif
222
223 return status;
224}
225
241INT fasp_pardiso_solve(dCSRmat* ptrA, dvector* b, dvector* u, Pardiso_data* pdata,
242 const SHORT prtlvl)
243{
244 INT status = FASP_SUCCESS;
245
246 MKL_INT n = ptrA->col;
247 MKL_INT* ia = ptrA->IA;
248 MKL_INT* ja = ptrA->JA;
249
250 REAL* a = ptrA->val;
251 REAL* f = b->val; /* RHS vector */
252 REAL* x = u->val; /* Solution vector */
253 MKL_INT nrhs = 1; /* Number of right hand sides */
254 MKL_INT idum; /* Integer dummy */
255 MKL_INT phase, error, msglvl; /* Auxiliary variables */
256
257 REAL start_time, end_time;
258 fasp_gettime(&start_time);
259
260 if (prtlvl > PRINT_MORE) {
261 msglvl = 1;
262 } else {
263 msglvl = 0; /* Do not print statistical information in file */
264 }
265
266 phase = 33; /* Back substitution and iterative refinement */
267 PARDISO(pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase, &n, a,
268 ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, f, x, &error);
269
270 if (error != 0) {
271 printf("### ERROR: %d, %s %d\n", error, __FUNCTION__, __LINE__);
272 printf("### ERROR: PARDISO solve failed!\n");
273 exit(ERROR_SOLVER_MISC);
274 }
275
276 if (prtlvl > PRINT_MIN) {
277 fasp_gettime(&end_time);
278 fasp_cputime("PARDISO solve", end_time - start_time);
279 }
280
281#if DEBUG_MODE
282 printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
283#endif
284
285 return status;
286}
287
297INT fasp_pardiso_free_internal_mem(Pardiso_data* pdata)
298{
299 INT status = FASP_SUCCESS;
300
301 MKL_INT* ia = NULL;
302 MKL_INT* ja = NULL;
303
304 double ddum; /* Double dummy */
305 MKL_INT idum; /* Integer dummy */
306 MKL_INT nrhs = 1; /* Number of right hand sides */
307 MKL_INT phase, error, msglvl; /* Auxiliary variables */
308
309#if DEBUG_MODE
310 printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
311#endif
312
313 phase = -1; /* Release internal memory */
314 PARDISO(pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase, &idum,
315 &ddum, ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, &ddum, &ddum, &error);
316
317#if DEBUG_MODE
318 printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
319#endif
320
321 return status;
322}
323
324#endif
325
326/*---------------------------------*/
327/*-- End of File --*/
328/*---------------------------------*/
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_pardiso
INT fasp_solver_pardiso(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Ax=b by PARDISO directly.
Definition: XtrPardiso.c:45
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
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_MORE
#define PRINT_MORE
Definition: fasp_const.h:76
ERROR_SOLVER_EXIT
#define ERROR_SOLVER_EXIT
Definition: fasp_const.h:49
PRINT_MIN
#define PRINT_MIN
Definition: fasp_const.h:74
Pardiso_data
Data for Intel MKL PARDISO interface.
Definition: fasp.h:625
Pardiso_data::pt
void * pt[64]
Internal solver memory pointer.
Definition: fasp.h:628
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