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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 700 - (hide annotations)
Thu Apr 6 00:13:40 2006 UTC (13 years, 10 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 jgs 150 /* $Id$ */
2    
3 dhawcroft 631
4     /*
5     ********************************************************************************
6 dhawcroft 633 * Copyright 2006 by ACcESS MNRF *
7 dhawcroft 631 * *
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 jgs 150 /**************************************************************/
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 gross 700 #include "Paso.h"
27 jgs 150 #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 gross 432 Paso_setError(TYPE_ERROR, "Paso_Solver_getJacobi: Jacobi preconditioner square block size.");
56 jgs 150 return NULL;
57     }
58     /* check matrix is square */
59     if (n_block>3) {
60 gross 432 Paso_setError(TYPE_ERROR, "Paso_Solver_getJacobi: Right now the Jacobi preconditioner supports block size less than 4 only");
61 jgs 150 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 gross 415 if (A_p->pattern->index[iPtr]==i) {
79 jgs 150 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 gross 415 if (A_p->pattern->index[iPtr]==i) {
94 jgs 150 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 gross 415 if (A_p->pattern->index[iPtr]==i) {
125 jgs 150 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