/[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 trunk/paso/src/Solver_ILU.c revision 1628 by phornby, Fri Jul 11 13:12:46 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 S11,S12,S13,S21,S22,S23,S31,S32,S33;
56    index_t iPtr,*index, *where_p;    register index_t i,iptr_main,iptr,iptr_ik,k,iptr_kj,j,iptr_ij,color,color2;
57    dim_t i,k;    double time0,time_color,time_fac;
58    Paso_SystemMatrix * schur=NULL;    /* allocations: */  
59    double A11,A12,A13,A21,A22,A23,A31,A32,A33,D,time0,time1,time2;    Paso_Solver_ILU* out=MEMALLOC(1,Paso_Solver_ILU);
60      if (Paso_checkPtr(out)) return NULL;
61      out->colorOf=MEMALLOC(n,index_t);
62      out->factors=MEMALLOC(A->len,double);
63      out->main_iptr=MEMALLOC(n,index_t);
64      out->pattern=Paso_Pattern_getReference(A->pattern);
65      out->n_block=n_block;
66      out->n=n;
67      if ( !(Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->factors)) ) {
68        time0=Paso_timer();
69        Paso_Pattern_color(A->pattern,&out->num_colors,out->colorOf);
70        time_color=Paso_timer()-time0;
71    
72        if (Paso_noError()) {
73           time0=Paso_timer();
74           /* find main diagonal and copy matrix values */
75           #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)
76           for (i = 0; i < n; ++i) {
77               for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
78                   iptr_main=A->pattern->ptr[0]-1;
79                    for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; iptr++) {
80                       if (A->pattern->index[iptr]==i) iptr_main=iptr;
81                       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];
82                   }
83                   out->main_iptr[i]=iptr_main;
84                   if (iptr_main==A->pattern->ptr[0]-1)
85                      Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: no main diagonal");
86               }
87           }
88           /* start factorization */
89        
90           #pragma omp barrier
91    /* identify independend set of rows/columns */         for (color=0;color<out->num_colors && Paso_noError();++color) {
92    mis_marker=TMPMEMALLOC(n,index_t);                if (n_block==1) {
93    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)
94    out=MEMALLOC(1,Paso_Solver_ILU);                   for (i = 0; i < n; ++i) {
95    out->ILU_of_Schur=NULL;                      if (out->colorOf[i]==color) {
96    out->inv_A_FF=NULL;                         for (color2=0;color2<color;++color2) {
97    out->A_FF_pivot=NULL;                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
98    out->A_FC=NULL;                               k=A->pattern->index[iptr_ik];                          
99    out->A_CF=NULL;                               if (out->colorOf[k]==color2) {
100    out->rows_in_F=NULL;                                  A11=out->factors[iptr_ik];
101    out->rows_in_C=NULL;                                  /* a_ij=a_ij-a_ik*a_kj */
102    out->mask_F=NULL;                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
103    out->mask_C=NULL;                                     j=A->pattern->index[iptr_kj];
104    out->x_F=NULL;                                     if (out->colorOf[j]>color2) {
105    out->b_F=NULL;                                        S11=out->factors[iptr_kj];
106    out->x_C=NULL;                                        for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
107    out->b_C=NULL;                                           if (j==A->pattern->index[iptr_ij]) {
108                                                out->factors[iptr_ij]-=A11*S11;
109    if ( !(Paso_checkPtr(mis_marker) || Paso_checkPtr(out) || Paso_checkPtr(counter) ) ) {                                              break;
110       /* identify independend set of rows/columns */                                           }
111       time0=Paso_timer();                                        }
112       Paso_SystemMatrixPattern_mis(A_p->pattern,mis_marker);                                     }
113       time2=Paso_timer()-time0;                                  }
114       if (Paso_noError()) {                               }
115          #pragma omp parallel for private(i) schedule(static)                            }
116          for (i = 0; i < n; ++i) counter[i]=mis_marker[i];                         }
117          out->n=n;                         iptr_main=out->main_iptr[i];
118          out->n_block=n_block;                         D=out->factors[iptr_main];
119          out->n_F=Paso_Util_cumsum(n,counter);                         if (ABS(D)>0.) {
120          out->mask_F=MEMALLOC(n,index_t);                            D=1./D;
121          out->rows_in_F=MEMALLOC(out->n_F,index_t);                            out->factors[iptr_main]=D;
122          out->inv_A_FF=MEMALLOC(n_block*n_block*out->n_F,double);                            /* a_ik=a_ii^{-1}*a_ik */
123          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) {
124          if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->inv_A_FF) || Paso_checkPtr(out->rows_in_F) ) ) {                               k=A->pattern->index[iptr_ik];
125             #pragma omp parallel                               if (out->colorOf[k]>color) {
126             {                                  A11=out->factors[iptr_ik];
127                /* creates an index for F from mask */                                  out->factors[iptr_ik]=A11*D;
128                #pragma omp for private(i) schedule(static)                               }                              
129                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];  
130                         } else {                         } else {
131                              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.");
132                         }                         }
133                      } else if (n_block==2) {                      }
134                         A11=A_p->val[iPtr*4];                   }
135                         A21=A_p->val[iPtr*4+1];                } else if (n_block==2) {
136                         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)
137                         A22=A_p->val[iPtr*4+3];                   for (i = 0; i < n; ++i) {
138                        if (out->colorOf[i]==color) {
139                           for (color2=0;color2<color;++color2) {
140                              for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
141                                 k=A->pattern->index[iptr_ik];                          
142                                 if (out->colorOf[k]==color2) {
143                                    A11=out->factors[iptr_ik*4  ];
144                                    A21=out->factors[iptr_ik*4+1];
145                                    A12=out->factors[iptr_ik*4+2];
146                                    A22=out->factors[iptr_ik*4+3];
147                                    /* a_ij=a_ij-a_ik*a_kj */
148                                    for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
149                                       j=A->pattern->index[iptr_kj];
150                                       if (out->colorOf[j]>color2) {
151                                          S11=out->factors[iptr_kj*4];
152                                          S21=out->factors[iptr_kj*4+1];
153                                          S12=out->factors[iptr_kj*4+2];
154                                          S22=out->factors[iptr_kj*4+3];
155                                          for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
156                                             if (j==A->pattern->index[iptr_ij]) {
157                                                out->factors[4*iptr_ij  ]-=A11*S11+A12*S21;
158                                                out->factors[4*iptr_ij+1]-=A21*S11+A22*S21;
159                                                out->factors[4*iptr_ij+2]-=A11*S12+A12*S22;
160                                                out->factors[4*iptr_ij+3]-=A21*S12+A22*S22;
161                                                break;
162                                             }
163                                          }
164                                       }
165                                    }
166                                 }
167                              }
168                           }
169                           iptr_main=out->main_iptr[i];
170                           A11=out->factors[iptr_main*4];
171                           A21=out->factors[iptr_main*4+1];
172                           A12=out->factors[iptr_main*4+2];
173                           A22=out->factors[iptr_main*4+3];
174                         D = A11*A22-A12*A21;                         D = A11*A22-A12*A21;
175                         if (ABS(D) > 0 ){                         if (ABS(D)>0.) {
176                              D=1./D;                            D=1./D;
177                              out->inv_A_FF[i*4]= A22*D;                            S11= A22*D;
178                              out->inv_A_FF[i*4+1]=-A21*D;                            S21=-A21*D;
179                              out->inv_A_FF[i*4+2]=-A12*D;                            S12=-A12*D;
180                              out->inv_A_FF[i*4+3]= A11*D;                            S22= A11*D;
181                              out->factors[iptr_main*4]  = S11;
182                              out->factors[iptr_main*4+1]= S21;
183                              out->factors[iptr_main*4+2]= S12;
184                              out->factors[iptr_main*4+3]= S22;
185                              /* a_ik=a_ii^{-1}*a_ik */
186                              for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
187                                 k=A->pattern->index[iptr_ik];
188                                 if (out->colorOf[k]>color) {
189                                    A11=out->factors[iptr_ik*4  ];
190                                    A21=out->factors[iptr_ik*4+1];
191                                    A12=out->factors[iptr_ik*4+2];
192                                    A22=out->factors[iptr_ik*4+3];
193                                    out->factors[4*iptr_ik  ]=S11*A11+S12*A21;
194                                    out->factors[4*iptr_ik+1]=S21*A11+S22*A21;
195                                    out->factors[4*iptr_ik+2]=S11*A12+S12*A22;
196                                    out->factors[4*iptr_ik+3]=S21*A12+S22*A22;
197                                 }                              
198                              }
199                         } else {                         } else {
200                              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.");
201                           }
202                        }
203                     }
204                  } else if (n_block==3) {
205                     #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)
206                     for (i = 0; i < n; ++i) {
207                        if (out->colorOf[i]==color) {
208                           for (color2=0;color2<color;++color2) {
209                              for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
210                                 k=A->pattern->index[iptr_ik];                          
211                                 if (out->colorOf[k]==color2) {
212                                    A11=out->factors[iptr_ik*9  ];
213                                    A21=out->factors[iptr_ik*9+1];
214                                    A31=out->factors[iptr_ik*9+2];
215                                    A12=out->factors[iptr_ik*9+3];
216                                    A22=out->factors[iptr_ik*9+4];
217                                    A32=out->factors[iptr_ik*9+5];
218                                    A13=out->factors[iptr_ik*9+6];
219                                    A23=out->factors[iptr_ik*9+7];
220                                    A33=out->factors[iptr_ik*9+8];
221                                    /* a_ij=a_ij-a_ik*a_kj */
222                                    for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
223                                       j=A->pattern->index[iptr_kj];
224                                       if (out->colorOf[j]>color2) {
225                                          S11=out->factors[iptr_kj*9  ];
226                                          S21=out->factors[iptr_kj*9+1];
227                                          S31=out->factors[iptr_kj*9+2];
228                                          S12=out->factors[iptr_kj*9+3];
229                                          S22=out->factors[iptr_kj*9+4];
230                                          S32=out->factors[iptr_kj*9+5];
231                                          S13=out->factors[iptr_kj*9+6];
232                                          S23=out->factors[iptr_kj*9+7];
233                                          S33=out->factors[iptr_kj*9+8];                                
234                                          for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
235                                             if (j==A->pattern->index[iptr_ij]) {
236                                                out->factors[iptr_ij*9  ]-=A11*S11+A12*S21+A13*S31;                            
237                                                out->factors[iptr_ij*9+1]-=A21*S11+A22*S21+A23*S31;
238                                                out->factors[iptr_ij*9+2]-=A31*S11+A32*S21+A33*S31;
239                                                out->factors[iptr_ij*9+3]-=A11*S12+A12*S22+A13*S32;                            
240                                                out->factors[iptr_ij*9+4]-=A21*S12+A22*S22+A23*S32;
241                                                out->factors[iptr_ij*9+5]-=A31*S12+A32*S22+A33*S32;
242                                                out->factors[iptr_ij*9+6]-=A11*S13+A12*S23+A13*S33;                            
243                                                out->factors[iptr_ij*9+7]-=A21*S13+A22*S23+A23*S33;
244                                                out->factors[iptr_ij*9+8]-=A31*S13+A32*S23+A33*S33;
245                                                break;
246                                             }
247                                          }
248                                       }
249                                    }
250                                 }
251                              }
252                         }                         }
253                      } else if (n_block==3) {                         iptr_main=out->main_iptr[i];
254                         A11=A_p->val[iPtr*9  ];                         A11=out->factors[iptr_main*9  ];
255                         A21=A_p->val[iPtr*9+1];                         A21=out->factors[iptr_main*9+1];
256                         A31=A_p->val[iPtr*9+2];                         A31=out->factors[iptr_main*9+2];
257                         A12=A_p->val[iPtr*9+3];                         A12=out->factors[iptr_main*9+3];
258                         A22=A_p->val[iPtr*9+4];                         A22=out->factors[iptr_main*9+4];
259                         A32=A_p->val[iPtr*9+5];                         A32=out->factors[iptr_main*9+5];
260                         A13=A_p->val[iPtr*9+6];                         A13=out->factors[iptr_main*9+6];
261                         A23=A_p->val[iPtr*9+7];                         A23=out->factors[iptr_main*9+7];
262                         A33=A_p->val[iPtr*9+8];                         A33=out->factors[iptr_main*9+8];
263                         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);
264                         if (ABS(D) > 0 ){                         if (ABS(D)>0.) {
265                              D=1./D;                            D=1./D;
266                              out->inv_A_FF[i*9  ]=(A22*A33-A23*A32)*D;                            S11=(A22*A33-A23*A32)*D;
267                              out->inv_A_FF[i*9+1]=(A31*A23-A21*A33)*D;                            S21=(A31*A23-A21*A33)*D;
268                              out->inv_A_FF[i*9+2]=(A21*A32-A31*A22)*D;                            S31=(A21*A32-A31*A22)*D;
269                              out->inv_A_FF[i*9+3]=(A13*A32-A12*A33)*D;                            S12=(A13*A32-A12*A33)*D;
270                              out->inv_A_FF[i*9+4]=(A11*A33-A31*A13)*D;                            S22=(A11*A33-A31*A13)*D;
271                              out->inv_A_FF[i*9+5]=(A12*A31-A11*A32)*D;                            S32=(A12*A31-A11*A32)*D;
272                              out->inv_A_FF[i*9+6]=(A12*A23-A13*A22)*D;                            S13=(A12*A23-A13*A22)*D;
273                              out->inv_A_FF[i*9+7]=(A13*A21-A11*A23)*D;                            S23=(A13*A21-A11*A23)*D;
274                              out->inv_A_FF[i*9+8]=(A11*A22-A12*A21)*D;                            S33=(A11*A22-A12*A21)*D;
275      
276                              out->factors[iptr_main*9  ]=S11;
277                              out->factors[iptr_main*9+1]=S21;
278                              out->factors[iptr_main*9+2]=S31;
279                              out->factors[iptr_main*9+3]=S12;
280                              out->factors[iptr_main*9+4]=S22;
281                              out->factors[iptr_main*9+5]=S32;
282                              out->factors[iptr_main*9+6]=S13;
283                              out->factors[iptr_main*9+7]=S23;
284                              out->factors[iptr_main*9+8]=S33;
285      
286                              /* a_ik=a_ii^{-1}*a_ik */
287                              for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
288                                 k=A->pattern->index[iptr_ik];
289                                 if (out->colorOf[k]>color) {
290                                    A11=out->factors[iptr_ik*9  ];
291                                    A21=out->factors[iptr_ik*9+1];
292                                    A31=out->factors[iptr_ik*9+2];
293                                    A12=out->factors[iptr_ik*9+3];
294                                    A22=out->factors[iptr_ik*9+4];
295                                    A32=out->factors[iptr_ik*9+5];
296                                    A13=out->factors[iptr_ik*9+6];
297                                    A23=out->factors[iptr_ik*9+7];
298                                    A33=out->factors[iptr_ik*9+8];
299                                    out->factors[iptr_ik*9  ]=S11*A11+S12*A21+S13*A31;                            
300                                    out->factors[iptr_ik*9+1]=S21*A11+S22*A21+S23*A31;
301                                    out->factors[iptr_ik*9+2]=S31*A11+S32*A21+S33*A31;
302                                    out->factors[iptr_ik*9+3]=S11*A12+S12*A22+S13*A32;                            
303                                    out->factors[iptr_ik*9+4]=S21*A12+S22*A22+S23*A32;
304                                    out->factors[iptr_ik*9+5]=S31*A12+S32*A22+S33*A32;
305                                    out->factors[iptr_ik*9+6]=S11*A13+S12*A23+S13*A33;                            
306                                    out->factors[iptr_ik*9+7]=S21*A13+S22*A23+S23*A33;
307                                    out->factors[iptr_ik*9+8]=S31*A13+S32*A23+S33*A33;
308                                 }                              
309                              }
310                         } else {                         } else {
311                              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.");
312                         }                         }
313                     }                      }
                 }  
               }  
            } /* 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 */  
                               }  
                             }  
                          }  
                      }  
314                   }                   }
315                }                } else {
316             }                   Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: block size greater than 3 is not supported.");
317          }                }      
318                  #pragma omp barrier
319           }
320           time_fac=Paso_timer()-time0;
321       }       }
322    }    }
   TMPMEMFREE(mis_marker);  
   TMPMEMFREE(counter);  
323    if (Paso_noError()) {    if (Paso_noError()) {
324        if (verbose) {        if (verbose) {
325           printf("ILU: %d unknowns eliminated. %d left.\n",out->n_F,n-out->n_F);           printf("ILU: %d color used \n",out->num_colors);
326           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);  
          }  
327       }       }
328       return out;       return out;
329    } else  {    } else  {
# Line 281  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 336  Paso_Solver_ILU* Paso_Solver_getILU(Paso
336    
337  /* apply ILU precondition b-> x                                /* apply ILU precondition b-> x                              
338    
339       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]  
340    
341   should be called within a parallel region                                                 should be called within a parallel region                                              
342   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 344  Paso_Solver_ILU* Paso_Solver_getILU(Paso
344  */  */
345    
346  void Paso_Solver_solveILU(Paso_Solver_ILU * ilu, double * x, double * b) {  void Paso_Solver_solveILU(Paso_Solver_ILU * ilu, double * x, double * b) {
347       dim_t i,k;       register dim_t i,k;
348         register index_t color,iptr_ik,iptr_main;
349         register double S1,S2,S3,R1,R2,R3;
350       dim_t n_block=ilu->n_block;       dim_t n_block=ilu->n_block;
351         dim_t n=ilu->n;
352            
353       if (ilu->n_C==0) {      
354          /* x=invA_FF*b  */       /* copy x into b*/
355          Paso_Solver_applyBlockDiagonalMatrix(n_block,ilu->n_F,ilu->inv_A_FF,ilu->A_FF_pivot,x,b);       #pragma omp parallel for private(i) schedule(static)
356       } else {       for (i=0;i<n*n_block;++i) x[i]=b[i];
357          /* b->[b_F,b_C]     */       /* forward substitution */
358          if (n_block==1) {       for (color=0;color<ilu->num_colors;++color) {
359             #pragma omp for private(i) schedule(static)             if (n_block==1) {
360             for (i=0;i<ilu->n_F;++i) ilu->b_F[i]=b[ilu->rows_in_F[i]];                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1,iptr_main)
361             #pragma omp for private(i) schedule(static)                for (i = 0; i < n; ++i) {
362             for (i=0;i<ilu->n_C;++i) ilu->b_C[i]=b[ilu->rows_in_C[i]];                     if (ilu->colorOf[i]==color) {
363          } else {                       /* x_i=x_i-a_ik*x_k */                    
364             #pragma omp for private(i,k) schedule(static)                       S1=x[i];
365             for (i=0;i<ilu->n_F;++i)                       for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
366                   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];                          
367             #pragma omp for private(i,k) schedule(static)                            if (ilu->colorOf[k]<color) {
368             for (i=0;i<ilu->n_C;++i)                               R1=x[k];                              
369                   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;
370          }                            }
371          /* x_F=invA_FF*b_F  */                       }
372          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];
373          /* b_C=b_C-A_CF*x_F */                       x[i]=ilu->factors[iptr_main]*S1;
374          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]];  
375                }                }
376             }             } else if (n_block==2) {
377          } else {                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2)
378             #pragma omp for private(i,k) schedule(static)                for (i = 0; i < n; ++i) {
379             for (i=0;i<ilu->n;++i) {                     if (ilu->colorOf[i]==color) {
380                   if (ilu->mask_C[i]>-1) {                       /* x_i=x_i-a_ik*x_k */
381                       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];
382                   } else {                       S2=x[2*i+1];
383                       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) {
384                              k=ilu->pattern->index[iptr_ik];                          
385                              if (ilu->colorOf[k]<color) {
386                                 R1=x[2*k];
387                                 R2=x[2*k+1];
388                                 S1-=ilu->factors[4*iptr_ik  ]*R1+ilu->factors[4*iptr_ik+2]*R2;
389                                 S2-=ilu->factors[4*iptr_ik+1]*R1+ilu->factors[4*iptr_ik+3]*R2;
390                              }
391                         }
392                         iptr_main=ilu->main_iptr[i];
393                         x[2*i  ]=ilu->factors[4*iptr_main  ]*S1+ilu->factors[4*iptr_main+2]*S2;
394                         x[2*i+1]=ilu->factors[4*iptr_main+1]*S1+ilu->factors[4*iptr_main+3]*S2;
395                       }
396    
397                  }
398               } else if (n_block==3) {
399                  #pragma omp parallel for schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3)
400                  for (i = 0; i < n; ++i) {
401                       if (ilu->colorOf[i]==color) {
402                         /* x_i=x_i-a_ik*x_k */
403                         S1=x[3*i];
404                         S2=x[3*i+1];
405                         S3=x[3*i+2];
406                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
407                              k=ilu->pattern->index[iptr_ik];                          
408                              if (ilu->colorOf[k]<color) {
409                                 R1=x[3*k];
410                                 R2=x[3*k+1];
411                                 R3=x[3*k+2];
412                                 S1-=ilu->factors[9*iptr_ik  ]*R1+ilu->factors[9*iptr_ik+3]*R2+ilu->factors[9*iptr_ik+6]*R3;
413                                 S2-=ilu->factors[9*iptr_ik+1]*R1+ilu->factors[9*iptr_ik+4]*R2+ilu->factors[9*iptr_ik+7]*R3;
414                                 S3-=ilu->factors[9*iptr_ik+2]*R1+ilu->factors[9*iptr_ik+5]*R2+ilu->factors[9*iptr_ik+8]*R3;
415                              }
416                         }
417                         iptr_main=ilu->main_iptr[i];
418                         x[3*i  ]=ilu->factors[9*iptr_main  ]*S1+ilu->factors[9*iptr_main+3]*S2+ilu->factors[9*iptr_main+6]*S3;
419                         x[3*i+1]=ilu->factors[9*iptr_main+1]*S1+ilu->factors[9*iptr_main+4]*S2+ilu->factors[9*iptr_main+7]*S3;
420                         x[3*i+2]=ilu->factors[9*iptr_main+2]*S1+ilu->factors[9*iptr_main+5]*S2+ilu->factors[9*iptr_main+8]*S3;
421                   }                   }
422                  }
423             }             }
424          }       }
425          /* all done */       /* backward substitution */
426         for (color=(ilu->num_colors)-1;color>-1;--color) {
427               if (n_block==1) {
428                  #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1)
429                  for (i = 0; i < n; ++i) {
430                       if (ilu->colorOf[i]==color) {
431                         /* x_i=x_i-a_ik*x_k */
432                         S1=x[i];
433                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
434                              k=ilu->pattern->index[iptr_ik];                          
435                              if (ilu->colorOf[k]>color) {
436                                 R1=x[k];
437                                 S1-=ilu->factors[iptr_ik]*R1;
438                              }
439                         }
440                         x[i]=S1;
441                       }
442                  }
443               } else if (n_block==2) {
444                  #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2)
445                  for (i = 0; i < n; ++i) {
446                       if (ilu->colorOf[i]==color) {
447                         /* x_i=x_i-a_ik*x_k */
448                         S1=x[2*i];
449                         S2=x[2*i+1];
450                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
451                              k=ilu->pattern->index[iptr_ik];                          
452                              if (ilu->colorOf[k]>color) {
453                                 R1=x[2*k];
454                                 R2=x[2*k+1];
455                                 S1-=ilu->factors[4*iptr_ik  ]*R1+ilu->factors[4*iptr_ik+2]*R2;
456                                 S2-=ilu->factors[4*iptr_ik+1]*R1+ilu->factors[4*iptr_ik+3]*R2;
457                              }
458                         }
459                         x[2*i]=S1;
460                         x[2*i+1]=S2;
461                       }
462                  }
463               } else if (n_block==3) {
464                  #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)
465                  for (i = 0; i < n; ++i) {
466                       if (ilu->colorOf[i]==color) {
467                         /* x_i=x_i-a_ik*x_k */
468                         S1=x[3*i  ];
469                         S2=x[3*i+1];
470                         S3=x[3*i+2];
471                         for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {
472                              k=ilu->pattern->index[iptr_ik];                          
473                              if (ilu->colorOf[k]>color) {
474                                 R1=x[3*k];
475                                 R2=x[3*k+1];
476                                 R3=x[3*k+2];
477                                 S1-=ilu->factors[9*iptr_ik  ]*R1+ilu->factors[9*iptr_ik+3]*R2+ilu->factors[9*iptr_ik+6]*R3;
478                                 S2-=ilu->factors[9*iptr_ik+1]*R1+ilu->factors[9*iptr_ik+4]*R2+ilu->factors[9*iptr_ik+7]*R3;
479                                 S3-=ilu->factors[9*iptr_ik+2]*R1+ilu->factors[9*iptr_ik+5]*R2+ilu->factors[9*iptr_ik+8]*R3;
480                              }
481                         }
482                         x[3*i]=S1;
483                         x[3*i+1]=S2;
484                         x[3*i+2]=S3;
485                       }
486                  }
487             }
488             #pragma omp barrier
489       }       }
490       return;       return;
491  }  }
   
492  /*  /*
493   * $Log$   * $Log$
494   * 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.1628

  ViewVC Help
Powered by ViewVC 1.1.26