3D EM based on Schur decomposition
Вы не можете выбрать более 25 тем Темы должны начинаться с буквы или цифры, могут содержать дефисы(-) и должны содержать не более 35 символов.

bicgstab.h 8.0KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274
  1. // ===========================================================================
  2. //
  3. // Filename: bicgstab.h
  4. //
  5. // Description:
  6. //
  7. // Version: 0.0
  8. // Created: 10/27/2009 03:15:06 PM
  9. // Revision: none
  10. // Compiler: g++ (c++)
  11. //
  12. // Author: Trevor Irons (ti)
  13. //
  14. // Organisation: Colorado School of Mines (CSM)
  15. // United States Geological Survey (USGS)
  16. //
  17. // Email: tirons@mines.edu, tirons@usgs.gov
  18. //
  19. // This program is free software: you can redistribute it and/or modify
  20. // it under the terms of the GNU General Public License as published by
  21. // the Free Software Foundation, either version 3 of the License, or
  22. // (at your option) any later version.
  23. //
  24. // This program is distributed in the hope that it will be useful,
  25. // but WITHOUT ANY WARRANTY; without even the implied warranty of
  26. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  27. // GNU General Public License for more details.
  28. //
  29. // You should have received a copy of the GNU General Public License
  30. // along with this program. If not, see <http://www.gnu.org/licenses/>.
  31. //
  32. // ===========================================================================
  33. #pragma once
  34. #ifndef BICGSTAB_INC
  35. #define BICGSTAB_INC
  36. #include <Eigen/Core>
  37. #include <Eigen/Sparse>
  38. #ifdef CHOLMODPRECONDITION
  39. #include <Eigen/CholmodSupport>
  40. #endif // CHOLMODPRECONDITION
  41. //#include "unsupported/Eigen/IterativeSolvers"
  42. //#include <unsupported/Eigen/SuperLUSupport>
  43. #include <iostream>
  44. #include <string>
  45. #include <complex>
  46. #include <fstream>
  47. #include <LemmaCore>
  48. #include "timer.h"
  49. using namespace Eigen;
  50. using namespace Lemma;
  51. typedef Eigen::SparseMatrix<Complex> SparseMat;
  52. // Computes implicit calculation
  53. template < typename Solver >
  54. inline VectorXcr implicitDCInvBPhi3 (
  55. const SparseMat& D,
  56. const Solver& solver,
  57. const Ref<VectorXcr const> ioms,
  58. const Ref<VectorXcr const> Phi,
  59. Real& tol, // not used
  60. int& max_it // not used
  61. ) {
  62. VectorXcr b = (ioms).asDiagonal() * (D.transpose()*Phi);
  63. VectorXcr y = solver.solve(b);
  64. //max_it = solver.iterations(); // actualy no need to pass this
  65. return D*y;
  66. //return y;
  67. }
  68. // // Simple extraction of indices in idx into reduceed array x1
  69. // void vmap( const Ref<VectorXcr const> x0, Ref<VectorXcr> x1, const std::vector<int>& idx ) {
  70. // for (unsigned int ii=0; ii<idx.size(); ++ii) {
  71. // x1(ii) = x0(idx[ii]);
  72. // }
  73. // }
  74. // Simple extraction of indices in idx into reduceed array x1
  75. VectorXcr vmap( const Ref<VectorXcr const> x0, const std::vector<int>& idx ) {
  76. VectorXcr x1 = VectorXcr::Zero( idx.size() );
  77. for (unsigned int ii=0; ii<idx.size(); ++ii) {
  78. x1(ii) = x0(idx[ii]);
  79. }
  80. return x1;
  81. }
  82. // reverse of above
  83. void ivmap( Ref<VectorXcr > x0, const Ref<VectorXcr const> x1, const std::vector<int>& idx ) {
  84. for (unsigned int ii=0; ii<idx.size(); ++ii) {
  85. x0(idx[ii]) = x1(ii);
  86. }
  87. }
  88. // On Input
  89. // A = Matrix
  90. // B = Right hand side
  91. // X = initial guess, and solution
  92. // maxit = maximum Number of iterations
  93. // tol = error tolerance
  94. // On Output
  95. // X real solution vector
  96. // errorn = Real error norm
  97. template < typename CSolver >
  98. int implicitbicgstab(//const SparseMat& D,
  99. //const SparseMat& C,
  100. const Ref< Eigen::SparseMatrix<Complex> const > D,
  101. const std::vector<int>& idx,
  102. const Ref< VectorXcr const > ioms,
  103. const Ref< VectorXcr const > rhs,
  104. Ref <VectorXcr> phi,
  105. CSolver& solver,
  106. int &max_it, const Real &tol, Real &errorn, int &iter_done, ofstream& logio) {
  107. logio << "using the preconditioned implicit solver" << std::endl;
  108. Complex omega, rho, rho_1, alpha, beta;
  109. Real tol2;
  110. int iter, max_it2, max_it1;
  111. // Look at reduced problem
  112. VectorXcr rhs2 = vmap(rhs, idx);
  113. VectorXcr phi2 = vmap(phi, idx);
  114. // Determine size of system and init vectors
  115. int n = idx.size(); // was phi.size();
  116. std::cout << "BiCGStab SIZES " << n << "\t" << phi.size() << "\t" << ioms.size() << std::endl;
  117. VectorXcr r(n);
  118. VectorXcr r_tld(n);
  119. VectorXcr p(n);
  120. VectorXcr s(n);
  121. VectorXcr v = VectorXcr::Zero(ioms.size());
  122. VectorXcr t = VectorXcr::Zero(ioms.size());
  123. // VectorXcr vm1 = VectorXcr::Zero(ioms.size());
  124. // VectorXcr tm1 = VectorXcr::Zero(ioms.size());
  125. // TODO, refigure for implicit large system
  126. // std::cout << "Start BiCGStab, memory needed: "
  127. // << (sizeof(Complex)*(9+2)*n/(1024.*1024*1024)) << " [Gb]\n";
  128. // Initialise
  129. iter_done = 0;
  130. Real eps = 1e-100;
  131. Real bnrm2 = rhs.norm();
  132. if (bnrm2 == 0) {
  133. phi.setConstant(0.0);
  134. errorn = 0;
  135. std::cerr << "Trivial case of Ax = b, where b is 0\n";
  136. return (0);
  137. }
  138. // If there is an initial guess
  139. if ( phi.norm() ) {
  140. tol2 = tol;
  141. max_it2 = 50000;
  142. //r = rhs - implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it2);
  143. //r = rhs - implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it2);
  144. r = rhs2 - vmap( implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it2), idx );
  145. } else {
  146. r = vmap(rhs, idx);
  147. }
  148. jsw_timer timer;
  149. errorn = r.norm() / bnrm2;
  150. omega = 1.;
  151. r_tld = r;
  152. Real errornold = 1e14;
  153. // Get down to business
  154. for (iter=0; iter<max_it; ++iter) {
  155. timer.begin();
  156. rho = r_tld.dot(r);
  157. if (abs(rho) < eps) {
  158. ivmap( phi, phi2, idx );
  159. return (0);
  160. }
  161. if (iter > 0) {
  162. beta = (rho/rho_1) * (alpha/omega);
  163. p = r.array() + beta*(p.array()-omega*v.array()).array();
  164. } else {
  165. p = r;
  166. }
  167. tol2 = tol;
  168. max_it2 = 500000;
  169. ivmap(phi, p, idx);
  170. v = vmap(implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it2), idx);
  171. alpha = rho / r_tld.dot(v);
  172. s = r.array() - alpha*v.array();
  173. errorn = s.norm()/bnrm2;
  174. if (errorn < tol && iter > 1) {
  175. phi2.array() += alpha*p.array();
  176. ivmap( phi, phi2, idx );
  177. return (0);
  178. }
  179. tol2 = tol;
  180. max_it1 = 500000;
  181. //t = implicitDCInvBPhi2(D, C, ioms, solver, s, tol2, max_it1);
  182. //t = implicitDCInvBPhi3(D, solver, ioms, s, tol2, max_it1);
  183. ivmap(phi, s, idx);
  184. t = vmap(implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it1), idx);
  185. omega = t.dot(s) / t.dot(t);
  186. r = s.array() - omega*t.array();
  187. errorn = r.norm() / bnrm2;
  188. iter_done = iter;
  189. if (errorn <= tol) {
  190. ivmap( phi, phi2, idx );
  191. return (0);
  192. }
  193. if (abs(omega) < eps) {
  194. ivmap( phi, phi2, idx );
  195. return (0);
  196. }
  197. rho_1 = rho;
  198. logio << "iteration " << std::setw(3) << iter
  199. << " errorn " << std::setw(6) << std::setprecision(4) << std::scientific << errorn
  200. //<< "\timplicit iterations " << std::setw(5) << max_it1+max_it2
  201. << " time " << std::setw(6) << std::fixed << std::setprecision(2) << timer.end()/60. << " [m] "
  202. << max_it2+max_it2 << " iterations " << std::endl;
  203. // Check to see how progress is going
  204. if (errornold - errorn < 0) {
  205. logio << "Irregular non-monotonic (negative) convergence. Recommend restart. \n";
  206. ivmap( phi, phi2, idx );
  207. return (2);
  208. }
  209. /*
  210. if (errornold - errorn < 1e-14) {
  211. logio << "not making any progress. Giving up\n";
  212. return (1);
  213. }
  214. */
  215. //std::cout << "|| p-s ||" << (alpha*p - omega*s).norm() << std::endl;
  216. // only update phi if good things are happening
  217. phi2.array() += alpha*p.array() + omega*s.array();
  218. errornold = errorn;
  219. }
  220. ivmap( phi, phi2, idx );
  221. return (0);
  222. }
  223. #endif // ----- #ifndef BICGSTAB_INC -----