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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3391 - (show annotations)
Thu Dec 2 02:25:53 2010 UTC (8 years, 8 months ago) by gross
File MIME type: text/plain
File size: 6521 byte(s)
seg fault fixed.
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 /* free any extra stuff possibly used by the MKL library */
35
36 void Paso_MKL_free(Paso_SparseMatrix* A) {
37 #ifdef MKL
38 index_t i;
39 if (A!=NULL) {
40
41 if ((A->solver_p!=NULL) && (A->solver_package==PASO_MKL) ) {
42 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
43 _INTEGER_t n = A->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 _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver_p);
53 _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;
54 for (i=0;i<64;++i) iparm[i]=0;
55
56 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
57 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
58 iparm, &msglvl,&ddum, &ddum, &error);
59 MEMFREE(A->solver_p);
60 A->solver_p=NULL;
61 if (error != MKL_ERROR_NO) Esys_setError(TYPE_ERROR,"memory release in PARDISO library failed.");
62 }
63 }
64 #endif
65 }
66
67 /* call the solver: */
68
69 void Paso_MKL(Paso_SparseMatrix* A,
70 double* out,
71 double* in,
72 index_t reordering,
73 dim_t numRefinements,
74 bool_t verbose)
75 {
76
77 #ifdef MKL
78 double time0;
79 index_t i;
80
81 if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
82 Esys_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
83 return;
84 }
85
86 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
87 _INTEGER_t n = A->numRows;
88 _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
89 _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
90 _INTEGER_t msglvl=0; /* message level */
91 _INTEGER_t nrhs=1; /* number of right hand sides */
92 _INTEGER_t idum; /* dummy integer */
93 _DOUBLE_PRECISION_t ddum; /* dummy float */
94 _INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
95
96 _INTEGER_t error=MKL_ERROR_NO; /* error code */
97 _INTEGER_t iparm[64]; /* parameters */
98 _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver_p);
99 /* set iparm */
100 for (i=0;i<64;++i) iparm[i]=0;
101 iparm[0] = 1; /* no default settings*/
102 switch (reordering) {
103 case PASO_MINIMUM_FILL_IN:
104 iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
105 break;
106 default:
107 iparm[1]=MKL_REORDERING_NESTED_DISSECTION;
108 break;
109 }
110 #ifdef _OPENMP
111 iparm[2] =omp_get_max_threads();
112 #else
113 iparm[2] = 1;
114 #endif
115 iparm[5] = 0; /* store solution into output array */
116 iparm[7] = numRefinements; /* maximum number of refinements */
117 iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */
118 iparm[10] = 1; /* rescaling the matrix before factorization started */
119 iparm[17] =0; /* =-1 report number of non-zeroes */
120 iparm[18] =0; /* =-1 report flops */
121
122 if (pt==NULL) {
123 /* allocate address pointer */
124 pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
125 if (Esys_checkPtr(pt)) return;
126 A->solver_p=(void*) pt;
127 A->solver_package=PASO_MKL;
128 for (i=0;i<64;++i) pt[i]=NULL;
129 time0=Esys_timer();
130 /* symbolic factorization */
131 phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
132 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
133 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
134 iparm, &msglvl, in, out, &error);
135 if (error != MKL_ERROR_NO) {
136 if (verbose) printf("MKL: symbolic factorization factorization failed.\n");
137 Esys_setError(VALUE_ERROR,"symbolic factorization in PARDISO library failed.");
138 Paso_MKL_free(A);
139 } else {
140 /* LDU factorization */
141 phase = MKL_PHASE_FACTORIZATION;
142 PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
143 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
144 iparm, &msglvl, in, out, &error);
145 if (error != MKL_ERROR_NO) {
146 if (verbose) printf("MKL: LDU factorization failed.\n");
147 Esys_setError(ZERO_DIVISION_ERROR,"factorization in PARDISO library failed. Most likely the matrix is singular.");
148 Paso_MKL_free(A);
149 }
150 if (verbose) printf("MKL: LDU factorization completed (time = %e).\n",Esys_timer()-time0);
151 }
152 }
153 /* forward backward substitution\ */
154 if (Esys_noError()) {
155 time0=Esys_timer();
156 phase = MKL_PHASE_SOLVE;
157 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
158 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
159 iparm, &msglvl, in, out, &error);
160 if (verbose) printf("MKL: solve completed.\n");
161 if (error != MKL_ERROR_NO) {
162 if (verbose) printf("MKL: forward/backward substitution failed.\n");
163 Esys_setError(VALUE_ERROR,"forward/backward substitution in PARDISO library failed. Most likely the matrix is singular.");
164 } else {
165 if (verbose) printf("MKL: forward/backward substitution completed (time = %e).\n",Esys_timer()-time0);
166 }
167 }
168 #else
169 Esys_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
170 #endif
171 }

  ViewVC Help
Powered by ViewVC 1.1.26