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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (hide annotations)
Thu Jan 28 02:03:15 2010 UTC (10 years, 1 month ago) by jfenwick
File MIME type: text/plain
File size: 11858 byte(s)
Don't panic.
Updating copyright stamps

1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 dhawcroft 631
14 ksteube 1811
15 gross 425 /**************************************************************/
16    
17     /* Paso: interface to the Intel MKL library */
18    
19     /**************************************************************/
20    
21     /* Copyrights by ACcESS Australia 2006 */
22 jfenwick 2625 /* Author: l.gross@uq.edu.au */
23 gross 425
24     /**************************************************************/
25    
26     #include "Paso.h"
27     #include "MKL.h"
28     #ifdef _OPENMP
29     #include <omp.h>
30     #endif
31    
32 gross 2479
33 gross 425 /**************************************************************/
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 ksteube 1312 _INTEGER_t n = A->mainBlock->numRows;
44 gross 425 _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 ksteube 1312 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
57 gross 425 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 gross 584 Paso_Options* options,
69     Paso_Performance* pp) {
70 gross 425 #ifdef MKL
71 gross 2479 double time0, time1;
72 gross 425 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 gross 2479 options->converged=FALSE;
79 gross 584 Performance_startMonitor(pp,PERFORMANCE_ALL);
80 gross 425 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
81     if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
82 ksteube 1312 _INTEGER_t n = A->mainBlock->numRows;
83 gross 425 _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 artak 2507 iparm[2] =omp_get_max_threads();
107 gross 425 #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 ksteube 1312 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
129 gross 425 iparm, &msglvl, in, out, &error);
130     if (error != MKL_ERROR_NO) {
131 gross 2748 if (options->verbose) printf("MKL: symbolic factorization factorization failed.\n");
132 gross 425 Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
133     Paso_MKL_free(A);
134 gross 584 } else {
135     /* LDU factorization */
136     phase = MKL_PHASE_FACTORIZATION;
137     PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
138 ksteube 1312 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
139 gross 425 iparm, &msglvl, in, out, &error);
140 gross 584 if (error != MKL_ERROR_NO) {
141 gross 2748 if (options->verbose) printf("MKL: LDU factorization failed.\n");
142 gross 727 Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
143 gross 425 Paso_MKL_free(A);
144 gross 584 }
145 artak 2488 if (options->verbose) printf("MKL: LDU factorization completed.\n");
146 gross 425 }
147 gross 2479 options->set_up_time=Paso_timer()-time0;
148     } else {
149     options->set_up_time=0;
150 gross 425 }
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 ksteube 1312 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
157 gross 425 iparm, &msglvl, in, out, &error);
158 artak 2488 if (options->verbose) printf("MKL: solve completed.\n");
159 gross 425 if (error != MKL_ERROR_NO) {
160 gross 2748 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 gross 2479 } else {
163 gross 2748 if (options->verbose) printf("MKL: forward/backward substitution completed.\n");
164 gross 2479 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 gross 425 }
170 gross 2479 options->time=Paso_timer()-time0 + options->set_up_time;
171 gross 425 }
172 gross 584 Performance_stopMonitor(pp,PERFORMANCE_ALL);
173 gross 425 #else
174     Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
175     #endif
176     }
177 artak 2488
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     if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
214     Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
215     return;
216     }
217     _INTEGER_t mtype = MKL_MTYPE_UNSYM;
218     if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
219     _INTEGER_t n = A->numRows;
220     _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
221     _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
222     _INTEGER_t msglvl=0; /* message level */
223     _INTEGER_t nrhs=1; /* number of right hand sides */
224     _INTEGER_t idum; /* dummy integer */
225     _DOUBLE_PRECISION_t ddum; /* dummy float */
226     _INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
227    
228     _INTEGER_t error=MKL_ERROR_NO; /* error code */
229     _INTEGER_t iparm[64]; /* parameters */
230     _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver);
231     /* set iparm */
232     for (i=0;i<64;++i) iparm[i]=0;
233     iparm[0] = 1; /* no default settings*/
234 artak 2507 iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
235 artak 2488 #ifdef _OPENMP
236     iparm[2] = omp_get_max_threads();
237     #else
238     iparm[2] = 1;
239     #endif
240    
241     iparm[5] = 0; /* store solution into output array */
242     iparm[7] = 2; /* maximum number of refinements */
243     iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */
244     iparm[10] = 1; /* rescaling the matrix before factorization started */
245     iparm[17] =0; /* =-1 report number of non-zeroes */
246     iparm[18] =0; /* =-1 report flops */
247    
248    
249     if (pt==NULL) {
250     /* allocate address pointer */
251     pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
252     if (Paso_checkPtr(pt)) return;
253     A->solver=(void*) pt;
254     for (i=0;i<64;++i) pt[i]=NULL;
255     /* symbolic factorization */
256     phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
257     PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
258     &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
259     iparm, &msglvl, in, out, &error);
260     if (error != MKL_ERROR_NO) {
261     Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
262     Paso_MKL_free1(A);
263     } else {
264     /* LDU factorization */
265     phase = MKL_PHASE_FACTORIZATION;
266     PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
267     &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
268     iparm, &msglvl, in, out, &error);
269     if (error != MKL_ERROR_NO) {
270     Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
271     Paso_MKL_free1(A);
272     }
273     if (verbose) printf("MKL: LDU factorization completed.\n");
274     }
275     }
276     /* forward backward substitution\ */
277     if (Paso_noError()) {
278     phase = MKL_PHASE_SOLVE;
279     PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
280     &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
281     iparm, &msglvl, in, out, &error);
282     if (verbose) printf("MKL: solve completed.\n");
283     if (error != MKL_ERROR_NO) {
284     Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
285     }
286     }
287     #else
288     Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
289     #endif
290     }
291 gross 425 /*
292     * $Log$
293     *
294     */

  ViewVC Help
Powered by ViewVC 1.1.26