/[escript]/trunk/esys2/finley/src/finleyC/Solvers/Solver_jacobi.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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