/[escript]/trunk/paso/src/Solver_jacobi.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 700 - (show annotations)
Thu Apr 6 00:13:40 2006 UTC (13 years, 4 months ago) by gross
File MIME type: text/plain
File size: 7071 byte(s)
A few changes in the build mechanism and the file structure so scons can build release tar files:

  * paso/src/Solver has been moved to paso/src 
  * all test_.py are now run_.py files and are assumed to be passing python tests. they can run by 
    scons py_tests and are part of the release test set
  * escript/py_src/test_ are moved to escript/test/python and are installed in to the build directory 
    (rather then the PYTHONPATH).
  * all py files in test/python which don't start with run_ or test_ are now 'local_py_tests'. they are installed i
    by not run automatically.
  * CppUnitTest is now treated as a escript module (against previous decisions).
  * scons realse builds nor tar/zip files with relvant source code (src and tests in seperate files)

the python tests don't pass yet due to path problems.


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

  ViewVC Help
Powered by ViewVC 1.1.26