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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 806 - (show annotations)
Thu Aug 10 11:58:52 2006 UTC (13 years ago) by gross
File MIME type: text/plain
File size: 3969 byte(s)
Interface to the direct solver library UMLPACK is no implemented.


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
58 if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
59 Paso_setError(TYPE_ERROR,"Paso_UMFPACK: UMFPACK requires CSR format with index offset 1 and block size 1.");
60 return;
61 }
62 Performance_startMonitor(pp,PERFORMANCE_ALL);
63 Paso_UMFPACK_Handler* pt = (Paso_UMFPACK_Handler *)(A->solver);
64 umfpack_di_defaults(control);
65
66 if (pt==NULL) {
67 int n = A->num_rows;
68 pt=MEMALLOC(1,Paso_UMFPACK_Handler);
69 if (Paso_checkPtr(pt)) return;
70 A->solver=(void*) pt;
71 time0=Paso_timer();
72 // call LDU symbolic factorization: //
73 error=umfpack_di_symbolic(n,n,A->pattern->ptr,A->pattern->index,A->val,&(pt->symbolic),control,info);
74 if (error != UMFPACK_OK) {
75 Paso_setError(VALUE_ERROR,"symbolic factorization failed.");
76 Paso_UMFPACK_free(A);
77 } else {
78 // call LDU factorization: //
79 error= umfpack_di_numeric(A->pattern->ptr,A->pattern->index,A->val,pt->symbolic,&(pt->numeric),control,info);
80 if (error != UMFPACK_OK) {
81 Paso_setError(ZERO_DIVISION_ERROR,"factorization failed. Most likely the matrix is singular.");
82 Paso_UMFPACK_free(A);
83 }
84 if (options->verbose) printf("timing UMFPACK: LDU factorization: %.4e sec.\n",Paso_timer()-time0);
85 }
86 }
87 if (Paso_noError()) {
88 time0=Paso_timer();
89 // call forward backward substitution: //
90 control[UMFPACK_IRSTEP]=2; /* number of refinement steps */
91 error=umfpack_di_solve(UMFPACK_A,A->pattern->ptr,A->pattern->index,A->val,out,in,pt->numeric,control,info);
92 if (options->verbose) printf("timing UMFPACK: solve: %.4e sec\n",Paso_timer()-time0);
93 if (error != UMFPACK_OK) {
94 Paso_setError(VALUE_ERROR,"forward/backward substition failed. Most likely the matrix is singular.");
95 }
96 }
97 Performance_stopMonitor(pp,PERFORMANCE_ALL);
98 #else
99 Paso_setError(SYSTEM_ERROR,"Paso_UMFPACK:UMFPACK is not avialble.");
100 #endif
101
102 }
103 /*
104 * $Log$
105 *
106 */

  ViewVC Help
Powered by ViewVC 1.1.26