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

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

Parent Directory Parent Directory | Revision Log Revision Log


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


1 gross 806 /* $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