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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 727 - (hide annotations)
Fri May 12 06:31:06 2006 UTC (13 years, 4 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 gross 425 /* $Id$ */
2    
3 dhawcroft 631 /*
4     ********************************************************************************
5 dhawcroft 633 * Copyright 2006 by ACcESS MNRF *
6 dhawcroft 631 * *
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 gross 425 /**************************************************************/
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 gross 584 Paso_Options* options,
67     Paso_Performance* pp) {
68 gross 425 #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 gross 584 Performance_startMonitor(pp,PERFORMANCE_ALL);
77 gross 425 _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 gross 584 } else {
131     /* LDU factorization */
132     phase = MKL_PHASE_FACTORIZATION;
133     PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
134 gross 425 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
135     iparm, &msglvl, in, out, &error);
136 gross 584 if (error != MKL_ERROR_NO) {
137 gross 727 Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
138 gross 425 Paso_MKL_free(A);
139 gross 584 }
140     if (options->verbose) printf("timing MKL: LDU factorization: %.4e sec.\n",Paso_timer()-time0);
141 gross 425 }
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 gross 727 Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
153 gross 425 }
154     }
155 gross 584 Performance_stopMonitor(pp,PERFORMANCE_ALL);
156 gross 425 #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