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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2748 - (show annotations)
Tue Nov 17 07:32:59 2009 UTC (9 years, 8 months ago) by gross
File MIME type: text/plain
File size: 4227 byte(s)
Macro elements are implemented now. VTK writer for macro elements still needs testing.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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_SystemMatrix* A) {
38 Paso_UMFPACK_Handler* pt =NULL;
39 if (A->solver!=NULL) {
40 pt=(Paso_UMFPACK_Handler*)(A->solver);
41 Paso_UMFPACK1_free((Paso_UMFPACK_Handler*)A->solver);
42 A->solver=NULL;
43 }
44 }
45 void Paso_UMFPACK1_free(Paso_UMFPACK_Handler* pt) {
46 if (pt!=NULL) {
47 #ifdef UMFPACK
48 umfpack_di_free_symbolic(&(pt->symbolic));
49 umfpack_di_free_numeric(&(pt->numeric));
50 #endif
51 MEMFREE(pt);
52 }
53 }
54
55
56 /* call the solver: */
57
58 void Paso_UMFPACK(Paso_SystemMatrix* A,
59 double* out,
60 double* in,
61 Paso_Options* options,
62 Paso_Performance* pp) {
63 #ifdef UMFPACK
64 double time0;
65 Paso_UMFPACK_Handler* pt = NULL;
66 options->converged=FALSE;
67
68 if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
69 Paso_setError(TYPE_ERROR,"Paso_UMFPACK: UMFPACK requires CSR format with index offset 1 and block size 1.");
70 return;
71 }
72 Performance_startMonitor(pp,PERFORMANCE_ALL);
73 pt = (Paso_UMFPACK_Handler *)(A->solver);
74
75 time0=Paso_timer();
76 Paso_UMFPACK1(&pt, A->mainBlock, out, in, 2);
77 options->set_up_time=0;
78 options->time=Paso_timer()-time0;
79 if (!Paso_noError()) {
80 Paso_UMFPACK_free(A);
81 } else {
82 if (options->verbose) printf("UMFPACK: solve completed.\n");
83 A->solver=(void*) pt;
84 options->converged=TRUE;
85 options->residual_norm=0;
86 options->num_iter=1;
87 options->num_level=0;
88 options->num_inner_iter=0;
89 }
90 Performance_stopMonitor(pp,PERFORMANCE_ALL);
91 }
92
93 void Paso_UMFPACK1(Paso_UMFPACK_Handler** pt, Paso_SparseMatrix* A, double* out, double* in, const int refines) {
94 double control[UMFPACK_CONTROL], info[UMFPACK_INFO];
95 int error = UMFPACK_OK;
96 umfpack_di_defaults(control);
97
98 if (*pt==NULL) {
99 int n = A->numRows;
100 *pt=(MEMALLOC(1,Paso_UMFPACK_Handler));
101 if (Paso_checkPtr(*pt)) return;
102 /* call LDU symbolic factorization: */
103 error=umfpack_di_symbolic(n,n,A->pattern->ptr,A->pattern->index,A->val,&((*pt)->symbolic),control,info);
104 if (error != UMFPACK_OK) {
105 Paso_setError(VALUE_ERROR,"symbolic factorization failed.");
106 return;
107 } else {
108 /* call LDU factorization: */
109 error= umfpack_di_numeric(A->pattern->ptr,A->pattern->index,A->val,(*pt)->symbolic,&((*pt)->numeric),control,info);
110 if (error != UMFPACK_OK) {
111 Paso_setError(ZERO_DIVISION_ERROR,"factorization failed. Most likely the matrix is singular.");
112 return;
113 }
114 }
115 }
116 if (Paso_noError()) {
117 /* call forward backward substitution: */
118 control[UMFPACK_IRSTEP]=refines; /* number of refinement steps */
119 error=umfpack_di_solve(UMFPACK_A,A->pattern->ptr,A->pattern->index,A->val,out,in,(*pt)->numeric,control,info);
120 if (error != UMFPACK_OK) {
121 Paso_setError(VALUE_ERROR,"forward/backward substition failed. Most likely the matrix is singular.");
122 return;
123 }
124 }
125 #else
126 Paso_setError(SYSTEM_ERROR,"Paso_UMFPACK:UMFPACK is not avialble.");
127 #endif
128
129 }

  ViewVC Help
Powered by ViewVC 1.1.26