/[escript]/trunk/paso/src/Solver_ILU.c
ViewVC logotype

Diff of /trunk/paso/src/Solver_ILU.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/paso/src/Solvers/Solver_ILU.c revision 431 by gross, Fri Jan 13 05:07:10 2006 UTC trunk/paso/src/Solver_ILU.c revision 1028 by gross, Wed Mar 14 00:15:24 2007 UTC
# Line 1  Line 1 
1  /* $Id$ */  /* $Id$ */
2    
3    
4    /*
5    ********************************************************************************
6    *               Copyright   2006 by ACcESS MNRF                                *
7    *                                                                              *
8    *                 http://www.access.edu.au                                     *
9    *           Primary Business: Queensland, Australia                            *
10    *     Licensed under the Open Software License version 3.0             *
11    *        http://www.opensource.org/licenses/osl-3.0.php                        *
12    ********************************************************************************
13    */
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /* Paso: ILU preconditioner with reordering                 */  /* Paso: ILU preconditioner with reordering                 */
# Line 13  Line 25 
25    
26  #include "Paso.h"  #include "Paso.h"
27  #include "Solver.h"  #include "Solver.h"
28  #include "Util.h"  #include "PasoUtil.h"
29    
30  /**************************************************************/  /**************************************************************/
31    
# Line 37  void Paso_Solver_ILU_free(Paso_Solver_IL Line 49  void Paso_Solver_ILU_free(Paso_Solver_IL
49  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SystemMatrix * A,bool_t verbose) {  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SystemMatrix * A,bool_t verbose) {
50    dim_t n=A->num_rows;    dim_t n=A->num_rows;
51    dim_t n_block=A->row_block_size;    dim_t n_block=A->row_block_size;
52    index_t num_colors=0;    index_t num_colors=0, *mis_marker=NULL;
53    register double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;    register double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;
54    register double mainA11,mainA12,mainA13,mainA21,mainA22,mainA23,mainA31,mainA32,mainA33;    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;
# Line 46  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 58  Paso_Solver_ILU* Paso_Solver_getILU(Paso
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;
61    index_t* mis_marker=TMPMEMALLOC(n,index_t);    mis_marker=TMPMEMALLOC(n,index_t);
62    out->colorOf=MEMALLOC(n,index_t);    out->colorOf=MEMALLOC(n,index_t);
63    out->factors=MEMALLOC(A->len,double);    out->factors=MEMALLOC(A->len,double);
64    out->main_iptr=MEMALLOC(n,double);    out->main_iptr=MEMALLOC(n,index_t);
65    out->pattern=Paso_SystemMatrixPattern_reference(A->pattern);    out->pattern=Paso_SystemMatrixPattern_reference(A->pattern);
66    out->n_block=n_block;    out->n_block=n_block;
67    out->n=n;    out->n=n;
# Line 73  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 85  Paso_Solver_ILU* Paso_Solver_getILU(Paso
85      time_color=Paso_timer()-time0;      time_color=Paso_timer()-time0;
86      time0=Paso_timer();      time0=Paso_timer();
87      /* find main diagonal and copy matrix values */      /* find main diagonal and copy matrix values */
88      #pragma omp for private(i) schedule(static) local(i,iptr,iptr_main,k)      #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)
89      for (i = 0; i < n; ++i) {      for (i = 0; i < n; ++i) {
90          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {          for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
91              iptr_main=A->pattern->ptr[0]-1;              iptr_main=A->pattern->ptr[0]-1;
# Line 88  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 100  Paso_Solver_ILU* Paso_Solver_getILU(Paso
100      }      }
101      /* start factorization */      /* start factorization */
102    
103        #pragma omp barrier
104      for (color=0;color<out->num_colors && Paso_noError();++color) {      for (color=0;color<out->num_colors && Paso_noError();++color) {
105             if (n_block==1) {             if (n_block==1) {
106                #pragma omp for private(i) 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)
107                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
108                   if (out->colorOf[i]==color) {                   if (out->colorOf[i]==color) {
109                      for (color2=0;color2<color;++color2) {                      for (color2=0;color2<color;++color2) {
# Line 133  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 146  Paso_Solver_ILU* Paso_Solver_getILU(Paso
146                   }                   }
147                }                }
148             } else if (n_block==2) {             } else if (n_block==2) {
149                #pragma omp for private(i) 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)
150                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
151                   if (out->colorOf[i]==color) {                   if (out->colorOf[i]==color) {
152                      for (color2=0;color2<color;++color2) {                      for (color2=0;color2<color;++color2) {
# Line 202  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 215  Paso_Solver_ILU* Paso_Solver_getILU(Paso
215                   }                   }
216                }                }
217             } else if (n_block==3) {             } else if (n_block==3) {
218                #pragma omp for private(i) 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)
219                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
220                   if (out->colorOf[i]==color) {                   if (out->colorOf[i]==color) {
221                      for (color2=0;color2<color;++color2) {                      for (color2=0;color2<color;++color2) {
# Line 315  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 328  Paso_Solver_ILU* Paso_Solver_getILU(Paso
328             } else {             } else {
329                Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: block size greater than 3 is not supported.");                Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: block size greater than 3 is not supported.");
330             }                   }      
331               #pragma omp barrier
332          }          }
333          time_fac=Paso_timer()-time0;          time_fac=Paso_timer()-time0;
334    }    }
# Line 357  void Paso_Solver_solveILU(Paso_Solver_IL Line 371  void Paso_Solver_solveILU(Paso_Solver_IL
371       /* forward substitution */       /* forward substitution */
372       for (color=0;color<ilu->num_colors;++color) {       for (color=0;color<ilu->num_colors;++color) {
373             if (n_block==1) {             if (n_block==1) {
374                #pragma omp for private(i) schedule(static) private(i,iptr_ik,k,S1,R1)                #pragma omp for schedule(static) private(i,iptr_ik,k,S1,R1,iptr_main)
375                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
376                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
377                       /* x_i=x_i-a_ik*x_k */                                           /* x_i=x_i-a_ik*x_k */                    
# Line 374  void Paso_Solver_solveILU(Paso_Solver_IL Line 388  void Paso_Solver_solveILU(Paso_Solver_IL
388                     }                     }
389                }                }
390             } else if (n_block==2) {             } else if (n_block==2) {
391                #pragma omp for private(i) schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2)                #pragma omp for schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2)
392                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
393                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
394                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 396  void Paso_Solver_solveILU(Paso_Solver_IL Line 410  void Paso_Solver_solveILU(Paso_Solver_IL
410    
411                }                }
412             } else if (n_block==3) {             } else if (n_block==3) {
413                #pragma omp for private(i) schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3)                #pragma omp for schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3)
414                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
415                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
416                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 426  void Paso_Solver_solveILU(Paso_Solver_IL Line 440  void Paso_Solver_solveILU(Paso_Solver_IL
440       /* backward substitution */       /* backward substitution */
441       for (color=(ilu->num_colors)-1;color>-1;--color) {       for (color=(ilu->num_colors)-1;color>-1;--color) {
442             if (n_block==1) {             if (n_block==1) {
443                #pragma omp for private(i) schedule(static) private(i,iptr_ik,k,S1)                #pragma omp for schedule(static) private(i,iptr_ik,k,S1,R1)
444                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
445                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
446                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 442  void Paso_Solver_solveILU(Paso_Solver_IL Line 456  void Paso_Solver_solveILU(Paso_Solver_IL
456                     }                     }
457                }                }
458             } else if (n_block==2) {             } else if (n_block==2) {
459                #pragma omp for private(i) schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2)                #pragma omp for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2)
460                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
461                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
462                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 462  void Paso_Solver_solveILU(Paso_Solver_IL Line 476  void Paso_Solver_solveILU(Paso_Solver_IL
476                     }                     }
477                }                }
478             } else if (n_block==3) {             } else if (n_block==3) {
479                #pragma omp for private(i) schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)                #pragma omp for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)
480                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
481                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
482                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */

Legend:
Removed from v.431  
changed lines
  Added in v.1028

  ViewVC Help
Powered by ViewVC 1.1.26