paso/src/MKL.c


/*******************************************************
*
* Copyright (c) 2003-2010 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: l.gross@uq.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_SparseMatrix* A) {
#ifdef MKL
      index_t i;
      if ((A->solver_p!=NULL) && (A->solver_package=PASO_MKL) ) {
          _INTEGER_t mtype = MKL_MTYPE_UNSYM;
          _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_p), &maxfct, &mnum, &mtype, &phase,
                   &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
                   iparm, &msglvl,&ddum, &ddum, &error);
          MEMFREE(A->solver_p);
          if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in PARDISO library failed.");
     }
#endif
}

/*  call the solver: */

void Paso_MKL(Paso_SparseMatrix* A,
              double* out,
              double* in,
          index_t reordering,
          dim_t numRefinements,
          bool_t verbose)
{         

#ifdef MKL
     double time0;
     index_t i;

     if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
        Esys_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
        return;
     }
     
     _INTEGER_t mtype = MKL_MTYPE_UNSYM;
     _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_p);
     /* set iparm */
     for (i=0;i<64;++i) iparm[i]=0;
     iparm[0] = 1; /* no default settings*/
     switch (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] = numRefinements; /* 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 (Esys_checkPtr(pt)) return;
        A->solver_p=(void*) pt;
    A->solver_package=PASO_MKL;
        for (i=0;i<64;++i) pt[i]=NULL;
        time0=Esys_timer();
        /* 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) {
             if (verbose) printf("MKL: symbolic factorization factorization failed.\n");
             Paso_setError(VALUE_ERROR,"symbolic factorization in PARDISO library failed.");
             Paso_MKL_free(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) {
             if (verbose) printf("MKL: LDU factorization failed.\n");
             Paso_setError(ZERO_DIVISION_ERROR,"factorization in PARDISO library failed. Most likely the matrix is singular.");
             Paso_MKL_free(A);
           }
           if (verbose) printf("MKL: LDU factorization completed (time = %e).\n",Paso_timer()-time0);
        }

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