paso/src/UMFPACK.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 UMFPACK library */

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

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

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

#include "Paso.h"
#include "UMFPACK.h"

#ifdef _OPENMP
#include <omp.h>
#endif

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

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

void Paso_UMFPACK_free(Paso_SparseMatrix* A) {
     Paso_UMFPACK_Handler* pt =NULL;
      if ( (A->solver_p!=NULL) && (A->solver_package == PASO_UMFPACK) ) {
           pt=(Paso_UMFPACK_Handler*)(A->solver_p);
       #ifdef UMFPACK
       umfpack_di_free_symbolic(&(pt->symbolic));
       umfpack_di_free_numeric(&(pt->numeric));
       #endif
       MEMFREE(pt);
           A->solver_p=NULL;
     }
}


/*  call the solver: */

void Paso_UMFPACK(Paso_SparseMatrix* A,
                          double* out,
                          double* in,
                  dim_t numRefinements,
                  bool_t verbose) 
{
#ifdef UMFPACK
     double time0;
     Paso_UMFPACK_Handler* pt = NULL;
     double control[UMFPACK_CONTROL], info[UMFPACK_INFO];
     int error = UMFPACK_OK;
     
     if (! (A->type & (MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC)) ) {
        Esys_setError(TYPE_ERROR,"Paso_UMFPACK: UMFPACK requires CSC format with index offset 1 and block size 1.");
        return;
     }
     
     pt = (Paso_UMFPACK_Handler *)(A->solver_p);
     umfpack_di_defaults(control);
     
     if (pt==NULL) {
    int n = A->numRows;
    pt=(MEMALLOC(1,Paso_UMFPACK_Handler));
    time0=Esys_timer();
    if (Esys_checkPtr(pt)) return;
    A->solver_p = (void*) pt;
    /* call LDU symbolic factorization: */
    error=umfpack_di_symbolic(n,n,A->pattern->ptr,A->pattern->index,A->val,&(pt->symbolic),control,info);
    if (error != UMFPACK_OK) {
       if (verbose) printf("UMFPACK: symbolic factorization factorization failed.\n");
       Esys_setError(VALUE_ERROR,"UMFPACK: symbolic factorization failed.");
       return;
    } else {
       /* call LDU factorization: */
       error= umfpack_di_numeric(A->pattern->ptr,A->pattern->index,A->val,pt->symbolic,&(pt->numeric),control,info);
       if (error != UMFPACK_OK) {
          if (verbose) printf("UMFPACK: LDU factorization failed.\n");
          Esys_setError(ZERO_DIVISION_ERROR,"factorization failed. Most likely the matrix is singular.");
          return;
       }
       if (verbose) printf("UMFPACK: LDU factorization completed (time = %e).\n",Esys_timer()-time0);
    }
     }
     if (Esys_noError())  {
    /* call forward backward substitution: */
    control[UMFPACK_IRSTEP]=numRefinements; /* number of refinement steps */
    time0=Esys_timer();
    error=umfpack_di_solve(UMFPACK_A,A->pattern->ptr,A->pattern->index,A->val,out,in,pt->numeric,control,info);
    if (error != UMFPACK_OK) {
       if (verbose) printf("UMFPACK: forward/backward substitution failed.\n");
       Esys_setError(VALUE_ERROR,"forward/backward substition failed. Most likely the matrix is singular.");
       return;
    }
    if (verbose) printf("UMFPACK: forward/backward substitution completed (time = %e).\n",Esys_timer()-time0);
     }
     #else
     Esys_setError(SYSTEM_ERROR,"Paso_UMFPACK:UMFPACK 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 UMFPACK 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 "UMFPACK.h"
28
29	#ifdef _OPENMP
30	#include <omp.h>
31	#endif
32
33	/**************************************************************/
34
35	/* free any extra stuff possibly used by the UMFPACK library */
36
37	void Paso_UMFPACK_free(Paso_SparseMatrix* A) {
38	Paso_UMFPACK_Handler* pt =NULL;
39	if ( (A->solver_p!=NULL) && (A->solver_package == PASO_UMFPACK) ) {
40	pt=(Paso_UMFPACK_Handler*)(A->solver_p);
41	#ifdef UMFPACK
42	umfpack_di_free_symbolic(&(pt->symbolic));
43	umfpack_di_free_numeric(&(pt->numeric));
44	#endif
45	MEMFREE(pt);
46	A->solver_p=NULL;
47	}
48	}
49
50
51
52	/* call the solver: */
53
54	void Paso_UMFPACK(Paso_SparseMatrix* A,
55	double* out,
56	double* in,
57	dim_t numRefinements,
58	bool_t verbose)
59	{
60	#ifdef UMFPACK
61	double time0;
62	Paso_UMFPACK_Handler* pt = NULL;
63	double control[UMFPACK_CONTROL], info[UMFPACK_INFO];
64	int error = UMFPACK_OK;
65
66	if (! (A->type & (MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC)) ) {
67	Esys_setError(TYPE_ERROR,"Paso_UMFPACK: UMFPACK requires CSC format with index offset 1 and block size 1.");
68	return;
69	}
70
71	pt = (Paso_UMFPACK_Handler *)(A->solver_p);
72	umfpack_di_defaults(control);
73
74	if (pt==NULL) {
75	int n = A->numRows;
76	pt=(MEMALLOC(1,Paso_UMFPACK_Handler));
77	time0=Esys_timer();
78	if (Esys_checkPtr(pt)) return;
79	A->solver_p = (void*) pt;
80	/* call LDU symbolic factorization: */
81	error=umfpack_di_symbolic(n,n,A->pattern->ptr,A->pattern->index,A->val,&(pt->symbolic),control,info);
82	if (error != UMFPACK_OK) {
83	if (verbose) printf("UMFPACK: symbolic factorization factorization failed.\n");
84	Esys_setError(VALUE_ERROR,"UMFPACK: symbolic factorization failed.");
85	return;
86	} else {
87	/* call LDU factorization: */
88	error= umfpack_di_numeric(A->pattern->ptr,A->pattern->index,A->val,pt->symbolic,&(pt->numeric),control,info);
89	if (error != UMFPACK_OK) {
90	if (verbose) printf("UMFPACK: LDU factorization failed.\n");
91	Esys_setError(ZERO_DIVISION_ERROR,"factorization failed. Most likely the matrix is singular.");
92	return;
93	}
94	if (verbose) printf("UMFPACK: LDU factorization completed (time = %e).\n",Esys_timer()-time0);
95	}
96	}
97	if (Esys_noError()) {
98	/* call forward backward substitution: */
99	control[UMFPACK_IRSTEP]=numRefinements; /* number of refinement steps */
100	time0=Esys_timer();
101	error=umfpack_di_solve(UMFPACK_A,A->pattern->ptr,A->pattern->index,A->val,out,in,pt->numeric,control,info);
102	if (error != UMFPACK_OK) {
103	if (verbose) printf("UMFPACK: forward/backward substitution failed.\n");
104	Esys_setError(VALUE_ERROR,"forward/backward substition failed. Most likely the matrix is singular.");
105	return;
106	}
107	if (verbose) printf("UMFPACK: forward/backward substitution completed (time = %e).\n",Esys_timer()-time0);
108	}
109	#else
110	Esys_setError(SYSTEM_ERROR,"Paso_UMFPACK:UMFPACK is not avialble.");
111	#endif
112	}