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

Contents of /trunk/paso/src/GMRES2.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1804 - (show annotations)
Wed Sep 24 07:52:19 2008 UTC (11 years, 2 months ago) by gross
File MIME type: text/plain
File size: 6600 byte(s)
a robister version of the upwinding scheme
1 /* $Id:$ */
2
3 /*******************************************************
4 *
5 * Copyright 2008 by University of Queensland
6 *
7 * http://esscc.uq.edu.au
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13 /*
14 * Purpose
15 * =======
16 *
17 * GMRES2 solves the non-linear system f0+J_0*d=0
18 * where f0=F(x0), J_0 is the jacobian of F at x0.
19 *
20 * Convergence test: norm(f0+J_0*d)<=tolerance*norm(f0)
21 *
22 * Arguments
23 * =========
24 *
25 * Paso_Function * F evaluation of F (including any preconditioner)
26 *
27 * x0 (input) current point
28 *
29 * f0 (input) function F at current point x0
30 *
31 * d (output) solution of f0+J0*d=0 with accuracy tolerance
32 *
33 * iter (input/output)
34 * On input, the maximum num_iterations to be performed.
35 * On output, actual number of num_iterations performed.
36 * tolerance (input/output)
37 * On input, the allowable convergence measure for norm(f0+J0*d)/norm(f0)
38 * On output, the final value for norm(f0+J0*d)
39 * return value
40 *
41 * =SOLVER_NO_ERROR: Successful exit. approximate solution returned.
42 * =SOLVER_MAXNUM_ITER_REACHED
43 * =SOLVER_INPUT_ERROR Illegal parameter:
44 * =SOLVER_BREAKDOWN: bad luck!
45 * =SOLVER_MEMORY_ERROR : no memory available
46 *
47 * ==============================================================
48 */
49
50 #include "Common.h"
51 #include "Solver.h"
52 #include "PasoUtil.h"
53 #ifdef _OPENMP
54 #include <omp.h>
55 #endif
56
57 err_t Paso_Solver_GMRES2(
58 Paso_Function * F,
59 const double* f0,
60 const double* x0,
61 double * x,
62 dim_t *iter,
63 double* tolerance,
64 Paso_Performance* pp)
65 {
66 double static RENORMALIZATION_CONST=0.001;
67 dim_t l=(*iter)+1, iter_max=*iter,k=0,i,j;
68 const dim_t n=F->n;
69 double rel_tol=*tolerance, abs_tol, normf0, normv, normv2, hh, hr, nu, norm_of_residual=0.;
70 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
71 double *h=NULL, **v=NULL, *c=NULL,*s=NULL,*g=NULL, *work=NULL;
72 err_t Status=SOLVER_NO_ERROR;
73
74 /* Test the input parameters. */
75 if (n < 0 || iter_max<=0 || l<1 || rel_tol<0) {
76 return SOLVER_INPUT_ERROR;
77 }
78 Paso_zeroes(n,x);
79 /*
80 * allocate memory:
81 */
82 h=TMPMEMALLOC(l*l,double);
83 v=TMPMEMALLOC(l,double*);
84 c=TMPMEMALLOC(l,double);
85 s=TMPMEMALLOC(l,double);
86 g=TMPMEMALLOC(l,double);
87 work=TMPMEMALLOC(n,double);
88
89 if (h==NULL || v ==NULL || c== NULL || s == NULL || g== NULL || work==NULL) {
90 Status=SOLVER_MEMORY_ERROR;
91 }
92 if ( Status ==SOLVER_NO_ERROR ) {
93 for (i=0;i<iter_max;i++) v[i]=NULL;
94 /*
95 * the show begins:
96 */
97 normf0=Paso_l2(n,f0,F->mpi_info);
98 k=0;
99 convergeFlag=(ABS(normf0)<=0);
100 if (! convergeFlag) {
101 abs_tol=rel_tol*normf0;
102 v[0]=TMPMEMALLOC(n,double);
103 if (v[0]==NULL) {
104 Status=SOLVER_MEMORY_ERROR;
105 } else {
106 Paso_zeroes(n,v[0]);
107 Paso_Update(n,1.,v[0],-1./normf0,f0);
108 g[0]=normf0;
109 }
110 while( (! (breakFlag || maxIterFlag || convergeFlag)) && (Status ==SOLVER_NO_ERROR) ) {
111 k++;
112 v[k]=TMPMEMALLOC(n,double);
113 if (v[k]==NULL) {
114 Status=SOLVER_MEMORY_ERROR;
115 } else {
116 /*
117 * call directional derivative function
118 */
119 Paso_FunctionDerivative(v[k],v[k-1],F,f0,x0,work,TRUE,pp);
120 normv=Paso_l2(n,v[k],F->mpi_info);
121 /*
122 * Modified Gram-Schmidt
123 */
124 for (j=0;j<k;j++){
125 hh=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);
126 Paso_Update(n,1.,v[k],(-hh),v[j]);
127 h[INDEX2(j,k-1,l)]=hh;
128 }
129 normv2=Paso_l2(n,v[k],F->mpi_info);
130 h[INDEX2(k,k-1,l)]=normv2;
131 /*
132 * reorthogonalize
133 */
134 if (! (normv + RENORMALIZATION_CONST*normv2 > normv)) {
135 for (j=0;j<k;j++){
136 hr=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);
137 h[INDEX2(j,k-1,l)]+=hr;
138 Paso_Update(n,1.,v[k],(-hr),v[j]);
139 }
140 normv2=Paso_l2(n,v[k],F->mpi_info);
141 h[INDEX2(k,k-1,l)]=normv2;
142 }
143 /*
144 * watch out for happy breakdown
145 */
146 if(normv2 > 0.) {
147 Paso_Update(n,1./normv2,v[k],0.,v[k]);
148 }
149 /*
150 * Form and store the information for the new Givens rotation
151 */
152 ApplyGivensRotations(k,&h[INDEX2(0,k-1,l)],c,s);
153 /*
154 * Don't divide by zero if solution has been found
155 */
156 g[k]=0;
157 nu=sqrt(h[INDEX2(k-1,k-1,l)]*h[INDEX2(k-1,k-1,l)]+h[INDEX2(k,k-1,l)]*h[INDEX2(k,k-1,l)]);
158 if (nu>0) {
159 c[k-1]=h[INDEX2(k-1,k-1,l)]/nu;
160 s[k-1]=-h[INDEX2(k,k-1,l)]/nu;
161 h[INDEX2(k-1,k-1,l)]=c[k-1]*h[INDEX2(k-1,k-1,l)]-s[k-1]*h[INDEX2(k,k-1,l)];
162 h[INDEX2(k,k-1,l)]=0;
163 ApplyGivensRotations(2,&(g[k-1]),&(c[k-1]),&(s[k-1]));
164 }
165 norm_of_residual=fabs(g[k]);
166 maxIterFlag = (k>=iter_max);
167 convergeFlag = (fabs(g[k]) <= abs_tol);
168 printf("GMRES2 step %d: error %e (tol=%e)\n",k,fabs(g[k]),abs_tol);
169 }
170 }
171 }
172 /*
173 * all done and ready for the forward substitution:
174 */
175 for (i=k-1;i>=0;--i) {
176 for (j=i+1;j<k;j++) {
177 g[i]-=h[INDEX2(i,j,l)]*g[j];
178 }
179 g[i]/=h[INDEX2(i,i,l)];
180 Paso_Update(n,1.,x,g[i],v[i]);
181 }
182 }
183 /*
184 * clean up:
185 */
186 if ( v !=NULL) {
187 for (i=0;i<iter_max;i++) TMPMEMFREE(v[i]);
188 }
189 TMPMEMFREE(h);
190 TMPMEMFREE(v);
191 TMPMEMFREE(c);
192 TMPMEMFREE(s);
193 TMPMEMFREE(g);
194 TMPMEMFREE(work);
195 *iter=k;
196 *tolerance=norm_of_residual;
197 return Status;
198 }
199

  ViewVC Help
Powered by ViewVC 1.1.26