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

Revision 633 - (show annotations)
Thu Mar 23 05:37:00 2006 UTC (13 years, 10 months ago) by dhawcroft
File MIME type: text/plain
File size: 7071 byte(s)

 1 /* \$Id\$ */ 2 3 4 /* 5 ******************************************************************************** 6 * Copyright 2006 by ACcESS MNRF * 7 * * 8 * http://www.access.edu.au * 9 * Primary Business: Queensland, Australia * 10 * Licensed under the Open Software License version 3.0 * 11 * http://www.opensource.org/licenses/osl-3.0.php * 12 ******************************************************************************** 13 */ 14 15 /**************************************************************/ 16 17 /* Paso: jacobi preconditioner: */ 18 19 /**************************************************************/ 20 21 /* Copyrights by ACcESS Australia 2003, 2004, 2005 */ 22 /* Author: gross@access.edu.au */ 23 24 /**************************************************************/ 25 26 #include "Paso.h" 27 #include "Solver.h" 28 29 /**************************************************************/ 30 31 /* frees the Jacobi preconditioner */ 32 33 void Paso_Solver_Jacobi_free(Paso_Solver_Jacobi * in) { 34 if (in!=NULL) { 35 MEMFREE(in->values); 36 MEMFREE(in->pivot); 37 MEMFREE(in); 38 } 39 } 40 41 /**************************************************************/ 42 43 /* Jacobi precondioner set up */ 44 45 Paso_Solver_Jacobi* Paso_Solver_getJacobi(Paso_SystemMatrix * A_p) { 46 Paso_Solver_Jacobi* out=NULL; 47 dim_t n = A_p->num_cols; 48 dim_t n_block=A_p->row_block_size; 49 dim_t block_size=A_p->block_size; 50 dim_t i; 51 index_t iPtr; 52 double A11,A12,A13,A21,A22,A23,A31,A32,A33,D; 53 /* check matrix is square */ 54 if (A_p->col_block_size !=A_p->row_block_size) { 55 Paso_setError(TYPE_ERROR, "Paso_Solver_getJacobi: Jacobi preconditioner square block size."); 56 return NULL; 57 } 58 /* check matrix is square */ 59 if (n_block>3) { 60 Paso_setError(TYPE_ERROR, "Paso_Solver_getJacobi: Right now the Jacobi preconditioner supports block size less than 4 only"); 61 return NULL; 62 } 63 /* allocate vector to hold main diagonal entries: */ 64 out=MEMALLOC(1,Paso_Solver_Jacobi); 65 if (! Paso_checkPtr(out)) { 66 /* allocate vector to hold main diagonal entries: */ 67 out->n_block=n_block; 68 out->n=n; 69 out->values = MEMALLOC( ((size_t) n) * ((size_t) block_size),double); 70 out->pivot = NULL; /* later use */ 71 if (! (Paso_checkPtr(out->values))) { 72 if (n_block==1) { 73 #pragma omp parallel for private(i, iPtr) schedule(static) 74 for (i = 0; i < A_p->pattern->n_ptr; i++) { 75 out->values[i]=1.; 76 /* find main diagonal */ 77 for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) { 78 if (A_p->pattern->index[iPtr]==i) { 79 if (ABS(A_p->val[iPtr])>0.) out->values[i]=1./A_p->val[iPtr]; 80 break; 81 } 82 } 83 } 84 } else if (n_block==2) { 85 #pragma omp parallel for private(i, iPtr, A11,A12,A21,A22,D) schedule(static) 86 for (i = 0; i < A_p->pattern->n_ptr; i++) { 87 out->values[i*4+0]= 1.; 88 out->values[i*4+1]= 0.; 89 out->values[i*4+2]= 0.; 90 out->values[i*4+3]= 1.; 91 /* find main diagonal */ 92 for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) { 93 if (A_p->pattern->index[iPtr]==i) { 94 A11=A_p->val[iPtr*4]; 95 A12=A_p->val[iPtr*4+2]; 96 A21=A_p->val[iPtr*4+1]; 97 A22=A_p->val[iPtr*4+3]; 98 D = A11*A22-A12*A21; 99 if (ABS(D) > 0 ){ 100 D=1./D; 101 out->values[i*4]= A22*D; 102 out->values[i*4+1]=-A21*D; 103 out->values[i*4+2]=-A12*D; 104 out->values[i*4+3]= A11*D; 105 } 106 break; 107 } 108 } 109 } 110 } else if (n_block==3) { 111 #pragma omp parallel for private(i, iPtr,A11,A12,A13,A21,A22,A23,A31,A32,A33,D) schedule(static) 112 for (i = 0; i < A_p->pattern->n_ptr; i++) { 113 out->values[i*9 ]=1.; 114 out->values[i*9+1]=0.; 115 out->values[i*9+2]=0.; 116 out->values[i*9+3]=0.; 117 out->values[i*9+4]=1.; 118 out->values[i*9+5]=0.; 119 out->values[i*9+6]=0.; 120 out->values[i*9+7]=0.; 121 out->values[i*9+8]=1.; 122 /* find main diagonal */ 123 for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) { 124 if (A_p->pattern->index[iPtr]==i) { 125 A11=A_p->val[iPtr*9 ]; 126 A21=A_p->val[iPtr*9+1]; 127 A31=A_p->val[iPtr*9+2]; 128 A12=A_p->val[iPtr*9+3]; 129 A22=A_p->val[iPtr*9+4]; 130 A32=A_p->val[iPtr*9+5]; 131 A13=A_p->val[iPtr*9+6]; 132 A23=A_p->val[iPtr*9+7]; 133 A33=A_p->val[iPtr*9+8]; 134 D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22); 135 if (ABS(D) > 0 ){ 136 D=1./D; 137 out->values[i*9 ]=(A22*A33-A23*A32)*D; 138 out->values[i*9+1]=(A31*A23-A21*A33)*D; 139 out->values[i*9+2]=(A21*A32-A31*A22)*D; 140 out->values[i*9+3]=(A13*A32-A12*A33)*D; 141 out->values[i*9+4]=(A11*A33-A31*A13)*D; 142 out->values[i*9+5]=(A12*A31-A11*A32)*D; 143 out->values[i*9+6]=(A12*A23-A13*A22)*D; 144 out->values[i*9+7]=(A13*A21-A11*A23)*D; 145 out->values[i*9+8]=(A11*A22-A12*A21)*D; 146 } 147 break; 148 } 149 } 150 } 151 } 152 } 153 } 154 if (Paso_noError()) { 155 return out; 156 } else { 157 Paso_Solver_Jacobi_free(out); 158 return NULL; 159 } 160 } 161 /**************************************************************/ 162 163 /* applies Jacobi preconditioner */ 164 165 /* should be called within a parallel region */ 166 /* barrier synconization should be performed to make sure that the input vector available */ 167 168 void Paso_Solver_solveJacobi(Paso_Solver_Jacobi * prec, double * x, double * b) { 169 Paso_Solver_applyBlockDiagonalMatrix(prec->n_block,prec->n,prec->values,prec->pivot,x,b); 170 return; 171 } 172 173 /* 174 * \$Log\$ 175 * Revision 1.2 2005/09/15 03:44:40 jgs 176 * Merge of development branch dev-02 back to main trunk on 2005-09-15 177 * 178 * Revision 1.1.2.1 2005/09/05 06:29:50 gross 179 * These files have been extracted from finley to define a stand alone libray for iterative 180 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but 181 * has not been tested yet. 182 * 183 * 184 */

## Properties

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