1 |
|
2 |
/******************************************************* |
3 |
* |
4 |
* Copyright (c) 2003-2008 by University of Queensland |
5 |
* Earth Systems Science Computational Center (ESSCC) |
6 |
* http://www.uq.edu.au/esscc |
7 |
* |
8 |
* Primary Business: Queensland, Australia |
9 |
* Licensed under the Open Software License version 3.0 |
10 |
* http://www.opensource.org/licenses/osl-3.0.php |
11 |
* |
12 |
*******************************************************/ |
13 |
|
14 |
|
15 |
/**************************************************************/ |
16 |
|
17 |
/* Paso: interface to the Intel MKL library */ |
18 |
|
19 |
/**************************************************************/ |
20 |
|
21 |
/* Copyrights by ACcESS Australia 2006 */ |
22 |
/* Author: gross@@access.edu.au */ |
23 |
|
24 |
/**************************************************************/ |
25 |
|
26 |
#include "Paso.h" |
27 |
#include "MKL.h" |
28 |
#ifdef _OPENMP |
29 |
#include <omp.h> |
30 |
#endif |
31 |
|
32 |
|
33 |
/**************************************************************/ |
34 |
|
35 |
/* free any extra stuff possibly used by the MKL library */ |
36 |
|
37 |
void Paso_MKL_free(Paso_SystemMatrix* A) { |
38 |
#ifdef MKL |
39 |
index_t i; |
40 |
if (A->solver!=NULL) { |
41 |
_INTEGER_t mtype = MKL_MTYPE_UNSYM; |
42 |
if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM; |
43 |
_INTEGER_t n = A->mainBlock->numRows; |
44 |
_INTEGER_t maxfct=1; /* number of factorizations on the same pattern */ |
45 |
_INTEGER_t mnum =1; /* factoriztion to be handeled in this call */ |
46 |
_INTEGER_t msglvl=0; /* message level */ |
47 |
_INTEGER_t nrhs=1; /* number of right hand sides */ |
48 |
_INTEGER_t idum; /* dummy integer */ |
49 |
_DOUBLE_PRECISION_t ddum; /* dummy float */ |
50 |
_INTEGER_t error=MKL_ERROR_NO; /* error code */ |
51 |
_INTEGER_t iparm[64]; /* parameters */ |
52 |
for (i=0;i<64;++i) iparm[i]=0; |
53 |
|
54 |
_INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY; |
55 |
PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase, |
56 |
&n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs, |
57 |
iparm, &msglvl,&ddum, &ddum, &error); |
58 |
MEMFREE(A->solver); |
59 |
if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in paradiso library failed."); |
60 |
} |
61 |
#endif |
62 |
} |
63 |
/* call the solver: */ |
64 |
|
65 |
void Paso_MKL(Paso_SystemMatrix* A, |
66 |
double* out, |
67 |
double* in, |
68 |
Paso_Options* options, |
69 |
Paso_Performance* pp) { |
70 |
#ifdef MKL |
71 |
double time0, time1; |
72 |
index_t i; |
73 |
|
74 |
if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) { |
75 |
Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1."); |
76 |
return; |
77 |
} |
78 |
options->converged=FALSE; |
79 |
Performance_startMonitor(pp,PERFORMANCE_ALL); |
80 |
_INTEGER_t mtype = MKL_MTYPE_UNSYM; |
81 |
if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM; |
82 |
_INTEGER_t n = A->mainBlock->numRows; |
83 |
_INTEGER_t maxfct=1; /* number of factorizations on the same pattern */ |
84 |
_INTEGER_t mnum =1; /* factoriztion to be handeled in this call */ |
85 |
_INTEGER_t msglvl=0; /* message level */ |
86 |
_INTEGER_t nrhs=1; /* number of right hand sides */ |
87 |
_INTEGER_t idum; /* dummy integer */ |
88 |
_DOUBLE_PRECISION_t ddum; /* dummy float */ |
89 |
_INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION; |
90 |
|
91 |
_INTEGER_t error=MKL_ERROR_NO; /* error code */ |
92 |
_INTEGER_t iparm[64]; /* parameters */ |
93 |
_MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver); |
94 |
/* set iparm */ |
95 |
for (i=0;i<64;++i) iparm[i]=0; |
96 |
iparm[0] = 1; /* no default settings*/ |
97 |
switch (options->reordering) { |
98 |
case PASO_MINIMUM_FILL_IN: |
99 |
iparm[1]=MKL_REORDERING_MINIMUM_DEGREE; |
100 |
break; |
101 |
default: |
102 |
iparm[1]=MKL_REORDERING_NESTED_DISSECTION; |
103 |
break; |
104 |
} |
105 |
#ifdef _OPENMP |
106 |
iparm[2] =omp_get_max_threads(); |
107 |
#else |
108 |
iparm[2] = 1; |
109 |
#endif |
110 |
iparm[5] = 0; /* store solution into output array */ |
111 |
iparm[7] = 2; /* maximum number of refinements */ |
112 |
iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */ |
113 |
iparm[10] = 1; /* rescaling the matrix before factorization started */ |
114 |
iparm[17] =0; /* =-1 report number of non-zeroes */ |
115 |
iparm[18] =0; /* =-1 report flops */ |
116 |
|
117 |
|
118 |
if (pt==NULL) { |
119 |
/* allocate address pointer */ |
120 |
pt=MEMALLOC(64,_MKL_DSS_HANDLE_t); |
121 |
if (Paso_checkPtr(pt)) return; |
122 |
A->solver=(void*) pt; |
123 |
for (i=0;i<64;++i) pt[i]=NULL; |
124 |
time0=Paso_timer(); |
125 |
/* symbolic factorization */ |
126 |
phase = MKL_PHASE_SYMBOLIC_FACTORIZATION; |
127 |
PARDISO (pt, &maxfct, &mnum, &mtype, &phase, |
128 |
&n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs, |
129 |
iparm, &msglvl, in, out, &error); |
130 |
if (error != MKL_ERROR_NO) { |
131 |
Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed."); |
132 |
Paso_MKL_free(A); |
133 |
} else { |
134 |
/* LDU factorization */ |
135 |
phase = MKL_PHASE_FACTORIZATION; |
136 |
PARDISO(pt, &maxfct, &mnum, &mtype, &phase, |
137 |
&n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs, |
138 |
iparm, &msglvl, in, out, &error); |
139 |
if (error != MKL_ERROR_NO) { |
140 |
Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular."); |
141 |
Paso_MKL_free(A); |
142 |
} |
143 |
if (options->verbose) printf("MKL: LDU factorization completed.\n"); |
144 |
} |
145 |
options->set_up_time=Paso_timer()-time0; |
146 |
} else { |
147 |
options->set_up_time=0; |
148 |
} |
149 |
/* forward backward substitution\ */ |
150 |
if (Paso_noError()) { |
151 |
time0=Paso_timer(); |
152 |
phase = MKL_PHASE_SOLVE; |
153 |
PARDISO (pt, &maxfct, &mnum, &mtype, &phase, |
154 |
&n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs, |
155 |
iparm, &msglvl, in, out, &error); |
156 |
if (options->verbose) printf("MKL: solve completed.\n"); |
157 |
if (error != MKL_ERROR_NO) { |
158 |
Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular."); |
159 |
} else { |
160 |
options->residual_norm=0.; |
161 |
options->num_iter=0; |
162 |
options->num_level=0; |
163 |
options->num_inner_iter=0; |
164 |
options->converged=TRUE; |
165 |
} |
166 |
options->time=Paso_timer()-time0 + options->set_up_time; |
167 |
} |
168 |
Performance_stopMonitor(pp,PERFORMANCE_ALL); |
169 |
#else |
170 |
Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble."); |
171 |
#endif |
172 |
} |
173 |
|
174 |
void Paso_MKL_free1(Paso_SparseMatrix* A) { |
175 |
#ifdef MKL |
176 |
index_t i; |
177 |
if (A->solver!=NULL) { |
178 |
_INTEGER_t mtype = MKL_MTYPE_UNSYM; |
179 |
if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM; |
180 |
_INTEGER_t n = A->numRows; |
181 |
_INTEGER_t maxfct=1; /* number of factorizations on the same pattern */ |
182 |
_INTEGER_t mnum =1; /* factoriztion to be handeled in this call */ |
183 |
_INTEGER_t msglvl=0; /* message level */ |
184 |
_INTEGER_t nrhs=1; /* number of right hand sides */ |
185 |
_INTEGER_t idum; /* dummy integer */ |
186 |
_DOUBLE_PRECISION_t ddum; /* dummy float */ |
187 |
_INTEGER_t error=MKL_ERROR_NO; /* error code */ |
188 |
_INTEGER_t iparm[64]; /* parameters */ |
189 |
for (i=0;i<64;++i) iparm[i]=0; |
190 |
|
191 |
_INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY; |
192 |
PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase, |
193 |
&n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs, |
194 |
iparm, &msglvl,&ddum, &ddum, &error); |
195 |
MEMFREE(A->solver); |
196 |
if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in paradiso library failed."); |
197 |
} |
198 |
#endif |
199 |
} |
200 |
/* call the solver: */ |
201 |
|
202 |
void Paso_MKL1(Paso_SparseMatrix* A, |
203 |
double* out, |
204 |
double* in, |
205 |
bool_t verbose) { |
206 |
#ifdef MKL |
207 |
index_t i; |
208 |
|
209 |
if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) { |
210 |
Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1."); |
211 |
return; |
212 |
} |
213 |
_INTEGER_t mtype = MKL_MTYPE_UNSYM; |
214 |
if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM; |
215 |
_INTEGER_t n = A->numRows; |
216 |
_INTEGER_t maxfct=1; /* number of factorizations on the same pattern */ |
217 |
_INTEGER_t mnum =1; /* factoriztion to be handeled in this call */ |
218 |
_INTEGER_t msglvl=0; /* message level */ |
219 |
_INTEGER_t nrhs=1; /* number of right hand sides */ |
220 |
_INTEGER_t idum; /* dummy integer */ |
221 |
_DOUBLE_PRECISION_t ddum; /* dummy float */ |
222 |
_INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION; |
223 |
|
224 |
_INTEGER_t error=MKL_ERROR_NO; /* error code */ |
225 |
_INTEGER_t iparm[64]; /* parameters */ |
226 |
_MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver); |
227 |
/* set iparm */ |
228 |
for (i=0;i<64;++i) iparm[i]=0; |
229 |
iparm[0] = 1; /* no default settings*/ |
230 |
iparm[1]=MKL_REORDERING_MINIMUM_DEGREE; |
231 |
#ifdef _OPENMP |
232 |
iparm[2] = omp_get_max_threads(); |
233 |
#else |
234 |
iparm[2] = 1; |
235 |
#endif |
236 |
|
237 |
iparm[5] = 0; /* store solution into output array */ |
238 |
iparm[7] = 2; /* maximum number of refinements */ |
239 |
iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */ |
240 |
iparm[10] = 1; /* rescaling the matrix before factorization started */ |
241 |
iparm[17] =0; /* =-1 report number of non-zeroes */ |
242 |
iparm[18] =0; /* =-1 report flops */ |
243 |
|
244 |
|
245 |
if (pt==NULL) { |
246 |
/* allocate address pointer */ |
247 |
pt=MEMALLOC(64,_MKL_DSS_HANDLE_t); |
248 |
if (Paso_checkPtr(pt)) return; |
249 |
A->solver=(void*) pt; |
250 |
for (i=0;i<64;++i) pt[i]=NULL; |
251 |
/* symbolic factorization */ |
252 |
phase = MKL_PHASE_SYMBOLIC_FACTORIZATION; |
253 |
PARDISO (pt, &maxfct, &mnum, &mtype, &phase, |
254 |
&n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs, |
255 |
iparm, &msglvl, in, out, &error); |
256 |
if (error != MKL_ERROR_NO) { |
257 |
Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed."); |
258 |
Paso_MKL_free1(A); |
259 |
} else { |
260 |
/* LDU factorization */ |
261 |
phase = MKL_PHASE_FACTORIZATION; |
262 |
PARDISO(pt, &maxfct, &mnum, &mtype, &phase, |
263 |
&n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs, |
264 |
iparm, &msglvl, in, out, &error); |
265 |
if (error != MKL_ERROR_NO) { |
266 |
Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular."); |
267 |
Paso_MKL_free1(A); |
268 |
} |
269 |
if (verbose) printf("MKL: LDU factorization completed.\n"); |
270 |
} |
271 |
} |
272 |
/* forward backward substitution\ */ |
273 |
if (Paso_noError()) { |
274 |
phase = MKL_PHASE_SOLVE; |
275 |
PARDISO (pt, &maxfct, &mnum, &mtype, &phase, |
276 |
&n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs, |
277 |
iparm, &msglvl, in, out, &error); |
278 |
if (verbose) printf("MKL: solve completed.\n"); |
279 |
if (error != MKL_ERROR_NO) { |
280 |
Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular."); |
281 |
} |
282 |
} |
283 |
#else |
284 |
Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble."); |
285 |
#endif |
286 |
} |
287 |
/* |
288 |
* $Log$ |
289 |
* |
290 |
*/ |