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

Revision 115 - (show annotations)
Fri Mar 4 07:12:47 2005 UTC (14 years, 1 month ago) by jgs
File MIME type: text/plain
File size: 6095 byte(s)
```*** empty log message ***

```
 1 /* \$Id\$ */ 2 3 /**************************************************************/ 4 5 /* Finley: jacobi preconditioner: */ 6 7 /**************************************************************/ 8 9 /* Copyrights by ACcESS Australia 2003, 2004, 2005 */ 10 /* Author: gross@access.edu.au */ 11 12 /**************************************************************/ 13 14 #include "Solver.h" 15 16 /**************************************************************/ 17 18 /* frees the Jacobi preconditioner */ 19 20 void Finley_Solver_Jacobi_free(Finley_Solver_Jacobi * in) { 21 if (in!=NULL) { 22 MEMFREE(in->values); 23 MEMFREE(in->pivot); 24 MEMFREE(in); 25 } 26 } 27 28 /**************************************************************/ 29 30 /* Jacobi precondioner set up */ 31 32 Finley_Solver_Jacobi* Finley_Solver_getJacobi(Finley_SystemMatrix * A_p) { 33 Finley_Solver_Jacobi* out=NULL; 34 int n = A_p->num_cols; 35 int n_block=A_p->row_block_size; 36 int block_size=A_p->block_size; 37 int i,iPtr; 38 double A11,A12,A13,A21,A22,A23,A31,A32,A33,D; 39 /* check matrix is square */ 40 if (A_p->col_block_size !=A_p->row_block_size) { 41 Finley_ErrorCode = TYPE_ERROR; 42 sprintf(Finley_ErrorMsg, "Jacobi preconditioner square block size."); 43 return NULL; 44 } 45 /* check matrix is square */ 46 if (n_block>3) { 47 Finley_ErrorCode = TYPE_ERROR; 48 sprintf(Finley_ErrorMsg, "Right now the Jacobi preconditioner supports block size less than 4 only"); 49 return NULL; 50 } 51 /* allocate vector to hold main diagonal entries: */ 52 out=MEMALLOC(1,Finley_Solver_Jacobi); 53 if (! Finley_checkPtr(out)) { 54 /* allocate vector to hold main diagonal entries: */ 55 out->n_block=n_block; 56 out->n=n; 57 out->values = MEMALLOC( ((size_t) n) * ((size_t) block_size),double); 58 out->pivot = NULL; /* later use */ 59 if (! (Finley_checkPtr(out->values))) { 60 if (n_block==1) { 61 #pragma omp parallel for private(i, iPtr) schedule(static) 62 for (i = 0; i < A_p->pattern->n_ptr; i++) { 63 out->values[i]=1.; 64 /* find main diagonal */ 65 for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) { 66 if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) { 67 if (ABS(A_p->val[iPtr])>0.) out->values[i]=1./A_p->val[iPtr]; 68 break; 69 } 70 } 71 } 72 } else if (n_block==2) { 73 #pragma omp parallel for private(i, iPtr, A11,A12,A21,A22,D) schedule(static) 74 for (i = 0; i < A_p->pattern->n_ptr; i++) { 75 out->values[i*4+0]= 1.; 76 out->values[i*4+1]= 0.; 77 out->values[i*4+2]= 0.; 78 out->values[i*4+3]= 1.; 79 /* find main diagonal */ 80 for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) { 81 if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) { 82 A11=A_p->val[iPtr*4]; 83 A12=A_p->val[iPtr*4+2]; 84 A21=A_p->val[iPtr*4+1]; 85 A22=A_p->val[iPtr*4+3]; 86 D = A11*A22-A12*A21; 87 if (ABS(D) > 0 ){ 88 D=1./D; 89 out->values[i*4]= A22*D; 90 out->values[i*4+1]=-A21*D; 91 out->values[i*4+2]=-A12*D; 92 out->values[i*4+3]= A11*D; 93 } 94 break; 95 } 96 } 97 } 98 } else if (n_block==3) { 99 #pragma omp parallel for private(i, iPtr,A11,A12,A13,A21,A22,A23,A31,A32,A33,D) schedule(static) 100 for (i = 0; i < A_p->pattern->n_ptr; i++) { 101 out->values[i*9 ]=1.; 102 out->values[i*9+1]=0.; 103 out->values[i*9+2]=0.; 104 out->values[i*9+3]=0.; 105 out->values[i*9+4]=1.; 106 out->values[i*9+5]=0.; 107 out->values[i*9+6]=0.; 108 out->values[i*9+7]=0.; 109 out->values[i*9+8]=1.; 110 /* find main diagonal */ 111 for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) { 112 if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) { 113 A11=A_p->val[iPtr*9 ]; 114 A21=A_p->val[iPtr*9+1]; 115 A31=A_p->val[iPtr*9+2]; 116 A12=A_p->val[iPtr*9+3]; 117 A22=A_p->val[iPtr*9+4]; 118 A32=A_p->val[iPtr*9+5]; 119 A13=A_p->val[iPtr*9+6]; 120 A23=A_p->val[iPtr*9+7]; 121 A33=A_p->val[iPtr*9+8]; 122 D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22); 123 if (ABS(D) > 0 ){ 124 D=1./D; 125 out->values[i*9 ]=(A22*A33-A23*A32)*D; 126 out->values[i*9+1]=(A31*A23-A21*A33)*D; 127 out->values[i*9+2]=(A21*A32-A31*A22)*D; 128 out->values[i*9+3]=(A13*A32-A12*A33)*D; 129 out->values[i*9+4]=(A11*A33-A31*A13)*D; 130 out->values[i*9+5]=(A12*A31-A11*A32)*D; 131 out->values[i*9+6]=(A12*A23-A13*A22)*D; 132 out->values[i*9+7]=(A13*A21-A11*A23)*D; 133 out->values[i*9+8]=(A11*A22-A12*A21)*D; 134 } 135 break; 136 } 137 } 138 } 139 } 140 } 141 } 142 if (Finley_ErrorCode == NO_ERROR) { 143 return out; 144 } else { 145 Finley_Solver_Jacobi_free(out); 146 return NULL; 147 } 148 } 149 /**************************************************************/ 150 151 /* applies Jacobi preconditioner */ 152 153 /* should be called within a parallel region */ 154 /* barrier synconization should be performed to make sure that the input vector available */ 155 156 void Finley_Solver_solveJacobi(Finley_Solver_Jacobi * prec, double * x, double * b) { 157 Finley_Solver_applyBlockDiagonalMatrix(prec->n_block,prec->n,prec->values,prec->pivot,x,b); 158 return; 159 }

## Properties

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

 ViewVC Help Powered by ViewVC 1.1.26