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!=NULL) {

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