# Contents of /trunk/paso/src/Solvers/Solver_jacobi.c

Revision 415 - (show annotations)
Wed Jan 4 05:37:33 2006 UTC (15 years, 7 months ago) by gross
File MIME type: text/plain
File size: 6399 byte(s)
```a better way of representing the matrix format type is introduced. this is needed for the Paradiso and UMFPACK interface
```
 1 /* \$Id\$ */ 2 3 /**************************************************************/ 4 5 /* Paso: jacobi preconditioner: */ 6 7 /**************************************************************/ 8 9 /* Copyrights by ACcESS Australia 2003, 2004, 2005 */ 10 /* Author: gross@access.edu.au */ 11 12 /**************************************************************/ 13 14 #include "Paso.h" 15 #include "Solver.h" 16 17 /**************************************************************/ 18 19 /* frees the Jacobi preconditioner */ 20 21 void Paso_Solver_Jacobi_free(Paso_Solver_Jacobi * in) { 22 if (in!=NULL) { 23 MEMFREE(in->values); 24 MEMFREE(in->pivot); 25 MEMFREE(in); 26 } 27 } 28 29 /**************************************************************/ 30 31 /* Jacobi precondioner set up */ 32 33 Paso_Solver_Jacobi* Paso_Solver_getJacobi(Paso_SystemMatrix * A_p) { 34 Paso_Solver_Jacobi* out=NULL; 35 dim_t n = A_p->num_cols; 36 dim_t n_block=A_p->row_block_size; 37 dim_t block_size=A_p->block_size; 38 dim_t i; 39 index_t iPtr; 40 double A11,A12,A13,A21,A22,A23,A31,A32,A33,D; 41 /* check matrix is square */ 42 if (A_p->col_block_size !=A_p->row_block_size) { 43 Paso_setError(TYPE_ERROR, "__FILE__: Jacobi preconditioner square block size."); 44 return NULL; 45 } 46 /* check matrix is square */ 47 if (n_block>3) { 48 Paso_setError(TYPE_ERROR, "__FILE__: 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,Paso_Solver_Jacobi); 53 if (! Paso_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 (! (Paso_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) { 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) { 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) { 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 (Paso_noError()) { 143 return out; 144 } else { 145 Paso_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 Paso_Solver_solveJacobi(Paso_Solver_Jacobi * prec, double * x, double * b) { 157 Paso_Solver_applyBlockDiagonalMatrix(prec->n_block,prec->n,prec->values,prec->pivot,x,b); 158 return; 159 } 160 161 /* 162 * \$Log\$ 163 * Revision 1.2 2005/09/15 03:44:40 jgs 164 * Merge of development branch dev-02 back to main trunk on 2005-09-15 165 * 166 * Revision 1.1.2.1 2005/09/05 06:29:50 gross 167 * These files have been extracted from finley to define a stand alone libray for iterative 168 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but 169 * has not been tested yet. 170 * 171 * 172 */

## Properties

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