/[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 1388 - (hide annotations)
Fri Jan 11 07:45:58 2008 UTC (11 years, 10 months ago) by trankine
File MIME type: text/plain
File size: 6469 byte(s)
And get the *(&(*&(* name right
1 ksteube 1312
2 jgs 150 /* $Id$ */
3    
4 ksteube 1312 /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15 dhawcroft 631
16 jgs 150 /**************************************************************/
17    
18     /* Paso: jacobi preconditioner: */
19    
20     /**************************************************************/
21    
22     /* Copyrights by ACcESS Australia 2003, 2004, 2005 */
23     /* Author: gross@access.edu.au */
24    
25     /**************************************************************/
26    
27 gross 700 #include "Paso.h"
28 jgs 150 #include "Solver.h"
29    
30     /**************************************************************/
31    
32     /* frees the Jacobi preconditioner */
33    
34     void Paso_Solver_Jacobi_free(Paso_Solver_Jacobi * in) {
35     if (in!=NULL) {
36     MEMFREE(in->values);
37     MEMFREE(in->pivot);
38     MEMFREE(in);
39     }
40     }
41    
42     /**************************************************************/
43    
44     /* Jacobi precondioner set up */
45    
46 ksteube 1312 Paso_Solver_Jacobi* Paso_Solver_getJacobi(Paso_SparseMatrix * A_p) {
47 jgs 150 Paso_Solver_Jacobi* out=NULL;
48 ksteube 1312 dim_t n = A_p->numCols;
49 jgs 150 dim_t n_block=A_p->row_block_size;
50     dim_t block_size=A_p->block_size;
51     dim_t i;
52     index_t iPtr;
53     double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;
54     /* check matrix is square */
55     if (A_p->col_block_size !=A_p->row_block_size) {
56 gross 432 Paso_setError(TYPE_ERROR, "Paso_Solver_getJacobi: Jacobi preconditioner square block size.");
57 jgs 150 return NULL;
58     }
59     /* check matrix is square */
60     if (n_block>3) {
61 gross 432 Paso_setError(TYPE_ERROR, "Paso_Solver_getJacobi: Right now the Jacobi preconditioner supports block size less than 4 only");
62 jgs 150 return NULL;
63     }
64     /* allocate vector to hold main diagonal entries: */
65     out=MEMALLOC(1,Paso_Solver_Jacobi);
66     if (! Paso_checkPtr(out)) {
67     /* allocate vector to hold main diagonal entries: */
68     out->n_block=n_block;
69     out->n=n;
70     out->values = MEMALLOC( ((size_t) n) * ((size_t) block_size),double);
71     out->pivot = NULL; /* later use */
72     if (! (Paso_checkPtr(out->values))) {
73     if (n_block==1) {
74     #pragma omp parallel for private(i, iPtr) schedule(static)
75 ksteube 1312 for (i = 0; i < A_p->pattern->numOutput; i++) {
76 jgs 150 out->values[i]=1.;
77     /* find main diagonal */
78     for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {
79 gross 415 if (A_p->pattern->index[iPtr]==i) {
80 jgs 150 if (ABS(A_p->val[iPtr])>0.) out->values[i]=1./A_p->val[iPtr];
81     break;
82     }
83     }
84     }
85     } else if (n_block==2) {
86     #pragma omp parallel for private(i, iPtr, A11,A12,A21,A22,D) schedule(static)
87 ksteube 1312 for (i = 0; i < A_p->pattern->numOutput; i++) {
88 jgs 150 out->values[i*4+0]= 1.;
89     out->values[i*4+1]= 0.;
90     out->values[i*4+2]= 0.;
91     out->values[i*4+3]= 1.;
92     /* find main diagonal */
93     for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {
94 gross 415 if (A_p->pattern->index[iPtr]==i) {
95 jgs 150 A11=A_p->val[iPtr*4];
96     A12=A_p->val[iPtr*4+2];
97     A21=A_p->val[iPtr*4+1];
98     A22=A_p->val[iPtr*4+3];
99     D = A11*A22-A12*A21;
100     if (ABS(D) > 0 ){
101     D=1./D;
102     out->values[i*4]= A22*D;
103     out->values[i*4+1]=-A21*D;
104     out->values[i*4+2]=-A12*D;
105     out->values[i*4+3]= A11*D;
106     }
107     break;
108     }
109     }
110     }
111     } else if (n_block==3) {
112     #pragma omp parallel for private(i, iPtr,A11,A12,A13,A21,A22,A23,A31,A32,A33,D) schedule(static)
113 ksteube 1312 for (i = 0; i < A_p->pattern->numOutput; i++) {
114 jgs 150 out->values[i*9 ]=1.;
115     out->values[i*9+1]=0.;
116     out->values[i*9+2]=0.;
117     out->values[i*9+3]=0.;
118     out->values[i*9+4]=1.;
119     out->values[i*9+5]=0.;
120     out->values[i*9+6]=0.;
121     out->values[i*9+7]=0.;
122     out->values[i*9+8]=1.;
123     /* find main diagonal */
124     for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {
125 gross 415 if (A_p->pattern->index[iPtr]==i) {
126 jgs 150 A11=A_p->val[iPtr*9 ];
127     A21=A_p->val[iPtr*9+1];
128     A31=A_p->val[iPtr*9+2];
129     A12=A_p->val[iPtr*9+3];
130     A22=A_p->val[iPtr*9+4];
131     A32=A_p->val[iPtr*9+5];
132     A13=A_p->val[iPtr*9+6];
133     A23=A_p->val[iPtr*9+7];
134     A33=A_p->val[iPtr*9+8];
135     D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
136     if (ABS(D) > 0 ){
137     D=1./D;
138     out->values[i*9 ]=(A22*A33-A23*A32)*D;
139     out->values[i*9+1]=(A31*A23-A21*A33)*D;
140     out->values[i*9+2]=(A21*A32-A31*A22)*D;
141     out->values[i*9+3]=(A13*A32-A12*A33)*D;
142     out->values[i*9+4]=(A11*A33-A31*A13)*D;
143     out->values[i*9+5]=(A12*A31-A11*A32)*D;
144     out->values[i*9+6]=(A12*A23-A13*A22)*D;
145     out->values[i*9+7]=(A13*A21-A11*A23)*D;
146     out->values[i*9+8]=(A11*A22-A12*A21)*D;
147     }
148     break;
149     }
150     }
151     }
152     }
153     }
154     }
155     if (Paso_noError()) {
156     return out;
157     } else {
158     Paso_Solver_Jacobi_free(out);
159     return NULL;
160     }
161     }
162     /**************************************************************/
163    
164     /* applies Jacobi preconditioner */
165    
166     /* should be called within a parallel region */
167     /* barrier synconization should be performed to make sure that the input vector available */
168    
169     void Paso_Solver_solveJacobi(Paso_Solver_Jacobi * prec, double * x, double * b) {
170     Paso_Solver_applyBlockDiagonalMatrix(prec->n_block,prec->n,prec->values,prec->pivot,x,b);
171     return;
172     }
173    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26