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

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

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

trunk/esys2/paso/src/Solvers/Solver_ILU.c revision 150 by jgs, Thu Sep 15 03:44:45 2005 UTC trunk/paso/src/Solvers/Solver_ILU.c revision 432 by gross, Fri Jan 13 07:38:54 2006 UTC
# Line 21  Line 21 
21    
22  void Paso_Solver_ILU_free(Paso_Solver_ILU * in) {  void Paso_Solver_ILU_free(Paso_Solver_ILU * in) {
23       if (in!=NULL) {       if (in!=NULL) {
24          Paso_Solver_ILU_free(in->ILU_of_Schur);          MEMFREE(in->colorOf);
25          MEMFREE(in->inv_A_FF);          MEMFREE(in->factors);
26          MEMFREE(in->A_FF_pivot);          MEMFREE(in->main_iptr);  
27          Paso_SystemMatrix_dealloc(in->A_FC);          Paso_SystemMatrixPattern_dealloc(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);  
28          MEMFREE(in);          MEMFREE(in);
29       }       }
30  }  }
31    
32  /**************************************************************/  /**************************************************************/
33    
34  /*   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  
35    
36  */  */
37  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SystemMatrix * A_p,bool_t verbose) {  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SystemMatrix * A,bool_t verbose) {
38    Paso_Solver_ILU* out=NULL;    dim_t n=A->num_rows;
39    dim_t n=A_p->num_rows;    dim_t n_block=A->row_block_size;
40    dim_t n_block=A_p->row_block_size;    index_t num_colors=0;
41    index_t* mis_marker=NULL;      register double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;
42    index_t* counter=NULL;      register double mainA11,mainA12,mainA13,mainA21,mainA22,mainA23,mainA31,mainA32,mainA33;
43    index_t iPtr,*index, *where_p;    register double S11,S12,S13,S21,S22,S23,S31,S32,S33;
44    dim_t i,k;    register index_t i,iptr_main,iptr,iptr_ik,k,iptr_kj,j,iptr_ij,color,color2;
45    Paso_SystemMatrix * schur=NULL;    double time0,time_color,time_fac;
46    double A11,A12,A13,A21,A22,A23,A31,A32,A33,D,time0,time1,time2;    /* allocations: */  
47        Paso_Solver_ILU* out=MEMALLOC(1,Paso_Solver_ILU);
48      if (Paso_checkPtr(out)) return NULL;
49    /* identify independend set of rows/columns */    index_t* mis_marker=TMPMEMALLOC(n,index_t);
50    mis_marker=TMPMEMALLOC(n,index_t);    out->colorOf=MEMALLOC(n,index_t);
51    counter=TMPMEMALLOC(n,index_t);    out->factors=MEMALLOC(A->len,double);
52    out=MEMALLOC(1,Paso_Solver_ILU);    out->main_iptr=MEMALLOC(n,index_t);
53    out->ILU_of_Schur=NULL;    out->pattern=Paso_SystemMatrixPattern_reference(A->pattern);
54    out->inv_A_FF=NULL;    out->n_block=n_block;
55    out->A_FF_pivot=NULL;    out->n=n;
56    out->A_FC=NULL;    time0=Paso_timer();
57    out->A_CF=NULL;    if ( !(Paso_checkPtr(mis_marker) ||  
58    out->rows_in_F=NULL;           Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->factors)) ) {
59    out->rows_in_C=NULL;      /* get coloring */
60    out->mask_F=NULL;      index_t num_colors=0;
61    out->mask_C=NULL;      #pragma omp parallel for private(i) schedule(static)
62    out->x_F=NULL;      for (i = 0; i < n; ++i) out->colorOf[i]=-1;
63    out->b_F=NULL;      while (Paso_Util_isAny(n,out->colorOf,-1) && Paso_noError()) {
64    out->x_C=NULL;         #pragma omp parallel for private(i) schedule(static)
65    out->b_C=NULL;         for (i = 0; i < n; ++i) mis_marker[i]=out->colorOf[i];
66           Paso_SystemMatrixPattern_mis(A->pattern,mis_marker);
67    if ( !(Paso_checkPtr(mis_marker) || Paso_checkPtr(out) || Paso_checkPtr(counter) ) ) {  
68       /* identify independend set of rows/columns */         #pragma omp parallel for private(i) schedule(static)
69       time0=Paso_timer();         for (i = 0; i < n; ++i) if (mis_marker[i]) out->colorOf[i]=num_colors;
70       Paso_SystemMatrixPattern_mis(A_p->pattern,mis_marker);         ++num_colors;
71       time2=Paso_timer()-time0;      }
72       if (Paso_noError()) {      out->num_colors=num_colors;
73          #pragma omp parallel for private(i) schedule(static)      time_color=Paso_timer()-time0;
74          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];      time0=Paso_timer();
75          out->n=n;      /* find main diagonal and copy matrix values */
76          out->n_block=n_block;      #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)
77          out->n_F=Paso_Util_cumsum(n,counter);      for (i = 0; i < n; ++i) {
78          out->mask_F=MEMALLOC(n,index_t);          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
79          out->rows_in_F=MEMALLOC(out->n_F,index_t);              iptr_main=A->pattern->ptr[0]-1;
80          out->inv_A_FF=MEMALLOC(n_block*n_block*out->n_F,double);              for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; iptr++) {
81          out->A_FF_pivot=NULL; /* later use for block size>3 */                  if (A->pattern->index[iptr]==i) iptr_main=iptr;
82          if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->inv_A_FF) || Paso_checkPtr(out->rows_in_F) ) ) {                  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             #pragma omp parallel              }
84             {              out->main_iptr[i]=iptr_main;
85                /* creates an index for F from mask */              if (iptr_main==A->pattern->ptr[0]-1)
86                #pragma omp for private(i) schedule(static)                 Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: no main diagonal");
87                for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1;          }
88                #pragma omp for private(i) schedule(static)      }
89        /* start factorization */
90    
91        #pragma omp barrier
92        for (color=0;color<out->num_colors && Paso_noError();++color) {
93               if (n_block==1) {
94                  #pragma omp parallel for schedule(static) private(i,color2,iptr_ik,k,iptr_kj,S11,j,iptr_ij,A11,iptr_main,D)
95                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
96                   if  (mis_marker[i]) {                   if (out->colorOf[i]==color) {
97                          out->rows_in_F[counter[i]]=i;                      for (color2=0;color2<color;++color2) {
98                          out->mask_F[i]=counter[i];                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
99                   } else {                            k=A->pattern->index[iptr_ik];                          
100                          out->mask_F[i]=-1;                            if (out->colorOf[k]==color2) {
101                                 A11=out->factors[iptr_ik];
102                                 /* a_ij=a_ij-a_ik*a_kj */
103                                 for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
104                                    j=A->pattern->index[iptr_kj];
105                                    if (out->colorOf[j]>color2) {
106                                       S11=out->factors[iptr_kj];
107                                       for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
108                                          if (j==A->pattern->index[iptr_ij]) {
109                                             out->factors[iptr_ij]-=A11*S11;
110                                             break;
111                                          }
112                                       }
113                                    }
114                                 }
115                              }
116                           }
117                        }
118                        iptr_main=out->main_iptr[i];
119                        D=out->factors[iptr_main];
120                        if (ABS(D)>0.) {
121                           D=1./D;
122                           out->factors[iptr_main]=D;
123                           /* a_ik=a_ii^{-1}*a_ik */
124                           for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
125                              k=A->pattern->index[iptr_ik];
126                              if (out->colorOf[k]>color) {
127                                 A11=out->factors[iptr_ik];
128                                 out->factors[iptr_ik]=A11*D;
129                              }                              
130                           }
131                        } else {
132                             Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
133                        }
134                   }                   }
135                }                }
136                #pragma omp for private(i, where_p,iPtr,A11,A12,A13,A21,A22,A23,A31,A32,A33,D,index) schedule(static)             } else if (n_block==2) {
137                for (i = 0; i < out->n_F; i++) {                #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                  /* find main diagonal */                for (i = 0; i < n; ++i) {
139                  iPtr=A_p->pattern->ptr[out->rows_in_F[i]];                   if (out->colorOf[i]==color) {
140                  index=&(A_p->pattern->index[iPtr]);                      for (color2=0;color2<color;++color2) {
141                  where_p=(index_t*)bsearch(&out->rows_in_F[i],                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
142                                          index,                            k=A->pattern->index[iptr_ik];                          
143                                          A_p->pattern->ptr[out->rows_in_F[i] + 1]-A_p->pattern->ptr[out->rows_in_F[i]],                            if (out->colorOf[k]==color2) {
144                                          sizeof(index_t),                               A11=out->factors[iptr_ik*4  ];
145                                          Paso_comparIndex);                               A21=out->factors[iptr_ik*4+1];
146                  if (where_p==NULL) {                               A12=out->factors[iptr_ik*4+2];
147                      Paso_setError(VALUE_ERROR, "__FILE__: main diagonal element missing.");                               A22=out->factors[iptr_ik*4+3];
148                  } else {                               /* a_ij=a_ij-a_ik*a_kj */
149                      iPtr+=(index_t)(where_p-index);                               for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
150                      /* get inverse of A_FF block: */                                  j=A->pattern->index[iptr_kj];
151                      if (n_block==1) {                                  if (out->colorOf[j]>color2) {
152                         if (ABS(A_p->val[iPtr])>0.) {                                     S11=out->factors[iptr_kj*4];
153                              out->inv_A_FF[i]=1./A_p->val[iPtr];                                     S21=out->factors[iptr_kj*4+1];
154                         } else {                                     S12=out->factors[iptr_kj*4+2];
155                              Paso_setError(ZERO_DIVISION_ERROR, "__FILE__: Break-down in ILU decomposition: non-regular main diagonal block.");                                     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                      } else if (n_block==2) {                                        if (j==A->pattern->index[iptr_ij]) {
158                         A11=A_p->val[iPtr*4];                                           out->factors[4*iptr_ij  ]-=A11*S11+A12*S21;
159                         A21=A_p->val[iPtr*4+1];                                           out->factors[4*iptr_ij+1]-=A21*S11+A22*S21;
160                         A12=A_p->val[iPtr*4+2];                                           out->factors[4*iptr_ij+2]-=A11*S12+A12*S22;
161                         A22=A_p->val[iPtr*4+3];                                           out->factors[4*iptr_ij+3]-=A21*S12+A22*S22;
162                         D = A11*A22-A12*A21;                                           break;
163                         if (ABS(D) > 0 ){                                        }
164                              D=1./D;                                     }
165                              out->inv_A_FF[i*4]= A22*D;                                  }
166                              out->inv_A_FF[i*4+1]=-A21*D;                               }
167                              out->inv_A_FF[i*4+2]=-A12*D;                            }
                             out->inv_A_FF[i*4+3]= A11*D;  
                        } else {  
                             Paso_setError(ZERO_DIVISION_ERROR, "__FILE__:Break-down in ILU decomposition: non-regular main diagonal block.");  
168                         }                         }
169                      } else if (n_block==3) {                      }
170                         A11=A_p->val[iPtr*9  ];                      iptr_main=out->main_iptr[i];
171                         A21=A_p->val[iPtr*9+1];                      A11=out->factors[iptr_main*4];
172                         A31=A_p->val[iPtr*9+2];                      A21=out->factors[iptr_main*4+1];
173                         A12=A_p->val[iPtr*9+3];                      A12=out->factors[iptr_main*4+2];
174                         A22=A_p->val[iPtr*9+4];                      A22=out->factors[iptr_main*4+3];
175                         A32=A_p->val[iPtr*9+5];                      D = A11*A22-A12*A21;
176                         A13=A_p->val[iPtr*9+6];                      if (ABS(D)>0.) {
177                         A23=A_p->val[iPtr*9+7];                         D=1./D;
178                         A33=A_p->val[iPtr*9+8];                         S11= A22*D;
179                         D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);                         S21=-A21*D;
180                         if (ABS(D) > 0 ){                         S12=-A12*D;
181                              D=1./D;                         S22= A11*D;
182                              out->inv_A_FF[i*9  ]=(A22*A33-A23*A32)*D;                         out->factors[iptr_main*4]  = S11;
183                              out->inv_A_FF[i*9+1]=(A31*A23-A21*A33)*D;                         out->factors[iptr_main*4+1]= S21;
184                              out->inv_A_FF[i*9+2]=(A21*A32-A31*A22)*D;                         out->factors[iptr_main*4+2]= S12;
185                              out->inv_A_FF[i*9+3]=(A13*A32-A12*A33)*D;                         out->factors[iptr_main*4+3]= S22;
186                              out->inv_A_FF[i*9+4]=(A11*A33-A31*A13)*D;                         /* a_ik=a_ii^{-1}*a_ik */
187                              out->inv_A_FF[i*9+5]=(A12*A31-A11*A32)*D;                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
188                              out->inv_A_FF[i*9+6]=(A12*A23-A13*A22)*D;                            k=A->pattern->index[iptr_ik];
189                              out->inv_A_FF[i*9+7]=(A13*A21-A11*A23)*D;                            if (out->colorOf[k]>color) {
190                              out->inv_A_FF[i*9+8]=(A11*A22-A12*A21)*D;                               A11=out->factors[iptr_ik*4  ];
191                         } else {                               A21=out->factors[iptr_ik*4+1];
192                              Paso_setError(ZERO_DIVISION_ERROR, "__FILE__:Break-down in ILU decomposition: non-regular main diagonal block.");                               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 {
201                  }                           Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
202                        }
203                     }
204                }                }
205             } /* end parallel region */             } 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             if( Paso_noError()) {                for (i = 0; i < n; ++i) {
208                /* if there are no nodes in the coarse level there is no more work to do */                   if (out->colorOf[i]==color) {
209                out->n_C=n-out->n_F;                      for (color2=0;color2<color;++color2) {
210                if (out->n_C>0) {                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
211                     out->rows_in_C=MEMALLOC(out->n_C,index_t);                            k=A->pattern->index[iptr_ik];                          
212                     out->mask_C=MEMALLOC(n,index_t);                            if (out->colorOf[k]==color2) {
213                     if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) {                               A11=out->factors[iptr_ik*9  ];
214                         /* creates an index for C from mask */                               A21=out->factors[iptr_ik*9+1];
215                         #pragma omp parallel for private(i) schedule(static)                               A31=out->factors[iptr_ik*9+2];
216                         for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];                               A12=out->factors[iptr_ik*9+3];
217                         Paso_Util_cumsum(n,counter);                               A22=out->factors[iptr_ik*9+4];
218                         #pragma omp parallel                               A32=out->factors[iptr_ik*9+5];
219                         {                               A13=out->factors[iptr_ik*9+6];
220                            #pragma omp for private(i) schedule(static)                               A23=out->factors[iptr_ik*9+7];
221                            for (i = 0; i < out->n_C; ++i) out->rows_in_C[i]=-1;                               A33=out->factors[iptr_ik*9+8];
222                            #pragma omp for private(i) schedule(static)                               /* a_ij=a_ij-a_ik*a_kj */
223                            for (i = 0; i < n; ++i) {                               for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
224                               if  (! mis_marker[i]) {                                  j=A->pattern->index[iptr_kj];
225                                  out->rows_in_C[counter[i]]=i;                                  if (out->colorOf[j]>color2) {
226                                  out->mask_C[i]=counter[i];                                     S11=out->factors[iptr_kj*9  ];
227                               } else {                                     S21=out->factors[iptr_kj*9+1];
228                                  out->mask_C[i]=-1;                                     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                        } /* end parallel region */                         }
253                        /* get A_CF block: */                      }
254                        out->A_CF=Paso_SystemMatrix_getSubmatrix(A_p,out->n_C,out->rows_in_C,out->mask_F);                      iptr_main=out->main_iptr[i];
255                        if (Paso_noError()) {                      A11=out->factors[iptr_main*9  ];
256                           /* get A_FC block: */                      A21=out->factors[iptr_main*9+1];
257                           out->A_FC=Paso_SystemMatrix_getSubmatrix(A_p,out->n_F,out->rows_in_F,out->mask_C);                      A31=out->factors[iptr_main*9+2];
258                           /* get A_FF block: */                      A12=out->factors[iptr_main*9+3];
259                           if (Paso_noError()) {                      A22=out->factors[iptr_main*9+4];
260                              schur=Paso_SystemMatrix_getSubmatrix(A_p,out->n_C,out->rows_in_C,out->mask_C);                      A32=out->factors[iptr_main*9+5];
261                              time0=Paso_timer()-time0;                      A13=out->factors[iptr_main*9+6];
262                              if (Paso_noError()) {                      A23=out->factors[iptr_main*9+7];
263                                  time1=Paso_timer();                      A33=out->factors[iptr_main*9+8];
264                                  /* update A_CC block to get Schur complement and then apply ILU to it */                      D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
265                                  Paso_Solver_updateIncompleteSchurComplement(schur,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);                      if (ABS(D)>0.) {
266                                  time1=Paso_timer()-time1;                         D=1./D;
267                                  out->ILU_of_Schur=Paso_Solver_getILU(schur,verbose);                         S11=(A22*A33-A23*A32)*D;
268                                  Paso_SystemMatrix_dealloc(schur);                         S21=(A31*A23-A21*A33)*D;
269                              }                         S31=(A21*A32-A31*A22)*D;
270                              /* allocate work arrays for ILU application */                         S12=(A13*A32-A12*A33)*D;
271                              if (Paso_noError()) {                         S22=(A11*A33-A31*A13)*D;
272                                out->x_F=MEMALLOC(n_block*out->n_F,double);                         S32=(A12*A31-A11*A32)*D;
273                                out->b_F=MEMALLOC(n_block*out->n_F,double);                         S13=(A12*A23-A13*A22)*D;
274                                out->x_C=MEMALLOC(n_block*out->n_C,double);                         S23=(A13*A21-A11*A23)*D;
275                                out->b_C=MEMALLOC(n_block*out->n_C,double);                         S33=(A11*A22-A12*A21)*D;
276                                if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {  
277                                    #pragma omp parallel                         out->factors[iptr_main*9  ]=S11;
278                                    {                         out->factors[iptr_main*9+1]=S21;
279                                      #pragma omp for private(i,k) schedule(static)                         out->factors[iptr_main*9+2]=S31;
280                                      for (i = 0; i < out->n_F; ++i) {                         out->factors[iptr_main*9+3]=S12;
281                                            for (k=0; k<n_block;++k) {                         out->factors[iptr_main*9+4]=S22;
282                                               out->x_F[i*n_block+k]=0.;                         out->factors[iptr_main*9+5]=S32;
283                                               out->b_F[i*n_block+k]=0.;                         out->factors[iptr_main*9+6]=S13;
284                                            }                         out->factors[iptr_main*9+7]=S23;
285                                      }                         out->factors[iptr_main*9+8]=S33;
286                                      #pragma omp for private(i,k) schedule(static)  
287                                      for (i = 0; i < out->n_C; ++i) {                         /* a_ik=a_ii^{-1}*a_ik */
288                                          for (k=0; k<n_block;++k) {                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
289                                            out->x_C[i*n_block+k]=0.;                            k=A->pattern->index[iptr_ik];
290                                            out->b_C[i*n_block+k]=0.;                            if (out->colorOf[k]>color) {
291                                          }                               A11=out->factors[iptr_ik*9  ];
292                                      }                               A21=out->factors[iptr_ik*9+1];
293                                    } /* end parallel region */                               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 {
312                             Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
313                        }
314                   }                   }
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    
324    TMPMEMFREE(mis_marker);    TMPMEMFREE(mis_marker);
   TMPMEMFREE(counter);  
325    if (Paso_noError()) {    if (Paso_noError()) {
326        if (verbose) {        if (verbose) {
327           printf("ILU: %d unknowns eliminated. %d left.\n",out->n_F,n-out->n_F);           printf("ILU: %d color used \n",out->num_colors);
328           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);  
          }  
329       }       }
330       return out;       return out;
331    } else  {    } else  {
# Line 281  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 338  Paso_Solver_ILU* Paso_Solver_getILU(Paso
338    
339  /* apply ILU precondition b-> x                                /* apply ILU precondition b-> x                              
340    
341       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]  
342    
343   should be called within a parallel region                                                 should be called within a parallel region                                              
344   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 346  Paso_Solver_ILU* Paso_Solver_getILU(Paso
346  */  */
347    
348  void Paso_Solver_solveILU(Paso_Solver_ILU * ilu, double * x, double * b) {  void Paso_Solver_solveILU(Paso_Solver_ILU * ilu, double * x, double * b) {
349       dim_t i,k;       register dim_t i,k;
350         register index_t color,iptr_ik,iptr_main;
351         register double S1,S2,S3,R1,R2,R3;
352       dim_t n_block=ilu->n_block;       dim_t n_block=ilu->n_block;
353         dim_t n=ilu->n;
354            
355       if (ilu->n_C==0) {      
356          /* x=invA_FF*b  */       /* copy x into b*/
357          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)
358       } else {       for (i=0;i<n*n_block;++i) x[i]=b[i];
359          /* b->[b_F,b_C]     */       /* forward substitution */
360          if (n_block==1) {       for (color=0;color<ilu->num_colors;++color) {
361             #pragma omp for private(i) schedule(static)             if (n_block==1) {
362             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)
363             #pragma omp for private(i) schedule(static)                for (i = 0; i < n; ++i) {
364             for (i=0;i<ilu->n_C;++i) ilu->b_C[i]=b[ilu->rows_in_C[i]];                     if (ilu->colorOf[i]==color) {
365          } else {                       /* x_i=x_i-a_ik*x_k */                    
366             #pragma omp for private(i,k) schedule(static)                       S1=x[i];
367             for (i=0;i<ilu->n_F;++i)                       for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
368                   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];                          
369             #pragma omp for private(i,k) schedule(static)                            if (ilu->colorOf[k]<color) {
370             for (i=0;i<ilu->n_C;++i)                               R1=x[k];                              
371                   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;
372          }                            }
373          /* x_F=invA_FF*b_F  */                       }
374          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];
375          /* b_C=b_C-A_CF*x_F */                       x[i]=ilu->factors[iptr_main]*S1;
376          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]];  
377                }                }
378             }             } else if (n_block==2) {
379          } else {                #pragma omp for schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2)
380             #pragma omp for private(i,k) schedule(static)                for (i = 0; i < n; ++i) {
381             for (i=0;i<ilu->n;++i) {                     if (ilu->colorOf[i]==color) {
382                   if (ilu->mask_C[i]>-1) {                       /* x_i=x_i-a_ik*x_k */
383                       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];
384                   } else {                       S2=x[2*i+1];
385                       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) {
386                              k=ilu->pattern->index[iptr_ik];                          
387                              if (ilu->colorOf[k]<color) {
388                                 R1=x[2*k];
389                                 R2=x[2*k+1];
390                                 S1-=ilu->factors[4*iptr_ik  ]*R1+ilu->factors[4*iptr_ik+2]*R2;
391                                 S2-=ilu->factors[4*iptr_ik+1]*R1+ilu->factors[4*iptr_ik+3]*R2;
392                              }
393                         }
394                         iptr_main=ilu->main_iptr[i];
395                         x[2*i  ]=ilu->factors[4*iptr_main  ]*S1+ilu->factors[4*iptr_main+2]*S2;
396                         x[2*i+1]=ilu->factors[4*iptr_main+1]*S1+ilu->factors[4*iptr_main+3]*S2;
397                       }
398    
399                  }
400               } else if (n_block==3) {
401                  #pragma omp for schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3)
402                  for (i = 0; i < n; ++i) {
403                       if (ilu->colorOf[i]==color) {
404                         /* x_i=x_i-a_ik*x_k */
405                         S1=x[3*i];
406                         S2=x[3*i+1];
407                         S3=x[3*i+2];
408                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
409                              k=ilu->pattern->index[iptr_ik];                          
410                              if (ilu->colorOf[k]<color) {
411                                 R1=x[3*k];
412                                 R2=x[3*k+1];
413                                 R3=x[3*k+2];
414                                 S1-=ilu->factors[9*iptr_ik  ]*R1+ilu->factors[9*iptr_ik+3]*R2+ilu->factors[9*iptr_ik+6]*R3;
415                                 S2-=ilu->factors[9*iptr_ik+1]*R1+ilu->factors[9*iptr_ik+4]*R2+ilu->factors[9*iptr_ik+7]*R3;
416                                 S3-=ilu->factors[9*iptr_ik+2]*R1+ilu->factors[9*iptr_ik+5]*R2+ilu->factors[9*iptr_ik+8]*R3;
417                              }
418                         }
419                         iptr_main=ilu->main_iptr[i];
420                         x[3*i  ]=ilu->factors[9*iptr_main  ]*S1+ilu->factors[9*iptr_main+3]*S2+ilu->factors[9*iptr_main+6]*S3;
421                         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;
422                         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;
423                   }                   }
424                  }
425             }             }
426          }             #pragma omp barrier
427          /* all done */       }
428         /* backward substitution */
429         for (color=(ilu->num_colors)-1;color>-1;--color) {
430               if (n_block==1) {
431                  #pragma omp for schedule(static) private(i,iptr_ik,k,S1,R1)
432                  for (i = 0; i < n; ++i) {
433                       if (ilu->colorOf[i]==color) {
434                         /* x_i=x_i-a_ik*x_k */
435                         S1=x[i];
436                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
437                              k=ilu->pattern->index[iptr_ik];                          
438                              if (ilu->colorOf[k]>color) {
439                                 R1=x[k];
440                                 S1-=ilu->factors[iptr_ik]*R1;
441                              }
442                         }
443                         x[i]=S1;
444                       }
445                  }
446               } else if (n_block==2) {
447                  #pragma omp for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2)
448                  for (i = 0; i < n; ++i) {
449                       if (ilu->colorOf[i]==color) {
450                         /* x_i=x_i-a_ik*x_k */
451                         S1=x[2*i];
452                         S2=x[2*i+1];
453                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
454                              k=ilu->pattern->index[iptr_ik];                          
455                              if (ilu->colorOf[k]>color) {
456                                 R1=x[2*k];
457                                 R2=x[2*k+1];
458                                 S1-=ilu->factors[4*iptr_ik  ]*R1+ilu->factors[4*iptr_ik+2]*R2;
459                                 S2-=ilu->factors[4*iptr_ik+1]*R1+ilu->factors[4*iptr_ik+3]*R2;
460                              }
461                         }
462                         x[2*i]=S1;
463                         x[2*i+1]=S2;
464                       }
465                  }
466               } else if (n_block==3) {
467                  #pragma omp for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)
468                  for (i = 0; i < n; ++i) {
469                       if (ilu->colorOf[i]==color) {
470                         /* x_i=x_i-a_ik*x_k */
471                         S1=x[3*i  ];
472                         S2=x[3*i+1];
473                         S3=x[3*i+2];
474                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
475                              k=ilu->pattern->index[iptr_ik];                          
476                              if (ilu->colorOf[k]>color) {
477                                 R1=x[3*k];
478                                 R2=x[3*k+1];
479                                 R3=x[3*k+2];
480                                 S1-=ilu->factors[9*iptr_ik  ]*R1+ilu->factors[9*iptr_ik+3]*R2+ilu->factors[9*iptr_ik+6]*R3;
481                                 S2-=ilu->factors[9*iptr_ik+1]*R1+ilu->factors[9*iptr_ik+4]*R2+ilu->factors[9*iptr_ik+7]*R3;
482                                 S3-=ilu->factors[9*iptr_ik+2]*R1+ilu->factors[9*iptr_ik+5]*R2+ilu->factors[9*iptr_ik+8]*R3;
483                              }
484                         }
485                         x[3*i]=S1;
486                         x[3*i+1]=S2;
487                         x[3*i+2]=S3;
488                       }
489                  }
490             }
491             #pragma omp barrier
492       }       }
493       return;       return;
494  }  }
   
495  /*  /*
496   * $Log$   * $Log$
497   * Revision 1.2  2005/09/15 03:44:40  jgs   * Revision 1.2  2005/09/15 03:44:40  jgs

Legend:
Removed from v.150  
changed lines
  Added in v.432

  ViewVC Help
Powered by ViewVC 1.1.26