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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 6 months ago) by ksteube
Original Path: trunk/paso/src/Solver_ILU.c
File MIME type: text/plain
File size: 25731 byte(s)
Copyright updated in all files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26