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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 425 - (show annotations)
Tue Jan 10 04:10:39 2006 UTC (13 years, 7 months ago) by gross
File MIME type: text/plain
File size: 5644 byte(s)
The sparse solver can be called by paso now. 

the building has been change to reduce some code redundancy:
now all scons SCscripts are importing scons/esys_options.py which
imports platform specific settings. 



1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Paso: interface to the Intel MKL library */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia 2006 */
10 /* Author: gross@@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Paso.h"
15 #include "MKL.h"
16 #ifdef _OPENMP
17 #include <omp.h>
18 #endif
19
20 /**************************************************************/
21
22 /* free any extra stuff possibly used by the MKL library */
23
24 void Paso_MKL_free(Paso_SystemMatrix* A) {
25 #ifdef MKL
26 index_t i;
27 if (A->solver!=NULL) {
28 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
29 if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
30 _INTEGER_t n = A->num_rows;
31 _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
32 _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
33 _INTEGER_t msglvl=0; /* message level */
34 _INTEGER_t nrhs=1; /* number of right hand sides */
35 _INTEGER_t idum; /* dummy integer */
36 _DOUBLE_PRECISION_t ddum; /* dummy float */
37 _INTEGER_t error=MKL_ERROR_NO; /* error code */
38 _INTEGER_t iparm[64]; /* parameters */
39 for (i=0;i<64;++i) iparm[i]=0;
40
41 _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;
42 PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase,
43 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
44 iparm, &msglvl,&ddum, &ddum, &error);
45 MEMFREE(A->solver);
46 if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in paradiso library failed.");
47 }
48 #endif
49 }
50 /* call the solver: */
51
52 void Paso_MKL(Paso_SystemMatrix* A,
53 double* out,
54 double* in,
55 Paso_Options* options) {
56 #ifdef MKL
57 double time0;
58 index_t i;
59
60 if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
61 Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
62 return;
63 }
64 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
65 if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
66 _INTEGER_t n = A->num_rows;
67 _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
68 _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
69 _INTEGER_t msglvl=0; /* message level */
70 _INTEGER_t nrhs=1; /* number of right hand sides */
71 _INTEGER_t idum; /* dummy integer */
72 _DOUBLE_PRECISION_t ddum; /* dummy float */
73 _INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
74
75 _INTEGER_t error=MKL_ERROR_NO; /* error code */
76 _INTEGER_t iparm[64]; /* parameters */
77 _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver);
78 /* set iparm */
79 for (i=0;i<64;++i) iparm[i]=0;
80 iparm[0] = 1; /* no default settings*/
81 switch (options->reordering) {
82 case PASO_MINIMUM_FILL_IN:
83 iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
84 break;
85 default:
86 iparm[1]=MKL_REORDERING_NESTED_DISSECTION;
87 break;
88 }
89 #ifdef _OPENMP
90 iparm[2] = omp_get_max_threads();
91 #else
92 iparm[2] = 1;
93 #endif
94 iparm[5] = 0; /* store solution into output array */
95 iparm[7] = 2; /* maximum number of refinements */
96 iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */
97 iparm[10] = 1; /* rescaling the matrix before factorization started */
98 iparm[17] =0; /* =-1 report number of non-zeroes */
99 iparm[18] =0; /* =-1 report flops */
100
101
102 if (pt==NULL) {
103 /* allocate address pointer */
104 pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
105 if (Paso_checkPtr(pt)) return;
106 A->solver=(void*) pt;
107 for (i=0;i<64;++i) pt[i]=NULL;
108 time0=Paso_timer();
109 /* symbolic factorization */
110 phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
111 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
112 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
113 iparm, &msglvl, in, out, &error);
114 if (error != MKL_ERROR_NO) {
115 Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
116 Paso_MKL_free(A);
117 return;
118 }
119 /* LDU factorization */
120 phase = MKL_PHASE_FACTORIZATION;
121 PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
122 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
123 iparm, &msglvl, in, out, &error);
124 if (error != MKL_ERROR_NO) {
125 Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed.");
126 Paso_MKL_free(A);
127 return;
128 }
129 if (options->verbose) printf("timing MKL: LDU factorization: %.4e sec.\n",Paso_timer()-time0);
130 }
131 /* forward backward substitution\ */
132 if (Paso_noError()) {
133 time0=Paso_timer();
134 phase = MKL_PHASE_SOLVE;
135 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
136 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
137 iparm, &msglvl, in, out, &error);
138 if (options->verbose) printf("timing MKL: solve: %.4e sec\n",Paso_timer()-time0);
139 if (error != MKL_ERROR_NO) {
140 Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed.");
141 return;
142 }
143 }
144 #else
145 Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
146 #endif
147 }
148 /*
149 * $Log$
150 *
151 */

  ViewVC Help
Powered by ViewVC 1.1.26