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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months, 1 week ago) by jfenwick
File size: 5117 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/paso/src/GMRES2.cpp:5567-5588 /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 /branches/trilinos_from_5897/paso/src/GMRES2.cpp:5898-6118 /release/3.0/paso/src/GMRES2.cpp:2591-2601 /release/4.0/paso/src/GMRES2.cpp:5380-5406 /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