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

  ViewVC Help
Powered by ViewVC 1.1.26