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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4873 - (hide annotations)
Wed Apr 16 06:38:51 2014 UTC (5 years, 6 months ago) by caltinay
File size: 4982 byte(s)
whitespace only changes.

1 gross 1476
2 jfenwick 3981 /*****************************************************************************
3 gross 1476 *
4 jfenwick 4657 * Copyright (c) 2003-2014 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 gross 1476 *
7 ksteube 1811 * Primary Business: Queensland, Australia
8     * Licensed under the Open Software License version 3.0
9     * http://www.opensource.org/licenses/osl-3.0.php
10 gross 1476 *
11 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 * Development 2012-2013 by School of Earth Sciences
13     * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 *
15     *****************************************************************************/
16 gross 1476
17 ksteube 1811
18 gross 1476 #include "Solver.h"
19 gross 1639 #include "PasoUtil.h"
20 gross 1476
21 caltinay 4856 namespace paso {
22    
23     err_t Solver_GMRES2(Function* F, const double* f0, const double* x0,
24     double* dx, dim_t* iter, double* tolerance,
25 caltinay 4873 Performance* pp)
26 gross 1476 {
27 caltinay 4856 static double RENORMALIZATION_CONST=0.001;
28     const dim_t l=(*iter)+1, iter_max=*iter;
29     dim_t k=0,i,j;
30     const dim_t n=F->getLen();
31     const double rel_tol=*tolerance;
32     double abs_tol, normf0, normv, normv2, hh, hr, nu, norm_of_residual=0.;
33 caltinay 4873 bool breakFlag=false, maxIterFlag=false, convergeFlag=false;
34 gross 1476
35 caltinay 4856 if (n < 0 || iter_max<=0 || l<1 || rel_tol<0) {
36     return SOLVER_INPUT_ERROR;
37     }
38 gross 1476
39 caltinay 4856 err_t Status=SOLVER_NO_ERROR;
40 gross 3005
41 caltinay 4856 double* h = new double[l*l];
42     double** v = new double*[l];
43     double* c = new double[l];
44     double* s = new double[l];
45     double* g = new double[l];
46     double* work = new double[n];
47 gross 3005
48 caltinay 4856 for (i=0; i<iter_max; i++)
49     v[i]=NULL;
50 gross 3005
51 caltinay 4856 util::zeroes(n,dx);
52    
53 caltinay 4873 /*
54 caltinay 4856 * the show begins:
55     */
56     normf0 = util::l2(n,f0,F->mpi_info);
57     k = 0;
58     convergeFlag = (ABS(normf0)<=0);
59     if (!convergeFlag) {
60     abs_tol = rel_tol*normf0;
61     printf("GMRES2 initial residual norm %e (rel. tol=%e)\n",normf0,rel_tol);
62     v[0] = new double[n];
63     util::zeroes(n, v[0]);
64     util::update(n, 1., v[0], -1./normf0, f0); // v = -1./normf0*f0
65     g[0] = normf0;
66     while (!breakFlag && !maxIterFlag && !convergeFlag && Status==SOLVER_NO_ERROR) {
67     k++;
68     v[k]=new double[n];
69     /*
70     * call directional derivative function
71     */
72     F->derivative(v[k], v[k-1], f0, x0, work, pp);
73     normv=util::l2(n,v[k],F->mpi_info);
74     /*
75     * Modified Gram-Schmidt
76     */
77     for (j=0; j<k; j++){
78     hh = util::innerProduct(n,v[j],v[k],F->mpi_info);
79     util::update(n,1.,v[k],(-hh),v[j]); // v[k]-hh*v[j]
80 caltinay 4873 h[INDEX2(j,k-1,l)]=hh;
81 caltinay 4856 //printf("%d : %d = %e\n",k,j,hh);
82     }
83     normv2=util::l2(n,v[k],F->mpi_info);
84     h[INDEX2(k,k-1,l)]=normv2;
85     /*
86     * reorthogonalize
87     */
88     if (!(normv + RENORMALIZATION_CONST*normv2 > normv)) {
89     // printf("GMRES2: renormalization!");
90     for (j=0; j<k; j++) {
91     hr = util::innerProduct(n,v[j],v[k],F->mpi_info);
92    
93     h[INDEX2(j,k-1,l)]+=hr;
94     util::update(n,1.,v[k],(-hr),v[j]);
95     }
96     normv2=util::l2(n,v[k],F->mpi_info);
97     h[INDEX2(k,k-1,l)]=normv2;
98     }
99 caltinay 4873 /*
100 caltinay 4856 * watch out for happy breakdown
101     */
102     if (normv2 > 0.) {
103     util::update(n,1./normv2,v[k],0.,v[k]); /* normalize v[k] */
104 caltinay 4873 }
105 caltinay 4856 /*
106     * Form and store the information for the new Givens rotation
107     */
108     util::applyGivensRotations(k,&h[INDEX2(0,k-1,l)],c,s);
109    
110     /*
111     * Don't divide by zero if solution has been found
112     */
113     g[k] = 0;
114     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)]);
115     if (nu > 0) {
116     c[k-1]= h[INDEX2(k-1,k-1,l)]/nu;
117     s[k-1]=-h[INDEX2(k,k-1,l)]/nu;
118     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)];
119     h[INDEX2(k,k-1,l)]=0;
120     util::applyGivensRotations(2,&(g[k-1]),&(c[k-1]),&(s[k-1]));
121     }
122     norm_of_residual = fabs(g[k]);
123     maxIterFlag = (k >= iter_max);
124     convergeFlag = (norm_of_residual <= abs_tol);
125     printf("GMRES2 step %d: residual %e (abs. tol=%e)\n", k, fabs(g[k]), abs_tol);
126     }
127     }
128    
129     // all done and ready for the forward substitution:
130     for (i=k-1;i>=0;--i) {
131     for (j=i+1;j<k;j++) {
132     g[i]-=h[INDEX2(i,j,l)]*g[j];
133     }
134     g[i] /= h[INDEX2(i,i,l)];
135     util::update(n, 1., dx, g[i], v[i]); // dx = dx+g[i]*v[i]
136     }
137     if ( v != NULL) {
138     for (i=0; i<iter_max; i++)
139     delete[] v[i];
140     }
141     delete[] h;
142     delete[] v;
143     delete[] c;
144     delete[] s;
145     delete[] g;
146     delete[] work;
147     *iter=k;
148     *tolerance=norm_of_residual;
149     return Status;
150 gross 1476 }
151    
152 caltinay 4856 } // namespace paso
153    

Properties

Name Value
svn:mergeinfo /branches/amg_from_3530/paso/src/GMRES2.cpp:3531-3826 /branches/lapack2681/paso/src/GMRES2.cpp:2682-2741 /branches/pasowrap/paso/src/GMRES2.cpp:3661-3674 /branches/py3_attempt2/paso/src/GMRES2.cpp:3871-3891 /branches/restext/paso/src/GMRES2.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/GMRES2.cpp:3669-3791 /branches/stage3.0/paso/src/GMRES2.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/GMRES2.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/GMRES2.cpp:3517-3974 /release/3.0/paso/src/GMRES2.cpp:2591-2601 /trunk/paso/src/GMRES2.cpp:4257-4344 /trunk/ripley/test/python/paso/src/GMRES2.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26