/[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 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 5 months ago) by jgs
Original Path: trunk/esys2/paso/src/Solvers/Solver_jacobi.c
File MIME type: text/plain
File size: 6438 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

1 jgs 150 /* $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+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 (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

  ViewVC Help
Powered by ViewVC 1.1.26