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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2548 - (hide annotations)
Mon Jul 20 06:20:06 2009 UTC (10 years, 3 months ago) by jfenwick
Original Path: trunk/paso/src/GMRES2.c
File MIME type: text/plain
File size: 5384 byte(s)
Updating copyright notices
1 gross 1476
2     /*******************************************************
3     *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7 gross 1476 *
8 ksteube 1811 * 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 gross 1476 *
12 ksteube 1811 *******************************************************/
13 gross 1476
14 ksteube 1811
15 gross 1476 #include "Common.h"
16     #include "Solver.h"
17 gross 1639 #include "PasoUtil.h"
18 gross 1476 #ifdef _OPENMP
19     #include <omp.h>
20     #endif
21    
22 gross 1639 err_t Paso_Solver_GMRES2(
23 gross 1476 Paso_Function * F,
24     const double* f0,
25     const double* x0,
26     double * x,
27     dim_t *iter,
28     double* tolerance,
29     Paso_Performance* pp)
30     {
31     double static RENORMALIZATION_CONST=0.001;
32 gross 1639 dim_t l=(*iter)+1, iter_max=*iter,k=0,i,j;
33     const dim_t n=F->n;
34 gross 1476 double rel_tol=*tolerance, abs_tol, normf0, normv, normv2, hh, hr, nu, norm_of_residual=0.;
35     bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
36     double *h=NULL, **v=NULL, *c=NULL,*s=NULL,*g=NULL, *work=NULL;
37     err_t Status=SOLVER_NO_ERROR;
38    
39     /* Test the input parameters. */
40     if (n < 0 || iter_max<=0 || l<1 || rel_tol<0) {
41     return SOLVER_INPUT_ERROR;
42     }
43     Paso_zeroes(n,x);
44     /*
45     * allocate memory:
46     */
47     h=TMPMEMALLOC(l*l,double);
48 gross 1804 v=TMPMEMALLOC(l,double*);
49 gross 1476 c=TMPMEMALLOC(l,double);
50     s=TMPMEMALLOC(l,double);
51     g=TMPMEMALLOC(l,double);
52     work=TMPMEMALLOC(n,double);
53    
54     if (h==NULL || v ==NULL || c== NULL || s == NULL || g== NULL || work==NULL) {
55     Status=SOLVER_MEMORY_ERROR;
56     }
57     if ( Status ==SOLVER_NO_ERROR ) {
58 gross 1639 for (i=0;i<iter_max;i++) v[i]=NULL;
59 gross 1476 /*
60     * the show begins:
61     */
62     normf0=Paso_l2(n,f0,F->mpi_info);
63     k=0;
64     convergeFlag=(ABS(normf0)<=0);
65     if (! convergeFlag) {
66     abs_tol=rel_tol*normf0;
67 gross 1639 v[0]=TMPMEMALLOC(n,double);
68     if (v[0]==NULL) {
69     Status=SOLVER_MEMORY_ERROR;
70     } else {
71     Paso_zeroes(n,v[0]);
72     Paso_Update(n,1.,v[0],-1./normf0,f0);
73     g[0]=normf0;
74     }
75 gross 1798 while( (! (breakFlag || maxIterFlag || convergeFlag)) && (Status ==SOLVER_NO_ERROR) ) {
76 gross 1476 k++;
77 gross 1639 v[k]=TMPMEMALLOC(n,double);
78     if (v[k]==NULL) {
79     Status=SOLVER_MEMORY_ERROR;
80     } else {
81     /*
82     * call directional derivative function
83     */
84 gross 1804 Paso_FunctionDerivative(v[k],v[k-1],F,f0,x0,work,TRUE,pp);
85 gross 1639 normv=Paso_l2(n,v[k],F->mpi_info);
86     /*
87     * Modified Gram-Schmidt
88     */
89     for (j=0;j<k;j++){
90     hh=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);
91     Paso_Update(n,1.,v[k],(-hh),v[j]);
92     h[INDEX2(j,k-1,l)]=hh;
93     }
94     normv2=Paso_l2(n,v[k],F->mpi_info);
95     h[INDEX2(k,k-1,l)]=normv2;
96     /*
97     * reorthogonalize
98     */
99     if (! (normv + RENORMALIZATION_CONST*normv2 > normv)) {
100     for (j=0;j<k;j++){
101     hr=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);
102     h[INDEX2(j,k-1,l)]+=hr;
103     Paso_Update(n,1.,v[k],(-hr),v[j]);
104     }
105     normv2=Paso_l2(n,v[k],F->mpi_info);
106     h[INDEX2(k,k-1,l)]=normv2;
107     }
108     /*
109     * watch out for happy breakdown
110     */
111     if(normv2 > 0.) {
112     Paso_Update(n,1./normv2,v[k],0.,v[k]);
113     }
114     /*
115     * Form and store the information for the new Givens rotation
116     */
117     ApplyGivensRotations(k,&h[INDEX2(0,k-1,l)],c,s);
118     /*
119     * Don't divide by zero if solution has been found
120     */
121     g[k]=0;
122     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)]);
123     if (nu>0) {
124     c[k-1]=h[INDEX2(k-1,k-1,l)]/nu;
125 gross 1798 s[k-1]=-h[INDEX2(k,k-1,l)]/nu;
126     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)];
127 gross 1639 h[INDEX2(k,k-1,l)]=0;
128     ApplyGivensRotations(2,&(g[k-1]),&(c[k-1]),&(s[k-1]));
129     }
130     norm_of_residual=fabs(g[k]);
131     maxIterFlag = (k>=iter_max);
132     convergeFlag = (fabs(g[k]) <= abs_tol);
133     printf("GMRES2 step %d: error %e (tol=%e)\n",k,fabs(g[k]),abs_tol);
134     }
135 gross 1476 }
136     }
137     /*
138     * all done and ready for the forward substitution:
139     */
140 gross 1798 for (i=k-1;i>=0;--i) {
141     for (j=i+1;j<k;j++) {
142     g[i]-=h[INDEX2(i,j,l)]*g[j];
143 gross 1476 }
144     g[i]/=h[INDEX2(i,i,l)];
145 gross 1477 Paso_Update(n,1.,x,g[i],v[i]);
146 gross 1476 }
147     }
148     /*
149     * clean up:
150     */
151     if ( v !=NULL) {
152 gross 1639 for (i=0;i<iter_max;i++) TMPMEMFREE(v[i]);
153 gross 1476 }
154     TMPMEMFREE(h);
155     TMPMEMFREE(v);
156     TMPMEMFREE(c);
157     TMPMEMFREE(s);
158     TMPMEMFREE(g);
159     TMPMEMFREE(work);
160     *iter=k;
161     *tolerance=norm_of_residual;
162     return Status;
163     }
164    

  ViewVC Help
Powered by ViewVC 1.1.26