/[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 2608 by jfenwick, Tue Aug 18 01:25:18 2009 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 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 50  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_SparseMatrix * A,bool_t verbose) {  Paso_Solver_ILU* Paso_Solver_getILU(Paso_SparseMatrix * A,bool_t verbose) {
50    dim_t n=A->numRows;    dim_t n=A->numRows;
51    dim_t n_block=A->row_block_size;    dim_t n_block=A->row_block_size;
   index_t num_colors=0, *mis_marker=NULL;  
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,iptr_ik,k,iptr_kj,j,iptr_ij,color,color2;
55    double time0,time_color,time_fac;    double time0=0,time_color=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;
# Line 75  Paso_Solver_ILU* Paso_Solver_getILU(Paso Line 72  Paso_Solver_ILU* Paso_Solver_getILU(Paso
72         /* find main diagonal and copy matrix values */         /* find main diagonal and copy matrix values */
73         #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)         #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)
74         for (i = 0; i < n; ++i) {         for (i = 0; i < n; ++i) {
            for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {  
75                 iptr_main=A->pattern->ptr[0]-1;                 iptr_main=A->pattern->ptr[0]-1;
76                  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++) {
77                     if (A->pattern->index[iptr]==i) iptr_main=iptr;                     if (A->pattern->index[iptr]==i) iptr_main=iptr;
78                     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];                     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];
79                 }                 }
80                 out->main_iptr[i]=iptr_main;                 out->main_iptr[i]=iptr_main;
81                 if (iptr_main==A->pattern->ptr[0]-1)                 if (iptr_main==A->pattern->ptr[0]-1)  {
82                    Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: no main diagonal");                    Paso_setError(VALUE_ERROR, "Paso_Solver_getILU: no main diagonal");
83             }                 }
84         }         }
85         /* start factorization */         /* start factorization */
86        
# Line 353  void Paso_Solver_solveILU(Paso_Solver_IL Line 349  void Paso_Solver_solveILU(Paso_Solver_IL
349            
350            
351       /* copy x into b*/       /* copy x into b*/
352       #pragma omp for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
353       for (i=0;i<n*n_block;++i) x[i]=b[i];       for (i=0;i<n*n_block;++i) x[i]=b[i];
354       /* forward substitution */       /* forward substitution */
355       for (color=0;color<ilu->num_colors;++color) {       for (color=0;color<ilu->num_colors;++color) {
356             if (n_block==1) {             if (n_block==1) {
357                #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)
358                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
359                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
360                       /* x_i=x_i-a_ik*x_k */                                           /* x_i=x_i-a_ik*x_k */                    
# Line 375  void Paso_Solver_solveILU(Paso_Solver_IL Line 371  void Paso_Solver_solveILU(Paso_Solver_IL
371                     }                     }
372                }                }
373             } else if (n_block==2) {             } else if (n_block==2) {
374                #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)
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 397  void Paso_Solver_solveILU(Paso_Solver_IL Line 393  void Paso_Solver_solveILU(Paso_Solver_IL
393    
394                }                }
395             } else if (n_block==3) {             } else if (n_block==3) {
396                #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)
397                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
398                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
399                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 422  void Paso_Solver_solveILU(Paso_Solver_IL Line 418  void Paso_Solver_solveILU(Paso_Solver_IL
418                   }                   }
419                }                }
420             }             }
            #pragma omp barrier  
421       }       }
422       /* backward substitution */       /* backward substitution */
423       for (color=(ilu->num_colors)-1;color>-1;--color) {       for (color=(ilu->num_colors)-1;color>-1;--color) {
424             if (n_block==1) {             if (n_block==1) {
425                #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)
426                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
427                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
428                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 443  void Paso_Solver_solveILU(Paso_Solver_IL Line 438  void Paso_Solver_solveILU(Paso_Solver_IL
438                     }                     }
439                }                }
440             } else if (n_block==2) {             } else if (n_block==2) {
441                #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)
442                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
443                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
444                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 463  void Paso_Solver_solveILU(Paso_Solver_IL Line 458  void Paso_Solver_solveILU(Paso_Solver_IL
458                     }                     }
459                }                }
460             } else if (n_block==3) {             } else if (n_block==3) {
461                #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)
462                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
463                     if (ilu->colorOf[i]==color) {                     if (ilu->colorOf[i]==color) {
464                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */

Legend:
Removed from v.1387  
changed lines
  Added in v.2608

  ViewVC Help
Powered by ViewVC 1.1.26