/[escript]/trunk/esys2/finley/src/finleyC/Solvers/Solver_ILU0.c
ViewVC logotype

Contents of /trunk/esys2/finley/src/finleyC/Solvers/Solver_ILU0.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (14 years, 4 months ago) by jgs
File MIME type: text/plain
File size: 7094 byte(s)
*** empty log message ***

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Finley: ILU(0) preconditioner: */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia 2003/04 */
10 /* Author: gross@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Finley.h"
15 #include "System.h"
16 #include "Solver.h"
17
18 void Finley_Solver_setILU0_1(int,maybelong*,maybelong*,double*,maybelong,maybelong*,maybelong*);
19 void Finley_Solver_solveILU0_1(int,double*,maybelong*,maybelong*,double*,maybelong,maybelong*,maybelong*);
20
21
22 /**************************************************************/
23
24 /* inverse preconditioner setup */
25
26 void Finley_Solver_setILU0(Finley_SystemMatrix * A_p) {
27 #if ITERATIVE_SOLVER == NO_LIB
28 Finley_Solver_Preconditioner* prec=(Finley_Solver_Preconditioner*) (A_p->iterative);
29 maybelong n=A_p->num_rows;
30 /* allocate arrays */
31 prec->mainDiag =MEMALLOC(n,maybelong);
32 prec->color = MEMALLOC(n,maybelong);
33 prec->values = MEMALLOC((A_p->len)*(size_t) n,double);
34 if (! (Finley_checkPtr(prec->mainDiag) || Finley_checkPtr(prec->color) || Finley_checkPtr(prec->values))) {
35 /* find the main diagonal: */
36 Finley_Solver_getMainDiagonal(A_p,prec->mainDiag);
37 /* color the rows : */
38 Finley_Solver_coloring(A_p->pattern,&(prec->numColors),prec->color);
39 if (Finley_ErrorCode==NO_ERROR) {
40 /* allocate vector to hold main diagonal entries: */
41 Finley_SystemMatrix_copy(A_p,prec->values);
42 /* allocate vector to hold main diagonal entries: */
43 if (A_p->row_block_size > 1) {
44 Finley_ErrorCode = TYPE_ERROR;
45 sprintf(Finley_ErrorMsg, "ILU preconditioner requires block size =1 .");
46 return;
47 } else {
48 Finley_Solver_setILU0_1(n,A_p->pattern->ptr,A_p->pattern->index,prec->values,prec->numColors,prec->color,prec->mainDiag);
49 }
50 }
51 }
52 if (Finley_ErrorCode!=NO_ERROR) Finley_Preconditioner_free(A_p->iterative);
53 return;
54 #endif
55 }
56
57 /**************************************************************/
58
59 /* inverse preconditioner solution using raw pointers */
60
61 /* should be called within a parallel region */
62 /* barrier synconization should be performed to make sure that the input vector available */
63
64 void Finley_Solver_solveILU0(Finley_SystemMatrix * A_p, double * x, double * b) {
65 #if ITERATIVE_SOLVER == NO_LIB
66 Finley_Solver_Preconditioner* prec=(Finley_Solver_Preconditioner*) (A_p->iterative);
67 maybelong n=A_p->num_rows;
68 maybelong i;
69 #pragma omp for private(i)
70 for (i=0;i<n*A_p->row_block_size;i++) x[i]=b[i];
71 if (A_p->row_block_size > 1) {
72 /* Finley_Solver_solveILU0_blk(n,x,A_p->pattern->ptr,A_p->pattern->index,prec->values,prec->numColors,prec->color,prec->mainDiag); */
73 } else {
74 Finley_Solver_solveILU0_1(n,x,A_p->pattern->ptr,A_p->pattern->index,prec->values,prec->numColors,prec->color,prec->mainDiag);
75 }
76 return;
77 #endif
78 }
79 /**************************************************************/
80
81 /* applies ILU preconditioner on x for block size 1 (in place) : */
82
83 void Finley_Solver_solveILU0_1(int n,double* x, maybelong* ptr,maybelong* index,double* values,
84 maybelong numColors,maybelong* color, maybelong* mainDiag) {
85 #if ITERATIVE_SOLVER == NO_LIB
86 int iColor,i,k,ik;
87 /* forward substitution */
88 iColor=0;
89 while (iColor<numColors) {
90 #pragma omp for private(i,k,ik)
91 for (i=0;i<n;i++) {
92 if (color[i]==iColor) {
93 for (ik=ptr[i]; ik < ptr[i+1]; ik++) {
94 k=index[ik];
95 if (color[k]<color[i]) x[i]-=values[ik]*x[k];
96 } /* end of ik-loop */
97 }
98 } /* end of i-loop */
99 iColor++;
100 }
101 /* backwards substitution */
102 iColor=numColors;
103 while (0<iColor) {
104 iColor--;
105 #pragma omp for private(i,k,ik)
106 for (i=0;i<n;i++) {
107 if (color[i]==iColor) {
108 for (ik=ptr[i]; ik < ptr[i+1]; ik++) {
109 k=index[ik];
110 if (color[k]>color[i]) x[i]-=values[ik]*x[k];
111 } /* end of ik-loop */
112 x[i]*=values[mainDiag[i]];
113 }
114 }
115 }
116 #endif
117 }
118 /**************************************************************/
119
120 /* ILU preconditioner for block size 1 : */
121
122 void Finley_Solver_setILU0_1(int n,maybelong* ptr,maybelong* index,double* values,
123 maybelong numColors,maybelong* color, maybelong* mainDiag) {
124 #if ITERATIVE_SOLVER == NO_LIB
125 double value_ik;
126 int iColor,i,k,j,ij,kj,ik,kColor;
127 #pragma omp parallel private(iColor)
128 {
129 for (iColor=0;iColor<numColors;iColor++) {
130 #pragma omp for private(i,k,j,ij,kj,ik,kColor,value_ik)
131 for (i=0;i<n;i++) {
132 if (color[i]==iColor) {
133 /* printf("processing row = %d\n",i); */
134 for (kColor=0;kColor<color[i];kColor++) {
135 /* run through the row elements a_i,k considering rows have been eleminated only */
136 for (ik=ptr[i]; ik < ptr[i+1]; ik++) {
137 k=index[ik];
138 if (color[k]==kColor) {
139 /* up-date the row entry (i,k) */
140 value_ik=values[ik]*values[mainDiag[k]];
141 values[ik]=value_ik;
142 /* printf(" working on element (%d,%d):%e\n",i,k,values[ik]); */
143 /* run through the entries a_k,j in row k considering rows have not been eleminated only */
144 for (kj=ptr[k];kj<ptr[k+1]; kj++) {
145 j=index[kj];
146 if (color[j]>kColor) {
147 /* is (i,j) in the shape of the matrix? */
148 /* this means that index[ij]=j */
149 for (ij=ptr[i];ij<ptr[i+1];ij++) {
150 if (index[ij]==j) {
151 values[ij]-=value_ik*values[kj];
152 /* printf(" updating (%d,%d) by (%d,%d): %e\n",i,j,k,j,values[ij]); */
153 break;
154 }
155 }
156 }
157 } /* end of kj-loop */
158 }
159 } /* end of ik-loop */
160 }
161 /* invert the main diagonal element */
162 if (ABS(values[mainDiag[i]])>0.) {
163 values[mainDiag[i]]=1./values[mainDiag[i]];
164 } else {
165 Finley_ErrorCode = ZERO_DIVISION_ERROR;
166 sprintf(Finley_ErrorMsg, "ILU factorization failed in row %d.",i);
167 }
168 /* printf("finalizing row %d : entry %e\n",i,values[mainDiag[i]]); */
169 } /* end of color if */
170 } /* end of row loop */
171 } /* end of color loop */
172 }
173 #endif
174 }
175

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26