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

Annotation of /trunk/paso/src/GMRES2.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1639 - (hide annotations)
Mon Jul 14 08:55:25 2008 UTC (11 years, 3 months ago) by gross
Original Path: trunk/paso/src/GMRES2.c
File MIME type: text/plain
File size: 6578 byte(s)


1 gross 1476 /* $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 gross 1639 * GMRES2 solves the non-linear system f0+J_0*d=0
18 gross 1476 * 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 gross 1639 #include "PasoUtil.h"
53 gross 1476 #ifdef _OPENMP
54     #include <omp.h>
55     #endif
56    
57 gross 1639 err_t Paso_Solver_GMRES2(
58 gross 1476 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 gross 1639 dim_t l=(*iter)+1, iter_max=*iter,k=0,i,j;
68     const dim_t n=F->n;
69 gross 1476 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(iter_max,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 gross 1639 for (i=0;i<iter_max;i++) v[i]=NULL;
94 gross 1476 /*
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 gross 1639 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 gross 1476 k++;
112 gross 1639 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);
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[k,k-1]/nu;
161     h[INDEX2(k-1,k-1,l)]=c[k-1]*h[INDEX2(k-1,k-1,l)]-s[k-1]*h[k,k-1];
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 gross 1476 }
171     }
172     /*
173     * all done and ready for the forward substitution:
174     */
175    
176     for (i=0;i<k;i++) {
177     for (j=0;j<i;j++) {
178     g[i]-=h[INDEX2(j,i,l)]*g[j];
179     }
180     g[i]/=h[INDEX2(i,i,l)];
181 gross 1477 Paso_Update(n,1.,x,g[i],v[i]);
182 gross 1476 }
183     }
184     /*
185     * clean up:
186     */
187     if ( v !=NULL) {
188 gross 1639 for (i=0;i<iter_max;i++) TMPMEMFREE(v[i]);
189 gross 1476 }
190     TMPMEMFREE(h);
191     TMPMEMFREE(v);
192     TMPMEMFREE(c);
193     TMPMEMFREE(s);
194     TMPMEMFREE(g);
195     TMPMEMFREE(work);
196     *iter=k;
197     *tolerance=norm_of_residual;
198     return Status;
199     }
200    

  ViewVC Help
Powered by ViewVC 1.1.26