paso/src/MKL.c


/*******************************************************
*
* Copyright (c) 2003-2008 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


/**************************************************************/

/* Paso: interface to the Intel MKL library */

/**************************************************************/

/* Copyrights by ACcESS Australia 2006 */
/* Author: gross@@access.edu.au */

/**************************************************************/

#include "Paso.h"
#include "MKL.h"
#ifdef _OPENMP
#include <omp.h>
#endif


/**************************************************************/

/*  free any extra stuff possibly used by the MKL library */

void Paso_MKL_free(Paso_SystemMatrix* A) {
#ifdef MKL
      index_t i;
      if (A->solver!=NULL) {
           _INTEGER_t mtype = MKL_MTYPE_UNSYM;
           if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
          _INTEGER_t n = A->mainBlock->numRows;
          _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
          _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
          _INTEGER_t msglvl=0; /* message level */
          _INTEGER_t nrhs=1; /* number of right hand sides */
          _INTEGER_t idum; /* dummy integer */
          _DOUBLE_PRECISION_t ddum; /* dummy float */
          _INTEGER_t error=MKL_ERROR_NO;  /* error code */
          _INTEGER_t iparm[64]; /* parameters */
          for (i=0;i<64;++i) iparm[i]=0;

          _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;
          PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase,
                   &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
                   iparm, &msglvl,&ddum, &ddum, &error);
          MEMFREE(A->solver);
          if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in paradiso library failed.");
     }
#endif
}
/*  call the solver: */

void Paso_MKL(Paso_SystemMatrix* A,
                          double* out,
                          double* in,
                          Paso_Options* options,
                          Paso_Performance* pp) {
#ifdef MKL
     double time0, time1;
     index_t i;

     if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
        Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
        return;
     }
     options->converged=FALSE;
     Performance_startMonitor(pp,PERFORMANCE_ALL);
     _INTEGER_t mtype = MKL_MTYPE_UNSYM;
     if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
     _INTEGER_t n = A->mainBlock->numRows;
     _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
     _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
     _INTEGER_t msglvl=0; /* message level */
     _INTEGER_t nrhs=1; /* number of right hand sides */
     _INTEGER_t idum; /* dummy integer */
     _DOUBLE_PRECISION_t ddum; /* dummy float */
     _INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;

     _INTEGER_t error=MKL_ERROR_NO;  /* error code */
     _INTEGER_t iparm[64]; /* parameters */
     _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver);
     /* set iparm */
     for (i=0;i<64;++i) iparm[i]=0;
     iparm[0] = 1; /* no default settings*/
     switch (options->reordering) {
            case PASO_MINIMUM_FILL_IN:
               iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
               break;
            default:
                iparm[1]=MKL_REORDERING_NESTED_DISSECTION;
                break;
     }
     #ifdef _OPENMP
     iparm[2] =omp_get_max_threads();
     #else
     iparm[2] = 1;
     #endif
     iparm[5] = 0; /* store solution into output array */
     iparm[7] = 2; /* maximum number of refinements */
     iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */
     iparm[10] = 1; /* rescaling the matrix before factorization started */
     iparm[17] =0; /* =-1 report number of non-zeroes */
     iparm[18] =0; /* =-1 report flops */


     if (pt==NULL) {
        /* allocate address pointer */
        pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
        if (Paso_checkPtr(pt)) return;
        A->solver=(void*) pt;
        for (i=0;i<64;++i) pt[i]=NULL;
        time0=Paso_timer();
        /* symbolic factorization */
        phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
        PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
                 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
                 iparm, &msglvl, in, out, &error);
        if (error != MKL_ERROR_NO) {
             Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
             Paso_MKL_free(A);
        } else {
           /* LDU factorization */
           phase = MKL_PHASE_FACTORIZATION;
           PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
                &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
                iparm, &msglvl, in, out, &error);
           if (error != MKL_ERROR_NO) {
             Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
             Paso_MKL_free(A);
           }
           if (options->verbose) printf("MKL: LDU factorization completed.\n");
        }
        options->set_up_time=Paso_timer()-time0;
     } else {
        options->set_up_time=0;
     }
     /* forward backward substitution\ */
     if (Paso_noError())  {
        time0=Paso_timer();
        phase = MKL_PHASE_SOLVE;
        PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
                 &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
                 iparm, &msglvl, in, out, &error);
        if (options->verbose) printf("MKL: solve completed.\n");
        if (error != MKL_ERROR_NO) {
              Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
        } else {
            options->residual_norm=0.;
            options->num_iter=0;
            options->num_level=0;
            options->num_inner_iter=0;
            options->converged=TRUE;
        }
        options->time=Paso_timer()-time0 + options->set_up_time;
     }
     Performance_stopMonitor(pp,PERFORMANCE_ALL);
#else
    Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
#endif
}

void Paso_MKL_free1(Paso_SparseMatrix* A) {
#ifdef MKL
      index_t i;
      if (A->solver!=NULL) {
           _INTEGER_t mtype = MKL_MTYPE_UNSYM;
           if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
          _INTEGER_t n = A->numRows;
          _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
          _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
          _INTEGER_t msglvl=0; /* message level */
          _INTEGER_t nrhs=1; /* number of right hand sides */
          _INTEGER_t idum; /* dummy integer */
          _DOUBLE_PRECISION_t ddum; /* dummy float */
          _INTEGER_t error=MKL_ERROR_NO;  /* error code */
          _INTEGER_t iparm[64]; /* parameters */
          for (i=0;i<64;++i) iparm[i]=0;

          _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;
          PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase,
                   &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
                   iparm, &msglvl,&ddum, &ddum, &error);
          MEMFREE(A->solver);
          if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in paradiso library failed.");
     }
#endif
}
/*  call the solver: */

void Paso_MKL1(Paso_SparseMatrix* A,
                          double* out,
                          double* in,
                          bool_t verbose) {
#ifdef MKL
     index_t i;

     if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
        Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
        return;
     }
     _INTEGER_t mtype = MKL_MTYPE_UNSYM;
     if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
     _INTEGER_t n = A->numRows;
     _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
     _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
     _INTEGER_t msglvl=0; /* message level */
     _INTEGER_t nrhs=1; /* number of right hand sides */
     _INTEGER_t idum; /* dummy integer */
     _DOUBLE_PRECISION_t ddum; /* dummy float */
     _INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;

     _INTEGER_t error=MKL_ERROR_NO;  /* error code */
     _INTEGER_t iparm[64]; /* parameters */
     _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver);
     /* set iparm */
     for (i=0;i<64;++i) iparm[i]=0;
     iparm[0] = 1; /* no default settings*/
     iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
     #ifdef _OPENMP
     iparm[2] = omp_get_max_threads();
     #else
     iparm[2] = 1;
     #endif
     
     iparm[5] = 0; /* store solution into output array */
     iparm[7] = 2; /* maximum number of refinements */
     iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */
     iparm[10] = 1; /* rescaling the matrix before factorization started */
     iparm[17] =0; /* =-1 report number of non-zeroes */
     iparm[18] =0; /* =-1 report flops */


     if (pt==NULL) {
        /* allocate address pointer */
        pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
        if (Paso_checkPtr(pt)) return;
        A->solver=(void*) pt;
        for (i=0;i<64;++i) pt[i]=NULL;
        /* symbolic factorization */
        phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
        PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
                 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
                 iparm, &msglvl, in, out, &error);
        if (error != MKL_ERROR_NO) {
             Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
             Paso_MKL_free1(A);
        } else {
           /* LDU factorization */
           phase = MKL_PHASE_FACTORIZATION;
           PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
                &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
                iparm, &msglvl, in, out, &error);
           if (error != MKL_ERROR_NO) {
             Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
             Paso_MKL_free1(A);
           }
           if (verbose) printf("MKL: LDU factorization completed.\n");
        }
     }
     /* forward backward substitution\ */
     if (Paso_noError())  {
        phase = MKL_PHASE_SOLVE;
        PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
                 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
                 iparm, &msglvl, in, out, &error);
        if (verbose) printf("MKL: solve completed.\n");
        if (error != MKL_ERROR_NO) {
              Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
        }
     }
#else
    Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
#endif
}
/*
 * $Log$
 *
 */
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
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, time1;
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	options->converged=FALSE;
79	Performance_startMonitor(pp,PERFORMANCE_ALL);
80	_INTEGER_t mtype = MKL_MTYPE_UNSYM;
81	if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
82	_INTEGER_t n = A->mainBlock->numRows;
83	_INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
84	_INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
85	_INTEGER_t msglvl=0; /* message level */
86	_INTEGER_t nrhs=1; /* number of right hand sides */
87	_INTEGER_t idum; /* dummy integer */
88	_DOUBLE_PRECISION_t ddum; /* dummy float */
89	_INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
90
91	_INTEGER_t error=MKL_ERROR_NO; /* error code */
92	_INTEGER_t iparm[64]; /* parameters */
93	_MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver);
94	/* set iparm */
95	for (i=0;i<64;++i) iparm[i]=0;
96	iparm[0] = 1; /* no default settings*/
97	switch (options->reordering) {
98	case PASO_MINIMUM_FILL_IN:
99	iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
100	break;
101	default:
102	iparm[1]=MKL_REORDERING_NESTED_DISSECTION;
103	break;
104	}
105	#ifdef _OPENMP
106	iparm[2] =omp_get_max_threads();
107	#else
108	iparm[2] = 1;
109	#endif
110	iparm[5] = 0; /* store solution into output array */
111	iparm[7] = 2; /* maximum number of refinements */
112	iparm[9] = 13; /* 10*(-iparm[9]) preturbation of pivot elements /
113	iparm[10] = 1; /* rescaling the matrix before factorization started */
114	iparm[17] =0; /* =-1 report number of non-zeroes */
115	iparm[18] =0; /* =-1 report flops */
116
117
118	if (pt==NULL) {
119	/* allocate address pointer */
120	pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
121	if (Paso_checkPtr(pt)) return;
122	A->solver=(void*) pt;
123	for (i=0;i<64;++i) pt[i]=NULL;
124	time0=Paso_timer();
125	/* symbolic factorization */
126	phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
127	PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
128	&n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
129	iparm, &msglvl, in, out, &error);
130	if (error != MKL_ERROR_NO) {
131	Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
132	Paso_MKL_free(A);
133	} else {
134	/* LDU factorization */
135	phase = MKL_PHASE_FACTORIZATION;
136	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
137	&n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
138	iparm, &msglvl, in, out, &error);
139	if (error != MKL_ERROR_NO) {
140	Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
141	Paso_MKL_free(A);
142	}
143	if (options->verbose) printf("MKL: LDU factorization completed.\n");
144	}
145	options->set_up_time=Paso_timer()-time0;
146	} else {
147	options->set_up_time=0;
148	}
149	/* forward backward substitution\ */
150	if (Paso_noError()) {
151	time0=Paso_timer();
152	phase = MKL_PHASE_SOLVE;
153	PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
154	&n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
155	iparm, &msglvl, in, out, &error);
156	if (options->verbose) printf("MKL: solve completed.\n");
157	if (error != MKL_ERROR_NO) {
158	Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
159	} else {
160	options->residual_norm=0.;
161	options->num_iter=0;
162	options->num_level=0;
163	options->num_inner_iter=0;
164	options->converged=TRUE;
165	}
166	options->time=Paso_timer()-time0 + options->set_up_time;
167	}
168	Performance_stopMonitor(pp,PERFORMANCE_ALL);
169	#else
170	Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
171	#endif
172	}
173
174	void Paso_MKL_free1(Paso_SparseMatrix* A) {
175	#ifdef MKL
176	index_t i;
177	if (A->solver!=NULL) {
178	_INTEGER_t mtype = MKL_MTYPE_UNSYM;
179	if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
180	_INTEGER_t n = A->numRows;
181	_INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
182	_INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
183	_INTEGER_t msglvl=0; /* message level */
184	_INTEGER_t nrhs=1; /* number of right hand sides */
185	_INTEGER_t idum; /* dummy integer */
186	_DOUBLE_PRECISION_t ddum; /* dummy float */
187	_INTEGER_t error=MKL_ERROR_NO; /* error code */
188	_INTEGER_t iparm[64]; /* parameters */
189	for (i=0;i<64;++i) iparm[i]=0;
190
191	_INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;
192	PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase,
193	&n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
194	iparm, &msglvl,&ddum, &ddum, &error);
195	MEMFREE(A->solver);
196	if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in paradiso library failed.");
197	}
198	#endif
199	}
200	/* call the solver: */
201
202	void Paso_MKL1(Paso_SparseMatrix* A,
203	double* out,
204	double* in,
205	bool_t verbose) {
206	#ifdef MKL
207	index_t i;
208
209	if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
210	Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
211	return;
212	}
213	_INTEGER_t mtype = MKL_MTYPE_UNSYM;
214	if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
215	_INTEGER_t n = A->numRows;
216	_INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
217	_INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
218	_INTEGER_t msglvl=0; /* message level */
219	_INTEGER_t nrhs=1; /* number of right hand sides */
220	_INTEGER_t idum; /* dummy integer */
221	_DOUBLE_PRECISION_t ddum; /* dummy float */
222	_INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
223
224	_INTEGER_t error=MKL_ERROR_NO; /* error code */
225	_INTEGER_t iparm[64]; /* parameters */
226	_MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver);
227	/* set iparm */
228	for (i=0;i<64;++i) iparm[i]=0;
229	iparm[0] = 1; /* no default settings*/
230	iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
231	#ifdef _OPENMP
232	iparm[2] = omp_get_max_threads();
233	#else
234	iparm[2] = 1;
235	#endif
236
237	iparm[5] = 0; /* store solution into output array */
238	iparm[7] = 2; /* maximum number of refinements */
239	iparm[9] = 13; /* 10*(-iparm[9]) preturbation of pivot elements /
240	iparm[10] = 1; /* rescaling the matrix before factorization started */
241	iparm[17] =0; /* =-1 report number of non-zeroes */
242	iparm[18] =0; /* =-1 report flops */
243
244
245	if (pt==NULL) {
246	/* allocate address pointer */
247	pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
248	if (Paso_checkPtr(pt)) return;
249	A->solver=(void*) pt;
250	for (i=0;i<64;++i) pt[i]=NULL;
251	/* symbolic factorization */
252	phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
253	PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
254	&n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
255	iparm, &msglvl, in, out, &error);
256	if (error != MKL_ERROR_NO) {
257	Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
258	Paso_MKL_free1(A);
259	} else {
260	/* LDU factorization */
261	phase = MKL_PHASE_FACTORIZATION;
262	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
263	&n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
264	iparm, &msglvl, in, out, &error);
265	if (error != MKL_ERROR_NO) {
266	Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
267	Paso_MKL_free1(A);
268	}
269	if (verbose) printf("MKL: LDU factorization completed.\n");
270	}
271	}
272	/* forward backward substitution\ */
273	if (Paso_noError()) {
274	phase = MKL_PHASE_SOLVE;
275	PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
276	&n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
277	iparm, &msglvl, in, out, &error);
278	if (verbose) printf("MKL: solve completed.\n");
279	if (error != MKL_ERROR_NO) {
280	Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
281	}
282	}
283	#else
284	Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
285	#endif
286	}
287	/*
288	* $Log$
289	*
290	*/