# Contents of /trunk/esys2/finley/src/finleyC/Solvers/Solver_jacobi.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: 4654 byte(s)
```*** empty log message ***

```
 1 /* \$Id\$ */ 2 3 /**************************************************************/ 4 5 /* Finley: jacobi 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 #include "Util.h" 18 #include "Common.h" 19 20 /**************************************************************/ 21 22 /* inverse preconditioner setup */ 23 24 void Finley_Solver_setJacobi(Finley_SystemMatrix * A_p) { 25 #if ITERATIVE_SOLVER == NO_LIB 26 Finley_Solver_Preconditioner* prec=(Finley_Solver_Preconditioner*) (A_p->iterative); 27 int n = A_p->num_cols; 28 int n_block=A_p->row_block_size; 29 int block_size=A_p->block_size; 30 int i,iPtr,info; 31 /* check matrix is square */ 32 if (A_p->col_block_size !=A_p->row_block_size) { 33 Finley_ErrorCode = TYPE_ERROR; 34 sprintf(Finley_ErrorMsg, "Jacobi preconditioner square block size."); 35 return; 36 } 37 /* check matrix is square */ 38 if (n_block>3) { 39 Finley_ErrorCode = TYPE_ERROR; 40 sprintf(Finley_ErrorMsg, "Right now the Jacobi preconditioner supports block size less than 4 only"); 41 return; 42 } 43 /* allocate vector to hold main diagonal entries: */ 44 prec->values = MEMALLOC( ((size_t) n) * ((size_t) block_size),double); 45 prec->pivot = MEMALLOC( ((size_t) n) * ((size_t) n_block),int); 46 if (! (Finley_checkPtr(prec->values) || Finley_checkPtr(prec->pivot)) ) { 47 48 if (n_block==1) { 49 #pragma omp parallel for private(i, iPtr) schedule(static) 50 for (i = 0; i < A_p->pattern->n_ptr; i++) { 51 /* find main diagonal */ 52 for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) { 53 if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) { 54 if (ABS(A_p->val[iPtr])>0.) { 55 prec->values[i]=1./A_p->val[iPtr]; 56 } else { 57 prec->values[i]=1.; 58 } 59 break; 60 } 61 } 62 } 63 } else { 64 #pragma omp parallel for private(i, iPtr,info) schedule(static) 65 for (i = 0; i < A_p->pattern->n_ptr; i++) { 66 /* find main diagonal */ 67 for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) { 68 if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) { 69 info=Finley_Util_SmallMatLU(n_block,&A_p->val[iPtr*block_size],&prec->values[i*block_size],&prec->pivot[i*n_block]); 70 if (info>0) { 71 Finley_ErrorCode = ZERO_DIVISION_ERROR; 72 sprintf(Finley_ErrorMsg, "non-regular main diagonal block in row %d",i); 73 } 74 break; 75 } 76 } 77 } 78 } 79 } 80 #endif 81 } 82 83 /**************************************************************/ 84 85 /* inverse preconditioner solution using raw pointers */ 86 87 /* should be called within a parallel region */ 88 /* barrier synconization should be performed to make sure that the input vector available */ 89 90 void Finley_Solver_solveJacobi(Finley_SystemMatrix * A_p, double * x, double * b) { 91 #if ITERATIVE_SOLVER == NO_LIB 92 Finley_Solver_Preconditioner* prec=(Finley_Solver_Preconditioner*) (A_p->iterative); 93 int n_block=A_p->row_block_size; 94 int block_size=A_p->block_size; 95 int i; 96 97 if (n_block==1) { 98 #pragma omp for private(i) schedule(static) 99 for (i = 0; i < A_p->pattern->n_ptr; i++) { 100 x[i] = prec->values[i] * b[i]; 101 } 102 } else { 103 #pragma omp for private(i) schedule(static) 104 for (i = 0; i < A_p->pattern->n_ptr; i++) { 105 Finley_Util_SmallMatForwardBackwardSolve(n_block,1,&prec->values[i*block_size],&prec->pivot[i*n_block],&x[i*n_block],&b[i*n_block]); 106 } 107 } 108 return; 109 #endif 110 } 111 112 /* 113 * \$Log\$ 114 * Revision 1.2 2004/12/14 05:39:32 jgs 115 * *** empty log message *** 116 * 117 * Revision 1.1.1.1.2.3 2004/12/07 10:12:06 gross 118 * GMRES added 119 * 120 * Revision 1.1.1.1.2.2 2004/11/24 01:37:17 gross 121 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now 122 * 123 * Revision 1.1.1.1.2.1 2004/11/12 06:58:21 gross 124 * 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 125 * 126 * Revision 1.1.1.1 2004/10/26 06:53:58 jgs 127 * initial import of project esys2 128 * 129 * Revision 1.1 2004/07/02 04:21:14 gross 130 * Finley C code has been included 131 * 132 * 133 */

## Properties

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