/[escript]/branches/doubleplusgood/paso/src/UMFPACK.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/UMFPACK.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4324 - (show annotations)
Wed Mar 20 00:55:44 2013 UTC (5 years, 11 months ago) by jfenwick
File size: 6671 byte(s)
continues
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /************************************************************************************/
18
19 /* Paso: interface to the Intel UMFPACK library */
20
21 /************************************************************************************/
22
23 /* Copyrights by ACcESS Australia 2006 */
24 /* Author: l.gross@uq.edu.au */
25
26 /************************************************************************************/
27
28 #include "Paso.h"
29 #include "UMFPACK.h"
30
31 #ifdef _OPENMP
32 #include <omp.h>
33 #endif
34
35 /************************************************************************************/
36
37 /* free any extra stuff possibly used by the UMFPACK library */
38
39 void Paso_UMFPACK_free(Paso_SparseMatrix* A) {
40 Paso_UMFPACK_Handler* pt =NULL;
41 if ( (A->solver_p!=NULL) && (A->solver_package == PASO_UMFPACK) ) {
42 pt=(Paso_UMFPACK_Handler*)(A->solver_p);
43 #ifdef UMFPACK
44 umfpack_di_free_symbolic(&(pt->symbolic));
45 umfpack_di_free_numeric(&(pt->numeric));
46 #endif
47 delete pt;
48 A->solver_p=NULL;
49 }
50 }
51
52
53
54 /* call the solver: */
55
56 void Paso_UMFPACK(Paso_SparseMatrix* A,
57 double* out,
58 double* in,
59 dim_t numRefinements,
60 bool_t verbose)
61 {
62 #ifdef UMFPACK
63 double time0;
64 Paso_UMFPACK_Handler* pt = NULL;
65 double control[UMFPACK_CONTROL], info[UMFPACK_INFO];
66 int error = UMFPACK_OK;
67 if (! ( (A->type & MATRIX_FORMAT_BLK1) && (A->type & MATRIX_FORMAT_CSC)) ) {
68 Esys_setError(TYPE_ERROR,"Paso_UMFPACK: UMFPACK requires CSC format with index offset 1 and block size 1.");
69 return;
70 }
71
72 pt = (Paso_UMFPACK_Handler *)(A->solver_p);
73 umfpack_di_defaults(control);
74
75 if (pt==NULL) {
76 int n = A->numRows;
77 pt=new Paso_UMFPACK_Handler;
78 time0=Esys_timer();
79 if (Esys_checkPtr(pt)) return;
80 A->solver_p = (void*) pt;
81 A->solver_package=PASO_UMFPACK;
82 /* call LDU symbolic factorization: */
83 error=umfpack_di_symbolic(n,n,A->pattern->ptr,A->pattern->index,A->val,&(pt->symbolic),control,info);
84 if (error == UMFPACK_ERROR_out_of_memory ) {
85 if (verbose) printf("UMFPACK: symbolic factorization failed because of memory overlow.\n");
86 Esys_setError(MEMORY_ERROR,"UMFPACK: symbolic factorization failed because of memory overlow.");
87 return;
88 } else if (error == UMFPACK_WARNING_singular_matrix ) {
89 if (verbose) printf("UMFPACK: symbolic factorization failed because of singular matrix.\n");
90 Esys_setError(ZERO_DIVISION_ERROR,"UMFPACK: symbolic factorization failed because of singular matrix.");
91 return;
92 } else if ( (error == UMFPACK_WARNING_determinant_underflow ) || ( error == UMFPACK_WARNING_determinant_overflow ) ) {
93 if (verbose) printf("UMFPACK: symbolic factorization failed because of under/overflow.\n");
94 Esys_setError(FLOATING_POINT_ERROR,"UMFPACK: symbolic factorization failed because of under/overflow.");
95 return;
96 } else if (error != UMFPACK_OK) {
97 if (verbose) printf("UMFPACK: symbolic factorization failed. Error code = %d.\n",error);
98 Esys_setError(SYSTEM_ERROR,"UMFPACK: symbolic factorization failed.");
99 return;
100 } else {
101 /* call LDU factorization: */
102 error= umfpack_di_numeric(A->pattern->ptr,A->pattern->index,A->val,pt->symbolic,&(pt->numeric),control,info);
103 if (error == UMFPACK_ERROR_out_of_memory ) {
104 if (verbose) printf("UMFPACK: LDU factorization failed because of memory overlow.\n");
105 Esys_setError(MEMORY_ERROR,"UMFPACK: LDU factorization failed because of memory overlow.");
106 return;
107 } else if (error == UMFPACK_WARNING_singular_matrix ) {
108 if (verbose) printf("UMFPACK: LDU factorization failed because of singular matrix.\n");
109 Esys_setError(ZERO_DIVISION_ERROR,"UMFPACK: LDU factorization failed because of singular matrix.");
110 return;
111 } else if ( (error == UMFPACK_WARNING_determinant_underflow ) || ( error == UMFPACK_WARNING_determinant_overflow ) ) {
112 if (verbose) printf("UMFPACK: symbolic factorization failed because of under/overflow.\n");
113 Esys_setError(FLOATING_POINT_ERROR,"UMFPACK: symbolic factorization failed because of under/overflow.");
114 return;
115 } else if (error != UMFPACK_OK) {
116 if (verbose) printf("UMFPACK: LDU factorization failed. Error code = %d.\n",error);
117 Esys_setError(SYSTEM_ERROR,"UMFPACK: factorization failed.");
118 return;
119 }
120 if (verbose) printf("UMFPACK: LDU factorization completed (time = %e).\n",Esys_timer()-time0);
121 }
122 }
123 if (Esys_noError()) {
124 /* call forward backward substitution: */
125 control[UMFPACK_IRSTEP]=numRefinements; /* number of refinement steps */
126 time0=Esys_timer();
127 error=umfpack_di_solve(UMFPACK_A,A->pattern->ptr,A->pattern->index,A->val,out,in,pt->numeric,control,info);
128 if (error == UMFPACK_ERROR_out_of_memory ) {
129 if (verbose) printf("UMFPACK: forward/backward substitution failed because of memory overlow.\n");
130 Esys_setError(MEMORY_ERROR,"UMFPACK: forward/backward substitution failed because of memory overlow.");
131 return;
132 } else if (error == UMFPACK_WARNING_singular_matrix ) {
133 if (verbose) printf("UMFPACK: forward/backward substitution because of singular matrix.\n");
134 Esys_setError(ZERO_DIVISION_ERROR,"UMFPACK: forward/backward substitution failed because of singular matrix.");
135 return;
136 } else if ( (error == UMFPACK_WARNING_determinant_underflow ) || ( error == UMFPACK_WARNING_determinant_overflow ) ) {
137 if (verbose) printf("UMFPACK: forward/backward substitution failed because of under/overflow.\n");
138 Esys_setError(FLOATING_POINT_ERROR,"UMFPACK: forward/backward substitution failed because of under/overflow.");
139 return;
140 } else if (error != UMFPACK_OK) {
141 if (verbose) printf("UMFPACK: forward/backward substitution failed. Error code = %d.\n", error);
142 Esys_setError(SYSTEM_ERROR,"UMFPACK: forward/backward substitution failed.");
143 return;
144 }
145 if (verbose) printf("UMFPACK: forward/backward substitution completed (time = %e).\n",Esys_timer()-time0);
146 }
147 #else
148 Esys_setError(SYSTEM_ERROR,"Paso_UMFPACK: UMFPACK is not available.");
149 #endif
150 }
151
152
153
154

  ViewVC Help
Powered by ViewVC 1.1.26