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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6000 - (show annotations)
Tue Mar 1 00:24:43 2016 UTC (3 years, 1 month ago) by caltinay
File size: 5126 byte(s)
a few more include rearrangements.

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

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 /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