/[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 483 by jgs, Thu Feb 2 02:10:15 2006 UTC trunk/paso/src/ILU.c revision 3120 by gross, Mon Aug 30 10:48:00 2010 UTC
# Line 1  Line 1 
1  /* $Id$ */  
2    /*******************************************************
3    *
4    * Copyright (c) 2003-2010 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 7  Line 19 
19  /**************************************************************/  /**************************************************************/
20    
21  /* Copyrights by ACcESS Australia 2003,2004,2005              */  /* Copyrights by ACcESS Australia 2003,2004,2005              */
22  /* Author: gross@access.edu.au                                */  /* Author: Lutz Gross, l.gross@uq.edu.au                                */
23    
24  /**************************************************************/  /**************************************************************/
25    
26  #include "Paso.h"  #include "Paso.h"
27  #include "Solver.h"  #include "Preconditioner.h"
28  #include "PasoUtil.h"  #include "PasoUtil.h"
29    
30  /**************************************************************/  /**************************************************************/
# 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) {
         MEMFREE(in->colorOf);  
36          MEMFREE(in->factors);          MEMFREE(in->factors);
         MEMFREE(in->main_iptr);    
         Paso_SystemMatrixPattern_dealloc(in->pattern);    
37          MEMFREE(in);          MEMFREE(in);
38       }       }
39  }  }
# Line 34  void Paso_Solver_ILU_free(Paso_Solver_IL Line 43  void Paso_Solver_ILU_free(Paso_Solver_IL
43  /*   constructs the incomplete block factorization of  /*   constructs the incomplete block factorization of
44    
45  */  */
46  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SystemMatrix * A,bool_t verbose) {  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SparseMatrix * A,bool_t verbose) {
47    dim_t n=A->num_rows;    const dim_t n=A->numRows;
48    dim_t n_block=A->row_block_size;    const dim_t n_block=A->row_block_size;
49    index_t num_colors=0;    const index_t* colorOf = Paso_Pattern_borrowColoringPointer(A->pattern);
50      const dim_t num_colors = Paso_Pattern_getNumColors(A->pattern);
51      const index_t *ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A);
52    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;  
53    register double S11,S12,S13,S21,S22,S23,S31,S32,S33;    register double S11,S12,S13,S21,S22,S23,S31,S32,S33;
54    register index_t i,iptr_main,iptr,iptr_ik,k,iptr_kj,j,iptr_ij,color,color2;    register index_t i,iptr_main,iptr_ik,k,iptr_kj,j,iptr_ij,color,color2, iptr;
55    double time0,time_color,time_fac;    double time0=0,time_fac=0;
56    /* allocations: */      /* allocations: */  
57    Paso_Solver_ILU* out=MEMALLOC(1,Paso_Solver_ILU);    Paso_Solver_ILU* out=MEMALLOC(1,Paso_Solver_ILU);
58    if (Paso_checkPtr(out)) return NULL;    if (Paso_checkPtr(out)) return NULL;
   index_t* mis_marker=TMPMEMALLOC(n,index_t);  
   out->colorOf=MEMALLOC(n,index_t);  
59    out->factors=MEMALLOC(A->len,double);    out->factors=MEMALLOC(A->len,double);
60    out->main_iptr=MEMALLOC(n,index_t);    
61    out->pattern=Paso_SystemMatrixPattern_reference(A->pattern);    if ( ! Paso_checkPtr(out->factors)  ) {
   out->n_block=n_block;  
   out->n=n;  
   time0=Paso_timer();  
   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;  
     time0=Paso_timer();  
     /* find main diagonal and copy matrix values */  
     #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)  
     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 */  
62    
63      #pragma omp barrier         time0=Paso_timer();
64      for (color=0;color<out->num_colors && Paso_noError();++color) {  
65             if (n_block==1) {         #pragma omp parallel for schedule(static) private(i,iptr,k)
66                #pragma omp parallel for schedule(static) private(i,color2,iptr_ik,k,iptr_kj,S11,j,iptr_ij,A11,iptr_main,D)         for (i = 0; i < n; ++i) {
67                for (i = 0; i < n; ++i) {                 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; iptr++) {
68                   if (out->colorOf[i]==color) {                       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];
69                      for (color2=0;color2<color;++color2) {                 }
70                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {         }    
71                            k=A->pattern->index[iptr_ik];                                   /* start factorization */
72                            if (out->colorOf[k]==color2) {         for (color=0;color<num_colors && Paso_noError();++color) {
73                               A11=out->factors[iptr_ik];                if (n_block==1) {
74                               /* a_ij=a_ij-a_ik*a_kj */                   #pragma omp parallel for schedule(static) private(i,color2,iptr_ik,k,iptr_kj,S11,j,iptr_ij,A11,iptr_main,D)
75                               for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {                   for (i = 0; i < n; ++i) {
76                                  j=A->pattern->index[iptr_kj];                      if (colorOf[i]==color) {
77                                  if (out->colorOf[j]>color2) {                         for (color2=0;color2<color;++color2) {
78                                     S11=out->factors[iptr_kj];                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
79                                     for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {                               k=A->pattern->index[iptr_ik];                          
80                                        if (j==A->pattern->index[iptr_ij]) {                               if (colorOf[k]==color2) {
81                                           out->factors[iptr_ij]-=A11*S11;                                  A11=out->factors[iptr_ik];
82                                           break;                                  /* a_ij=a_ij-a_ik*a_kj */
83                                    for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
84                                       j=A->pattern->index[iptr_kj];
85                                       if (colorOf[j]>color2) {
86                                          S11=out->factors[iptr_kj];
87                                          for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
88                                             if (j==A->pattern->index[iptr_ij]) {
89                                                out->factors[iptr_ij]-=A11*S11;
90                                                break;
91                                             }
92                                        }                                        }
93                                     }                                     }
94                                  }                                  }
95                               }                               }
96                            }                            }
97                         }                         }
98                      }                         iptr_main=ptr_main[i];
99                      iptr_main=out->main_iptr[i];                         D=out->factors[iptr_main];
100                      D=out->factors[iptr_main];                         if (ABS(D)>0.) {
101                      if (ABS(D)>0.) {                            D=1./D;
102                         D=1./D;                            out->factors[iptr_main]=D;
103                         out->factors[iptr_main]=D;                            /* a_ik=a_ii^{-1}*a_ik */
104                         /* a_ik=a_ii^{-1}*a_ik */                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
105                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                               k=A->pattern->index[iptr_ik];
106                            k=A->pattern->index[iptr_ik];                               if (colorOf[k]>color) {
107                            if (out->colorOf[k]>color) {                                  A11=out->factors[iptr_ik];
108                               A11=out->factors[iptr_ik];                                  out->factors[iptr_ik]=A11*D;
109                               out->factors[iptr_ik]=A11*D;                               }                              
110                            }                                                          }
111                           } else {
112                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
113                         }                         }
                     } else {  
                          Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");  
114                      }                      }
115                   }                   }
116                }                } else if (n_block==2) {
117             } 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)
118                #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) {
119                for (i = 0; i < n; ++i) {                      if (colorOf[i]==color) {
120                   if (out->colorOf[i]==color) {                         for (color2=0;color2<color;++color2) {
121                      for (color2=0;color2<color;++color2) {                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
122                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                               k=A->pattern->index[iptr_ik];                          
123                            k=A->pattern->index[iptr_ik];                                                         if (colorOf[k]==color2) {
124                            if (out->colorOf[k]==color2) {                                  A11=out->factors[iptr_ik*4  ];
125                               A11=out->factors[iptr_ik*4  ];                                  A21=out->factors[iptr_ik*4+1];
126                               A21=out->factors[iptr_ik*4+1];                                  A12=out->factors[iptr_ik*4+2];
127                               A12=out->factors[iptr_ik*4+2];                                  A22=out->factors[iptr_ik*4+3];
128                               A22=out->factors[iptr_ik*4+3];                                  /* a_ij=a_ij-a_ik*a_kj */
129                               /* a_ij=a_ij-a_ik*a_kj */                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
130                               for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {                                     j=A->pattern->index[iptr_kj];
131                                  j=A->pattern->index[iptr_kj];                                     if (colorOf[j]>color2) {
132                                  if (out->colorOf[j]>color2) {                                        S11=out->factors[iptr_kj*4];
133                                     S11=out->factors[iptr_kj*4];                                        S21=out->factors[iptr_kj*4+1];
134                                     S21=out->factors[iptr_kj*4+1];                                        S12=out->factors[iptr_kj*4+2];
135                                     S12=out->factors[iptr_kj*4+2];                                        S22=out->factors[iptr_kj*4+3];
136                                     S22=out->factors[iptr_kj*4+3];                                        for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
137                                     for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {                                           if (j==A->pattern->index[iptr_ij]) {
138                                        if (j==A->pattern->index[iptr_ij]) {                                              out->factors[4*iptr_ij  ]-=A11*S11+A12*S21;
139                                           out->factors[4*iptr_ij  ]-=A11*S11+A12*S21;                                              out->factors[4*iptr_ij+1]-=A21*S11+A22*S21;
140                                           out->factors[4*iptr_ij+1]-=A21*S11+A22*S21;                                              out->factors[4*iptr_ij+2]-=A11*S12+A12*S22;
141                                           out->factors[4*iptr_ij+2]-=A11*S12+A12*S22;                                              out->factors[4*iptr_ij+3]-=A21*S12+A22*S22;
142                                           out->factors[4*iptr_ij+3]-=A21*S12+A22*S22;                                              break;
143                                           break;                                           }
144                                        }                                        }
145                                     }                                     }
146                                  }                                  }
147                               }                               }
148                            }                            }
149                         }                         }
150                      }                         iptr_main=ptr_main[i];
151                      iptr_main=out->main_iptr[i];                         A11=out->factors[iptr_main*4];
152                      A11=out->factors[iptr_main*4];                         A21=out->factors[iptr_main*4+1];
153                      A21=out->factors[iptr_main*4+1];                         A12=out->factors[iptr_main*4+2];
154                      A12=out->factors[iptr_main*4+2];                         A22=out->factors[iptr_main*4+3];
155                      A22=out->factors[iptr_main*4+3];                         D = A11*A22-A12*A21;
156                      D = A11*A22-A12*A21;                         if (ABS(D)>0.) {
157                      if (ABS(D)>0.) {                            D=1./D;
158                         D=1./D;                            S11= A22*D;
159                         S11= A22*D;                            S21=-A21*D;
160                         S21=-A21*D;                            S12=-A12*D;
161                         S12=-A12*D;                            S22= A11*D;
162                         S22= A11*D;                            out->factors[iptr_main*4]  = S11;
163                         out->factors[iptr_main*4]  = S11;                            out->factors[iptr_main*4+1]= S21;
164                         out->factors[iptr_main*4+1]= S21;                            out->factors[iptr_main*4+2]= S12;
165                         out->factors[iptr_main*4+2]= S12;                            out->factors[iptr_main*4+3]= S22;
166                         out->factors[iptr_main*4+3]= S22;                            /* a_ik=a_ii^{-1}*a_ik */
167                         /* a_ik=a_ii^{-1}*a_ik */                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
168                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                               k=A->pattern->index[iptr_ik];
169                            k=A->pattern->index[iptr_ik];                               if (colorOf[k]>color) {
170                            if (out->colorOf[k]>color) {                                  A11=out->factors[iptr_ik*4  ];
171                               A11=out->factors[iptr_ik*4  ];                                  A21=out->factors[iptr_ik*4+1];
172                               A21=out->factors[iptr_ik*4+1];                                  A12=out->factors[iptr_ik*4+2];
173                               A12=out->factors[iptr_ik*4+2];                                  A22=out->factors[iptr_ik*4+3];
174                               A22=out->factors[iptr_ik*4+3];                                  out->factors[4*iptr_ik  ]=S11*A11+S12*A21;
175                               out->factors[4*iptr_ik  ]=S11*A11+S12*A21;                                  out->factors[4*iptr_ik+1]=S21*A11+S22*A21;
176                               out->factors[4*iptr_ik+1]=S21*A11+S22*A21;                                  out->factors[4*iptr_ik+2]=S11*A12+S12*A22;
177                               out->factors[4*iptr_ik+2]=S11*A12+S12*A22;                                  out->factors[4*iptr_ik+3]=S21*A12+S22*A22;
178                               out->factors[4*iptr_ik+3]=S21*A12+S22*A22;                               }                              
179                            }                                                          }
180                           } else {
181                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
182                         }                         }
                     } else {  
                          Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");  
183                      }                      }
184                   }                   }
185                }                } else if (n_block==3) {
186             } 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)
187                #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) {
188                for (i = 0; i < n; ++i) {                      if (colorOf[i]==color) {
189                   if (out->colorOf[i]==color) {                         for (color2=0;color2<color;++color2) {
190                      for (color2=0;color2<color;++color2) {                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
191                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                               k=A->pattern->index[iptr_ik];                          
192                            k=A->pattern->index[iptr_ik];                                                         if (colorOf[k]==color2) {
193                            if (out->colorOf[k]==color2) {                                  A11=out->factors[iptr_ik*9  ];
194                               A11=out->factors[iptr_ik*9  ];                                  A21=out->factors[iptr_ik*9+1];
195                               A21=out->factors[iptr_ik*9+1];                                  A31=out->factors[iptr_ik*9+2];
196                               A31=out->factors[iptr_ik*9+2];                                  A12=out->factors[iptr_ik*9+3];
197                               A12=out->factors[iptr_ik*9+3];                                  A22=out->factors[iptr_ik*9+4];
198                               A22=out->factors[iptr_ik*9+4];                                  A32=out->factors[iptr_ik*9+5];
199                               A32=out->factors[iptr_ik*9+5];                                  A13=out->factors[iptr_ik*9+6];
200                               A13=out->factors[iptr_ik*9+6];                                  A23=out->factors[iptr_ik*9+7];
201                               A23=out->factors[iptr_ik*9+7];                                  A33=out->factors[iptr_ik*9+8];
202                               A33=out->factors[iptr_ik*9+8];                                  /* a_ij=a_ij-a_ik*a_kj */
203                               /* a_ij=a_ij-a_ik*a_kj */                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
204                               for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {                                     j=A->pattern->index[iptr_kj];
205                                  j=A->pattern->index[iptr_kj];                                     if (colorOf[j]>color2) {
206                                  if (out->colorOf[j]>color2) {                                        S11=out->factors[iptr_kj*9  ];
207                                     S11=out->factors[iptr_kj*9  ];                                        S21=out->factors[iptr_kj*9+1];
208                                     S21=out->factors[iptr_kj*9+1];                                        S31=out->factors[iptr_kj*9+2];
209                                     S31=out->factors[iptr_kj*9+2];                                        S12=out->factors[iptr_kj*9+3];
210                                     S12=out->factors[iptr_kj*9+3];                                        S22=out->factors[iptr_kj*9+4];
211                                     S22=out->factors[iptr_kj*9+4];                                        S32=out->factors[iptr_kj*9+5];
212                                     S32=out->factors[iptr_kj*9+5];                                        S13=out->factors[iptr_kj*9+6];
213                                     S13=out->factors[iptr_kj*9+6];                                        S23=out->factors[iptr_kj*9+7];
214                                     S23=out->factors[iptr_kj*9+7];                                        S33=out->factors[iptr_kj*9+8];                                
215                                     S33=out->factors[iptr_kj*9+8];                                                                        for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
216                                     for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {                                           if (j==A->pattern->index[iptr_ij]) {
217                                        if (j==A->pattern->index[iptr_ij]) {                                              out->factors[iptr_ij*9  ]-=A11*S11+A12*S21+A13*S31;                            
218                                           out->factors[iptr_ij*9  ]-=A11*S11+A12*S21+A13*S31;                                                                          out->factors[iptr_ij*9+1]-=A21*S11+A22*S21+A23*S31;
219                                           out->factors[iptr_ij*9+1]-=A21*S11+A22*S21+A23*S31;                                              out->factors[iptr_ij*9+2]-=A31*S11+A32*S21+A33*S31;
220                                           out->factors[iptr_ij*9+2]-=A31*S11+A32*S21+A33*S31;                                              out->factors[iptr_ij*9+3]-=A11*S12+A12*S22+A13*S32;                            
221                                           out->factors[iptr_ij*9+3]-=A11*S12+A12*S22+A13*S32;                                                                          out->factors[iptr_ij*9+4]-=A21*S12+A22*S22+A23*S32;
222                                           out->factors[iptr_ij*9+4]-=A21*S12+A22*S22+A23*S32;                                              out->factors[iptr_ij*9+5]-=A31*S12+A32*S22+A33*S32;
223                                           out->factors[iptr_ij*9+5]-=A31*S12+A32*S22+A33*S32;                                              out->factors[iptr_ij*9+6]-=A11*S13+A12*S23+A13*S33;                            
224                                           out->factors[iptr_ij*9+6]-=A11*S13+A12*S23+A13*S33;                                                                          out->factors[iptr_ij*9+7]-=A21*S13+A22*S23+A23*S33;
225                                           out->factors[iptr_ij*9+7]-=A21*S13+A22*S23+A23*S33;                                              out->factors[iptr_ij*9+8]-=A31*S13+A32*S23+A33*S33;
226                                           out->factors[iptr_ij*9+8]-=A31*S13+A32*S23+A33*S33;                                              break;
227                                           break;                                           }
228                                        }                                        }
229                                     }                                     }
230                                  }                                  }
231                               }                               }
232                            }                            }
233                         }                         }
234                      }                         iptr_main=ptr_main[i];
235                      iptr_main=out->main_iptr[i];                         A11=out->factors[iptr_main*9  ];
236                      A11=out->factors[iptr_main*9  ];                         A21=out->factors[iptr_main*9+1];
237                      A21=out->factors[iptr_main*9+1];                         A31=out->factors[iptr_main*9+2];
238                      A31=out->factors[iptr_main*9+2];                         A12=out->factors[iptr_main*9+3];
239                      A12=out->factors[iptr_main*9+3];                         A22=out->factors[iptr_main*9+4];
240                      A22=out->factors[iptr_main*9+4];                         A32=out->factors[iptr_main*9+5];
241                      A32=out->factors[iptr_main*9+5];                         A13=out->factors[iptr_main*9+6];
242                      A13=out->factors[iptr_main*9+6];                         A23=out->factors[iptr_main*9+7];
243                      A23=out->factors[iptr_main*9+7];                         A33=out->factors[iptr_main*9+8];
244                      A33=out->factors[iptr_main*9+8];                         D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
245                      D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);                         if (ABS(D)>0.) {
246                      if (ABS(D)>0.) {                            D=1./D;
247                         D=1./D;                            S11=(A22*A33-A23*A32)*D;
248                         S11=(A22*A33-A23*A32)*D;                            S21=(A31*A23-A21*A33)*D;
249                         S21=(A31*A23-A21*A33)*D;                            S31=(A21*A32-A31*A22)*D;
250                         S31=(A21*A32-A31*A22)*D;                            S12=(A13*A32-A12*A33)*D;
251                         S12=(A13*A32-A12*A33)*D;                            S22=(A11*A33-A31*A13)*D;
252                         S22=(A11*A33-A31*A13)*D;                            S32=(A12*A31-A11*A32)*D;
253                         S32=(A12*A31-A11*A32)*D;                            S13=(A12*A23-A13*A22)*D;
254                         S13=(A12*A23-A13*A22)*D;                            S23=(A13*A21-A11*A23)*D;
255                         S23=(A13*A21-A11*A23)*D;                            S33=(A11*A22-A12*A21)*D;
256                         S33=(A11*A22-A12*A21)*D;    
257                              out->factors[iptr_main*9  ]=S11;
258                         out->factors[iptr_main*9  ]=S11;                            out->factors[iptr_main*9+1]=S21;
259                         out->factors[iptr_main*9+1]=S21;                            out->factors[iptr_main*9+2]=S31;
260                         out->factors[iptr_main*9+2]=S31;                            out->factors[iptr_main*9+3]=S12;
261                         out->factors[iptr_main*9+3]=S12;                            out->factors[iptr_main*9+4]=S22;
262                         out->factors[iptr_main*9+4]=S22;                            out->factors[iptr_main*9+5]=S32;
263                         out->factors[iptr_main*9+5]=S32;                            out->factors[iptr_main*9+6]=S13;
264                         out->factors[iptr_main*9+6]=S13;                            out->factors[iptr_main*9+7]=S23;
265                         out->factors[iptr_main*9+7]=S23;                            out->factors[iptr_main*9+8]=S33;
266                         out->factors[iptr_main*9+8]=S33;    
267                              /* a_ik=a_ii^{-1}*a_ik */
268                         /* a_ik=a_ii^{-1}*a_ik */                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
269                         for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                               k=A->pattern->index[iptr_ik];
270                            k=A->pattern->index[iptr_ik];                               if (colorOf[k]>color) {
271                            if (out->colorOf[k]>color) {                                  A11=out->factors[iptr_ik*9  ];
272                               A11=out->factors[iptr_ik*9  ];                                  A21=out->factors[iptr_ik*9+1];
273                               A21=out->factors[iptr_ik*9+1];                                  A31=out->factors[iptr_ik*9+2];
274                               A31=out->factors[iptr_ik*9+2];                                  A12=out->factors[iptr_ik*9+3];
275                               A12=out->factors[iptr_ik*9+3];                                  A22=out->factors[iptr_ik*9+4];
276                               A22=out->factors[iptr_ik*9+4];                                  A32=out->factors[iptr_ik*9+5];
277                               A32=out->factors[iptr_ik*9+5];                                  A13=out->factors[iptr_ik*9+6];
278                               A13=out->factors[iptr_ik*9+6];                                  A23=out->factors[iptr_ik*9+7];
279                               A23=out->factors[iptr_ik*9+7];                                  A33=out->factors[iptr_ik*9+8];
280                               A33=out->factors[iptr_ik*9+8];                                  out->factors[iptr_ik*9  ]=S11*A11+S12*A21+S13*A31;                            
281                               out->factors[iptr_ik*9  ]=S11*A11+S12*A21+S13*A31;                                                              out->factors[iptr_ik*9+1]=S21*A11+S22*A21+S23*A31;
282                               out->factors[iptr_ik*9+1]=S21*A11+S22*A21+S23*A31;                                  out->factors[iptr_ik*9+2]=S31*A11+S32*A21+S33*A31;
283                               out->factors[iptr_ik*9+2]=S31*A11+S32*A21+S33*A31;                                  out->factors[iptr_ik*9+3]=S11*A12+S12*A22+S13*A32;                            
284                               out->factors[iptr_ik*9+3]=S11*A12+S12*A22+S13*A32;                                                              out->factors[iptr_ik*9+4]=S21*A12+S22*A22+S23*A32;
285                               out->factors[iptr_ik*9+4]=S21*A12+S22*A22+S23*A32;                                  out->factors[iptr_ik*9+5]=S31*A12+S32*A22+S33*A32;
286                               out->factors[iptr_ik*9+5]=S31*A12+S32*A22+S33*A32;                                  out->factors[iptr_ik*9+6]=S11*A13+S12*A23+S13*A33;                            
287                               out->factors[iptr_ik*9+6]=S11*A13+S12*A23+S13*A33;                                                              out->factors[iptr_ik*9+7]=S21*A13+S22*A23+S23*A33;
288                               out->factors[iptr_ik*9+7]=S21*A13+S22*A23+S23*A33;                                  out->factors[iptr_ik*9+8]=S31*A13+S32*A23+S33*A33;
289                               out->factors[iptr_ik*9+8]=S31*A13+S32*A23+S33*A33;                               }                              
290                            }                                                          }
291                           } else {
292                                Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");
293                         }                         }
                     } else {  
                          Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getILU: non-regular main diagonal block.");  
294                      }                      }
295                   }                   }
296                }                } else {
297             } else {                   Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: block size greater than 3 is not supported.");
298                Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: block size greater than 3 is not supported.");                }      
299             }                      #pragma omp barrier
300             #pragma omp barrier         }
301          }         time_fac=Paso_timer()-time0;
         time_fac=Paso_timer()-time0;  
302    }    }
   
   TMPMEMFREE(mis_marker);  
303    if (Paso_noError()) {    if (Paso_noError()) {
304        if (verbose) {        if (verbose) printf("timing: ILU: coloring/elemination : %e sec\n",time_fac);
          printf("ILU: %d color used \n",out->num_colors);  
          printf("timing: ILU: coloring/elemination : %e/%e\n",time_color,time_fac);  
      }  
305       return out;       return out;
306    } else  {    } else  {
307       Paso_Solver_ILU_free(out);       Paso_Solver_ILU_free(out);
# Line 345  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 320  Paso_Solver_ILU* Paso_Solver_getILU(Paso
320    
321  */  */
322    
323  void Paso_Solver_solveILU(Paso_Solver_ILU * ilu, double * x, double * b) {  void Paso_Solver_solveILU(Paso_SparseMatrix * A, Paso_Solver_ILU * ilu, double * x, const double * b) {
324       register dim_t i,k;       register dim_t i,k;
325       register index_t color,iptr_ik,iptr_main;       register index_t color,iptr_ik,iptr_main;
326       register double S1,S2,S3,R1,R2,R3;       register double S1,S2,S3,R1,R2,R3;
327       dim_t n_block=ilu->n_block;       const dim_t n=A->numRows;
328       dim_t n=ilu->n;       const dim_t n_block=A->row_block_size;
329         const index_t* colorOf = Paso_Pattern_borrowColoringPointer(A->pattern);
330         const dim_t num_colors = Paso_Pattern_getNumColors(A->pattern);
331         const index_t *ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A);
332            
333            
334       /* copy x into b*/       /* copy x into b*/
335       #pragma omp for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
336       for (i=0;i<n*n_block;++i) x[i]=b[i];       for (i=0;i<n*n_block;++i) x[i]=b[i];
337       /* forward substitution */       /* forward substitution */
338       for (color=0;color<ilu->num_colors;++color) {       for (color=0;color<num_colors;++color) {
339             if (n_block==1) {             if (n_block==1) {
340                #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)
341                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
342                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
343                       /* x_i=x_i-a_ik*x_k */                                           /* x_i=x_i-a_ik*x_k */                    
344                       S1=x[i];                       S1=x[i];
345                       for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {                       for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
346                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
347                            if (ilu->colorOf[k]<color) {                            if (colorOf[k]<color) {
348                               R1=x[k];                                                             R1=x[k];                              
349                               S1-=ilu->factors[iptr_ik]*R1;                               S1-=ilu->factors[iptr_ik]*R1;
350                            }                            }
351                       }                       }
352                       iptr_main=ilu->main_iptr[i];                       iptr_main=ptr_main[i];
353                       x[i]=ilu->factors[iptr_main]*S1;                       x[i]=ilu->factors[iptr_main]*S1;
354                     }                     }
355                }                }
356             } else if (n_block==2) {             } else if (n_block==2) {
357                #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)
358                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
359                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
360                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
361                       S1=x[2*i];                       S1=x[2*i];
362                       S2=x[2*i+1];                       S2=x[2*i+1];
363                       for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {                       for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
364                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
365                            if (ilu->colorOf[k]<color) {                            if (colorOf[k]<color) {
366                               R1=x[2*k];                               R1=x[2*k];
367                               R2=x[2*k+1];                               R2=x[2*k+1];
368                               S1-=ilu->factors[4*iptr_ik  ]*R1+ilu->factors[4*iptr_ik+2]*R2;                               S1-=ilu->factors[4*iptr_ik  ]*R1+ilu->factors[4*iptr_ik+2]*R2;
369                               S2-=ilu->factors[4*iptr_ik+1]*R1+ilu->factors[4*iptr_ik+3]*R2;                               S2-=ilu->factors[4*iptr_ik+1]*R1+ilu->factors[4*iptr_ik+3]*R2;
370                            }                            }
371                       }                       }
372                       iptr_main=ilu->main_iptr[i];                       iptr_main=ptr_main[i];
373                       x[2*i  ]=ilu->factors[4*iptr_main  ]*S1+ilu->factors[4*iptr_main+2]*S2;                       x[2*i  ]=ilu->factors[4*iptr_main  ]*S1+ilu->factors[4*iptr_main+2]*S2;
374                       x[2*i+1]=ilu->factors[4*iptr_main+1]*S1+ilu->factors[4*iptr_main+3]*S2;                       x[2*i+1]=ilu->factors[4*iptr_main+1]*S1+ilu->factors[4*iptr_main+3]*S2;
375                     }                     }
376    
377                }                }
378             } else if (n_block==3) {             } else if (n_block==3) {
379                #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)
380                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
381                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
382                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
383                       S1=x[3*i];                       S1=x[3*i];
384                       S2=x[3*i+1];                       S2=x[3*i+1];
385                       S3=x[3*i+2];                       S3=x[3*i+2];
386                       for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {                       for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
387                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
388                            if (ilu->colorOf[k]<color) {                            if (colorOf[k]<color) {
389                               R1=x[3*k];                               R1=x[3*k];
390                               R2=x[3*k+1];                               R2=x[3*k+1];
391                               R3=x[3*k+2];                               R3=x[3*k+2];
# Line 416  void Paso_Solver_solveILU(Paso_Solver_IL Line 394  void Paso_Solver_solveILU(Paso_Solver_IL
394                               S3-=ilu->factors[9*iptr_ik+2]*R1+ilu->factors[9*iptr_ik+5]*R2+ilu->factors[9*iptr_ik+8]*R3;                               S3-=ilu->factors[9*iptr_ik+2]*R1+ilu->factors[9*iptr_ik+5]*R2+ilu->factors[9*iptr_ik+8]*R3;
395                            }                            }
396                       }                       }
397                       iptr_main=ilu->main_iptr[i];                       iptr_main=ptr_main[i];
398                       x[3*i  ]=ilu->factors[9*iptr_main  ]*S1+ilu->factors[9*iptr_main+3]*S2+ilu->factors[9*iptr_main+6]*S3;                       x[3*i  ]=ilu->factors[9*iptr_main  ]*S1+ilu->factors[9*iptr_main+3]*S2+ilu->factors[9*iptr_main+6]*S3;
399                       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;                       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;
400                       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;                       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;
401                   }                   }
402                }                }
403             }             }
            #pragma omp barrier  
404       }       }
405       /* backward substitution */       /* backward substitution */
406       for (color=(ilu->num_colors)-1;color>-1;--color) {       for (color=(num_colors)-1;color>-1;--color) {
407             if (n_block==1) {             if (n_block==1) {
408                #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)
409                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
410                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
411                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
412                       S1=x[i];                       S1=x[i];
413                       for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {                       for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
414                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
415                            if (ilu->colorOf[k]>color) {                            if (colorOf[k]>color) {
416                               R1=x[k];                               R1=x[k];
417                               S1-=ilu->factors[iptr_ik]*R1;                               S1-=ilu->factors[iptr_ik]*R1;
418                            }                            }
# Line 444  void Paso_Solver_solveILU(Paso_Solver_IL Line 421  void Paso_Solver_solveILU(Paso_Solver_IL
421                     }                     }
422                }                }
423             } else if (n_block==2) {             } else if (n_block==2) {
424                #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)
425                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
426                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
427                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
428                       S1=x[2*i];                       S1=x[2*i];
429                       S2=x[2*i+1];                       S2=x[2*i+1];
430                       for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {                       for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
431                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
432                            if (ilu->colorOf[k]>color) {                            if (colorOf[k]>color) {
433                               R1=x[2*k];                               R1=x[2*k];
434                               R2=x[2*k+1];                               R2=x[2*k+1];
435                               S1-=ilu->factors[4*iptr_ik  ]*R1+ilu->factors[4*iptr_ik+2]*R2;                               S1-=ilu->factors[4*iptr_ik  ]*R1+ilu->factors[4*iptr_ik+2]*R2;
# Line 464  void Paso_Solver_solveILU(Paso_Solver_IL Line 441  void Paso_Solver_solveILU(Paso_Solver_IL
441                     }                     }
442                }                }
443             } else if (n_block==3) {             } else if (n_block==3) {
444                #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)
445                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
446                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
447                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
448                       S1=x[3*i  ];                       S1=x[3*i  ];
449                       S2=x[3*i+1];                       S2=x[3*i+1];
450                       S3=x[3*i+2];                       S3=x[3*i+2];
451                       for (iptr_ik=ilu->pattern->ptr[i];iptr_ik<ilu->pattern->ptr[i+1]; ++iptr_ik) {                       for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
452                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
453                            if (ilu->colorOf[k]>color) {                            if (colorOf[k]>color) {
454                               R1=x[3*k];                               R1=x[3*k];
455                               R2=x[3*k+1];                               R2=x[3*k+1];
456                               R3=x[3*k+2];                               R3=x[3*k+2];
# Line 492  void Paso_Solver_solveILU(Paso_Solver_IL Line 469  void Paso_Solver_solveILU(Paso_Solver_IL
469       }       }
470       return;       return;
471  }  }
472  /*  
  * $Log$  
  * Revision 1.2  2005/09/15 03:44:40  jgs  
  * Merge of development branch dev-02 back to main trunk on 2005-09-15  
  *  
  * Revision 1.1.2.1  2005/09/05 06:29:50  gross  
  * These files have been extracted from finley to define a stand alone libray for iterative  
  * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but  
  * has not been tested yet.  
  *  
  *  
  */  

Legend:
Removed from v.483  
changed lines
  Added in v.3120

  ViewVC Help
Powered by ViewVC 1.1.26