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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 727 - (show annotations)
Fri May 12 06:31:06 2006 UTC (13 years, 3 months ago) by gross
File MIME type: text/plain
File size: 6481 byte(s)
In case of an error in paso a seg fault occured. this was caused by
inproper deallocation of memory in this case. this problem is fixed now.


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

  ViewVC Help
Powered by ViewVC 1.1.26