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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3010 - (show annotations)
Tue Apr 27 05:10:46 2010 UTC (12 years, 3 months ago) by artak
File MIME type: text/plain
File size: 5786 byte(s)
Preparation for AMG on Systems without cut using frobenius norm
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 #define MAX_BLOCK_SIZE 3
34
35 /**************************************************************/
36
37 /* free any extra stuff possibly used by the UMFPACK library */
38
39 void Paso_UMFPACK_free(Paso_SystemMatrix* A) {
40 Paso_UMFPACK_Handler* pt =NULL;
41 if (A->solver!=NULL) {
42 pt=(Paso_UMFPACK_Handler*)(A->solver);
43 Paso_UMFPACK1_free((Paso_UMFPACK_Handler*)A->solver);
44 A->solver=NULL;
45 }
46 }
47 void Paso_UMFPACK1_free(Paso_UMFPACK_Handler* pt) {
48 if (pt!=NULL) {
49 #ifdef UMFPACK
50 umfpack_di_free_symbolic(&(pt->symbolic));
51 umfpack_di_free_numeric(&(pt->numeric));
52 #endif
53 MEMFREE(pt);
54 }
55 }
56
57
58 /* call the solver: */
59
60 void Paso_UMFPACK(Paso_SystemMatrix* A,
61 double* out,
62 double* in,
63 Paso_Options* options,
64 Paso_Performance* pp) {
65 #ifdef UMFPACK
66 double time0;
67 Paso_UMFPACK_Handler* pt = NULL;
68 options->converged=FALSE;
69
70 if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
71 Paso_setError(TYPE_ERROR,"Paso_UMFPACK: UMFPACK requires CSR format with index offset 1 and block size 1.");
72 return;
73 }
74 Performance_startMonitor(pp,PERFORMANCE_ALL);
75 pt = (Paso_UMFPACK_Handler *)(A->solver);
76
77 time0=Paso_timer();
78 Paso_UMFPACK1(&pt, A->mainBlock, out, in, 2);
79 options->set_up_time=0;
80 options->time=Paso_timer()-time0;
81 if (!Paso_noError()) {
82 Paso_UMFPACK_free(A);
83 } else {
84 if (options->verbose) printf("UMFPACK: solve completed.\n");
85 A->solver=(void*) pt;
86 options->converged=TRUE;
87 options->residual_norm=0;
88 options->num_iter=1;
89 options->num_level=0;
90 options->num_inner_iter=0;
91 }
92 Performance_stopMonitor(pp,PERFORMANCE_ALL);
93 }
94
95 void Paso_UMFPACK1(Paso_UMFPACK_Handler** pt, Paso_SparseMatrix* A, double* out, double* in, const int refines) {
96 double control[UMFPACK_CONTROL], info[UMFPACK_INFO];
97 int error = UMFPACK_OK;
98 dim_t i,j;
99 /*****/
100 double **xx;
101 double **bb;
102 Paso_SparseMatrix *block[MAX_BLOCK_SIZE];
103
104 xx=MEMALLOC(A->row_block_size,double*);
105 bb=MEMALLOC(A->row_block_size,double*);
106 if (Paso_checkPtr(xx) || Paso_checkPtr(bb)) return;
107 /****/
108
109 umfpack_di_defaults(control);
110
111 for (i=0;i<A->row_block_size;i++) {
112 xx[i]=MEMALLOC(A->numRows,double);
113 bb[i]=MEMALLOC(A->numRows,double);
114 if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
115 }
116
117 /*#pragma omp parallel for private(i,j) schedule(static)*/
118 for (i=0;i<A->numRows;i++) {
119 for (j=0;j<A->row_block_size;j++) {
120 bb[j][i]=in[A->row_block_size*i+j];
121 xx[j][i]=0;
122 }
123 }
124
125 for (i=0;i<MAX_BLOCK_SIZE;++i) {
126 block[i]=NULL;
127 }
128
129 for (i=0;i<A->row_block_size;++i) {
130 block[i]=Paso_SparseMatrix_getBlock(A,i+1);
131
132 if (*pt==NULL) {
133 int n = block[i]->numRows;
134 *pt=(MEMALLOC(1,Paso_UMFPACK_Handler));
135 if (Paso_checkPtr(*pt)) return;
136 /* call LDU symbolic factorization: */
137 error=umfpack_di_symbolic(n,n,block[i]->pattern->ptr,block[i]->pattern->index,block[i]->val,&((*pt)->symbolic),control,info);
138 if (error != UMFPACK_OK) {
139 Paso_setError(VALUE_ERROR,"symbolic factorization failed.");
140 return;
141 } else {
142 /* call LDU factorization: */
143 error= umfpack_di_numeric(block[i]->pattern->ptr,block[i]->pattern->index,block[i]->val,(*pt)->symbolic,&((*pt)->numeric),control,info);
144 if (error != UMFPACK_OK) {
145 Paso_setError(ZERO_DIVISION_ERROR,"factorization failed. Most likely the matrix is singular.");
146 return;
147 }
148 }
149 }
150
151 if (Paso_noError()) {
152 /* call forward backward substitution: */
153 control[UMFPACK_IRSTEP]=refines; /* number of refinement steps */
154 error=umfpack_di_solve(UMFPACK_A,block[i]->pattern->ptr,block[i]->pattern->index,block[i]->val,xx[i],bb[i],(*pt)->numeric,control,info);
155 if (error != UMFPACK_OK) {
156 Paso_setError(VALUE_ERROR,"forward/backward substition failed. Most likely the matrix is singular.");
157 return;
158 }
159 }
160
161 }
162
163 /*#pragma omp parallel for private(i,j) schedule(static)*/
164 for (i=0;i<A->numRows;i++) {
165 for (j=0;j<A->row_block_size;j++) {
166 out[A->row_block_size*i+j]=xx[j][i];
167 }
168 }
169
170 for (i=0;i<A->row_block_size;i++) {
171 MEMFREE(xx[i]);
172 MEMFREE(bb[i]);
173 Paso_SparseMatrix_free(block[i]);
174 }
175
176 #else
177 Paso_setError(SYSTEM_ERROR,"Paso_UMFPACK:UMFPACK is not avialble.");
178 #endif
179
180 }

  ViewVC Help
Powered by ViewVC 1.1.26