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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3314 - (show annotations)
Tue Oct 26 22:22:20 2010 UTC (8 years, 11 months ago) by gross
File MIME type: text/plain
File size: 3726 byte(s)
output removed/leaned up
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) && (A->type & 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 }

  ViewVC Help
Powered by ViewVC 1.1.26