/[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

temp/paso/src/Solver_ILU.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC trunk/paso/src/Solver_ILU.c revision 3094 by gross, Fri Aug 13 08:38:06 2010 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2010 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 20  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    
# Line 34  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_Pattern_free(in->pattern);    
37          MEMFREE(in);          MEMFREE(in);
38       }       }
39  }  }
# Line 48  void Paso_Solver_ILU_free(Paso_Solver_IL Line 44  void Paso_Solver_ILU_free(Paso_Solver_IL
44    
45  */  */
46  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SparseMatrix * A,bool_t verbose) {  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SparseMatrix * A,bool_t verbose) {
47    dim_t n=A->numRows;    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, *mis_marker=NULL;    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;
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;
   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_Pattern_getReference(A->pattern);    if ( ! Paso_checkPtr(out->factors)  ) {
   out->n_block=n_block;  
   out->n=n;  
   if ( !(Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->factors)) ) {  
     time0=Paso_timer();  
     Paso_Pattern_color(A->pattern,&out->num_colors,out->colorOf);  
     time_color=Paso_timer()-time0;  
62    
     if (Paso_noError()) {  
63         time0=Paso_timer();         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");  
            }  
        }  
64         /* start factorization */         /* start factorization */
65        
66         #pragma omp barrier         #pragma omp barrier
67         for (color=0;color<out->num_colors && Paso_noError();++color) {         for (color=0;color<num_colors && Paso_noError();++color) {
68                if (n_block==1) {                if (n_block==1) {
69                   #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,color2,iptr_ik,k,iptr_kj,S11,j,iptr_ij,A11,iptr_main,D)
70                   for (i = 0; i < n; ++i) {                   for (i = 0; i < n; ++i) {
71                      if (out->colorOf[i]==color) {                      if (colorOf[i]==color) {
72                         for (color2=0;color2<color;++color2) {                         for (color2=0;color2<color;++color2) {
73                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
74                               k=A->pattern->index[iptr_ik];                                                         k=A->pattern->index[iptr_ik];                          
75                               if (out->colorOf[k]==color2) {                               if (colorOf[k]==color2) {
76                                  A11=out->factors[iptr_ik];                                  A11=out->factors[iptr_ik];
77                                  /* a_ij=a_ij-a_ik*a_kj */                                  /* a_ij=a_ij-a_ik*a_kj */
78                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
79                                     j=A->pattern->index[iptr_kj];                                     j=A->pattern->index[iptr_kj];
80                                     if (out->colorOf[j]>color2) {                                     if (colorOf[j]>color2) {
81                                        S11=out->factors[iptr_kj];                                        S11=out->factors[iptr_kj];
82                                        for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {                                        for (iptr_ij=A->pattern->ptr[i];iptr_ij<A->pattern->ptr[i+1]; iptr_ij++) {
83                                           if (j==A->pattern->index[iptr_ij]) {                                           if (j==A->pattern->index[iptr_ij]) {
# Line 115  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 90  Paso_Solver_ILU* Paso_Solver_getILU(Paso
90                               }                               }
91                            }                            }
92                         }                         }
93                         iptr_main=out->main_iptr[i];                         iptr_main=ptr_main[i];
94                         D=out->factors[iptr_main];                         D=out->factors[iptr_main];
95                         if (ABS(D)>0.) {                         if (ABS(D)>0.) {
96                            D=1./D;                            D=1./D;
# Line 123  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 98  Paso_Solver_ILU* Paso_Solver_getILU(Paso
98                            /* a_ik=a_ii^{-1}*a_ik */                            /* a_ik=a_ii^{-1}*a_ik */
99                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
100                               k=A->pattern->index[iptr_ik];                               k=A->pattern->index[iptr_ik];
101                               if (out->colorOf[k]>color) {                               if (colorOf[k]>color) {
102                                  A11=out->factors[iptr_ik];                                  A11=out->factors[iptr_ik];
103                                  out->factors[iptr_ik]=A11*D;                                  out->factors[iptr_ik]=A11*D;
104                               }                                                             }                              
# Line 136  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 111  Paso_Solver_ILU* Paso_Solver_getILU(Paso
111                } else if (n_block==2) {                } else if (n_block==2) {
112                   #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)                   #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)
113                   for (i = 0; i < n; ++i) {                   for (i = 0; i < n; ++i) {
114                      if (out->colorOf[i]==color) {                      if (colorOf[i]==color) {
115                         for (color2=0;color2<color;++color2) {                         for (color2=0;color2<color;++color2) {
116                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
117                               k=A->pattern->index[iptr_ik];                                                         k=A->pattern->index[iptr_ik];                          
118                               if (out->colorOf[k]==color2) {                               if (colorOf[k]==color2) {
119                                  A11=out->factors[iptr_ik*4  ];                                  A11=out->factors[iptr_ik*4  ];
120                                  A21=out->factors[iptr_ik*4+1];                                  A21=out->factors[iptr_ik*4+1];
121                                  A12=out->factors[iptr_ik*4+2];                                  A12=out->factors[iptr_ik*4+2];
# Line 148  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 123  Paso_Solver_ILU* Paso_Solver_getILU(Paso
123                                  /* a_ij=a_ij-a_ik*a_kj */                                  /* a_ij=a_ij-a_ik*a_kj */
124                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
125                                     j=A->pattern->index[iptr_kj];                                     j=A->pattern->index[iptr_kj];
126                                     if (out->colorOf[j]>color2) {                                     if (colorOf[j]>color2) {
127                                        S11=out->factors[iptr_kj*4];                                        S11=out->factors[iptr_kj*4];
128                                        S21=out->factors[iptr_kj*4+1];                                        S21=out->factors[iptr_kj*4+1];
129                                        S12=out->factors[iptr_kj*4+2];                                        S12=out->factors[iptr_kj*4+2];
# Line 167  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 142  Paso_Solver_ILU* Paso_Solver_getILU(Paso
142                               }                               }
143                            }                            }
144                         }                         }
145                         iptr_main=out->main_iptr[i];                         iptr_main=ptr_main[i];
146                         A11=out->factors[iptr_main*4];                         A11=out->factors[iptr_main*4];
147                         A21=out->factors[iptr_main*4+1];                         A21=out->factors[iptr_main*4+1];
148                         A12=out->factors[iptr_main*4+2];                         A12=out->factors[iptr_main*4+2];
# Line 186  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 161  Paso_Solver_ILU* Paso_Solver_getILU(Paso
161                            /* a_ik=a_ii^{-1}*a_ik */                            /* a_ik=a_ii^{-1}*a_ik */
162                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
163                               k=A->pattern->index[iptr_ik];                               k=A->pattern->index[iptr_ik];
164                               if (out->colorOf[k]>color) {                               if (colorOf[k]>color) {
165                                  A11=out->factors[iptr_ik*4  ];                                  A11=out->factors[iptr_ik*4  ];
166                                  A21=out->factors[iptr_ik*4+1];                                  A21=out->factors[iptr_ik*4+1];
167                                  A12=out->factors[iptr_ik*4+2];                                  A12=out->factors[iptr_ik*4+2];
# Line 205  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 180  Paso_Solver_ILU* Paso_Solver_getILU(Paso
180                } else if (n_block==3) {                } else if (n_block==3) {
181                   #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)                   #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)
182                   for (i = 0; i < n; ++i) {                   for (i = 0; i < n; ++i) {
183                      if (out->colorOf[i]==color) {                      if (colorOf[i]==color) {
184                         for (color2=0;color2<color;++color2) {                         for (color2=0;color2<color;++color2) {
185                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
186                               k=A->pattern->index[iptr_ik];                                                         k=A->pattern->index[iptr_ik];                          
187                               if (out->colorOf[k]==color2) {                               if (colorOf[k]==color2) {
188                                  A11=out->factors[iptr_ik*9  ];                                  A11=out->factors[iptr_ik*9  ];
189                                  A21=out->factors[iptr_ik*9+1];                                  A21=out->factors[iptr_ik*9+1];
190                                  A31=out->factors[iptr_ik*9+2];                                  A31=out->factors[iptr_ik*9+2];
# Line 222  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 197  Paso_Solver_ILU* Paso_Solver_getILU(Paso
197                                  /* a_ij=a_ij-a_ik*a_kj */                                  /* a_ij=a_ij-a_ik*a_kj */
198                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {                                  for (iptr_kj=A->pattern->ptr[k];iptr_kj<A->pattern->ptr[k+1]; iptr_kj++) {
199                                     j=A->pattern->index[iptr_kj];                                     j=A->pattern->index[iptr_kj];
200                                     if (out->colorOf[j]>color2) {                                     if (colorOf[j]>color2) {
201                                        S11=out->factors[iptr_kj*9  ];                                        S11=out->factors[iptr_kj*9  ];
202                                        S21=out->factors[iptr_kj*9+1];                                        S21=out->factors[iptr_kj*9+1];
203                                        S31=out->factors[iptr_kj*9+2];                                        S31=out->factors[iptr_kj*9+2];
# Line 251  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 226  Paso_Solver_ILU* Paso_Solver_getILU(Paso
226                               }                               }
227                            }                            }
228                         }                         }
229                         iptr_main=out->main_iptr[i];                         iptr_main=ptr_main[i];
230                         A11=out->factors[iptr_main*9  ];                         A11=out->factors[iptr_main*9  ];
231                         A21=out->factors[iptr_main*9+1];                         A21=out->factors[iptr_main*9+1];
232                         A31=out->factors[iptr_main*9+2];                         A31=out->factors[iptr_main*9+2];
# Line 287  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 262  Paso_Solver_ILU* Paso_Solver_getILU(Paso
262                            /* a_ik=a_ii^{-1}*a_ik */                            /* a_ik=a_ii^{-1}*a_ik */
263                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {                            for (iptr_ik=A->pattern->ptr[i];iptr_ik<A->pattern->ptr[i+1]; ++iptr_ik) {
264                               k=A->pattern->index[iptr_ik];                               k=A->pattern->index[iptr_ik];
265                               if (out->colorOf[k]>color) {                               if (colorOf[k]>color) {
266                                  A11=out->factors[iptr_ik*9  ];                                  A11=out->factors[iptr_ik*9  ];
267                                  A21=out->factors[iptr_ik*9+1];                                  A21=out->factors[iptr_ik*9+1];
268                                  A31=out->factors[iptr_ik*9+2];                                  A31=out->factors[iptr_ik*9+2];
# Line 319  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 294  Paso_Solver_ILU* Paso_Solver_getILU(Paso
294                #pragma omp barrier                #pragma omp barrier
295         }         }
296         time_fac=Paso_timer()-time0;         time_fac=Paso_timer()-time0;
      }  
297    }    }
298    if (Paso_noError()) {    if (Paso_noError()) {
299        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);  
      }  
300       return out;       return out;
301    } else  {    } else  {
302       Paso_Solver_ILU_free(out);       Paso_Solver_ILU_free(out);
# Line 344  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 315  Paso_Solver_ILU* Paso_Solver_getILU(Paso
315    
316  */  */
317    
318  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) {
319       register dim_t i,k;       register dim_t i,k;
320       register index_t color,iptr_ik,iptr_main;       register index_t color,iptr_ik,iptr_main;
321       register double S1,S2,S3,R1,R2,R3;       register double S1,S2,S3,R1,R2,R3;
322       dim_t n_block=ilu->n_block;       const dim_t n=A->numRows;
323       dim_t n=ilu->n;       const dim_t n_block=A->row_block_size;
324         const index_t* colorOf = Paso_Pattern_borrowColoringPointer(A->pattern);
325         const dim_t num_colors = Paso_Pattern_getNumColors(A->pattern);
326         const index_t *ptr_main = Paso_SparseMatrix_borrowMainDiagonalPointer(A);
327            
328            
329       /* copy x into b*/       /* copy x into b*/
330       #pragma omp for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
331       for (i=0;i<n*n_block;++i) x[i]=b[i];       for (i=0;i<n*n_block;++i) x[i]=b[i];
332       /* forward substitution */       /* forward substitution */
333       for (color=0;color<ilu->num_colors;++color) {       for (color=0;color<num_colors;++color) {
334             if (n_block==1) {             if (n_block==1) {
335                #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)
336                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
337                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
338                       /* x_i=x_i-a_ik*x_k */                                           /* x_i=x_i-a_ik*x_k */                    
339                       S1=x[i];                       S1=x[i];
340                       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) {
341                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
342                            if (ilu->colorOf[k]<color) {                            if (colorOf[k]<color) {
343                               R1=x[k];                                                             R1=x[k];                              
344                               S1-=ilu->factors[iptr_ik]*R1;                               S1-=ilu->factors[iptr_ik]*R1;
345                            }                            }
346                       }                       }
347                       iptr_main=ilu->main_iptr[i];                       iptr_main=ptr_main[i];
348                       x[i]=ilu->factors[iptr_main]*S1;                       x[i]=ilu->factors[iptr_main]*S1;
349                     }                     }
350                }                }
351             } else if (n_block==2) {             } else if (n_block==2) {
352                #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)
353                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
354                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
355                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
356                       S1=x[2*i];                       S1=x[2*i];
357                       S2=x[2*i+1];                       S2=x[2*i+1];
358                       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) {
359                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
360                            if (ilu->colorOf[k]<color) {                            if (colorOf[k]<color) {
361                               R1=x[2*k];                               R1=x[2*k];
362                               R2=x[2*k+1];                               R2=x[2*k+1];
363                               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;
364                               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;
365                            }                            }
366                       }                       }
367                       iptr_main=ilu->main_iptr[i];                       iptr_main=ptr_main[i];
368                       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;
369                       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;
370                     }                     }
371    
372                }                }
373             } else if (n_block==3) {             } else if (n_block==3) {
374                #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)
375                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
376                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
377                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
378                       S1=x[3*i];                       S1=x[3*i];
379                       S2=x[3*i+1];                       S2=x[3*i+1];
380                       S3=x[3*i+2];                       S3=x[3*i+2];
381                       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) {
382                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
383                            if (ilu->colorOf[k]<color) {                            if (colorOf[k]<color) {
384                               R1=x[3*k];                               R1=x[3*k];
385                               R2=x[3*k+1];                               R2=x[3*k+1];
386                               R3=x[3*k+2];                               R3=x[3*k+2];
# Line 415  void Paso_Solver_solveILU(Paso_Solver_IL Line 389  void Paso_Solver_solveILU(Paso_Solver_IL
389                               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;
390                            }                            }
391                       }                       }
392                       iptr_main=ilu->main_iptr[i];                       iptr_main=ptr_main[i];
393                       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;
394                       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;
395                       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;
396                   }                   }
397                }                }
398             }             }
            #pragma omp barrier  
399       }       }
400       /* backward substitution */       /* backward substitution */
401       for (color=(ilu->num_colors)-1;color>-1;--color) {       for (color=(num_colors)-1;color>-1;--color) {
402             if (n_block==1) {             if (n_block==1) {
403                #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)
404                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
405                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
406                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
407                       S1=x[i];                       S1=x[i];
408                       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) {
409                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
410                            if (ilu->colorOf[k]>color) {                            if (colorOf[k]>color) {
411                               R1=x[k];                               R1=x[k];
412                               S1-=ilu->factors[iptr_ik]*R1;                               S1-=ilu->factors[iptr_ik]*R1;
413                            }                            }
# Line 443  void Paso_Solver_solveILU(Paso_Solver_IL Line 416  void Paso_Solver_solveILU(Paso_Solver_IL
416                     }                     }
417                }                }
418             } else if (n_block==2) {             } else if (n_block==2) {
419                #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)
420                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
421                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
422                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
423                       S1=x[2*i];                       S1=x[2*i];
424                       S2=x[2*i+1];                       S2=x[2*i+1];
425                       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) {
426                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
427                            if (ilu->colorOf[k]>color) {                            if (colorOf[k]>color) {
428                               R1=x[2*k];                               R1=x[2*k];
429                               R2=x[2*k+1];                               R2=x[2*k+1];
430                               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 463  void Paso_Solver_solveILU(Paso_Solver_IL Line 436  void Paso_Solver_solveILU(Paso_Solver_IL
436                     }                     }
437                }                }
438             } else if (n_block==3) {             } else if (n_block==3) {
439                #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)
440                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
441                     if (ilu->colorOf[i]==color) {                     if (colorOf[i]==color) {
442                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
443                       S1=x[3*i  ];                       S1=x[3*i  ];
444                       S2=x[3*i+1];                       S2=x[3*i+1];
445                       S3=x[3*i+2];                       S3=x[3*i+2];
446                       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) {
447                            k=ilu->pattern->index[iptr_ik];                                                      k=A->pattern->index[iptr_ik];                          
448                            if (ilu->colorOf[k]>color) {                            if (colorOf[k]>color) {
449                               R1=x[3*k];                               R1=x[3*k];
450                               R2=x[3*k+1];                               R2=x[3*k+1];
451                               R3=x[3*k+2];                               R3=x[3*k+2];
# Line 491  void Paso_Solver_solveILU(Paso_Solver_IL Line 464  void Paso_Solver_solveILU(Paso_Solver_IL
464       }       }
465       return;       return;
466  }  }
467  /*  
  * $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.1387  
changed lines
  Added in v.3094

  ViewVC Help
Powered by ViewVC 1.1.26