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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 584 - (hide annotations)
Thu Mar 9 23:03:38 2006 UTC (13 years, 6 months ago) by gross
File MIME type: text/plain
File size: 5764 byte(s)
eigenvalues: compiles and passes tests on altix now
1 gross 425 /* $Id$ */
2    
3     /**************************************************************/
4    
5     /* Paso: interface to the Intel MKL library */
6    
7     /**************************************************************/
8    
9     /* Copyrights by ACcESS Australia 2006 */
10     /* Author: gross@@access.edu.au */
11    
12     /**************************************************************/
13    
14     #include "Paso.h"
15     #include "MKL.h"
16     #ifdef _OPENMP
17     #include <omp.h>
18     #endif
19    
20     /**************************************************************/
21    
22     /* free any extra stuff possibly used by the MKL library */
23    
24     void Paso_MKL_free(Paso_SystemMatrix* A) {
25     #ifdef MKL
26     index_t i;
27     if (A->solver!=NULL) {
28     _INTEGER_t mtype = MKL_MTYPE_UNSYM;
29     if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
30     _INTEGER_t n = A->num_rows;
31     _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
32     _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
33     _INTEGER_t msglvl=0; /* message level */
34     _INTEGER_t nrhs=1; /* number of right hand sides */
35     _INTEGER_t idum; /* dummy integer */
36     _DOUBLE_PRECISION_t ddum; /* dummy float */
37     _INTEGER_t error=MKL_ERROR_NO; /* error code */
38     _INTEGER_t iparm[64]; /* parameters */
39     for (i=0;i<64;++i) iparm[i]=0;
40    
41     _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;
42     PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase,
43     &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
44     iparm, &msglvl,&ddum, &ddum, &error);
45     MEMFREE(A->solver);
46     if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in paradiso library failed.");
47     }
48     #endif
49     }
50     /* call the solver: */
51    
52     void Paso_MKL(Paso_SystemMatrix* A,
53     double* out,
54     double* in,
55 gross 584 Paso_Options* options,
56     Paso_Performance* pp) {
57 gross 425 #ifdef MKL
58     double time0;
59     index_t i;
60    
61     if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
62     Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
63     return;
64     }
65 gross 584 Performance_startMonitor(pp,PERFORMANCE_ALL);
66 gross 425 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
67     if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
68     _INTEGER_t n = A->num_rows;
69     _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
70     _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
71     _INTEGER_t msglvl=0; /* message level */
72     _INTEGER_t nrhs=1; /* number of right hand sides */
73     _INTEGER_t idum; /* dummy integer */
74     _DOUBLE_PRECISION_t ddum; /* dummy float */
75     _INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
76    
77     _INTEGER_t error=MKL_ERROR_NO; /* error code */
78     _INTEGER_t iparm[64]; /* parameters */
79     _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver);
80     /* set iparm */
81     for (i=0;i<64;++i) iparm[i]=0;
82     iparm[0] = 1; /* no default settings*/
83     switch (options->reordering) {
84     case PASO_MINIMUM_FILL_IN:
85     iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
86     break;
87     default:
88     iparm[1]=MKL_REORDERING_NESTED_DISSECTION;
89     break;
90     }
91     #ifdef _OPENMP
92     iparm[2] = omp_get_max_threads();
93     #else
94     iparm[2] = 1;
95     #endif
96     iparm[5] = 0; /* store solution into output array */
97     iparm[7] = 2; /* maximum number of refinements */
98     iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */
99     iparm[10] = 1; /* rescaling the matrix before factorization started */
100     iparm[17] =0; /* =-1 report number of non-zeroes */
101     iparm[18] =0; /* =-1 report flops */
102    
103    
104     if (pt==NULL) {
105     /* allocate address pointer */
106     pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
107     if (Paso_checkPtr(pt)) return;
108     A->solver=(void*) pt;
109     for (i=0;i<64;++i) pt[i]=NULL;
110     time0=Paso_timer();
111     /* symbolic factorization */
112     phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
113     PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
114     &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
115     iparm, &msglvl, in, out, &error);
116     if (error != MKL_ERROR_NO) {
117     Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
118     Paso_MKL_free(A);
119 gross 584 } else {
120     /* LDU factorization */
121     phase = MKL_PHASE_FACTORIZATION;
122     PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
123 gross 425 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
124     iparm, &msglvl, in, out, &error);
125 gross 584 if (error != MKL_ERROR_NO) {
126 gross 425 Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed.");
127     Paso_MKL_free(A);
128 gross 584 }
129     if (options->verbose) printf("timing MKL: LDU factorization: %.4e sec.\n",Paso_timer()-time0);
130 gross 425 }
131     }
132     /* forward backward substitution\ */
133     if (Paso_noError()) {
134     time0=Paso_timer();
135     phase = MKL_PHASE_SOLVE;
136     PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
137     &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
138     iparm, &msglvl, in, out, &error);
139     if (options->verbose) printf("timing MKL: solve: %.4e sec\n",Paso_timer()-time0);
140     if (error != MKL_ERROR_NO) {
141     Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed.");
142     }
143     }
144 gross 584 Performance_stopMonitor(pp,PERFORMANCE_ALL);
145 gross 425 #else
146     Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
147     #endif
148     }
149     /*
150     * $Log$
151     *
152     */

  ViewVC Help
Powered by ViewVC 1.1.26