/[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 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years, 1 month ago) by ksteube
File MIME type: text/plain
File size: 6469 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

1
2 /* $Id$ */
3
4 /*******************************************************
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
16 /**************************************************************/
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 #include "Paso.h"
28 #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 Paso_Solver_Jacobi* Paso_Solver_getJacobi(Paso_SparseMatrix * A_p) {
47 Paso_Solver_Jacobi* out=NULL;
48 dim_t n = A_p->numCols;
49 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 Paso_setError(TYPE_ERROR, "Paso_Solver_getJacobi: Jacobi preconditioner square block size.");
57 return NULL;
58 }
59 /* check matrix is square */
60 if (n_block>3) {
61 Paso_setError(TYPE_ERROR, "Paso_Solver_getJacobi: Right now the Jacobi preconditioner supports block size less than 4 only");
62 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 for (i = 0; i < A_p->pattern->numOutput; i++) {
76 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 if (A_p->pattern->index[iPtr]==i) {
80 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 for (i = 0; i < A_p->pattern->numOutput; i++) {
88 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 if (A_p->pattern->index[iPtr]==i) {
95 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 for (i = 0; i < A_p->pattern->numOutput; i++) {
114 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 if (A_p->pattern->index[iPtr]==i) {
126 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