/[escript]/trunk/esys2/finley/src/finleyC/Solvers/Solver_jacobi.c
ViewVC logotype

Contents of /trunk/esys2/finley/src/finleyC/Solvers/Solver_jacobi.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 97 - (show annotations)
Tue Dec 14 05:39:33 2004 UTC (15 years, 1 month ago) by jgs
File MIME type: text/plain
File size: 4654 byte(s)
*** empty log message ***

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Finley: jacobi preconditioner: */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia 2003/04 */
10 /* Author: gross@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Finley.h"
15 #include "System.h"
16 #include "Solver.h"
17 #include "Util.h"
18 #include "Common.h"
19
20 /**************************************************************/
21
22 /* inverse preconditioner setup */
23
24 void Finley_Solver_setJacobi(Finley_SystemMatrix * A_p) {
25 #if ITERATIVE_SOLVER == NO_LIB
26 Finley_Solver_Preconditioner* prec=(Finley_Solver_Preconditioner*) (A_p->iterative);
27 int n = A_p->num_cols;
28 int n_block=A_p->row_block_size;
29 int block_size=A_p->block_size;
30 int i,iPtr,info;
31 /* check matrix is square */
32 if (A_p->col_block_size !=A_p->row_block_size) {
33 Finley_ErrorCode = TYPE_ERROR;
34 sprintf(Finley_ErrorMsg, "Jacobi preconditioner square block size.");
35 return;
36 }
37 /* check matrix is square */
38 if (n_block>3) {
39 Finley_ErrorCode = TYPE_ERROR;
40 sprintf(Finley_ErrorMsg, "Right now the Jacobi preconditioner supports block size less than 4 only");
41 return;
42 }
43 /* allocate vector to hold main diagonal entries: */
44 prec->values = MEMALLOC( ((size_t) n) * ((size_t) block_size),double);
45 prec->pivot = MEMALLOC( ((size_t) n) * ((size_t) n_block),int);
46 if (! (Finley_checkPtr(prec->values) || Finley_checkPtr(prec->pivot)) ) {
47
48 if (n_block==1) {
49 #pragma omp parallel for private(i, iPtr) schedule(static)
50 for (i = 0; i < A_p->pattern->n_ptr; i++) {
51 /* find main diagonal */
52 for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {
53 if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) {
54 if (ABS(A_p->val[iPtr])>0.) {
55 prec->values[i]=1./A_p->val[iPtr];
56 } else {
57 prec->values[i]=1.;
58 }
59 break;
60 }
61 }
62 }
63 } else {
64 #pragma omp parallel for private(i, iPtr,info) schedule(static)
65 for (i = 0; i < A_p->pattern->n_ptr; i++) {
66 /* find main diagonal */
67 for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {
68 if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) {
69 info=Finley_Util_SmallMatLU(n_block,&A_p->val[iPtr*block_size],&prec->values[i*block_size],&prec->pivot[i*n_block]);
70 if (info>0) {
71 Finley_ErrorCode = ZERO_DIVISION_ERROR;
72 sprintf(Finley_ErrorMsg, "non-regular main diagonal block in row %d",i);
73 }
74 break;
75 }
76 }
77 }
78 }
79 }
80 #endif
81 }
82
83 /**************************************************************/
84
85 /* inverse preconditioner solution using raw pointers */
86
87 /* should be called within a parallel region */
88 /* barrier synconization should be performed to make sure that the input vector available */
89
90 void Finley_Solver_solveJacobi(Finley_SystemMatrix * A_p, double * x, double * b) {
91 #if ITERATIVE_SOLVER == NO_LIB
92 Finley_Solver_Preconditioner* prec=(Finley_Solver_Preconditioner*) (A_p->iterative);
93 int n_block=A_p->row_block_size;
94 int block_size=A_p->block_size;
95 int i;
96
97 if (n_block==1) {
98 #pragma omp for private(i) schedule(static)
99 for (i = 0; i < A_p->pattern->n_ptr; i++) {
100 x[i] = prec->values[i] * b[i];
101 }
102 } else {
103 #pragma omp for private(i) schedule(static)
104 for (i = 0; i < A_p->pattern->n_ptr; i++) {
105 Finley_Util_SmallMatForwardBackwardSolve(n_block,1,&prec->values[i*block_size],&prec->pivot[i*n_block],&x[i*n_block],&b[i*n_block]);
106 }
107 }
108 return;
109 #endif
110 }
111
112 /*
113 * $Log$
114 * Revision 1.2 2004/12/14 05:39:32 jgs
115 * *** empty log message ***
116 *
117 * Revision 1.1.1.1.2.3 2004/12/07 10:12:06 gross
118 * GMRES added
119 *
120 * Revision 1.1.1.1.2.2 2004/11/24 01:37:17 gross
121 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
122 *
123 * Revision 1.1.1.1.2.1 2004/11/12 06:58:21 gross
124 * a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry
125 *
126 * Revision 1.1.1.1 2004/10/26 06:53:58 jgs
127 * initial import of project esys2
128 *
129 * Revision 1.1 2004/07/02 04:21:14 gross
130 * Finley C code has been included
131 *
132 *
133 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26