/[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 631 by dhawcroft, Thu Mar 23 04:27:32 2006 UTC trunk/paso/src/Solver_ILU.c revision 1628 by phornby, Fri Jul 11 13:12:46 2008 UTC
# Line 1  Line 1 
 /* $Id$ */  
1    
2    /* $Id$ */
3    
4  /*  /*******************************************************
5  ********************************************************************************   *
6  *               Copyright 2006 by ACcESS MNRF                                *   *           Copyright 2003-2007 by ACceSS MNRF
7  *                                                                              *   *       Copyright 2007 by University of Queensland
8  *                 http://www.access.edu.au                                     *   *
9  *           Primary Business: Queensland, Australia                            *   *                http://esscc.uq.edu.au
10  *     Licensed under the Open Software License version 3.0             *   *        Primary Business: Queensland, Australia
11  *        http://www.opensource.org/licenses/osl-3.0.php                        *   *  Licensed under the Open Software License version 3.0
12  ********************************************************************************   *     http://www.opensource.org/licenses/osl-3.0.php
13  */   *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
# Line 36  void Paso_Solver_ILU_free(Paso_Solver_IL Line 37  void Paso_Solver_ILU_free(Paso_Solver_IL
37          MEMFREE(in->colorOf);          MEMFREE(in->colorOf);
38          MEMFREE(in->factors);          MEMFREE(in->factors);
39          MEMFREE(in->main_iptr);            MEMFREE(in->main_iptr);  
40          Paso_SystemMatrixPattern_dealloc(in->pattern);            Paso_Pattern_free(in->pattern);  
41          MEMFREE(in);          MEMFREE(in);
42       }       }
43  }  }
# Line 46  void Paso_Solver_ILU_free(Paso_Solver_IL Line 47  void Paso_Solver_ILU_free(Paso_Solver_IL
47  /*   constructs the incomplete block factorization of  /*   constructs the incomplete block factorization of
48    
49  */  */
50  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SystemMatrix * A,bool_t verbose) {  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SparseMatrix * A,bool_t verbose) {
51    dim_t n=A->num_rows;    dim_t n=A->numRows;
52    dim_t n_block=A->row_block_size;    dim_t n_block=A->row_block_size;
53    index_t num_colors=0;    index_t num_colors=0, *mis_marker=NULL;
54    register double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;    register double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;
   register double mainA11,mainA12,mainA13,mainA21,mainA22,mainA23,mainA31,mainA32,mainA33;  
55    register double S11,S12,S13,S21,S22,S23,S31,S32,S33;    register double S11,S12,S13,S21,S22,S23,S31,S32,S33;
56    register index_t i,iptr_main,iptr,iptr_ik,k,iptr_kj,j,iptr_ij,color,color2;    register index_t i,iptr_main,iptr,iptr_ik,k,iptr_kj,j,iptr_ij,color,color2;
57    double time0,time_color,time_fac;    double time0,time_color,time_fac;
58    /* allocations: */      /* allocations: */  
59    Paso_Solver_ILU* out=MEMALLOC(1,Paso_Solver_ILU);    Paso_Solver_ILU* out=MEMALLOC(1,Paso_Solver_ILU);
60    if (Paso_checkPtr(out)) return NULL;    if (Paso_checkPtr(out)) return NULL;
   index_t* mis_marker=TMPMEMALLOC(n,index_t);  
61    out->colorOf=MEMALLOC(n,index_t);    out->colorOf=MEMALLOC(n,index_t);
62    out->factors=MEMALLOC(A->len,double);    out->factors=MEMALLOC(A->len,double);
63    out->main_iptr=MEMALLOC(n,index_t);    out->main_iptr=MEMALLOC(n,index_t);
64    out->pattern=Paso_SystemMatrixPattern_reference(A->pattern);    out->pattern=Paso_Pattern_getReference(A->pattern);
65    out->n_block=n_block;    out->n_block=n_block;
66    out->n=n;    out->n=n;
67    time0=Paso_timer();    if ( !(Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->factors)) ) {
   if ( !(Paso_checkPtr(mis_marker) ||    
          Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->factors)) ) {  
     /* get coloring */  
     index_t num_colors=0;  
     #pragma omp parallel for private(i) schedule(static)  
     for (i = 0; i < n; ++i) out->colorOf[i]=-1;  
     while (Paso_Util_isAny(n,out->colorOf,-1) && Paso_noError()) {  
        #pragma omp parallel for private(i) schedule(static)  
        for (i = 0; i < n; ++i) mis_marker[i]=out->colorOf[i];  
        Paso_SystemMatrixPattern_mis(A->pattern,mis_marker);  
   
        #pragma omp parallel for private(i) schedule(static)  
        for (i = 0; i < n; ++i) if (mis_marker[i]) out->colorOf[i]=num_colors;  
        ++num_colors;  
     }  
     out->num_colors=num_colors;  
     time_color=Paso_timer()-time0;  
68      time0=Paso_timer();      time0=Paso_timer();
69      /* find main diagonal and copy matrix values */      Paso_Pattern_color(A->pattern,&out->num_colors,out->colorOf);
70      #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)      time_color=Paso_timer()-time0;
     for (i = 0; i < n; ++i) {  
         for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
             iptr_main=A->pattern->ptr[0]-1;  
             for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; iptr++) {  
                 if (A->pattern->index[iptr]==i) iptr_main=iptr;  
                 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];  
             }  
             out->main_iptr[i]=iptr_main;  
             if (iptr_main==A->pattern->ptr[0]-1)  
                Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: no main diagonal");  
         }  
     }  
     /* start factorization */  
71    
72      #pragma omp barrier      if (Paso_noError()) {
73      for (color=0;color<out->num_colors && Paso_noError();++color) {         time0=Paso_timer();
74             if (n_block==1) {         /* find main diagonal and copy matrix values */
75                #pragma omp parallel for schedule(static) private(i,color2,iptr_ik,k,iptr_kj,S11,j,iptr_ij,A11,iptr_main,D)         #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)
76                for (i = 0; i < n; ++i) {         for (i = 0; i < n; ++i) {
77                   if (out->colorOf[i]==color) {             for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
78                      for (color2=0;color2<color;++color2) {                 iptr_main=A->pattern->ptr[0]-1;
79                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                  for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; iptr++) {
80                            k=A->pattern->index[iptr_ik];                                               if (A->pattern->index[iptr]==i) iptr_main=iptr;
81                            if (out->colorOf[k]==color2) {                     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                               A11=out->factors[iptr_ik];                 }
83                               /* a_ij=a_ij-a_ik*a_kj */                 out->main_iptr[i]=iptr_main;
84                               for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {                 if (iptr_main==A->pattern->ptr[0]-1)
85                                  j=A->pattern->index[iptr_kj];                    Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: no main diagonal");
86                                  if (out->colorOf[j]>color2) {             }
87                                     S11=out->factors[iptr_kj];         }
88                                     for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {         /* start factorization */
89                                        if (j==A->pattern->index[iptr_ij]) {    
90                                           out->factors[iptr_ij]-=A11*S11;         #pragma omp barrier
91                                           break;         for (color=0;color<out->num_colors && Paso_noError();++color) {
92                  if (n_block==1) {
93                     #pragma omp parallel for schedule(static) private(i,color2,iptr_ik,k,iptr_kj,S11,j,iptr_ij,A11,iptr_main,D)
94                     for (i = 0; i < n; ++i) {
95                        if (out->colorOf[i]==color) {
96                           for (color2=0;color2<color;++color2) {
97                              for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
98                                 k=A->pattern->index[iptr_ik];                          
99                                 if (out->colorOf[k]==color2) {
100                                    A11=out->factors[iptr_ik];
101                                    /* a_ij=a_ij-a_ik*a_kj */
102                                    for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
103                                       j=A->pattern->index[iptr_kj];
104                                       if (out->colorOf[j]>color2) {
105                                          S11=out->factors[iptr_kj];
106                                          for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
107                                             if (j==A->pattern->index[iptr_ij]) {
108                                                out->factors[iptr_ij]-=A11*S11;
109                                                break;
110                                             }
111                                        }                                        }
112                                     }                                     }
113                                  }                                  }
114                               }                               }
115                            }                            }
116                         }                         }
117                      }                         iptr_main=out->main_iptr[i];
118                      iptr_main=out->main_iptr[i];                         D=out->factors[iptr_main];
119                      D=out->factors[iptr_main];                         if (ABS(D)>0.) {
120                      if (ABS(D)>0.) {                            D=1./D;
121                         D=1./D;                            out->factors[iptr_main]=D;
122                         out->factors[iptr_main]=D;                            /* a_ik=a_ii^{-1}*a_ik */
123                         /* a_ik=a_ii^{-1}*a_ik */                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
124                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                               k=A->pattern->index[iptr_ik];
125                            k=A->pattern->index[iptr_ik];                               if (out->colorOf[k]>color) {
126                            if (out->colorOf[k]>color) {                                  A11=out->factors[iptr_ik];
127                               A11=out->factors[iptr_ik];                                  out->factors[iptr_ik]=A11*D;
128                               out->factors[iptr_ik]=A11*D;                               }                              
129                            }                                                          }
130                           } else {
131                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
132                         }                         }
                     } else {  
                          Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");  
133                      }                      }
134                   }                   }
135                }                } else if (n_block==2) {
136             } else if (n_block==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                #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)                   for (i = 0; i < n; ++i) {
138                for (i = 0; i < n; ++i) {                      if (out->colorOf[i]==color) {
139                   if (out->colorOf[i]==color) {                         for (color2=0;color2<color;++color2) {
140                      for (color2=0;color2<color;++color2) {                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
141                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                               k=A->pattern->index[iptr_ik];                          
142                            k=A->pattern->index[iptr_ik];                                                         if (out->colorOf[k]==color2) {
143                            if (out->colorOf[k]==color2) {                                  A11=out->factors[iptr_ik*4  ];
144                               A11=out->factors[iptr_ik*4  ];                                  A21=out->factors[iptr_ik*4+1];
145                               A21=out->factors[iptr_ik*4+1];                                  A12=out->factors[iptr_ik*4+2];
146                               A12=out->factors[iptr_ik*4+2];                                  A22=out->factors[iptr_ik*4+3];
147                               A22=out->factors[iptr_ik*4+3];                                  /* a_ij=a_ij-a_ik*a_kj */
148                               /* a_ij=a_ij-a_ik*a_kj */                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
149                               for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {                                     j=A->pattern->index[iptr_kj];
150                                  j=A->pattern->index[iptr_kj];                                     if (out->colorOf[j]>color2) {
151                                  if (out->colorOf[j]>color2) {                                        S11=out->factors[iptr_kj*4];
152                                     S11=out->factors[iptr_kj*4];                                        S21=out->factors[iptr_kj*4+1];
153                                     S21=out->factors[iptr_kj*4+1];                                        S12=out->factors[iptr_kj*4+2];
154                                     S12=out->factors[iptr_kj*4+2];                                        S22=out->factors[iptr_kj*4+3];
155                                     S22=out->factors[iptr_kj*4+3];                                        for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
156                                     for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {                                           if (j==A->pattern->index[iptr_ij]) {
157                                        if (j==A->pattern->index[iptr_ij]) {                                              out->factors[4*iptr_ij  ]-=A11*S11+A12*S21;
158                                           out->factors[4*iptr_ij  ]-=A11*S11+A12*S21;                                              out->factors[4*iptr_ij+1]-=A21*S11+A22*S21;
159                                           out->factors[4*iptr_ij+1]-=A21*S11+A22*S21;                                              out->factors[4*iptr_ij+2]-=A11*S12+A12*S22;
160                                           out->factors[4*iptr_ij+2]-=A11*S12+A12*S22;                                              out->factors[4*iptr_ij+3]-=A21*S12+A22*S22;
161                                           out->factors[4*iptr_ij+3]-=A21*S12+A22*S22;                                              break;
162                                           break;                                           }
163                                        }                                        }
164                                     }                                     }
165                                  }                                  }
166                               }                               }
167                            }                            }
168                         }                         }
169                      }                         iptr_main=out->main_iptr[i];
170                      iptr_main=out->main_iptr[i];                         A11=out->factors[iptr_main*4];
171                      A11=out->factors[iptr_main*4];                         A21=out->factors[iptr_main*4+1];
172                      A21=out->factors[iptr_main*4+1];                         A12=out->factors[iptr_main*4+2];
173                      A12=out->factors[iptr_main*4+2];                         A22=out->factors[iptr_main*4+3];
174                      A22=out->factors[iptr_main*4+3];                         D = A11*A22-A12*A21;
175                      D = A11*A22-A12*A21;                         if (ABS(D)>0.) {
176                      if (ABS(D)>0.) {                            D=1./D;
177                         D=1./D;                            S11= A22*D;
178                         S11= A22*D;                            S21=-A21*D;
179                         S21=-A21*D;                            S12=-A12*D;
180                         S12=-A12*D;                            S22= A11*D;
181                         S22= A11*D;                            out->factors[iptr_main*4]  = S11;
182                         out->factors[iptr_main*4]  = S11;                            out->factors[iptr_main*4+1]= S21;
183                         out->factors[iptr_main*4+1]= S21;                            out->factors[iptr_main*4+2]= S12;
184                         out->factors[iptr_main*4+2]= S12;                            out->factors[iptr_main*4+3]= S22;
185                         out->factors[iptr_main*4+3]= S22;                            /* a_ik=a_ii^{-1}*a_ik */
186                         /* a_ik=a_ii^{-1}*a_ik */                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
187                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                               k=A->pattern->index[iptr_ik];
188                            k=A->pattern->index[iptr_ik];                               if (out->colorOf[k]>color) {
189                            if (out->colorOf[k]>color) {                                  A11=out->factors[iptr_ik*4  ];
190                               A11=out->factors[iptr_ik*4  ];                                  A21=out->factors[iptr_ik*4+1];
191                               A21=out->factors[iptr_ik*4+1];                                  A12=out->factors[iptr_ik*4+2];
192                               A12=out->factors[iptr_ik*4+2];                                  A22=out->factors[iptr_ik*4+3];
193                               A22=out->factors[iptr_ik*4+3];                                  out->factors[4*iptr_ik  ]=S11*A11+S12*A21;
194                               out->factors[4*iptr_ik  ]=S11*A11+S12*A21;                                  out->factors[4*iptr_ik+1]=S21*A11+S22*A21;
195                               out->factors[4*iptr_ik+1]=S21*A11+S22*A21;                                  out->factors[4*iptr_ik+2]=S11*A12+S12*A22;
196                               out->factors[4*iptr_ik+2]=S11*A12+S12*A22;                                  out->factors[4*iptr_ik+3]=S21*A12+S22*A22;
197                               out->factors[4*iptr_ik+3]=S21*A12+S22*A22;                               }                              
198                            }                                                          }
199                           } else {
200                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
201                         }                         }
                     } else {  
                          Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");  
202                      }                      }
203                   }                   }
204                }                } else if (n_block==3) {
205             } else if (n_block==3) {                   #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                #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)                   for (i = 0; i < n; ++i) {
207                for (i = 0; i < n; ++i) {                      if (out->colorOf[i]==color) {
208                   if (out->colorOf[i]==color) {                         for (color2=0;color2<color;++color2) {
209                      for (color2=0;color2<color;++color2) {                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
210                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                               k=A->pattern->index[iptr_ik];                          
211                            k=A->pattern->index[iptr_ik];                                                         if (out->colorOf[k]==color2) {
212                            if (out->colorOf[k]==color2) {                                  A11=out->factors[iptr_ik*9  ];
213                               A11=out->factors[iptr_ik*9  ];                                  A21=out->factors[iptr_ik*9+1];
214                               A21=out->factors[iptr_ik*9+1];                                  A31=out->factors[iptr_ik*9+2];
215                               A31=out->factors[iptr_ik*9+2];                                  A12=out->factors[iptr_ik*9+3];
216                               A12=out->factors[iptr_ik*9+3];                                  A22=out->factors[iptr_ik*9+4];
217                               A22=out->factors[iptr_ik*9+4];                                  A32=out->factors[iptr_ik*9+5];
218                               A32=out->factors[iptr_ik*9+5];                                  A13=out->factors[iptr_ik*9+6];
219                               A13=out->factors[iptr_ik*9+6];                                  A23=out->factors[iptr_ik*9+7];
220                               A23=out->factors[iptr_ik*9+7];                                  A33=out->factors[iptr_ik*9+8];
221                               A33=out->factors[iptr_ik*9+8];                                  /* a_ij=a_ij-a_ik*a_kj */
222                               /* a_ij=a_ij-a_ik*a_kj */                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
223                               for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {                                     j=A->pattern->index[iptr_kj];
224                                  j=A->pattern->index[iptr_kj];                                     if (out->colorOf[j]>color2) {
225                                  if (out->colorOf[j]>color2) {                                        S11=out->factors[iptr_kj*9  ];
226                                     S11=out->factors[iptr_kj*9  ];                                        S21=out->factors[iptr_kj*9+1];
227                                     S21=out->factors[iptr_kj*9+1];                                        S31=out->factors[iptr_kj*9+2];
228                                     S31=out->factors[iptr_kj*9+2];                                        S12=out->factors[iptr_kj*9+3];
229                                     S12=out->factors[iptr_kj*9+3];                                        S22=out->factors[iptr_kj*9+4];
230                                     S22=out->factors[iptr_kj*9+4];                                        S32=out->factors[iptr_kj*9+5];
231                                     S32=out->factors[iptr_kj*9+5];                                        S13=out->factors[iptr_kj*9+6];
232                                     S13=out->factors[iptr_kj*9+6];                                        S23=out->factors[iptr_kj*9+7];
233                                     S23=out->factors[iptr_kj*9+7];                                        S33=out->factors[iptr_kj*9+8];                                
234                                     S33=out->factors[iptr_kj*9+8];                                                                        for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
235                                     for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {                                           if (j==A->pattern->index[iptr_ij]) {
236                                        if (j==A->pattern->index[iptr_ij]) {                                              out->factors[iptr_ij*9  ]-=A11*S11+A12*S21+A13*S31;                            
237                                           out->factors[iptr_ij*9  ]-=A11*S11+A12*S21+A13*S31;                                                                          out->factors[iptr_ij*9+1]-=A21*S11+A22*S21+A23*S31;
238                                           out->factors[iptr_ij*9+1]-=A21*S11+A22*S21+A23*S31;                                              out->factors[iptr_ij*9+2]-=A31*S11+A32*S21+A33*S31;
239                                           out->factors[iptr_ij*9+2]-=A31*S11+A32*S21+A33*S31;                                              out->factors[iptr_ij*9+3]-=A11*S12+A12*S22+A13*S32;                            
240                                           out->factors[iptr_ij*9+3]-=A11*S12+A12*S22+A13*S32;                                                                          out->factors[iptr_ij*9+4]-=A21*S12+A22*S22+A23*S32;
241                                           out->factors[iptr_ij*9+4]-=A21*S12+A22*S22+A23*S32;                                              out->factors[iptr_ij*9+5]-=A31*S12+A32*S22+A33*S32;
242                                           out->factors[iptr_ij*9+5]-=A31*S12+A32*S22+A33*S32;                                              out->factors[iptr_ij*9+6]-=A11*S13+A12*S23+A13*S33;                            
243                                           out->factors[iptr_ij*9+6]-=A11*S13+A12*S23+A13*S33;                                                                          out->factors[iptr_ij*9+7]-=A21*S13+A22*S23+A23*S33;
244                                           out->factors[iptr_ij*9+7]-=A21*S13+A22*S23+A23*S33;                                              out->factors[iptr_ij*9+8]-=A31*S13+A32*S23+A33*S33;
245                                           out->factors[iptr_ij*9+8]-=A31*S13+A32*S23+A33*S33;                                              break;
246                                           break;                                           }
247                                        }                                        }
248                                     }                                     }
249                                  }                                  }
250                               }                               }
251                            }                            }
252                         }                         }
253                      }                         iptr_main=out->main_iptr[i];
254                      iptr_main=out->main_iptr[i];                         A11=out->factors[iptr_main*9  ];
255                      A11=out->factors[iptr_main*9  ];                         A21=out->factors[iptr_main*9+1];
256                      A21=out->factors[iptr_main*9+1];                         A31=out->factors[iptr_main*9+2];
257                      A31=out->factors[iptr_main*9+2];                         A12=out->factors[iptr_main*9+3];
258                      A12=out->factors[iptr_main*9+3];                         A22=out->factors[iptr_main*9+4];
259                      A22=out->factors[iptr_main*9+4];                         A32=out->factors[iptr_main*9+5];
260                      A32=out->factors[iptr_main*9+5];                         A13=out->factors[iptr_main*9+6];
261                      A13=out->factors[iptr_main*9+6];                         A23=out->factors[iptr_main*9+7];
262                      A23=out->factors[iptr_main*9+7];                         A33=out->factors[iptr_main*9+8];
263                      A33=out->factors[iptr_main*9+8];                         D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
264                      D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);                         if (ABS(D)>0.) {
265                      if (ABS(D)>0.) {                            D=1./D;
266                         D=1./D;                            S11=(A22*A33-A23*A32)*D;
267                         S11=(A22*A33-A23*A32)*D;                            S21=(A31*A23-A21*A33)*D;
268                         S21=(A31*A23-A21*A33)*D;                            S31=(A21*A32-A31*A22)*D;
269                         S31=(A21*A32-A31*A22)*D;                            S12=(A13*A32-A12*A33)*D;
270                         S12=(A13*A32-A12*A33)*D;                            S22=(A11*A33-A31*A13)*D;
271                         S22=(A11*A33-A31*A13)*D;                            S32=(A12*A31-A11*A32)*D;
272                         S32=(A12*A31-A11*A32)*D;                            S13=(A12*A23-A13*A22)*D;
273                         S13=(A12*A23-A13*A22)*D;                            S23=(A13*A21-A11*A23)*D;
274                         S23=(A13*A21-A11*A23)*D;                            S33=(A11*A22-A12*A21)*D;
275                         S33=(A11*A22-A12*A21)*D;    
276                              out->factors[iptr_main*9  ]=S11;
277                         out->factors[iptr_main*9  ]=S11;                            out->factors[iptr_main*9+1]=S21;
278                         out->factors[iptr_main*9+1]=S21;                            out->factors[iptr_main*9+2]=S31;
279                         out->factors[iptr_main*9+2]=S31;                            out->factors[iptr_main*9+3]=S12;
280                         out->factors[iptr_main*9+3]=S12;                            out->factors[iptr_main*9+4]=S22;
281                         out->factors[iptr_main*9+4]=S22;                            out->factors[iptr_main*9+5]=S32;
282                         out->factors[iptr_main*9+5]=S32;                            out->factors[iptr_main*9+6]=S13;
283                         out->factors[iptr_main*9+6]=S13;                            out->factors[iptr_main*9+7]=S23;
284                         out->factors[iptr_main*9+7]=S23;                            out->factors[iptr_main*9+8]=S33;
285                         out->factors[iptr_main*9+8]=S33;    
286                              /* a_ik=a_ii^{-1}*a_ik */
287                         /* a_ik=a_ii^{-1}*a_ik */                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
288                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                               k=A->pattern->index[iptr_ik];
289                            k=A->pattern->index[iptr_ik];                               if (out->colorOf[k]>color) {
290                            if (out->colorOf[k]>color) {                                  A11=out->factors[iptr_ik*9  ];
291                               A11=out->factors[iptr_ik*9  ];                                  A21=out->factors[iptr_ik*9+1];
292                               A21=out->factors[iptr_ik*9+1];                                  A31=out->factors[iptr_ik*9+2];
293                               A31=out->factors[iptr_ik*9+2];                                  A12=out->factors[iptr_ik*9+3];
294                               A12=out->factors[iptr_ik*9+3];                                  A22=out->factors[iptr_ik*9+4];
295                               A22=out->factors[iptr_ik*9+4];                                  A32=out->factors[iptr_ik*9+5];
296                               A32=out->factors[iptr_ik*9+5];                                  A13=out->factors[iptr_ik*9+6];
297                               A13=out->factors[iptr_ik*9+6];                                  A23=out->factors[iptr_ik*9+7];
298                               A23=out->factors[iptr_ik*9+7];                                  A33=out->factors[iptr_ik*9+8];
299                               A33=out->factors[iptr_ik*9+8];                                  out->factors[iptr_ik*9  ]=S11*A11+S12*A21+S13*A31;                            
300                               out->factors[iptr_ik*9  ]=S11*A11+S12*A21+S13*A31;                                                              out->factors[iptr_ik*9+1]=S21*A11+S22*A21+S23*A31;
301                               out->factors[iptr_ik*9+1]=S21*A11+S22*A21+S23*A31;                                  out->factors[iptr_ik*9+2]=S31*A11+S32*A21+S33*A31;
302                               out->factors[iptr_ik*9+2]=S31*A11+S32*A21+S33*A31;                                  out->factors[iptr_ik*9+3]=S11*A12+S12*A22+S13*A32;                            
303                               out->factors[iptr_ik*9+3]=S11*A12+S12*A22+S13*A32;                                                              out->factors[iptr_ik*9+4]=S21*A12+S22*A22+S23*A32;
304                               out->factors[iptr_ik*9+4]=S21*A12+S22*A22+S23*A32;                                  out->factors[iptr_ik*9+5]=S31*A12+S32*A22+S33*A32;
305                               out->factors[iptr_ik*9+5]=S31*A12+S32*A22+S33*A32;                                  out->factors[iptr_ik*9+6]=S11*A13+S12*A23+S13*A33;                            
306                               out->factors[iptr_ik*9+6]=S11*A13+S12*A23+S13*A33;                                                              out->factors[iptr_ik*9+7]=S21*A13+S22*A23+S23*A33;
307                               out->factors[iptr_ik*9+7]=S21*A13+S22*A23+S23*A33;                                  out->factors[iptr_ik*9+8]=S31*A13+S32*A23+S33*A33;
308                               out->factors[iptr_ik*9+8]=S31*A13+S32*A23+S33*A33;                               }                              
309                            }                                                          }
310                           } else {
311                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
312                         }                         }
                     } else {  
                          Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");  
313                      }                      }
314                   }                   }
315                }                } else {
316             } else {                   Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: block size greater than 3 is not supported.");
317                Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: block size greater than 3 is not supported.");                }      
318             }                      #pragma omp barrier
319             #pragma omp barrier         }
320          }         time_fac=Paso_timer()-time0;
321          time_fac=Paso_timer()-time0;       }
322    }    }
   
   TMPMEMFREE(mis_marker);  
323    if (Paso_noError()) {    if (Paso_noError()) {
324        if (verbose) {        if (verbose) {
325           printf("ILU: %d color used \n",out->num_colors);           printf("ILU: %d color used \n",out->num_colors);
# Line 366  void Paso_Solver_solveILU(Paso_Solver_IL Line 352  void Paso_Solver_solveILU(Paso_Solver_IL
352            
353            
354       /* copy x into b*/       /* copy x into b*/
355       #pragma omp for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
356       for (i=0;i<n*n_block;++i) x[i]=b[i];       for (i=0;i<n*n_block;++i) x[i]=b[i];
357       /* forward substitution */       /* forward substitution */
358       for (color=0;color<ilu->num_colors;++color) {       for (color=0;color<ilu->num_colors;++color) {
359             if (n_block==1) {             if (n_block==1) {
360                #pragma omp for schedule(static) private(i,iptr_ik,k,S1,R1,iptr_main)                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1,iptr_main)
361                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
362                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
363                       /* x_i=x_i-a_ik*x_k */                                           /* x_i=x_i-a_ik*x_k */                    
# Line 388  void Paso_Solver_solveILU(Paso_Solver_IL Line 374  void Paso_Solver_solveILU(Paso_Solver_IL
374                     }                     }
375                }                }
376             } else if (n_block==2) {             } else if (n_block==2) {
377                #pragma omp for schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2)                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2)
378                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
379                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
380                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 410  void Paso_Solver_solveILU(Paso_Solver_IL Line 396  void Paso_Solver_solveILU(Paso_Solver_IL
396    
397                }                }
398             } else if (n_block==3) {             } else if (n_block==3) {
399                #pragma omp for schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3)                #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) {                for (i = 0; i < n; ++i) {
401                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
402                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 435  void Paso_Solver_solveILU(Paso_Solver_IL Line 421  void Paso_Solver_solveILU(Paso_Solver_IL
421                   }                   }
422                }                }
423             }             }
            #pragma omp barrier  
424       }       }
425       /* backward substitution */       /* backward substitution */
426       for (color=(ilu->num_colors)-1;color>-1;--color) {       for (color=(ilu->num_colors)-1;color>-1;--color) {
427             if (n_block==1) {             if (n_block==1) {
428                #pragma omp for schedule(static) private(i,iptr_ik,k,S1,R1)                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1)
429                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
430                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
431                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 456  void Paso_Solver_solveILU(Paso_Solver_IL Line 441  void Paso_Solver_solveILU(Paso_Solver_IL
441                     }                     }
442                }                }
443             } else if (n_block==2) {             } else if (n_block==2) {
444                #pragma omp for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2)                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2)
445                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
446                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
447                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 476  void Paso_Solver_solveILU(Paso_Solver_IL Line 461  void Paso_Solver_solveILU(Paso_Solver_IL
461                     }                     }
462                }                }
463             } else if (n_block==3) {             } else if (n_block==3) {
464                #pragma omp for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)
465                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
466                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
467                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */

Legend:
Removed from v.631  
changed lines
  Added in v.1628

  ViewVC Help
Powered by ViewVC 1.1.26