# Contents of /trunk/esys2/finley/src/finleyC/Solvers/Solver_ILU0.c

Revision 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (14 years, 4 months ago) by jgs
File MIME type: text/plain
File size: 7094 byte(s)
*** empty log message ***

 1 /* \$Id\$ */ 2 3 /**************************************************************/ 4 5 /* Finley: ILU(0) preconditioner: */ 6 7 /**************************************************************/ 8 9 /* Copyrights by ACcESS Australia 2003/04 */ 10 /* Author: gross@access.edu.au */ 11 12 /**************************************************************/ 13 14 #include "Finley.h" 15 #include "System.h" 16 #include "Solver.h" 17 18 void Finley_Solver_setILU0_1(int,maybelong*,maybelong*,double*,maybelong,maybelong*,maybelong*); 19 void Finley_Solver_solveILU0_1(int,double*,maybelong*,maybelong*,double*,maybelong,maybelong*,maybelong*); 20 21 22 /**************************************************************/ 23 24 /* inverse preconditioner setup */ 25 26 void Finley_Solver_setILU0(Finley_SystemMatrix * A_p) { 27 #if ITERATIVE_SOLVER == NO_LIB 28 Finley_Solver_Preconditioner* prec=(Finley_Solver_Preconditioner*) (A_p->iterative); 29 maybelong n=A_p->num_rows; 30 /* allocate arrays */ 31 prec->mainDiag =MEMALLOC(n,maybelong); 32 prec->color = MEMALLOC(n,maybelong); 33 prec->values = MEMALLOC((A_p->len)*(size_t) n,double); 34 if (! (Finley_checkPtr(prec->mainDiag) || Finley_checkPtr(prec->color) || Finley_checkPtr(prec->values))) { 35 /* find the main diagonal: */ 36 Finley_Solver_getMainDiagonal(A_p,prec->mainDiag); 37 /* color the rows : */ 38 Finley_Solver_coloring(A_p->pattern,&(prec->numColors),prec->color); 39 if (Finley_ErrorCode==NO_ERROR) { 40 /* allocate vector to hold main diagonal entries: */ 41 Finley_SystemMatrix_copy(A_p,prec->values); 42 /* allocate vector to hold main diagonal entries: */ 43 if (A_p->row_block_size > 1) { 44 Finley_ErrorCode = TYPE_ERROR; 45 sprintf(Finley_ErrorMsg, "ILU preconditioner requires block size =1 ."); 46 return; 47 } else { 48 Finley_Solver_setILU0_1(n,A_p->pattern->ptr,A_p->pattern->index,prec->values,prec->numColors,prec->color,prec->mainDiag); 49 } 50 } 51 } 52 if (Finley_ErrorCode!=NO_ERROR) Finley_Preconditioner_free(A_p->iterative); 53 return; 54 #endif 55 } 56 57 /**************************************************************/ 58 59 /* inverse preconditioner solution using raw pointers */ 60 61 /* should be called within a parallel region */ 62 /* barrier synconization should be performed to make sure that the input vector available */ 63 64 void Finley_Solver_solveILU0(Finley_SystemMatrix * A_p, double * x, double * b) { 65 #if ITERATIVE_SOLVER == NO_LIB 66 Finley_Solver_Preconditioner* prec=(Finley_Solver_Preconditioner*) (A_p->iterative); 67 maybelong n=A_p->num_rows; 68 maybelong i; 69 #pragma omp for private(i) 70 for (i=0;irow_block_size;i++) x[i]=b[i]; 71 if (A_p->row_block_size > 1) { 72 /* Finley_Solver_solveILU0_blk(n,x,A_p->pattern->ptr,A_p->pattern->index,prec->values,prec->numColors,prec->color,prec->mainDiag); */ 73 } else { 74 Finley_Solver_solveILU0_1(n,x,A_p->pattern->ptr,A_p->pattern->index,prec->values,prec->numColors,prec->color,prec->mainDiag); 75 } 76 return; 77 #endif 78 } 79 /**************************************************************/ 80 81 /* applies ILU preconditioner on x for block size 1 (in place) : */ 82 83 void Finley_Solver_solveILU0_1(int n,double* x, maybelong* ptr,maybelong* index,double* values, 84 maybelong numColors,maybelong* color, maybelong* mainDiag) { 85 #if ITERATIVE_SOLVER == NO_LIB 86 int iColor,i,k,ik; 87 /* forward substitution */ 88 iColor=0; 89 while (iColorcolor[i]) x[i]-=values[ik]*x[k]; 111 } /* end of ik-loop */ 112 x[i]*=values[mainDiag[i]]; 113 } 114 } 115 } 116 #endif 117 } 118 /**************************************************************/ 119 120 /* ILU preconditioner for block size 1 : */ 121 122 void Finley_Solver_setILU0_1(int n,maybelong* ptr,maybelong* index,double* values, 123 maybelong numColors,maybelong* color, maybelong* mainDiag) { 124 #if ITERATIVE_SOLVER == NO_LIB 125 double value_ik; 126 int iColor,i,k,j,ij,kj,ik,kColor; 127 #pragma omp parallel private(iColor) 128 { 129 for (iColor=0;iColorkColor) { 147 /* is (i,j) in the shape of the matrix? */ 148 /* this means that index[ij]=j */ 149 for (ij=ptr[i];ij0.) { 163 values[mainDiag[i]]=1./values[mainDiag[i]]; 164 } else { 165 Finley_ErrorCode = ZERO_DIVISION_ERROR; 166 sprintf(Finley_ErrorMsg, "ILU factorization failed in row %d.",i); 167 } 168 /* printf("finalizing row %d : entry %e\n",i,values[mainDiag[i]]); */ 169 } /* end of color if */ 170 } /* end of row loop */ 171 } /* end of color loop */ 172 } 173 #endif 174 } 175

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision