/[escript]/trunk/paso/src/MKL.c
ViewVC logotype

Contents of /trunk/paso/src/MKL.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3010 - (show annotations)
Tue Apr 27 05:10:46 2010 UTC (12 years, 3 months ago) by artak
File MIME type: text/plain
File size: 13554 byte(s)
Preparation for AMG on Systems without cut using frobenius norm
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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: l.gross@uq.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 if (options->verbose) printf("MKL: symbolic factorization factorization failed.\n");
132 Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
133 Paso_MKL_free(A);
134 } else {
135 /* LDU factorization */
136 phase = MKL_PHASE_FACTORIZATION;
137 PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
138 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
139 iparm, &msglvl, in, out, &error);
140 if (error != MKL_ERROR_NO) {
141 if (options->verbose) printf("MKL: LDU factorization failed.\n");
142 Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
143 Paso_MKL_free(A);
144 }
145 if (options->verbose) printf("MKL: LDU factorization completed.\n");
146 }
147 options->set_up_time=Paso_timer()-time0;
148 } else {
149 options->set_up_time=0;
150 }
151 /* forward backward substitution\ */
152 if (Paso_noError()) {
153 time0=Paso_timer();
154 phase = MKL_PHASE_SOLVE;
155 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
156 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
157 iparm, &msglvl, in, out, &error);
158 if (options->verbose) printf("MKL: solve completed.\n");
159 if (error != MKL_ERROR_NO) {
160 if (options->verbose) printf("MKL: forward/backward substitution failed.\n");
161 Paso_setError(VALUE_ERROR,"forward/backward substitution in paradiso library failed. Most likely the matrix is singular.");
162 } else {
163 if (options->verbose) printf("MKL: forward/backward substitution completed.\n");
164 options->residual_norm=0.;
165 options->num_iter=0;
166 options->num_level=0;
167 options->num_inner_iter=0;
168 options->converged=TRUE;
169 }
170 options->time=Paso_timer()-time0 + options->set_up_time;
171 }
172 Performance_stopMonitor(pp,PERFORMANCE_ALL);
173 #else
174 Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
175 #endif
176 }
177
178 void Paso_MKL_free1(Paso_SparseMatrix* A) {
179 #ifdef MKL
180 index_t i;
181 if (A->solver!=NULL) {
182 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
183 if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
184 _INTEGER_t n = A->numRows;
185 _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
186 _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
187 _INTEGER_t msglvl=0; /* message level */
188 _INTEGER_t nrhs=1; /* number of right hand sides */
189 _INTEGER_t idum; /* dummy integer */
190 _DOUBLE_PRECISION_t ddum; /* dummy float */
191 _INTEGER_t error=MKL_ERROR_NO; /* error code */
192 _INTEGER_t iparm[64]; /* parameters */
193 for (i=0;i<64;++i) iparm[i]=0;
194
195 _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;
196 PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase,
197 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
198 iparm, &msglvl,&ddum, &ddum, &error);
199 MEMFREE(A->solver);
200 if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in paradiso library failed.");
201 }
202 #endif
203 }
204 /* call the solver: */
205
206 void Paso_MKL1(Paso_SparseMatrix* A,
207 double* out,
208 double* in,
209 bool_t verbose) {
210 #ifdef MKL
211 index_t i;
212
213 /***** FOR AMG **/
214 double **xx;
215 double **bb;
216 Paso_SparseMatrix *block[MAX_BLOCK_SIZE];
217
218 xx=MEMALLOC(A->row_block_size,double*);
219 bb=MEMALLOC(A->row_block_size,double*);
220 if (Paso_checkPtr(xx) || Paso_checkPtr(bb)) return;
221 /****/
222
223
224 if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
225 Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
226 return;
227 }
228 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
229 if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
230 _INTEGER_t n = A->numRows;
231 _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
232 _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
233 _INTEGER_t msglvl=0; /* message level */
234 _INTEGER_t nrhs=1; /* number of right hand sides */
235 _INTEGER_t idum; /* dummy integer */
236 _DOUBLE_PRECISION_t ddum; /* dummy float */
237 _INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
238
239 _INTEGER_t error=MKL_ERROR_NO; /* error code */
240 _INTEGER_t iparm[64]; /* parameters */
241 _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver);
242 /* set iparm */
243 for (i=0;i<64;++i) iparm[i]=0;
244 iparm[0] = 1; /* no default settings*/
245 iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
246 #ifdef _OPENMP
247 iparm[2] = omp_get_max_threads();
248 #else
249 iparm[2] = 1;
250 #endif
251
252 iparm[5] = 0; /* store solution into output array */
253 iparm[7] = 2; /* maximum number of refinements */
254 iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */
255 iparm[10] = 1; /* rescaling the matrix before factorization started */
256 iparm[17] =0; /* =-1 report number of non-zeroes */
257 iparm[18] =0; /* =-1 report flops */
258
259
260 /******* FOR AMG ****/
261 for (i=0;i<A->row_block_size;i++) {
262 xx[i]=MEMALLOC(A->numRows,double);
263 bb[i]=MEMALLOC(A->numRows,double);
264 if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
265 }
266
267 /*#pragma omp parallel for private(i,j) schedule(static)*/
268 for (i=0;i<A->numRows;i++) {
269 for (j=0;j<A->row_block_size;j++) {
270 bb[j][i]=in[A->row_block_size*i+j];
271 xx[j][i]=0;
272 }
273 }
274
275 for (i=0;i<MAX_BLOCK_SIZE;++i) {
276 block[i]=NULL;
277 }
278 /*****************/
279
280 for (i=0;i<A->row_block_size;++i) {
281 block[i]=Paso_SparseMatrix_getBlock(A,i+1);
282
283 if (pt==NULL) {
284 /* allocate address pointer */
285 pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
286 if (Paso_checkPtr(pt)) return;
287 block[i]->solver=(void*) pt;
288 for (i=0;i<64;++i) pt[i]=NULL;
289 /* symbolic factorization */
290 phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
291 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
292 &n, block[i]->val, block[i]->pattern->ptr, block[i]->pattern->index, &idum, &nrhs,
293 iparm, &msglvl, bb[i], xx[i], &error);
294 if (error != MKL_ERROR_NO) {
295 Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
296 Paso_MKL_free1(block[i]);
297 } else {
298 /* LDU factorization */
299 phase = MKL_PHASE_FACTORIZATION;
300 PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
301 &n, block[i]->val, block[i]->pattern->ptr, block[i]->pattern->index, &idum, &nrhs,
302 iparm, &msglvl, bb[i], xx[i], &error);
303 if (error != MKL_ERROR_NO) {
304 Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
305 Paso_MKL_free1(block[i]);
306 }
307 if (verbose) printf("MKL: LDU factorization completed.\n");
308 }
309 }
310 /* forward backward substitution\ */
311 if (Paso_noError()) {
312 phase = MKL_PHASE_SOLVE;
313 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
314 &n, block[i]->val, block[i]->pattern->ptr, block[i]->pattern->index, &idum, &nrhs,
315 iparm, &msglvl, bb[i], xx[i], &error);
316 if (verbose) printf("MKL: solve completed.\n");
317 if (error != MKL_ERROR_NO) {
318 Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
319 }
320 }
321 }
322
323 /***** FOR AMG ********/
324 /*#pragma omp parallel for private(i,j) schedule(static)*/
325 for (i=0;i<A->numRows;i++) {
326 for (j=0;j<A->row_block_size;j++) {
327 out[A->row_block_size*i+j]=xx[j][i];
328 }
329 }
330
331 for (i=0;i<A->row_block_size;i++) {
332 MEMFREE(xx[i]);
333 MEMFREE(bb[i]);
334 Paso_SparseMatrix_free(block[i]);
335 }
336 /*****************/
337 #else
338 Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
339 #endif
340 }
341 /*
342 * $Log$
343 *
344 */

  ViewVC Help
Powered by ViewVC 1.1.26