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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years ago) by ksteube
File MIME type: text/plain
File size: 6381 byte(s)
Copyright updated in all files

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

  ViewVC Help
Powered by ViewVC 1.1.26