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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1384 - (hide annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 8 months ago) by phornby
Original Path: temp_trunk_copy/paso/src/MKL.c
File MIME type: text/plain
File size: 6416 byte(s)
Make a temp copy of the trunk before checking in the windows changes


1 ksteube 1312
2 gross 425 /* $Id$ */
3    
4 ksteube 1312 /*******************************************************
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 dhawcroft 631
16 gross 425 /**************************************************************/
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 ksteube 1312 _INTEGER_t n = A->mainBlock->numRows;
44 gross 425 _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 ksteube 1312 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
57 gross 425 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 gross 584 Paso_Options* options,
69     Paso_Performance* pp) {
70 gross 425 #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 gross 584 Performance_startMonitor(pp,PERFORMANCE_ALL);
79 gross 425 _INTEGER_t mtype = MKL_MTYPE_UNSYM;
80     if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
81 ksteube 1312 _INTEGER_t n = A->mainBlock->numRows;
82 gross 425 _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 ksteube 1312 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
128 gross 425 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 gross 584 } else {
133     /* LDU factorization */
134     phase = MKL_PHASE_FACTORIZATION;
135     PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
136 ksteube 1312 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
137 gross 425 iparm, &msglvl, in, out, &error);
138 gross 584 if (error != MKL_ERROR_NO) {
139 gross 727 Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
140 gross 425 Paso_MKL_free(A);
141 gross 584 }
142     if (options->verbose) printf("timing MKL: LDU factorization: %.4e sec.\n",Paso_timer()-time0);
143 gross 425 }
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 ksteube 1312 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
151 gross 425 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 gross 727 Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
155 gross 425 }
156     }
157 gross 584 Performance_stopMonitor(pp,PERFORMANCE_ALL);
158 gross 425 #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