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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26