Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
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_types.h"
22#include "mkl_spblas.h"
23#endif
24
25/*---------------------------------*/
26/*-- Public Functions --*/
27/*---------------------------------*/
28
46 dvector *b,
47 dvector *u,
48 const SHORT prtlvl)
49{
50#if WITH_PARDISO
51
52 INT status = FASP_SUCCESS;
53
54 MKL_INT n = ptrA->col;
55 MKL_INT *ia = ptrA->IA;
56 MKL_INT *ja = ptrA->JA;
57 REAL *a = ptrA->val;
58
59 MKL_INT mtype = 11; /* Real unsymmetric matrix */
60 MKL_INT nrhs = 1; /* Number of right hand sides */
61 MKL_INT idum; /* Integer dummy */
62 MKL_INT iparm[64]; /* Pardiso control parameters */
63 MKL_INT maxfct, mnum, phase, error, msglvl; /* Auxiliary variables */
64
65 REAL * f = b->val; /* RHS vector */
66 REAL * x = u->val; /* Solution vector */
67 void *pt[64]; /* Internal solver memory pointer pt */
68 double ddum; /* Double dummy */
69
70#if DEBUG_MODE
71 printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
72 printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", m, n, nnz);
73#endif
74
75 REAL start_time, end_time;
76 fasp_gettime(&start_time);
77
78 PARDISOINIT(pt, &mtype, iparm); /* Initialize */
79 iparm[34] = 1; /* Use 0-based indexing */
80 maxfct = 1; /* Maximum number of numerical factorizations */
81 mnum = 1; /* Which factorization to use */
82 msglvl = 0; /* Do not print statistical information in file */
83 error = 0; /* Initialize error flag */
84
85 phase = 11; /* Reordering and symbolic factorization */
86 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
87 &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
88 if ( error != 0 ) {
89 printf ("### ERROR: Symbolic factorization failed %d!\n", error);
90 exit (1);
91 }
92
93 phase = 22; /* Numerical factorization */
94 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
95 &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
96 if ( error != 0 ) {
97 printf ("\n### ERROR: Numerical factorization failed %d!\n", error);
98 exit (2);
99 }
100
101 phase = 33; /* Back substitution and iterative refinement */
102 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
103 &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, f, x, &error);
104
105 if ( error != 0 ) {
106 printf ("\n### ERROR: Solution failed %d!\n", error);
107 exit (3);
108 }
109
110 if ( prtlvl > PRINT_MIN ) {
111 fasp_gettime(&end_time);
112 fasp_cputime("PARDISO solver", end_time - start_time);
113 }
114
115 phase = -1; /* Release internal memory */
116 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
117 &n, &ddum, ia, ja, &idum, &nrhs,
118 iparm, &msglvl, &ddum, &ddum, &error);
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: PARDISO is not available!\n");
129 return ERROR_SOLVER_EXIT;
130
131#endif
132
133}
134
135#if WITH_PARDISO
149INT fasp_pardiso_factorize (dCSRmat *ptrA,
150 Pardiso_data *pdata,
151 const SHORT prtlvl)
152{
153 INT status = FASP_SUCCESS;
154
155 MKL_INT n = ptrA->col;
156 MKL_INT *ia = ptrA->IA;
157 MKL_INT *ja = ptrA->JA;
158 REAL *a = ptrA->val;
159
160 double ddum; /* Double dummy */
161 MKL_INT nrhs = 1; /* Number of right hand sides */
162 MKL_INT idum; /* Integer dummy */
163 MKL_INT phase, error, msglvl; /* Auxiliary variables */
164
165#if DEBUG_MODE
166 printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
167 printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", m, n, nnz);
168#endif
169
170 pdata->mtype = 11; /* Real unsymmetric matrix */
171
172 PARDISOINIT(pdata->pt, &(pdata->mtype), pdata->iparm); /* Initialize */
173 pdata->iparm[34] = 1; /* Use 0-based indexing */
174
175 /* Numbers of processors, value of OMP_NUM_THREADS */
176#ifdef _OPENMP
177 pdata->iparm[2] = fasp_get_num_threads();
178#endif
179
180 REAL start_time, end_time;
181 fasp_gettime(&start_time);
182
183 pdata->maxfct = 1; /* Maximum number of numerical factorizations */
184 pdata->mnum = 1; /* Which factorization to use */
185 msglvl = 0; /* Do not print statistical information in file */
186 error = 0; /* Initialize error flag */
187
188 phase = 11; /* Reordering and symbolic factorization */
189 PARDISO (pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase, &n,
190 a, ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, &ddum, &ddum, &error);
191 if ( error != 0 ) {
192 printf ("### ERROR: Symbolic factorization failed %d!\n", error);
193 exit (1);
194 }
195
196 phase = 22; /* Numerical factorization */
197 PARDISO (pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase, &n,
198 a, ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, &ddum, &ddum, &error);
199
200 if ( error != 0 ) {
201 printf ("\n### ERROR: Numerical factorization failed %d!\n", error);
202 exit (2);
203 }
204
205 if ( prtlvl > PRINT_MIN ) {
206 fasp_gettime(&end_time);
207 fasp_cputime("PARDISO setup", end_time - start_time);
208 }
209
210#if DEBUG_MODE
211 printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
212#endif
213
214 return status;
215}
216
232INT fasp_pardiso_solve (dCSRmat *ptrA,
233 dvector *b,
234 dvector *u,
235 Pardiso_data *pdata,
236 const SHORT prtlvl)
237{
238 INT status = FASP_SUCCESS;
239
240 MKL_INT n = ptrA->col;
241 MKL_INT *ia = ptrA->IA;
242 MKL_INT *ja = ptrA->JA;
243
244 REAL *a = ptrA->val;
245 REAL * f = b->val; /* RHS vector */
246 REAL * x = u->val; /* Solution vector */
247 MKL_INT nrhs = 1; /* Number of right hand sides */
248 MKL_INT idum; /* Integer dummy */
249 MKL_INT phase, error, msglvl; /* Auxiliary variables */
250
251 REAL start_time, end_time;
252 fasp_gettime(&start_time);
253
254 msglvl = 0; /* Do not print statistical information in file */
255
256 phase = 33; /* Back substitution and iterative refinement */
257 PARDISO (pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase,
258 &n, a, ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, f, x, &error);
259
260 if ( error != 0 ) {
261 printf ("### ERROR: Solution failed %d!\n", error);
262 exit (3);
263 }
264
265 if ( prtlvl > PRINT_MIN ) {
266 fasp_gettime(&end_time);
267 fasp_cputime("PARDISO solve", end_time - start_time);
268 }
269
270#if DEBUG_MODE
271 printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
272#endif
273
274 return status;
275}
276
286INT fasp_pardiso_free_internal_mem (Pardiso_data *pdata)
287{
288 INT status = FASP_SUCCESS;
289
290 MKL_INT *ia = NULL;
291 MKL_INT *ja = NULL;
292
293 double ddum; /* Double dummy */
294 MKL_INT idum; /* Integer dummy */
295 MKL_INT nrhs = 1; /* Number of right hand sides */
296 MKL_INT phase, error, msglvl; /* Auxiliary variables */
297
298 msglvl = 0; /* Do not print statistical information in file */
299
300#if DEBUG_MODE
301 printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
302#endif
303
304 phase = -1; /* Release internal memory */
305 PARDISO (pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase,
306 &idum, &ddum, ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, &ddum,
307 &ddum, &error);
308
309#if DEBUG_MODE
310 printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
311#endif
312
313 return status;
314}
315
316#endif
317
318/*---------------------------------*/
319/*-- End of File --*/
320/*---------------------------------*/
void fasp_cputime(const char *message, const REAL cputime)
Print CPU walltime.
Definition: AuxMessage.c:179
void fasp_gettime(REAL *time)
Get system time.
Definition: AuxTiming.c:36
INT fasp_solver_pardiso(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Ax=b by PARDISO directly.
Definition: XtrPardiso.c:45
Main header file for the FASP project.
#define REAL
Definition: fasp.h:67
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define INT
Definition: fasp.h:64
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:19
#define ERROR_SOLVER_EXIT
Definition: fasp_const.h:49
#define PRINT_MIN
Definition: fasp_const.h:74
Data for Intel MKL PARDISO interface.
Definition: fasp.h:611
void * pt[64]
Internal solver memory pointer.
Definition: fasp.h:614
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:143
INT col
column of matrix A, n
Definition: fasp.h:149
REAL * val
nonzero entries of A
Definition: fasp.h:161
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:155
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:158
Vector with n entries of REAL type.
Definition: fasp.h:346
REAL * val
actual vector entries
Definition: fasp.h:352