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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1384 - (hide annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 9 months ago) by phornby
Original Path: temp_trunk_copy/paso/src/UMFPACK.c
File MIME type: text/plain
File size: 3873 byte(s)
Make a temp copy of the trunk before checking in the windows changes


1 ksteube 1312
2 gross 806 /* $Id$ */
3    
4 ksteube 1312 /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15 gross 806
16     /**************************************************************/
17    
18     /* Paso: interface to the Intel UMFPACK library */
19    
20     /**************************************************************/
21    
22     /* Copyrights by ACcESS Australia 2006 */
23     /* Author: gross@@access.edu.au */
24    
25     /**************************************************************/
26    
27     #include "Paso.h"
28     #include "UMFPACK.h"
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 ksteube 1312 Paso_UMFPACK_Handler* pt =NULL;
39 gross 806 #ifdef UMFPACK
40     if (A->solver!=NULL) {
41 ksteube 1312 pt=(Paso_UMFPACK_Handler*)(A->solver);
42 gross 806 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 gross 1149 Paso_UMFPACK_Handler* pt = NULL;
61 gross 806
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 gross 1149 pt = (Paso_UMFPACK_Handler *)(A->solver);
68 gross 806 umfpack_di_defaults(control);
69    
70     if (pt==NULL) {
71 ksteube 1312 int n = A->mainBlock->numRows;
72 gross 806 pt=MEMALLOC(1,Paso_UMFPACK_Handler);
73     if (Paso_checkPtr(pt)) return;
74     A->solver=(void*) pt;
75     time0=Paso_timer();
76 gross 1149 /* call LDU symbolic factorization: */
77 ksteube 1312 error=umfpack_di_symbolic(n,n,A->mainBlock->pattern->ptr,A->mainBlock->pattern->index,A->mainBlock->val,&(pt->symbolic),control,info);
78 gross 806 if (error != UMFPACK_OK) {
79     Paso_setError(VALUE_ERROR,"symbolic factorization failed.");
80     Paso_UMFPACK_free(A);
81     } else {
82 gross 1149 /* call LDU factorization: */
83 ksteube 1312 error= umfpack_di_numeric(A->mainBlock->pattern->ptr,A->mainBlock->pattern->index,A->mainBlock->val,pt->symbolic,&(pt->numeric),control,info);
84 gross 806 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 gross 1149 /* call forward backward substitution: */
94 gross 806 control[UMFPACK_IRSTEP]=2; /* number of refinement steps */
95 ksteube 1312 error=umfpack_di_solve(UMFPACK_A,A->mainBlock->pattern->ptr,A->mainBlock->pattern->index,A->mainBlock->val,out,in,pt->numeric,control,info);
96 gross 806 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     }

  ViewVC Help
Powered by ViewVC 1.1.26