paso/src/UMFPACK.c

/* $Id$ */

/*
********************************************************************************
*               Copyright   2006 by ACcESS MNRF                                *
*                                                                              * 
*                 http://www.access.edu.au                                     *
*           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: gross@@access.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_SystemMatrix* A) {
#ifdef UMFPACK
      if (A->solver!=NULL) {
           Paso_UMFPACK_Handler* pt=(Paso_UMFPACK_Handler*)(A->solver);
           umfpack_di_free_symbolic(&(pt->symbolic));
           umfpack_di_free_numeric(&(pt->numeric));
           MEMFREE(pt);
           A->solver=NULL;
     }
#endif
}
/*  call the solver: */

void Paso_UMFPACK(Paso_SystemMatrix* A,
                          double* out,
                          double* in,
                          Paso_Options* options,
                          Paso_Performance* pp) {
#ifdef UMFPACK
     double time0;
     double control[UMFPACK_CONTROL], info[UMFPACK_INFO];
     int error = UMFPACK_OK;
     Paso_UMFPACK_Handler* pt = NULL;

     if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
        Paso_setError(TYPE_ERROR,"Paso_UMFPACK: UMFPACK requires CSR format with index offset 1 and block size 1.");
        return;
     }
     Performance_startMonitor(pp,PERFORMANCE_ALL);
     pt = (Paso_UMFPACK_Handler *)(A->solver);
     umfpack_di_defaults(control);

     if (pt==NULL) {
        int n = A->num_rows;
        pt=MEMALLOC(1,Paso_UMFPACK_Handler);
        if (Paso_checkPtr(pt)) return;
        A->solver=(void*) pt;
        time0=Paso_timer();
        /* 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) {
             Paso_setError(VALUE_ERROR,"symbolic factorization failed.");
             Paso_UMFPACK_free(A);
        } 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) {
             Paso_setError(ZERO_DIVISION_ERROR,"factorization failed. Most likely the matrix is singular.");
             Paso_UMFPACK_free(A);
           }
           if (options->verbose) printf("timing UMFPACK: LDU factorization: %.4e sec.\n",Paso_timer()-time0);
        }
     }
     if (Paso_noError())  {
        time0=Paso_timer();
        /* call forward backward substitution: */
        control[UMFPACK_IRSTEP]=2; /* number of refinement steps */
        error=umfpack_di_solve(UMFPACK_A,A->pattern->ptr,A->pattern->index,A->val,out,in,pt->numeric,control,info);
        if (options->verbose) printf("timing UMFPACK: solve: %.4e sec\n",Paso_timer()-time0);
        if (error != UMFPACK_OK) {
              Paso_setError(VALUE_ERROR,"forward/backward substition failed. Most likely the matrix is singular.");
        }
     }
     Performance_stopMonitor(pp,PERFORMANCE_ALL);
#else
    Paso_setError(SYSTEM_ERROR,"Paso_UMFPACK:UMFPACK is not avialble.");
#endif

}
/*
 * $Log$
 *
 */
1	/* $Id$ */
2
3	/*
4	********************************************************************************
5	* Copyright 2006 by ACcESS MNRF *
6	* *
7	* http://www.access.edu.au *
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	/* Paso: interface to the Intel UMFPACK library */
17
18	/**************************************************************/
19
20	/* Copyrights by ACcESS Australia 2006 */
21	/* Author: gross@@access.edu.au */
22
23	/**************************************************************/
24
25	#include "Paso.h"
26	#include "UMFPACK.h"
27	#ifdef _OPENMP
28	#include <omp.h>
29	#endif
30
31	/**************************************************************/
32
33	/* free any extra stuff possibly used by the UMFPACK library */
34
35	void Paso_UMFPACK_free(Paso_SystemMatrix* A) {
36	#ifdef UMFPACK
37	if (A->solver!=NULL) {
38	Paso_UMFPACK_Handler* pt=(Paso_UMFPACK_Handler*)(A->solver);
39	umfpack_di_free_symbolic(&(pt->symbolic));
40	umfpack_di_free_numeric(&(pt->numeric));
41	MEMFREE(pt);
42	A->solver=NULL;
43	}
44	#endif
45	}
46	/* call the solver: */
47
48	void Paso_UMFPACK(Paso_SystemMatrix* A,
49	double* out,
50	double* in,
51	Paso_Options* options,
52	Paso_Performance* pp) {
53	#ifdef UMFPACK
54	double time0;
55	double control[UMFPACK_CONTROL], info[UMFPACK_INFO];
56	int error = UMFPACK_OK;
57	Paso_UMFPACK_Handler* pt = NULL;
58
59	if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
60	Paso_setError(TYPE_ERROR,"Paso_UMFPACK: UMFPACK requires CSR format with index offset 1 and block size 1.");
61	return;
62	}
63	Performance_startMonitor(pp,PERFORMANCE_ALL);
64	pt = (Paso_UMFPACK_Handler *)(A->solver);
65	umfpack_di_defaults(control);
66
67	if (pt==NULL) {
68	int n = A->num_rows;
69	pt=MEMALLOC(1,Paso_UMFPACK_Handler);
70	if (Paso_checkPtr(pt)) return;
71	A->solver=(void*) pt;
72	time0=Paso_timer();
73	/* call LDU symbolic factorization: */
74	error=umfpack_di_symbolic(n,n,A->pattern->ptr,A->pattern->index,A->val,&(pt->symbolic),control,info);
75	if (error != UMFPACK_OK) {
76	Paso_setError(VALUE_ERROR,"symbolic factorization failed.");
77	Paso_UMFPACK_free(A);
78	} else {
79	/* call LDU factorization: */
80	error= umfpack_di_numeric(A->pattern->ptr,A->pattern->index,A->val,pt->symbolic,&(pt->numeric),control,info);
81	if (error != UMFPACK_OK) {
82	Paso_setError(ZERO_DIVISION_ERROR,"factorization failed. Most likely the matrix is singular.");
83	Paso_UMFPACK_free(A);
84	}
85	if (options->verbose) printf("timing UMFPACK: LDU factorization: %.4e sec.\n",Paso_timer()-time0);
86	}
87	}
88	if (Paso_noError()) {
89	time0=Paso_timer();
90	/* call forward backward substitution: */
91	control[UMFPACK_IRSTEP]=2; /* number of refinement steps */
92	error=umfpack_di_solve(UMFPACK_A,A->pattern->ptr,A->pattern->index,A->val,out,in,pt->numeric,control,info);
93	if (options->verbose) printf("timing UMFPACK: solve: %.4e sec\n",Paso_timer()-time0);
94	if (error != UMFPACK_OK) {
95	Paso_setError(VALUE_ERROR,"forward/backward substition failed. Most likely the matrix is singular.");
96	}
97	}
98	Performance_stopMonitor(pp,PERFORMANCE_ALL);
99	#else
100	Paso_setError(SYSTEM_ERROR,"Paso_UMFPACK:UMFPACK is not avialble.");
101	#endif
102
103	}
104	/*
105	* $Log$
106	*
107	*/