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

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 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 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 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; impi_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; jmpi_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; jmpi_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

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