/[escript]/trunk/paso/src/UMFPACK.c
ViewVC logotype

Contents of /trunk/paso/src/UMFPACK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1981 - (show annotations)
Thu Nov 6 05:27:33 2008 UTC (10 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 5986 byte(s)
More warning removal.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 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: gross@@access.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_SystemMatrix* A) {
38 #ifdef UMFPACK
39 Paso_UMFPACK_Handler* pt =NULL;
40 if (A->solver!=NULL) {
41 pt=(Paso_UMFPACK_Handler*)(A->solver);
42 umfpack_di_free_symbolic(&(pt->symbolic));
43 umfpack_di_free_numeric(&(pt->numeric));
44 MEMFREE(pt);
45 A->solver=NULL;
46 }
47 #endif
48 }
49 /* call the solver: */
50
51 void Paso_UMFPACK(Paso_SystemMatrix* A,
52 double* out,
53 double* in,
54 Paso_Options* options,
55 Paso_Performance* pp) {
56 #ifdef UMFPACK
57 double time0;
58 double control[UMFPACK_CONTROL], info[UMFPACK_INFO];
59 int error = UMFPACK_OK;
60 Paso_UMFPACK_Handler* pt = NULL;
61
62 if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
63 Paso_setError(TYPE_ERROR,"Paso_UMFPACK: UMFPACK requires CSR format with index offset 1 and block size 1.");
64 return;
65 }
66 Performance_startMonitor(pp,PERFORMANCE_ALL);
67 pt = (Paso_UMFPACK_Handler *)(A->solver);
68 umfpack_di_defaults(control);
69
70 if (pt==NULL) {
71 int n = A->mainBlock->numRows;
72 pt=MEMALLOC(1,Paso_UMFPACK_Handler);
73 if (Paso_checkPtr(pt)) return;
74 A->solver=(void*) pt;
75 time0=Paso_timer();
76 /* call LDU symbolic factorization: */
77 error=umfpack_di_symbolic(n,n,A->mainBlock->pattern->ptr,A->mainBlock->pattern->index,A->mainBlock->val,&(pt->symbolic),control,info);
78 if (error != UMFPACK_OK) {
79 Paso_setError(VALUE_ERROR,"symbolic factorization failed.");
80 Paso_UMFPACK_free(A);
81 } else {
82 /* call LDU factorization: */
83 error= umfpack_di_numeric(A->mainBlock->pattern->ptr,A->mainBlock->pattern->index,A->mainBlock->val,pt->symbolic,&(pt->numeric),control,info);
84 if (error != UMFPACK_OK) {
85 Paso_setError(ZERO_DIVISION_ERROR,"factorization failed. Most likely the matrix is singular.");
86 Paso_UMFPACK_free(A);
87 }
88 if (options->verbose) printf("timing UMFPACK: LDU factorization: %.4e sec.\n",Paso_timer()-time0);
89 }
90 }
91 if (Paso_noError()) {
92 time0=Paso_timer();
93 /* call forward backward substitution: */
94 control[UMFPACK_IRSTEP]=2; /* number of refinement steps */
95 error=umfpack_di_solve(UMFPACK_A,A->mainBlock->pattern->ptr,A->mainBlock->pattern->index,A->mainBlock->val,out,in,pt->numeric,control,info);
96 if (options->verbose) printf("timing UMFPACK: solve: %.4e sec\n",Paso_timer()-time0);
97 if (error != UMFPACK_OK) {
98 Paso_setError(VALUE_ERROR,"forward/backward substition failed. Most likely the matrix is singular.");
99 }
100 }
101 Performance_stopMonitor(pp,PERFORMANCE_ALL);
102 #else
103 Paso_setError(SYSTEM_ERROR,"Paso_UMFPACK:UMFPACK is not avialble.");
104 #endif
105
106 }
107
108
109 void Paso_UMFPACK1(Paso_SparseMatrix* A,
110 double* out,
111 double* in,
112 bool_t verbose) {
113 #ifdef UMFPACK
114 double time0;
115 double control[UMFPACK_CONTROL], info[UMFPACK_INFO];
116 int error = UMFPACK_OK;
117 Paso_UMFPACK_Handler* pt = NULL;
118
119 if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
120 Paso_setError(TYPE_ERROR,"Paso_UMFPACK: UMFPACK requires CSR format with index offset 1 and block size 1.");
121 return;
122 }
123 umfpack_di_defaults(control);
124
125 if (pt==NULL) {
126 int n = A->numRows;
127 pt=MEMALLOC(1,Paso_UMFPACK_Handler);
128 if (Paso_checkPtr(pt)) return;
129 time0=Paso_timer();
130 /* call LDU symbolic factorization: */
131 error=umfpack_di_symbolic(n,n,A->pattern->ptr,A->pattern->index,A->val,&(pt->symbolic),control,info);
132 if (error != UMFPACK_OK) {
133 Paso_setError(VALUE_ERROR,"symbolic factorization failed.");
134 MEMFREE(pt);
135 } else {
136 /* call LDU factorization: */
137 error= umfpack_di_numeric(A->pattern->ptr,A->pattern->index,A->val,pt->symbolic,&(pt->numeric),control,info);
138 if (error != UMFPACK_OK) {
139 Paso_setError(ZERO_DIVISION_ERROR,"factorization failed. Most likely the matrix is singular.");
140 MEMFREE(pt);
141 }
142 if (verbose) printf("timing UMFPACK: LDU factorization: %.4e sec.\n",Paso_timer()-time0);
143 }
144 }
145 if (Paso_noError()) {
146 time0=Paso_timer();
147 /* call forward backward substitution: */
148 control[UMFPACK_IRSTEP]=2; /* number of refinement steps */
149 error=umfpack_di_solve(UMFPACK_A,A->pattern->ptr,A->pattern->index,A->val,out,in,pt->numeric,control,info);
150 if (verbose) printf("timing UMFPACK: solve: %.4e sec\n",Paso_timer()-time0);
151 if (error != UMFPACK_OK) {
152 Paso_setError(VALUE_ERROR,"forward/backward substition failed. Most likely the matrix is singular.");
153 }
154 }
155 #else
156 Paso_setError(SYSTEM_ERROR,"Paso_UMFPACK:UMFPACK is not avialble.");
157 #endif
158
159 }

  ViewVC Help
Powered by ViewVC 1.1.26