Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
interface_pardiso.c
Go to the documentation of this file.
1 
8 #include <time.h>
9 
10 #include "fasp.h"
11 #include "fasp_functs.h"
12 
13 #if WITH_PARDISO
14 #include "mkl_pardiso.h"
15 #include "mkl_types.h"
16 #include "mkl_spblas.h"
17 #endif
18 
19 /*---------------------------------*/
20 /*-- Public Functions --*/
21 /*---------------------------------*/
22 
39  dvector *b,
40  dvector *u,
41  const SHORT prtlvl)
42 {
43 #if WITH_PARDISO
44  INT status = FASP_SUCCESS;
45 
46  MKL_INT n = ptrA->col;
47  MKL_INT m = ptrA->row;
48  MKL_INT nnz = ptrA->nnz;
49  MKL_INT *ia = ptrA->IA;
50  MKL_INT *ja = ptrA->JA;
51  REAL *a = ptrA->val;
52 
53  MKL_INT mtype = 11; /* Real unsymmetric matrix */
54  MKL_INT nrhs = 1; /* Number of right hand sides. */
55  MKL_INT idum; /* Integer dummy. */
56  MKL_INT iparm[64]; /* Pardiso control parameters. */
57  MKL_INT maxfct, mnum, phase, error, msglvl; /* Auxiliary variables. */
58 
59  REAL * f = b->val; /* RHS vector. */
60  REAL * x = u->val; /* Solution vector. */
61  void *pt[64]; /* Internal solver memory pointer pt. */
62  double ddum; /* Double dummy */
63 
64 #if DEBUG_MODE
65  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
66  printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n", m, n, nnz);
67 #endif
68 
69  PARDISOINIT(pt, &mtype, iparm); /* Initialize. */
70  iparm[34] = 1; /* Use 0-based indexing. */
71  maxfct = 1; /* Maximum number of numerical factorizations. */
72  mnum = 1; /* Which factorization to use. */
73  msglvl = 0; /* Do not print statistical information in file */
74  error = 0; /* Initialize error flag */
75 
76  clock_t start_time = clock();
77 
78  phase = 11; /* Reordering and symbolic factorization. */
79  PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
80  &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
81  if ( error != 0 ) {
82  printf ("### ERROR: Symbolic factorization failed %d!\n", error);
83  exit (1);
84  }
85 
86  phase = 22; /* Numerical factorization. */
87  PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
88  &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
89  if ( error != 0 ) {
90  printf ("\n### ERROR: Numerical factorization failed %d!\n", error);
91  exit (2);
92  }
93 
94  phase = 33; /* Back substitution and iterative refinement. */
95  PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
96  &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, f, x, &error);
97 
98  if ( error != 0 ) {
99  printf ("\n### ERROR: Solution failed %d!\n", error);
100  exit (3);
101  }
102 
103  if ( prtlvl > PRINT_MIN ) {
104  clock_t end_time = clock();
105  double solve_time = (double)(end_time - start_time)/(double)(CLOCKS_PER_SEC);
106  printf("PARDISO costs %f seconds.\n", solve_time);
107  }
108 
109  phase = -1; /* Release internal memory. */
110  PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
111  &n, &ddum, ia, ja, &idum, &nrhs,
112  iparm, &msglvl, &ddum, &ddum, &error);
113 
114 #if DEBUG_MODE
115  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
116 #endif
117 
118  return status;
119 
120 #else
121 
122 #if DEBUG_MODE
123  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
124 #endif
125 
126  printf("### ERROR: PARDISO is not available!\n");
127  return ERROR_SOLVER_EXIT;
128 
129 #endif
130 
131 }
132 
133 #if WITH_PARDISO
134 
147 INT fasp_pardiso_factorize (dCSRmat *ptrA,
148  Pardiso_data *pdata,
149  const SHORT prtlvl)
150 {
151  INT status = FASP_SUCCESS;
152 
153  MKL_INT n = ptrA->col;
154  MKL_INT m = ptrA->row;
155  MKL_INT nnz = ptrA->nnz;
156  MKL_INT *ia = ptrA->IA;
157  MKL_INT *ja = ptrA->JA;
158  REAL *a = ptrA->val;
159 
160  pdata->mtype = 11; /* Real unsymmetric matrix */
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", m, n, nnz);
170 #endif
171 
172  PARDISOINIT(pdata->pt, &(pdata->mtype), pdata->iparm); /* Initialize. */
173  pdata->iparm[34] = 1; /* Use 0-based indexing. */
174  pdata->maxfct = 1; /* Maximum number of numerical factorizations. */
175  pdata->mnum = 1; /* Which factorization to use. */
176  msglvl = 0; /* Do not print statistical information in file */
177  error = 0; /* Initialize error flag */
178 
179  clock_t start_time = clock();
180 
181  phase = 11; /* Reordering and symbolic factorization. */
182  PARDISO (pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase,
183  &n, a, ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, &ddum, &ddum, &error);
184  if ( error != 0 ) {
185  printf ("\n### ERROR: Symbolic factorization failed %d!\n", error);
186  exit (1);
187  }
188 
189  phase = 22; /* Numerical factorization. */
190  PARDISO (pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase,
191  &n, a, ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, &ddum, &ddum, &error);
192  if ( error != 0 )
193  {
194  printf ("\n### ERROR: Numerical factorization failed %d!\n", error);
195  exit (2);
196  }
197 
198  if ( prtlvl > PRINT_MIN ) {
199  clock_t end_time = clock();
200  double solve_time = (double)(end_time - start_time)/(double)(CLOCKS_PER_SEC);
201  printf("PARDISO factoization costs %f seconds.\n", solve_time);
202  }
203 
204 #if DEBUG_MODE
205  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
206 #endif
207  return status;
208 }
209 
225 INT fasp_pardiso_solve (dCSRmat *ptrA,
226  dvector *b,
227  dvector *u,
228  Pardiso_data *pdata,
229  const SHORT prtlvl)
230 {
231  INT status = FASP_SUCCESS;
232 
233  MKL_INT n = ptrA->col;
234  MKL_INT m = ptrA->row;
235  MKL_INT nnz = ptrA->nnz;
236  MKL_INT *ia = ptrA->IA;
237  MKL_INT *ja = ptrA->JA;
238 
239  REAL *a = ptrA->val;
240  REAL * f = b->val; /* RHS vector. */
241  REAL * x = u->val; /* Solution vector. */
242  MKL_INT mtype = 11; /* Real unsymmetric matrix */
243  MKL_INT nrhs = 1; /* Number of right hand sides. */
244  MKL_INT idum; /* Integer dummy. */
245  MKL_INT phase, error, msglvl; /* Auxiliary variables. */
246 
247  msglvl = 0; /* Do not print statistical information in file */
248 
249  clock_t start_time = clock();
250 
251  phase = 33; /* Back substitution and iterative refinement. */
252  PARDISO (pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase,
253  &n, a, ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, f, x, &error);
254 
255  if ( error != 0 ) {
256  printf ("### ERROR: Solution failed %d!\n", error);
257  exit (3);
258  }
259 
260  if ( prtlvl > PRINT_MIN ) {
261  clock_t end_time = clock();
262  double solve_time = (double)(end_time - start_time)/(double)(CLOCKS_PER_SEC);
263  printf("PARDISO costs %f seconds.\n", solve_time);
264  }
265 
266 #if DEBUG_MODE
267  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
268 #endif
269 
270  return status;
271 }
272 
282 INT fasp_pardiso_free_internal_mem (Pardiso_data *pdata)
283 {
284  INT status = FASP_SUCCESS;
285 
286  MKL_INT n, m, nnz;
287  MKL_INT *ia = NULL;
288  MKL_INT *ja = NULL;
289  REAL *a = NULL;
290 
291  REAL * f = NULL; /* RHS vector. */
292  REAL * x = NULL; /* Solution vector. */
293  double ddum; /* Double dummy */
294 
295  MKL_INT idum; /* Integer dummy. */
296  MKL_INT nrhs = 1; /* Number of right hand sides. */
297  MKL_INT phase, error, msglvl; /* Auxiliary variables. */
298 
299  msglvl = 0; /* Do not print statistical information in file */
300 
301 #if DEBUG_MODE
302  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
303 #endif
304 
305  phase = -1; /* Release internal memory. */
306  PARDISO (pdata->pt, &(pdata->maxfct), &(pdata->mnum), &(pdata->mtype), &phase,
307  &idum, a, ia, ja, &idum, &nrhs, pdata->iparm, &msglvl, &ddum, &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 /*---------------------------------*/
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
#define REAL
Definition: fasp.h:67
REAL * val
nonzero entries of A
Definition: fasp.h:166
REAL * val
actual vector entries
Definition: fasp.h:348
Vector with n entries of REAL type.
Definition: fasp.h:342
#define INT
Definition: fasp.h:64
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
Parameters for Intel MKL PARDISO interface.
Definition: fasp.h:477
INT nnz
number of nonzero entries
Definition: fasp.h:157
INT col
column of matrix A, n
Definition: fasp.h:154
Main header file for FASP.
INT fasp_solver_pardiso(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Ax=b by PARDISO directly. Each row of A should be in ascending order w.r.t. column indices...
INT row
row number of matrix A, m
Definition: fasp.h:151
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
#define ERROR_SOLVER_EXIT
Definition: fasp_const.h:55
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160
#define PRINT_MIN
Definition: fasp_const.h:80