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

Revision 97 - (show annotations)
Tue Dec 14 05:39:33 2004 UTC (15 years, 1 month ago) by jgs
File MIME type: text/plain
File size: 7778 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 176 /* 177 * \$Log\$ 178 * Revision 1.2 2004/12/14 05:39:32 jgs 179 * *** empty log message *** 180 * 181 * Revision 1.1.1.1.2.2 2004/11/24 01:37:17 gross 182 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now 183 * 184 * Revision 1.1.1.1.2.1 2004/11/12 06:58:21 gross 185 * a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry 186 * 187 * Revision 1.1.1.1 2004/10/26 06:53:58 jgs 188 * initial import of project esys2 189 * 190 * Revision 1.1 2004/07/02 04:21:14 gross 191 * Finley C code has been included 192 * 193 * 194 */

## Properties

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