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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/paso/src/Solvers/Solver_ILU.c revision 155 by jgs, Wed Nov 9 02:02:19 2005 UTC temp_trunk_copy/paso/src/Solver_ILU.c revision 1384 by phornby, Fri Jan 11 02:29:38 2008 UTC
# Line 1  Line 1 
1    
2  /* $Id$ */  /* $Id$ */
3    
4    /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
18  /* Paso: ILU preconditioner with reordering                 */  /* Paso: ILU preconditioner with reordering                 */
# Line 13  Line 26 
26    
27  #include "Paso.h"  #include "Paso.h"
28  #include "Solver.h"  #include "Solver.h"
29  #include "Util.h"  #include "PasoUtil.h"
30    
31  /**************************************************************/  /**************************************************************/
32    
# Line 21  Line 34 
34    
35  void Paso_Solver_ILU_free(Paso_Solver_ILU * in) {  void Paso_Solver_ILU_free(Paso_Solver_ILU * in) {
36       if (in!=NULL) {       if (in!=NULL) {
37          Paso_Solver_ILU_free(in->ILU_of_Schur);          MEMFREE(in->colorOf);
38          MEMFREE(in->inv_A_FF);          MEMFREE(in->factors);
39          MEMFREE(in->A_FF_pivot);          MEMFREE(in->main_iptr);  
40          Paso_SystemMatrix_dealloc(in->A_FC);          Paso_Pattern_free(in->pattern);  
         Paso_SystemMatrix_dealloc(in->A_CF);  
         MEMFREE(in->rows_in_F);  
         MEMFREE(in->rows_in_C);  
         MEMFREE(in->mask_F);  
         MEMFREE(in->mask_C);  
         MEMFREE(in->x_F);  
         MEMFREE(in->b_F);  
         MEMFREE(in->x_C);  
         MEMFREE(in->b_C);  
41          MEMFREE(in);          MEMFREE(in);
42       }       }
43  }  }
44    
45  /**************************************************************/  /**************************************************************/
46    
47  /*   constructs the block-block factorization of  /*   constructs the incomplete block factorization of
   
         [ A_FF A_FC ]  
    A_p=  
         [ A_CF A_FF ]  
   
 to  
   
   [      I         0  ]  [ A_FF 0 ] [ I    invA_FF*A_FF ]    
   [ A_CF*invA_FF   I  ]  [   0  S ] [ 0          I      ]  
   
   
    where S=A_FF-ACF*invA_FF*A_FC within the shape of S  
   
    then ILU is applied to S again until S becomes empty  
48    
49  */  */
50  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SystemMatrix * A_p,bool_t verbose) {  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SparseMatrix * A,bool_t verbose) {
51    Paso_Solver_ILU* out=NULL;    dim_t n=A->numRows;
52    dim_t n=A_p->num_rows;    dim_t n_block=A->row_block_size;
53    dim_t n_block=A_p->row_block_size;    index_t num_colors=0, *mis_marker=NULL;
54    index_t* mis_marker=NULL;      register double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;
55    index_t* counter=NULL;      register double mainA11,mainA12,mainA13,mainA21,mainA22,mainA23,mainA31,mainA32,mainA33;
56    index_t iPtr,*index, *where_p;    register double S11,S12,S13,S21,S22,S23,S31,S32,S33;
57    dim_t i,k;    register index_t i,iptr_main,iptr,iptr_ik,k,iptr_kj,j,iptr_ij,color,color2;
58    Paso_SystemMatrix * schur=NULL;    double time0,time_color,time_fac;
59    double A11,A12,A13,A21,A22,A23,A31,A32,A33,D,time0,time1,time2;    /* allocations: */  
60      Paso_Solver_ILU* out=MEMALLOC(1,Paso_Solver_ILU);
61      if (Paso_checkPtr(out)) return NULL;
62      out->colorOf=MEMALLOC(n,index_t);
63      out->factors=MEMALLOC(A->len,double);
64      out->main_iptr=MEMALLOC(n,index_t);
65      out->pattern=Paso_Pattern_getReference(A->pattern);
66      out->n_block=n_block;
67      out->n=n;
68      if ( !(Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->factors)) ) {
69        time0=Paso_timer();
70        Paso_Pattern_color(A->pattern,&out->num_colors,out->colorOf);
71        time_color=Paso_timer()-time0;
72    
73        if (Paso_noError()) {
74           time0=Paso_timer();
75           /* find main diagonal and copy matrix values */
76           #pragma omp parallel for schedule(static) private(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           #pragma omp barrier
92    /* identify independend set of rows/columns */         for (color=0;color<out->num_colors && Paso_noError();++color) {
93    mis_marker=TMPMEMALLOC(n,index_t);                if (n_block==1) {
94    counter=TMPMEMALLOC(n,index_t);                   #pragma omp parallel for schedule(static) private(i,color2,iptr_ik,k,iptr_kj,S11,j,iptr_ij,A11,iptr_main,D)
95    out=MEMALLOC(1,Paso_Solver_ILU);                   for (i = 0; i < n; ++i) {
96    out->ILU_of_Schur=NULL;                      if (out->colorOf[i]==color) {
97    out->inv_A_FF=NULL;                         for (color2=0;color2<color;++color2) {
98    out->A_FF_pivot=NULL;                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
99    out->A_FC=NULL;                               k=A->pattern->index[iptr_ik];                          
100    out->A_CF=NULL;                               if (out->colorOf[k]==color2) {
101    out->rows_in_F=NULL;                                  A11=out->factors[iptr_ik];
102    out->rows_in_C=NULL;                                  /* a_ij=a_ij-a_ik*a_kj */
103    out->mask_F=NULL;                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
104    out->mask_C=NULL;                                     j=A->pattern->index[iptr_kj];
105    out->x_F=NULL;                                     if (out->colorOf[j]>color2) {
106    out->b_F=NULL;                                        S11=out->factors[iptr_kj];
107    out->x_C=NULL;                                        for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
108    out->b_C=NULL;                                           if (j==A->pattern->index[iptr_ij]) {
109                                                out->factors[iptr_ij]-=A11*S11;
110    if ( !(Paso_checkPtr(mis_marker) || Paso_checkPtr(out) || Paso_checkPtr(counter) ) ) {                                              break;
111       /* identify independend set of rows/columns */                                           }
112       time0=Paso_timer();                                        }
113       Paso_SystemMatrixPattern_mis(A_p->pattern,mis_marker);                                     }
114       time2=Paso_timer()-time0;                                  }
115       if (Paso_noError()) {                               }
116          #pragma omp parallel for private(i) schedule(static)                            }
117          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];                         }
118          out->n=n;                         iptr_main=out->main_iptr[i];
119          out->n_block=n_block;                         D=out->factors[iptr_main];
120          out->n_F=Paso_Util_cumsum(n,counter);                         if (ABS(D)>0.) {
121          out->mask_F=MEMALLOC(n,index_t);                            D=1./D;
122          out->rows_in_F=MEMALLOC(out->n_F,index_t);                            out->factors[iptr_main]=D;
123          out->inv_A_FF=MEMALLOC(n_block*n_block*out->n_F,double);                            /* a_ik=a_ii^{-1}*a_ik */
124          out->A_FF_pivot=NULL; /* later use for block size>3 */                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
125          if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->inv_A_FF) || Paso_checkPtr(out->rows_in_F) ) ) {                               k=A->pattern->index[iptr_ik];
126             #pragma omp parallel                               if (out->colorOf[k]>color) {
127             {                                  A11=out->factors[iptr_ik];
128                /* creates an index for F from mask */                                  out->factors[iptr_ik]=A11*D;
129                #pragma omp for private(i) schedule(static)                               }                              
130                for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1;                            }
               #pragma omp for private(i) schedule(static)  
               for (i = 0; i < n; ++i) {  
                  if  (mis_marker[i]) {  
                         out->rows_in_F[counter[i]]=i;  
                         out->mask_F[i]=counter[i];  
                  } else {  
                         out->mask_F[i]=-1;  
                  }  
               }  
               #pragma omp for private(i, where_p,iPtr,A11,A12,A13,A21,A22,A23,A31,A32,A33,D,index) schedule(static)  
               for (i = 0; i < out->n_F; i++) {  
                 /* find main diagonal */  
                 iPtr=A_p->pattern->ptr[out->rows_in_F[i]];  
                 index=&(A_p->pattern->index[iPtr]);  
                 where_p=(index_t*)bsearch(&out->rows_in_F[i],  
                                         index,  
                                         A_p->pattern->ptr[out->rows_in_F[i] + 1]-A_p->pattern->ptr[out->rows_in_F[i]],  
                                         sizeof(index_t),  
                                         Paso_comparIndex);  
                 if (where_p==NULL) {  
                     Paso_setError(VALUE_ERROR, "__FILE__: main diagonal element missing.");  
                 } else {  
                     iPtr+=(index_t)(where_p-index);  
                     /* get inverse of A_FF block: */  
                     if (n_block==1) {  
                        if (ABS(A_p->val[iPtr])>0.) {  
                             out->inv_A_FF[i]=1./A_p->val[iPtr];  
131                         } else {                         } else {
132                              Paso_setError(ZERO_DIVISION_ERROR, "__FILE__: Break-down in ILU decomposition: non-regular main diagonal block.");                              Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
133                         }                         }
134                      } else if (n_block==2) {                      }
135                         A11=A_p->val[iPtr*4];                   }
136                         A21=A_p->val[iPtr*4+1];                } else if (n_block==2) {
137                         A12=A_p->val[iPtr*4+2];                   #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)
138                         A22=A_p->val[iPtr*4+3];                   for (i = 0; i < n; ++i) {
139                        if (out->colorOf[i]==color) {
140                           for (color2=0;color2<color;++color2) {
141                              for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
142                                 k=A->pattern->index[iptr_ik];                          
143                                 if (out->colorOf[k]==color2) {
144                                    A11=out->factors[iptr_ik*4  ];
145                                    A21=out->factors[iptr_ik*4+1];
146                                    A12=out->factors[iptr_ik*4+2];
147                                    A22=out->factors[iptr_ik*4+3];
148                                    /* a_ij=a_ij-a_ik*a_kj */
149                                    for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
150                                       j=A->pattern->index[iptr_kj];
151                                       if (out->colorOf[j]>color2) {
152                                          S11=out->factors[iptr_kj*4];
153                                          S21=out->factors[iptr_kj*4+1];
154                                          S12=out->factors[iptr_kj*4+2];
155                                          S22=out->factors[iptr_kj*4+3];
156                                          for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
157                                             if (j==A->pattern->index[iptr_ij]) {
158                                                out->factors[4*iptr_ij  ]-=A11*S11+A12*S21;
159                                                out->factors[4*iptr_ij+1]-=A21*S11+A22*S21;
160                                                out->factors[4*iptr_ij+2]-=A11*S12+A12*S22;
161                                                out->factors[4*iptr_ij+3]-=A21*S12+A22*S22;
162                                                break;
163                                             }
164                                          }
165                                       }
166                                    }
167                                 }
168                              }
169                           }
170                           iptr_main=out->main_iptr[i];
171                           A11=out->factors[iptr_main*4];
172                           A21=out->factors[iptr_main*4+1];
173                           A12=out->factors[iptr_main*4+2];
174                           A22=out->factors[iptr_main*4+3];
175                         D = A11*A22-A12*A21;                         D = A11*A22-A12*A21;
176                         if (ABS(D) > 0 ){                         if (ABS(D)>0.) {
177                              D=1./D;                            D=1./D;
178                              out->inv_A_FF[i*4]= A22*D;                            S11= A22*D;
179                              out->inv_A_FF[i*4+1]=-A21*D;                            S21=-A21*D;
180                              out->inv_A_FF[i*4+2]=-A12*D;                            S12=-A12*D;
181                              out->inv_A_FF[i*4+3]= A11*D;                            S22= A11*D;
182                              out->factors[iptr_main*4]  = S11;
183                              out->factors[iptr_main*4+1]= S21;
184                              out->factors[iptr_main*4+2]= S12;
185                              out->factors[iptr_main*4+3]= S22;
186                              /* a_ik=a_ii^{-1}*a_ik */
187                              for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
188                                 k=A->pattern->index[iptr_ik];
189                                 if (out->colorOf[k]>color) {
190                                    A11=out->factors[iptr_ik*4  ];
191                                    A21=out->factors[iptr_ik*4+1];
192                                    A12=out->factors[iptr_ik*4+2];
193                                    A22=out->factors[iptr_ik*4+3];
194                                    out->factors[4*iptr_ik  ]=S11*A11+S12*A21;
195                                    out->factors[4*iptr_ik+1]=S21*A11+S22*A21;
196                                    out->factors[4*iptr_ik+2]=S11*A12+S12*A22;
197                                    out->factors[4*iptr_ik+3]=S21*A12+S22*A22;
198                                 }                              
199                              }
200                         } else {                         } else {
201                              Paso_setError(ZERO_DIVISION_ERROR, "__FILE__:Break-down in ILU decomposition: non-regular main diagonal block.");                              Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
202                           }
203                        }
204                     }
205                  } else if (n_block==3) {
206                     #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)
207                     for (i = 0; i < n; ++i) {
208                        if (out->colorOf[i]==color) {
209                           for (color2=0;color2<color;++color2) {
210                              for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
211                                 k=A->pattern->index[iptr_ik];                          
212                                 if (out->colorOf[k]==color2) {
213                                    A11=out->factors[iptr_ik*9  ];
214                                    A21=out->factors[iptr_ik*9+1];
215                                    A31=out->factors[iptr_ik*9+2];
216                                    A12=out->factors[iptr_ik*9+3];
217                                    A22=out->factors[iptr_ik*9+4];
218                                    A32=out->factors[iptr_ik*9+5];
219                                    A13=out->factors[iptr_ik*9+6];
220                                    A23=out->factors[iptr_ik*9+7];
221                                    A33=out->factors[iptr_ik*9+8];
222                                    /* a_ij=a_ij-a_ik*a_kj */
223                                    for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
224                                       j=A->pattern->index[iptr_kj];
225                                       if (out->colorOf[j]>color2) {
226                                          S11=out->factors[iptr_kj*9  ];
227                                          S21=out->factors[iptr_kj*9+1];
228                                          S31=out->factors[iptr_kj*9+2];
229                                          S12=out->factors[iptr_kj*9+3];
230                                          S22=out->factors[iptr_kj*9+4];
231                                          S32=out->factors[iptr_kj*9+5];
232                                          S13=out->factors[iptr_kj*9+6];
233                                          S23=out->factors[iptr_kj*9+7];
234                                          S33=out->factors[iptr_kj*9+8];                                
235                                          for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
236                                             if (j==A->pattern->index[iptr_ij]) {
237                                                out->factors[iptr_ij*9  ]-=A11*S11+A12*S21+A13*S31;                            
238                                                out->factors[iptr_ij*9+1]-=A21*S11+A22*S21+A23*S31;
239                                                out->factors[iptr_ij*9+2]-=A31*S11+A32*S21+A33*S31;
240                                                out->factors[iptr_ij*9+3]-=A11*S12+A12*S22+A13*S32;                            
241                                                out->factors[iptr_ij*9+4]-=A21*S12+A22*S22+A23*S32;
242                                                out->factors[iptr_ij*9+5]-=A31*S12+A32*S22+A33*S32;
243                                                out->factors[iptr_ij*9+6]-=A11*S13+A12*S23+A13*S33;                            
244                                                out->factors[iptr_ij*9+7]-=A21*S13+A22*S23+A23*S33;
245                                                out->factors[iptr_ij*9+8]-=A31*S13+A32*S23+A33*S33;
246                                                break;
247                                             }
248                                          }
249                                       }
250                                    }
251                                 }
252                              }
253                         }                         }
254                      } else if (n_block==3) {                         iptr_main=out->main_iptr[i];
255                         A11=A_p->val[iPtr*9  ];                         A11=out->factors[iptr_main*9  ];
256                         A21=A_p->val[iPtr*9+1];                         A21=out->factors[iptr_main*9+1];
257                         A31=A_p->val[iPtr*9+2];                         A31=out->factors[iptr_main*9+2];
258                         A12=A_p->val[iPtr*9+3];                         A12=out->factors[iptr_main*9+3];
259                         A22=A_p->val[iPtr*9+4];                         A22=out->factors[iptr_main*9+4];
260                         A32=A_p->val[iPtr*9+5];                         A32=out->factors[iptr_main*9+5];
261                         A13=A_p->val[iPtr*9+6];                         A13=out->factors[iptr_main*9+6];
262                         A23=A_p->val[iPtr*9+7];                         A23=out->factors[iptr_main*9+7];
263                         A33=A_p->val[iPtr*9+8];                         A33=out->factors[iptr_main*9+8];
264                         D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);                         D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
265                         if (ABS(D) > 0 ){                         if (ABS(D)>0.) {
266                              D=1./D;                            D=1./D;
267                              out->inv_A_FF[i*9  ]=(A22*A33-A23*A32)*D;                            S11=(A22*A33-A23*A32)*D;
268                              out->inv_A_FF[i*9+1]=(A31*A23-A21*A33)*D;                            S21=(A31*A23-A21*A33)*D;
269                              out->inv_A_FF[i*9+2]=(A21*A32-A31*A22)*D;                            S31=(A21*A32-A31*A22)*D;
270                              out->inv_A_FF[i*9+3]=(A13*A32-A12*A33)*D;                            S12=(A13*A32-A12*A33)*D;
271                              out->inv_A_FF[i*9+4]=(A11*A33-A31*A13)*D;                            S22=(A11*A33-A31*A13)*D;
272                              out->inv_A_FF[i*9+5]=(A12*A31-A11*A32)*D;                            S32=(A12*A31-A11*A32)*D;
273                              out->inv_A_FF[i*9+6]=(A12*A23-A13*A22)*D;                            S13=(A12*A23-A13*A22)*D;
274                              out->inv_A_FF[i*9+7]=(A13*A21-A11*A23)*D;                            S23=(A13*A21-A11*A23)*D;
275                              out->inv_A_FF[i*9+8]=(A11*A22-A12*A21)*D;                            S33=(A11*A22-A12*A21)*D;
276      
277                              out->factors[iptr_main*9  ]=S11;
278                              out->factors[iptr_main*9+1]=S21;
279                              out->factors[iptr_main*9+2]=S31;
280                              out->factors[iptr_main*9+3]=S12;
281                              out->factors[iptr_main*9+4]=S22;
282                              out->factors[iptr_main*9+5]=S32;
283                              out->factors[iptr_main*9+6]=S13;
284                              out->factors[iptr_main*9+7]=S23;
285                              out->factors[iptr_main*9+8]=S33;
286      
287                              /* a_ik=a_ii^{-1}*a_ik */
288                              for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
289                                 k=A->pattern->index[iptr_ik];
290                                 if (out->colorOf[k]>color) {
291                                    A11=out->factors[iptr_ik*9  ];
292                                    A21=out->factors[iptr_ik*9+1];
293                                    A31=out->factors[iptr_ik*9+2];
294                                    A12=out->factors[iptr_ik*9+3];
295                                    A22=out->factors[iptr_ik*9+4];
296                                    A32=out->factors[iptr_ik*9+5];
297                                    A13=out->factors[iptr_ik*9+6];
298                                    A23=out->factors[iptr_ik*9+7];
299                                    A33=out->factors[iptr_ik*9+8];
300                                    out->factors[iptr_ik*9  ]=S11*A11+S12*A21+S13*A31;                            
301                                    out->factors[iptr_ik*9+1]=S21*A11+S22*A21+S23*A31;
302                                    out->factors[iptr_ik*9+2]=S31*A11+S32*A21+S33*A31;
303                                    out->factors[iptr_ik*9+3]=S11*A12+S12*A22+S13*A32;                            
304                                    out->factors[iptr_ik*9+4]=S21*A12+S22*A22+S23*A32;
305                                    out->factors[iptr_ik*9+5]=S31*A12+S32*A22+S33*A32;
306                                    out->factors[iptr_ik*9+6]=S11*A13+S12*A23+S13*A33;                            
307                                    out->factors[iptr_ik*9+7]=S21*A13+S22*A23+S23*A33;
308                                    out->factors[iptr_ik*9+8]=S31*A13+S32*A23+S33*A33;
309                                 }                              
310                              }
311                         } else {                         } else {
312                              Paso_setError(ZERO_DIVISION_ERROR, "__FILE__:Break-down in ILU decomposition: non-regular main diagonal block.");                              Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
313                         }                         }
314                     }                      }
                 }  
               }  
            } /* end parallel region */  
   
            if( Paso_noError()) {  
               /* if there are no nodes in the coarse level there is no more work to do */  
               out->n_C=n-out->n_F;  
               if (out->n_C>0) {  
                    out->rows_in_C=MEMALLOC(out->n_C,index_t);  
                    out->mask_C=MEMALLOC(n,index_t);  
                    if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) {  
                        /* creates an index for C from mask */  
                        #pragma omp parallel for private(i) schedule(static)  
                        for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];  
                        Paso_Util_cumsum(n,counter);  
                        #pragma omp parallel  
                        {  
                           #pragma omp for private(i) schedule(static)  
                           for (i = 0; i < out->n_C; ++i) out->rows_in_C[i]=-1;  
                           #pragma omp for private(i) schedule(static)  
                           for (i = 0; i < n; ++i) {  
                              if  (! mis_marker[i]) {  
                                 out->rows_in_C[counter[i]]=i;  
                                 out->mask_C[i]=counter[i];  
                              } else {  
                                 out->mask_C[i]=-1;  
                              }  
                           }  
                       } /* end parallel region */  
                       /* get A_CF block: */  
                       out->A_CF=Paso_SystemMatrix_getSubmatrix(A_p,out->n_C,out->rows_in_C,out->mask_F);  
                       if (Paso_noError()) {  
                          /* get A_FC block: */  
                          out->A_FC=Paso_SystemMatrix_getSubmatrix(A_p,out->n_F,out->rows_in_F,out->mask_C);  
                          /* get A_FF block: */  
                          if (Paso_noError()) {  
                             schur=Paso_SystemMatrix_getSubmatrix(A_p,out->n_C,out->rows_in_C,out->mask_C);  
                             time0=Paso_timer()-time0;  
                             if (Paso_noError()) {  
                                 time1=Paso_timer();  
                                 /* update A_CC block to get Schur complement and then apply ILU to it */  
                                 Paso_Solver_updateIncompleteSchurComplement(schur,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);  
                                 time1=Paso_timer()-time1;  
                                 out->ILU_of_Schur=Paso_Solver_getILU(schur,verbose);  
                                 Paso_SystemMatrix_dealloc(schur);  
                             }  
                             /* allocate work arrays for ILU application */  
                             if (Paso_noError()) {  
                               out->x_F=MEMALLOC(n_block*out->n_F,double);  
                               out->b_F=MEMALLOC(n_block*out->n_F,double);  
                               out->x_C=MEMALLOC(n_block*out->n_C,double);  
                               out->b_C=MEMALLOC(n_block*out->n_C,double);  
                               if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {  
                                   #pragma omp parallel  
                                   {  
                                     #pragma omp for private(i,k) schedule(static)  
                                     for (i = 0; i < out->n_F; ++i) {  
                                           for (k=0; k<n_block;++k) {  
                                              out->x_F[i*n_block+k]=0.;  
                                              out->b_F[i*n_block+k]=0.;  
                                           }  
                                     }  
                                     #pragma omp for private(i,k) schedule(static)  
                                     for (i = 0; i < out->n_C; ++i) {  
                                         for (k=0; k<n_block;++k) {  
                                           out->x_C[i*n_block+k]=0.;  
                                           out->b_C[i*n_block+k]=0.;  
                                         }  
                                     }  
                                   } /* end parallel region */  
                               }  
                             }  
                          }  
                      }  
315                   }                   }
316                }                } else {
317             }                   Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: block size greater than 3 is not supported.");
318          }                }      
319                  #pragma omp barrier
320           }
321           time_fac=Paso_timer()-time0;
322       }       }
323    }    }
   TMPMEMFREE(mis_marker);  
   TMPMEMFREE(counter);  
324    if (Paso_noError()) {    if (Paso_noError()) {
325        if (verbose) {        if (verbose) {
326           printf("ILU: %d unknowns eliminated. %d left.\n",out->n_F,n-out->n_F);           printf("ILU: %d color used \n",out->num_colors);
327           if (out->n_C>0) {           printf("timing: ILU: coloring/elemination : %e/%e\n",time_color,time_fac);
             printf("timing: ILU: MIS/reordering/elemination : %e/%e/%e\n",time2,time0,time1);  
          } else {  
             printf("timing: ILU: MIS: %e\n",time2);  
          }  
328       }       }
329       return out;       return out;
330    } else  {    } else  {
# Line 281  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 337  Paso_Solver_ILU* Paso_Solver_getILU(Paso
337    
338  /* apply ILU precondition b-> x                                /* apply ILU precondition b-> x                              
339    
340       in fact it solves       in fact it solves LUx=b in the form x= U^{-1} L^{-1}b
   
   [      I         0  ]  [ A_FF 0 ] [ I    invA_FF*A_FF ]  [ x_F ]  = [b_F]  
   [ A_CF*invA_FF   I  ]  [   0  S ] [ 0          I      ]  [ x_C ]  = [b_C]  
   
  in the form  
   
    b->[b_F,b_C]  
    x_F=invA_FF*b_F  
    b_C=b_C-A_CF*x_F  
    x_C=ILU(b_C)  
    b_F=b_F-A_FC*x_C  
    x_F=invA_FF*b_F  
    x<-[x_F,x_C]  
341    
342   should be called within a parallel region                                                 should be called within a parallel region                                              
343   barrier synconization should be performed to make sure that the input vector available   barrier synconization should be performed to make sure that the input vector available
# Line 302  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 345  Paso_Solver_ILU* Paso_Solver_getILU(Paso
345  */  */
346    
347  void Paso_Solver_solveILU(Paso_Solver_ILU * ilu, double * x, double * b) {  void Paso_Solver_solveILU(Paso_Solver_ILU * ilu, double * x, double * b) {
348       dim_t i,k;       register dim_t i,k;
349         register index_t color,iptr_ik,iptr_main;
350         register double S1,S2,S3,R1,R2,R3;
351       dim_t n_block=ilu->n_block;       dim_t n_block=ilu->n_block;
352         dim_t n=ilu->n;
353            
354       if (ilu->n_C==0) {      
355          /* x=invA_FF*b  */       /* copy x into b*/
356          Paso_Solver_applyBlockDiagonalMatrix(n_block,ilu->n_F,ilu->inv_A_FF,ilu->A_FF_pivot,x,b);       #pragma omp for private(i) schedule(static)
357       } else {       for (i=0;i<n*n_block;++i) x[i]=b[i];
358          /* b->[b_F,b_C]     */       /* forward substitution */
359          if (n_block==1) {       for (color=0;color<ilu->num_colors;++color) {
360             #pragma omp for private(i) schedule(static)             if (n_block==1) {
361             for (i=0;i<ilu->n_F;++i) ilu->b_F[i]=b[ilu->rows_in_F[i]];                #pragma omp for schedule(static) private(i,iptr_ik,k,S1,R1,iptr_main)
362             #pragma omp for private(i) schedule(static)                for (i = 0; i < n; ++i) {
363             for (i=0;i<ilu->n_C;++i) ilu->b_C[i]=b[ilu->rows_in_C[i]];                     if (ilu->colorOf[i]==color) {
364          } else {                       /* x_i=x_i-a_ik*x_k */                    
365             #pragma omp for private(i,k) schedule(static)                       S1=x[i];
366             for (i=0;i<ilu->n_F;++i)                       for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
367                   for (k=0;k<n_block;k++) ilu->b_F[ilu->n_block*i+k]=b[n_block*ilu->rows_in_F[i]+k];                            k=ilu->pattern->index[iptr_ik];                          
368             #pragma omp for private(i,k) schedule(static)                            if (ilu->colorOf[k]<color) {
369             for (i=0;i<ilu->n_C;++i)                               R1=x[k];                              
370                   for (k=0;k<n_block;k++) ilu->b_C[ilu->n_block*i+k]=b[n_block*ilu->rows_in_C[i]+k];                               S1-=ilu->factors[iptr_ik]*R1;
371          }                            }
372          /* x_F=invA_FF*b_F  */                       }
373          Paso_Solver_applyBlockDiagonalMatrix(n_block,ilu->n_F,ilu->inv_A_FF,ilu->A_FF_pivot,ilu->x_F,ilu->b_F);                       iptr_main=ilu->main_iptr[i];
374          /* b_C=b_C-A_CF*x_F */                       x[i]=ilu->factors[iptr_main]*S1;
375          Paso_SystemMatrix_MatrixVector(-1.,ilu->A_CF,ilu->x_F,1.,ilu->b_C);                     }
         /* x_C=ILU(b_C)     */  
         Paso_Solver_solveILU(ilu->ILU_of_Schur,ilu->x_C,ilu->b_C);  
         /* b_F=b_F-A_FC*x_C */  
         Paso_SystemMatrix_MatrixVector(-1.,ilu->A_FC,ilu->x_C,1.,ilu->b_F);  
         /* x_F=invA_FF*b_F  */  
         Paso_Solver_applyBlockDiagonalMatrix(n_block,ilu->n_F,ilu->inv_A_FF,ilu->A_FF_pivot,ilu->x_F,ilu->b_F);  
         /* x<-[x_F,x_C]     */  
         if (n_block==1) {  
            #pragma omp for private(i) schedule(static)  
            for (i=0;i<ilu->n;++i) {  
               if (ilu->mask_C[i]>-1) {  
                   x[i]=ilu->x_C[ilu->mask_C[i]];  
               } else {  
                   x[i]=ilu->x_F[ilu->mask_F[i]];  
376                }                }
377             }             } else if (n_block==2) {
378          } else {                #pragma omp for schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2)
379             #pragma omp for private(i,k) schedule(static)                for (i = 0; i < n; ++i) {
380             for (i=0;i<ilu->n;++i) {                     if (ilu->colorOf[i]==color) {
381                   if (ilu->mask_C[i]>-1) {                       /* x_i=x_i-a_ik*x_k */
382                       for (k=0;k<n_block;k++) x[n_block*i+k]=ilu->x_C[n_block*ilu->mask_C[i]+k];                       S1=x[2*i];
383                   } else {                       S2=x[2*i+1];
384                       for (k=0;k<n_block;k++) x[n_block*i+k]=ilu->x_F[n_block*ilu->mask_F[i]+k];                       for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
385                              k=ilu->pattern->index[iptr_ik];                          
386                              if (ilu->colorOf[k]<color) {
387                                 R1=x[2*k];
388                                 R2=x[2*k+1];
389                                 S1-=ilu->factors[4*iptr_ik  ]*R1+ilu->factors[4*iptr_ik+2]*R2;
390                                 S2-=ilu->factors[4*iptr_ik+1]*R1+ilu->factors[4*iptr_ik+3]*R2;
391                              }
392                         }
393                         iptr_main=ilu->main_iptr[i];
394                         x[2*i  ]=ilu->factors[4*iptr_main  ]*S1+ilu->factors[4*iptr_main+2]*S2;
395                         x[2*i+1]=ilu->factors[4*iptr_main+1]*S1+ilu->factors[4*iptr_main+3]*S2;
396                       }
397    
398                  }
399               } else if (n_block==3) {
400                  #pragma omp for schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3)
401                  for (i = 0; i < n; ++i) {
402                       if (ilu->colorOf[i]==color) {
403                         /* x_i=x_i-a_ik*x_k */
404                         S1=x[3*i];
405                         S2=x[3*i+1];
406                         S3=x[3*i+2];
407                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
408                              k=ilu->pattern->index[iptr_ik];                          
409                              if (ilu->colorOf[k]<color) {
410                                 R1=x[3*k];
411                                 R2=x[3*k+1];
412                                 R3=x[3*k+2];
413                                 S1-=ilu->factors[9*iptr_ik  ]*R1+ilu->factors[9*iptr_ik+3]*R2+ilu->factors[9*iptr_ik+6]*R3;
414                                 S2-=ilu->factors[9*iptr_ik+1]*R1+ilu->factors[9*iptr_ik+4]*R2+ilu->factors[9*iptr_ik+7]*R3;
415                                 S3-=ilu->factors[9*iptr_ik+2]*R1+ilu->factors[9*iptr_ik+5]*R2+ilu->factors[9*iptr_ik+8]*R3;
416                              }
417                         }
418                         iptr_main=ilu->main_iptr[i];
419                         x[3*i  ]=ilu->factors[9*iptr_main  ]*S1+ilu->factors[9*iptr_main+3]*S2+ilu->factors[9*iptr_main+6]*S3;
420                         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;
421                         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;
422                   }                   }
423                  }
424             }             }
425          }             #pragma omp barrier
426          /* all done */       }
427         /* backward substitution */
428         for (color=(ilu->num_colors)-1;color>-1;--color) {
429               if (n_block==1) {
430                  #pragma omp for schedule(static) private(i,iptr_ik,k,S1,R1)
431                  for (i = 0; i < n; ++i) {
432                       if (ilu->colorOf[i]==color) {
433                         /* x_i=x_i-a_ik*x_k */
434                         S1=x[i];
435                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
436                              k=ilu->pattern->index[iptr_ik];                          
437                              if (ilu->colorOf[k]>color) {
438                                 R1=x[k];
439                                 S1-=ilu->factors[iptr_ik]*R1;
440                              }
441                         }
442                         x[i]=S1;
443                       }
444                  }
445               } else if (n_block==2) {
446                  #pragma omp for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2)
447                  for (i = 0; i < n; ++i) {
448                       if (ilu->colorOf[i]==color) {
449                         /* x_i=x_i-a_ik*x_k */
450                         S1=x[2*i];
451                         S2=x[2*i+1];
452                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
453                              k=ilu->pattern->index[iptr_ik];                          
454                              if (ilu->colorOf[k]>color) {
455                                 R1=x[2*k];
456                                 R2=x[2*k+1];
457                                 S1-=ilu->factors[4*iptr_ik  ]*R1+ilu->factors[4*iptr_ik+2]*R2;
458                                 S2-=ilu->factors[4*iptr_ik+1]*R1+ilu->factors[4*iptr_ik+3]*R2;
459                              }
460                         }
461                         x[2*i]=S1;
462                         x[2*i+1]=S2;
463                       }
464                  }
465               } else if (n_block==3) {
466                  #pragma omp for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)
467                  for (i = 0; i < n; ++i) {
468                       if (ilu->colorOf[i]==color) {
469                         /* x_i=x_i-a_ik*x_k */
470                         S1=x[3*i  ];
471                         S2=x[3*i+1];
472                         S3=x[3*i+2];
473                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
474                              k=ilu->pattern->index[iptr_ik];                          
475                              if (ilu->colorOf[k]>color) {
476                                 R1=x[3*k];
477                                 R2=x[3*k+1];
478                                 R3=x[3*k+2];
479                                 S1-=ilu->factors[9*iptr_ik  ]*R1+ilu->factors[9*iptr_ik+3]*R2+ilu->factors[9*iptr_ik+6]*R3;
480                                 S2-=ilu->factors[9*iptr_ik+1]*R1+ilu->factors[9*iptr_ik+4]*R2+ilu->factors[9*iptr_ik+7]*R3;
481                                 S3-=ilu->factors[9*iptr_ik+2]*R1+ilu->factors[9*iptr_ik+5]*R2+ilu->factors[9*iptr_ik+8]*R3;
482                              }
483                         }
484                         x[3*i]=S1;
485                         x[3*i+1]=S2;
486                         x[3*i+2]=S3;
487                       }
488                  }
489             }
490             #pragma omp barrier
491       }       }
492       return;       return;
493  }  }
   
494  /*  /*
495   * $Log$   * $Log$
496   * Revision 1.2  2005/09/15 03:44:40  jgs   * Revision 1.2  2005/09/15 03:44:40  jgs

Legend:
Removed from v.155  
changed lines
  Added in v.1384

  ViewVC Help
Powered by ViewVC 1.1.26