/[escript]/trunk/paso/src/Solvers/Solver_ILU.c
ViewVC logotype

Contents of /trunk/paso/src/Solvers/Solver_ILU.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 431 - (show annotations)
Fri Jan 13 05:07:10 2006 UTC (13 years, 10 months ago) by gross
File MIME type: text/plain
File size: 25279 byte(s)
new implementation of ILU0
1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Paso: ILU preconditioner with reordering */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia 2003,2004,2005 */
10 /* Author: gross@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Paso.h"
15 #include "Solver.h"
16 #include "Util.h"
17
18 /**************************************************************/
19
20 /* free all memory used by ILU */
21
22 void Paso_Solver_ILU_free(Paso_Solver_ILU * in) {
23 if (in!=NULL) {
24 MEMFREE(in->colorOf);
25 MEMFREE(in->factors);
26 MEMFREE(in->main_iptr);
27 Paso_SystemMatrixPattern_dealloc(in->pattern);
28 MEMFREE(in);
29 }
30 }
31
32 /**************************************************************/
33
34 /* constructs the incomplete block factorization of
35
36 */
37 Paso_Solver_ILU* Paso_Solver_getILU(Paso_SystemMatrix * A,bool_t verbose) {
38 dim_t n=A->num_rows;
39 dim_t n_block=A->row_block_size;
40 index_t num_colors=0;
41 register double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;
42 register double mainA11,mainA12,mainA13,mainA21,mainA22,mainA23,mainA31,mainA32,mainA33;
43 register double S11,S12,S13,S21,S22,S23,S31,S32,S33;
44 register index_t i,iptr_main,iptr,iptr_ik,k,iptr_kj,j,iptr_ij,color,color2;
45 double time0,time_color,time_fac;
46 /* allocations: */
47 Paso_Solver_ILU* out=MEMALLOC(1,Paso_Solver_ILU);
48 if (Paso_checkPtr(out)) return NULL;
49 index_t* mis_marker=TMPMEMALLOC(n,index_t);
50 out->colorOf=MEMALLOC(n,index_t);
51 out->factors=MEMALLOC(A->len,double);
52 out->main_iptr=MEMALLOC(n,double);
53 out->pattern=Paso_SystemMatrixPattern_reference(A->pattern);
54 out->n_block=n_block;
55 out->n=n;
56 time0=Paso_timer();
57 if ( !(Paso_checkPtr(mis_marker) ||
58 Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->factors)) ) {
59 /* get coloring */
60 index_t num_colors=0;
61 #pragma omp parallel for private(i) schedule(static)
62 for (i = 0; i < n; ++i) out->colorOf[i]=-1;
63 while (Paso_Util_isAny(n,out->colorOf,-1) && Paso_noError()) {
64 #pragma omp parallel for private(i) schedule(static)
65 for (i = 0; i < n; ++i) mis_marker[i]=out->colorOf[i];
66 Paso_SystemMatrixPattern_mis(A->pattern,mis_marker);
67
68 #pragma omp parallel for private(i) schedule(static)
69 for (i = 0; i < n; ++i) if (mis_marker[i]) out->colorOf[i]=num_colors;
70 ++num_colors;
71 }
72 out->num_colors=num_colors;
73 time_color=Paso_timer()-time0;
74 time0=Paso_timer();
75 /* find main diagonal and copy matrix values */
76 #pragma omp for private(i) schedule(static) local(i,iptr,iptr_main,k)
77 for (i = 0; i < n; ++i) {
78 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
79 iptr_main=A->pattern->ptr[0]-1;
80 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; iptr++) {
81 if (A->pattern->index[iptr]==i) iptr_main=iptr;
82 for (k=0;k<n_block*n_block;++k) out->factors[n_block*n_block*iptr+k]=A->val[n_block*n_block*iptr+k];
83 }
84 out->main_iptr[i]=iptr_main;
85 if (iptr_main==A->pattern->ptr[0]-1)
86 Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: no main diagonal");
87 }
88 }
89 /* start factorization */
90
91 for (color=0;color<out->num_colors && Paso_noError();++color) {
92 if (n_block==1) {
93 #pragma omp for private(i) schedule(static) private(i,color2,iptr_ik,k,iptr_kj,S11,j,iptr_ij,A11,iptr_main,D)
94 for (i = 0; i < n; ++i) {
95 if (out->colorOf[i]==color) {
96 for (color2=0;color2<color;++color2) {
97 for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
98 k=A->pattern->index[iptr_ik];
99 if (out->colorOf[k]==color2) {
100 A11=out->factors[iptr_ik];
101 /* a_ij=a_ij-a_ik*a_kj */
102 for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
103 j=A->pattern->index[iptr_kj];
104 if (out->colorOf[j]>color2) {
105 S11=out->factors[iptr_kj];
106 for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
107 if (j==A->pattern->index[iptr_ij]) {
108 out->factors[iptr_ij]-=A11*S11;
109 break;
110 }
111 }
112 }
113 }
114 }
115 }
116 }
117 iptr_main=out->main_iptr[i];
118 D=out->factors[iptr_main];
119 if (ABS(D)>0.) {
120 D=1./D;
121 out->factors[iptr_main]=D;
122 /* a_ik=a_ii^{-1}*a_ik */
123 for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
124 k=A->pattern->index[iptr_ik];
125 if (out->colorOf[k]>color) {
126 A11=out->factors[iptr_ik];
127 out->factors[iptr_ik]=A11*D;
128 }
129 }
130 } else {
131 Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
132 }
133 }
134 }
135 } else if (n_block==2) {
136 #pragma omp for private(i) schedule(static) private(i,color2,iptr_ik,k,iptr_kj,S11,S21,S12,S22,j,iptr_ij,A11,A21,A12,A22,iptr_main,D)
137 for (i = 0; i < n; ++i) {
138 if (out->colorOf[i]==color) {
139 for (color2=0;color2<color;++color2) {
140 for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
141 k=A->pattern->index[iptr_ik];
142 if (out->colorOf[k]==color2) {
143 A11=out->factors[iptr_ik*4 ];
144 A21=out->factors[iptr_ik*4+1];
145 A12=out->factors[iptr_ik*4+2];
146 A22=out->factors[iptr_ik*4+3];
147 /* a_ij=a_ij-a_ik*a_kj */
148 for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
149 j=A->pattern->index[iptr_kj];
150 if (out->colorOf[j]>color2) {
151 S11=out->factors[iptr_kj*4];
152 S21=out->factors[iptr_kj*4+1];
153 S12=out->factors[iptr_kj*4+2];
154 S22=out->factors[iptr_kj*4+3];
155 for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
156 if (j==A->pattern->index[iptr_ij]) {
157 out->factors[4*iptr_ij ]-=A11*S11+A12*S21;
158 out->factors[4*iptr_ij+1]-=A21*S11+A22*S21;
159 out->factors[4*iptr_ij+2]-=A11*S12+A12*S22;
160 out->factors[4*iptr_ij+3]-=A21*S12+A22*S22;
161 break;
162 }
163 }
164 }
165 }
166 }
167 }
168 }
169 iptr_main=out->main_iptr[i];
170 A11=out->factors[iptr_main*4];
171 A21=out->factors[iptr_main*4+1];
172 A12=out->factors[iptr_main*4+2];
173 A22=out->factors[iptr_main*4+3];
174 D = A11*A22-A12*A21;
175 if (ABS(D)>0.) {
176 D=1./D;
177 S11= A22*D;
178 S21=-A21*D;
179 S12=-A12*D;
180 S22= A11*D;
181 out->factors[iptr_main*4] = S11;
182 out->factors[iptr_main*4+1]= S21;
183 out->factors[iptr_main*4+2]= S12;
184 out->factors[iptr_main*4+3]= S22;
185 /* a_ik=a_ii^{-1}*a_ik */
186 for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
187 k=A->pattern->index[iptr_ik];
188 if (out->colorOf[k]>color) {
189 A11=out->factors[iptr_ik*4 ];
190 A21=out->factors[iptr_ik*4+1];
191 A12=out->factors[iptr_ik*4+2];
192 A22=out->factors[iptr_ik*4+3];
193 out->factors[4*iptr_ik ]=S11*A11+S12*A21;
194 out->factors[4*iptr_ik+1]=S21*A11+S22*A21;
195 out->factors[4*iptr_ik+2]=S11*A12+S12*A22;
196 out->factors[4*iptr_ik+3]=S21*A12+S22*A22;
197 }
198 }
199 } else {
200 Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
201 }
202 }
203 }
204 } else if (n_block==3) {
205 #pragma omp for private(i) schedule(static) private(i,color2,iptr_ik,k,iptr_kj,S11,S21,S31,S12,S22,S32,S13,S23,S33,j,iptr_ij,A11,A21,A31,A12,A22,A32,A13,A23,A33,iptr_main,D)
206 for (i = 0; i < n; ++i) {
207 if (out->colorOf[i]==color) {
208 for (color2=0;color2<color;++color2) {
209 for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
210 k=A->pattern->index[iptr_ik];
211 if (out->colorOf[k]==color2) {
212 A11=out->factors[iptr_ik*9 ];
213 A21=out->factors[iptr_ik*9+1];
214 A31=out->factors[iptr_ik*9+2];
215 A12=out->factors[iptr_ik*9+3];
216 A22=out->factors[iptr_ik*9+4];
217 A32=out->factors[iptr_ik*9+5];
218 A13=out->factors[iptr_ik*9+6];
219 A23=out->factors[iptr_ik*9+7];
220 A33=out->factors[iptr_ik*9+8];
221 /* a_ij=a_ij-a_ik*a_kj */
222 for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
223 j=A->pattern->index[iptr_kj];
224 if (out->colorOf[j]>color2) {
225 S11=out->factors[iptr_kj*9 ];
226 S21=out->factors[iptr_kj*9+1];
227 S31=out->factors[iptr_kj*9+2];
228 S12=out->factors[iptr_kj*9+3];
229 S22=out->factors[iptr_kj*9+4];
230 S32=out->factors[iptr_kj*9+5];
231 S13=out->factors[iptr_kj*9+6];
232 S23=out->factors[iptr_kj*9+7];
233 S33=out->factors[iptr_kj*9+8];
234 for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
235 if (j==A->pattern->index[iptr_ij]) {
236 out->factors[iptr_ij*9 ]-=A11*S11+A12*S21+A13*S31;
237 out->factors[iptr_ij*9+1]-=A21*S11+A22*S21+A23*S31;
238 out->factors[iptr_ij*9+2]-=A31*S11+A32*S21+A33*S31;
239 out->factors[iptr_ij*9+3]-=A11*S12+A12*S22+A13*S32;
240 out->factors[iptr_ij*9+4]-=A21*S12+A22*S22+A23*S32;
241 out->factors[iptr_ij*9+5]-=A31*S12+A32*S22+A33*S32;
242 out->factors[iptr_ij*9+6]-=A11*S13+A12*S23+A13*S33;
243 out->factors[iptr_ij*9+7]-=A21*S13+A22*S23+A23*S33;
244 out->factors[iptr_ij*9+8]-=A31*S13+A32*S23+A33*S33;
245 break;
246 }
247 }
248 }
249 }
250 }
251 }
252 }
253 iptr_main=out->main_iptr[i];
254 A11=out->factors[iptr_main*9 ];
255 A21=out->factors[iptr_main*9+1];
256 A31=out->factors[iptr_main*9+2];
257 A12=out->factors[iptr_main*9+3];
258 A22=out->factors[iptr_main*9+4];
259 A32=out->factors[iptr_main*9+5];
260 A13=out->factors[iptr_main*9+6];
261 A23=out->factors[iptr_main*9+7];
262 A33=out->factors[iptr_main*9+8];
263 D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
264 if (ABS(D)>0.) {
265 D=1./D;
266 S11=(A22*A33-A23*A32)*D;
267 S21=(A31*A23-A21*A33)*D;
268 S31=(A21*A32-A31*A22)*D;
269 S12=(A13*A32-A12*A33)*D;
270 S22=(A11*A33-A31*A13)*D;
271 S32=(A12*A31-A11*A32)*D;
272 S13=(A12*A23-A13*A22)*D;
273 S23=(A13*A21-A11*A23)*D;
274 S33=(A11*A22-A12*A21)*D;
275
276 out->factors[iptr_main*9 ]=S11;
277 out->factors[iptr_main*9+1]=S21;
278 out->factors[iptr_main*9+2]=S31;
279 out->factors[iptr_main*9+3]=S12;
280 out->factors[iptr_main*9+4]=S22;
281 out->factors[iptr_main*9+5]=S32;
282 out->factors[iptr_main*9+6]=S13;
283 out->factors[iptr_main*9+7]=S23;
284 out->factors[iptr_main*9+8]=S33;
285
286 /* a_ik=a_ii^{-1}*a_ik */
287 for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
288 k=A->pattern->index[iptr_ik];
289 if (out->colorOf[k]>color) {
290 A11=out->factors[iptr_ik*9 ];
291 A21=out->factors[iptr_ik*9+1];
292 A31=out->factors[iptr_ik*9+2];
293 A12=out->factors[iptr_ik*9+3];
294 A22=out->factors[iptr_ik*9+4];
295 A32=out->factors[iptr_ik*9+5];
296 A13=out->factors[iptr_ik*9+6];
297 A23=out->factors[iptr_ik*9+7];
298 A33=out->factors[iptr_ik*9+8];
299 out->factors[iptr_ik*9 ]=S11*A11+S12*A21+S13*A31;
300 out->factors[iptr_ik*9+1]=S21*A11+S22*A21+S23*A31;
301 out->factors[iptr_ik*9+2]=S31*A11+S32*A21+S33*A31;
302 out->factors[iptr_ik*9+3]=S11*A12+S12*A22+S13*A32;
303 out->factors[iptr_ik*9+4]=S21*A12+S22*A22+S23*A32;
304 out->factors[iptr_ik*9+5]=S31*A12+S32*A22+S33*A32;
305 out->factors[iptr_ik*9+6]=S11*A13+S12*A23+S13*A33;
306 out->factors[iptr_ik*9+7]=S21*A13+S22*A23+S23*A33;
307 out->factors[iptr_ik*9+8]=S31*A13+S32*A23+S33*A33;
308 }
309 }
310 } else {
311 Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
312 }
313 }
314 }
315 } else {
316 Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: block size greater than 3 is not supported.");
317 }
318 }
319 time_fac=Paso_timer()-time0;
320 }
321
322 TMPMEMFREE(mis_marker);
323 if (Paso_noError()) {
324 if (verbose) {
325 printf("ILU: %d color used \n",out->num_colors);
326 printf("timing: ILU: coloring/elemination : %e/%e\n",time_color,time_fac);
327 }
328 return out;
329 } else {
330 Paso_Solver_ILU_free(out);
331 return NULL;
332 }
333 }
334
335 /**************************************************************/
336
337 /* apply ILU precondition b-> x
338
339 in fact it solves LUx=b in the form x= U^{-1} L^{-1}b
340
341 should be called within a parallel region
342 barrier synconization should be performed to make sure that the input vector available
343
344 */
345
346 void Paso_Solver_solveILU(Paso_Solver_ILU * ilu, double * x, double * b) {
347 register dim_t i,k;
348 register index_t color,iptr_ik,iptr_main;
349 register double S1,S2,S3,R1,R2,R3;
350 dim_t n_block=ilu->n_block;
351 dim_t n=ilu->n;
352
353
354 /* copy x into b*/
355 #pragma omp for private(i) schedule(static)
356 for (i=0;i<n*n_block;++i) x[i]=b[i];
357 /* forward substitution */
358 for (color=0;color<ilu->num_colors;++color) {
359 if (n_block==1) {
360 #pragma omp for private(i) schedule(static) private(i,iptr_ik,k,S1,R1)
361 for (i = 0; i < n; ++i) {
362 if (ilu->colorOf[i]==color) {
363 /* x_i=x_i-a_ik*x_k */
364 S1=x[i];
365 for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
366 k=ilu->pattern->index[iptr_ik];
367 if (ilu->colorOf[k]<color) {
368 R1=x[k];
369 S1-=ilu->factors[iptr_ik]*R1;
370 }
371 }
372 iptr_main=ilu->main_iptr[i];
373 x[i]=ilu->factors[iptr_main]*S1;
374 }
375 }
376 } else if (n_block==2) {
377 #pragma omp for private(i) schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2)
378 for (i = 0; i < n; ++i) {
379 if (ilu->colorOf[i]==color) {
380 /* x_i=x_i-a_ik*x_k */
381 S1=x[2*i];
382 S2=x[2*i+1];
383 for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
384 k=ilu->pattern->index[iptr_ik];
385 if (ilu->colorOf[k]<color) {
386 R1=x[2*k];
387 R2=x[2*k+1];
388 S1-=ilu->factors[4*iptr_ik ]*R1+ilu->factors[4*iptr_ik+2]*R2;
389 S2-=ilu->factors[4*iptr_ik+1]*R1+ilu->factors[4*iptr_ik+3]*R2;
390 }
391 }
392 iptr_main=ilu->main_iptr[i];
393 x[2*i ]=ilu->factors[4*iptr_main ]*S1+ilu->factors[4*iptr_main+2]*S2;
394 x[2*i+1]=ilu->factors[4*iptr_main+1]*S1+ilu->factors[4*iptr_main+3]*S2;
395 }
396
397 }
398 } else if (n_block==3) {
399 #pragma omp for private(i) schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3)
400 for (i = 0; i < n; ++i) {
401 if (ilu->colorOf[i]==color) {
402 /* x_i=x_i-a_ik*x_k */
403 S1=x[3*i];
404 S2=x[3*i+1];
405 S3=x[3*i+2];
406 for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
407 k=ilu->pattern->index[iptr_ik];
408 if (ilu->colorOf[k]<color) {
409 R1=x[3*k];
410 R2=x[3*k+1];
411 R3=x[3*k+2];
412 S1-=ilu->factors[9*iptr_ik ]*R1+ilu->factors[9*iptr_ik+3]*R2+ilu->factors[9*iptr_ik+6]*R3;
413 S2-=ilu->factors[9*iptr_ik+1]*R1+ilu->factors[9*iptr_ik+4]*R2+ilu->factors[9*iptr_ik+7]*R3;
414 S3-=ilu->factors[9*iptr_ik+2]*R1+ilu->factors[9*iptr_ik+5]*R2+ilu->factors[9*iptr_ik+8]*R3;
415 }
416 }
417 iptr_main=ilu->main_iptr[i];
418 x[3*i ]=ilu->factors[9*iptr_main ]*S1+ilu->factors[9*iptr_main+3]*S2+ilu->factors[9*iptr_main+6]*S3;
419 x[3*i+1]=ilu->factors[9*iptr_main+1]*S1+ilu->factors[9*iptr_main+4]*S2+ilu->factors[9*iptr_main+7]*S3;
420 x[3*i+2]=ilu->factors[9*iptr_main+2]*S1+ilu->factors[9*iptr_main+5]*S2+ilu->factors[9*iptr_main+8]*S3;
421 }
422 }
423 }
424 #pragma omp barrier
425 }
426 /* backward substitution */
427 for (color=(ilu->num_colors)-1;color>-1;--color) {
428 if (n_block==1) {
429 #pragma omp for private(i) schedule(static) private(i,iptr_ik,k,S1)
430 for (i = 0; i < n; ++i) {
431 if (ilu->colorOf[i]==color) {
432 /* x_i=x_i-a_ik*x_k */
433 S1=x[i];
434 for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
435 k=ilu->pattern->index[iptr_ik];
436 if (ilu->colorOf[k]>color) {
437 R1=x[k];
438 S1-=ilu->factors[iptr_ik]*R1;
439 }
440 }
441 x[i]=S1;
442 }
443 }
444 } else if (n_block==2) {
445 #pragma omp for private(i) schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2)
446 for (i = 0; i < n; ++i) {
447 if (ilu->colorOf[i]==color) {
448 /* x_i=x_i-a_ik*x_k */
449 S1=x[2*i];
450 S2=x[2*i+1];
451 for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
452 k=ilu->pattern->index[iptr_ik];
453 if (ilu->colorOf[k]>color) {
454 R1=x[2*k];
455 R2=x[2*k+1];
456 S1-=ilu->factors[4*iptr_ik ]*R1+ilu->factors[4*iptr_ik+2]*R2;
457 S2-=ilu->factors[4*iptr_ik+1]*R1+ilu->factors[4*iptr_ik+3]*R2;
458 }
459 }
460 x[2*i]=S1;
461 x[2*i+1]=S2;
462 }
463 }
464 } else if (n_block==3) {
465 #pragma omp for private(i) schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)
466 for (i = 0; i < n; ++i) {
467 if (ilu->colorOf[i]==color) {
468 /* x_i=x_i-a_ik*x_k */
469 S1=x[3*i ];
470 S2=x[3*i+1];
471 S3=x[3*i+2];
472 for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
473 k=ilu->pattern->index[iptr_ik];
474 if (ilu->colorOf[k]>color) {
475 R1=x[3*k];
476 R2=x[3*k+1];
477 R3=x[3*k+2];
478 S1-=ilu->factors[9*iptr_ik ]*R1+ilu->factors[9*iptr_ik+3]*R2+ilu->factors[9*iptr_ik+6]*R3;
479 S2-=ilu->factors[9*iptr_ik+1]*R1+ilu->factors[9*iptr_ik+4]*R2+ilu->factors[9*iptr_ik+7]*R3;
480 S3-=ilu->factors[9*iptr_ik+2]*R1+ilu->factors[9*iptr_ik+5]*R2+ilu->factors[9*iptr_ik+8]*R3;
481 }
482 }
483 x[3*i]=S1;
484 x[3*i+1]=S2;
485 x[3*i+2]=S3;
486 }
487 }
488 }
489 #pragma omp barrier
490 }
491 return;
492 }
493 /*
494 * $Log$
495 * Revision 1.2 2005/09/15 03:44:40 jgs
496 * Merge of development branch dev-02 back to main trunk on 2005-09-15
497 *
498 * Revision 1.1.2.1 2005/09/05 06:29:50 gross
499 * These files have been extracted from finley to define a stand alone libray for iterative
500 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
501 * has not been tested yet.
502 *
503 *
504 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26