/[escript]/branches/doubleplusgood/paso/src/MKL.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/MKL.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4291 - (show annotations)
Thu Mar 7 08:57:58 2013 UTC (6 years, 1 month ago) by jfenwick
File size: 6713 byte(s)


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /************************************************************************************/
18
19 /* Paso: interface to the Intel MKL library */
20
21 /************************************************************************************/
22
23 /* Copyrights by ACcESS Australia 2006 */
24 /* Author: l.gross@uq.edu.au */
25
26 /************************************************************************************/
27
28 #include "Paso.h"
29 #include "MKL.h"
30 #ifdef _OPENMP
31 #include <omp.h>
32 #endif
33
34 /************************************************************************************/
35
36 /* free any extra stuff possibly used by the MKL library */
37
38 void Paso_MKL_free(Paso_SparseMatrix* A) {
39 #ifdef MKL
40 index_t i;
41 if (A!=NULL) {
42
43 if ((A->solver_p!=NULL) && (A->solver_package==PASO_MKL) ) {
44 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
45 _INTEGER_t n = A->numRows;
46 _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
47 _INTEGER_t mnum =1; /* factorization to be handled in this call */
48 _INTEGER_t msglvl=0; /* message level */
49 _INTEGER_t nrhs=1; /* number of right hand sides */
50 _INTEGER_t idum; /* dummy integer */
51 _DOUBLE_PRECISION_t ddum; /* dummy float */
52 _INTEGER_t error=MKL_ERROR_NO; /* error code */
53 _INTEGER_t iparm[64]; /* parameters */
54 _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver_p);
55 _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;
56 for (i=0;i<64;++i) iparm[i]=0;
57
58 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
59 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
60 iparm, &msglvl,&ddum, &ddum, &error);
61 delete[] A->solver_p;
62 A->solver_p=NULL;
63 if (error != MKL_ERROR_NO) Esys_setError(TYPE_ERROR,"memory release in PARDISO library failed.");
64 }
65 }
66 #endif
67 }
68
69 /* call the solver: */
70
71 void Paso_MKL(Paso_SparseMatrix* A,
72 double* out,
73 double* in,
74 index_t reordering,
75 dim_t numRefinements,
76 bool_t verbose)
77 {
78
79 #ifdef MKL
80 double time0;
81 index_t i;
82
83 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
84 _INTEGER_t n = A->numRows;
85 _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
86 _INTEGER_t mnum =1; /* factorization to be handled in this call */
87 _INTEGER_t msglvl=0; /* message level */
88 _INTEGER_t nrhs=1; /* number of right hand sides */
89 _INTEGER_t idum; /* dummy integer */
90 _DOUBLE_PRECISION_t ddum; /* dummy float */
91 _INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
92
93 _INTEGER_t error=MKL_ERROR_NO; /* error code */
94 _INTEGER_t iparm[64]; /* parameters */
95 _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver_p);
96
97 if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
98 Esys_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
99 return;
100 }
101
102 /* set iparm */
103 for (i=0;i<64;++i) iparm[i]=0;
104 iparm[0] = 1; /* no default settings*/
105 switch (reordering) {
106 case PASO_MINIMUM_FILL_IN:
107 iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
108 break;
109 default:
110 iparm[1]=MKL_REORDERING_NESTED_DISSECTION;
111 break;
112 }
113 #ifdef _OPENMP
114 iparm[2] =omp_get_max_threads();
115 #else
116 iparm[2] = 1;
117 #endif
118 iparm[5] = 0; /* store solution into output array */
119 iparm[7] = numRefinements; /* maximum number of refinements */
120 iparm[9] = 13; /* 10**(-iparm[9]) perturbation of pivot elements */
121 iparm[10] = 1; /* rescaling the matrix before factorization started */
122 iparm[17] =0; /* =-1 report number of non-zeroes */
123 iparm[18] =0; /* =-1 report flops */
124
125 if (pt==NULL) {
126 /* allocate address pointer */
127 pt=new _MKL_DSS_HANDLE_t[64];
128 if (Esys_checkPtr(pt)) return;
129 for (i=0;i<64;++i) pt[i]=NULL;
130 A->solver_p=(void*) pt;
131 A->solver_package=PASO_MKL;
132 /* symbolic factorization */
133 phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
134 time0=Esys_timer();
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 (error != MKL_ERROR_NO) {
139 if (verbose) printf("MKL: symbolic factorization failed.\n");
140 Esys_setError(VALUE_ERROR,"symbolic factorization in PARDISO library failed.");
141 Paso_MKL_free(A);
142 } else {
143 /* LDU factorization */
144 phase = MKL_PHASE_FACTORIZATION;
145 PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
146 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
147 iparm, &msglvl, in, out, &error);
148 if (error != MKL_ERROR_NO) {
149 if (verbose) printf("MKL: LDU factorization failed.\n");
150 Esys_setError(ZERO_DIVISION_ERROR,"factorization in PARDISO library failed. Most likely the matrix is singular.");
151 Paso_MKL_free(A);
152 }
153 if (verbose) printf("MKL: LDU factorization completed (time = %e).\n",Esys_timer()-time0);
154 }
155 }
156 /* forward backward substitution\ */
157 if (Esys_noError()) {
158 time0=Esys_timer();
159 phase = MKL_PHASE_SOLVE;
160 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
161 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
162 iparm, &msglvl, in, out, &error);
163 if (verbose) printf("MKL: solve completed.\n");
164 if (error != MKL_ERROR_NO) {
165 if (verbose) printf("MKL: forward/backward substitution failed.\n");
166 Esys_setError(VALUE_ERROR,"forward/backward substitution in PARDISO library failed. Most likely the matrix is singular.");
167 } else {
168 if (verbose) printf("MKL: forward/backward substitution completed (time = %e).\n",Esys_timer()-time0);
169 }
170 }
171 #else
172 Esys_setError(SYSTEM_ERROR,"Paso_MKL: MKL is not available.");
173 #endif
174 }
175

  ViewVC Help
Powered by ViewVC 1.1.26