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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1388 - (show annotations)
Fri Jan 11 07:45:58 2008 UTC (11 years, 7 months ago) by trankine
File MIME type: text/plain
File size: 6416 byte(s)
And get the *(&(*&(* name right
1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 /**************************************************************/
17
18 /* Paso: interface to the Intel MKL library */
19
20 /**************************************************************/
21
22 /* Copyrights by ACcESS Australia 2006 */
23 /* Author: gross@@access.edu.au */
24
25 /**************************************************************/
26
27 #include "Paso.h"
28 #include "MKL.h"
29 #ifdef _OPENMP
30 #include <omp.h>
31 #endif
32
33 /**************************************************************/
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 _INTEGER_t n = A->mainBlock->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 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 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
57 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 Paso_Options* options,
69 Paso_Performance* pp) {
70 #ifdef MKL
71 double time0;
72 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 Performance_startMonitor(pp,PERFORMANCE_ALL);
79 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
80 if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
81 _INTEGER_t n = A->mainBlock->numRows;
82 _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
83 _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
84 _INTEGER_t msglvl=0; /* message level */
85 _INTEGER_t nrhs=1; /* number of right hand sides */
86 _INTEGER_t idum; /* dummy integer */
87 _DOUBLE_PRECISION_t ddum; /* dummy float */
88 _INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
89
90 _INTEGER_t error=MKL_ERROR_NO; /* error code */
91 _INTEGER_t iparm[64]; /* parameters */
92 _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver);
93 /* set iparm */
94 for (i=0;i<64;++i) iparm[i]=0;
95 iparm[0] = 1; /* no default settings*/
96 switch (options->reordering) {
97 case PASO_MINIMUM_FILL_IN:
98 iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
99 break;
100 default:
101 iparm[1]=MKL_REORDERING_NESTED_DISSECTION;
102 break;
103 }
104 #ifdef _OPENMP
105 iparm[2] = omp_get_max_threads();
106 #else
107 iparm[2] = 1;
108 #endif
109 iparm[5] = 0; /* store solution into output array */
110 iparm[7] = 2; /* maximum number of refinements */
111 iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */
112 iparm[10] = 1; /* rescaling the matrix before factorization started */
113 iparm[17] =0; /* =-1 report number of non-zeroes */
114 iparm[18] =0; /* =-1 report flops */
115
116
117 if (pt==NULL) {
118 /* allocate address pointer */
119 pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
120 if (Paso_checkPtr(pt)) return;
121 A->solver=(void*) pt;
122 for (i=0;i<64;++i) pt[i]=NULL;
123 time0=Paso_timer();
124 /* symbolic factorization */
125 phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
126 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
127 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
128 iparm, &msglvl, in, out, &error);
129 if (error != MKL_ERROR_NO) {
130 Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
131 Paso_MKL_free(A);
132 } else {
133 /* LDU factorization */
134 phase = MKL_PHASE_FACTORIZATION;
135 PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
136 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
137 iparm, &msglvl, in, out, &error);
138 if (error != MKL_ERROR_NO) {
139 Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
140 Paso_MKL_free(A);
141 }
142 if (options->verbose) printf("timing MKL: LDU factorization: %.4e sec.\n",Paso_timer()-time0);
143 }
144 }
145 /* forward backward substitution\ */
146 if (Paso_noError()) {
147 time0=Paso_timer();
148 phase = MKL_PHASE_SOLVE;
149 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
150 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
151 iparm, &msglvl, in, out, &error);
152 if (options->verbose) printf("timing MKL: solve: %.4e sec\n",Paso_timer()-time0);
153 if (error != MKL_ERROR_NO) {
154 Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
155 }
156 }
157 Performance_stopMonitor(pp,PERFORMANCE_ALL);
158 #else
159 Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
160 #endif
161 }
162 /*
163 * $Log$
164 *
165 */

  ViewVC Help
Powered by ViewVC 1.1.26