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

EMSchur3D.h 26KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630
  1. /* This file is part of Lemma, a geophysical modelling and inversion API.
  2. * More information is available at http://lemmasoftware.org
  3. */
  4. /* This Source Code Form is subject to the terms of the Mozilla Public
  5. * License, v. 2.0. If a copy of the MPL was not distributed with this
  6. * file, You can obtain one at http://mozilla.org/MPL/2.0/.
  7. */
  8. /**
  9. * @file
  10. * @date 02/19/2015 01:10:39 PM
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email Trevor.Irons@xri-geo.com
  14. * @copyright Copyright (c) 2015, XRI Geophysics, LLC
  15. * @copyright Copyright (c) 2015, Trevor Irons
  16. * @copyright Copyright (c) 2011, Trevor Irons
  17. * @copyright Copyright (c) 2011, Colorado School of Mines
  18. */
  19. #ifndef EMSCHUR3D_INC
  20. #define EMSCHUR3D_INC
  21. #include "EMSchur3DBase.h"
  22. #include "bicgstab.h"
  23. //#include "CSymSimplicialCholesky.h"
  24. namespace Lemma {
  25. /**
  26. \brief Templated concrete classes of EMSChur3DBase.
  27. \details
  28. */
  29. template < class Solver >
  30. class EMSchur3D : public EMSchur3DBase {
  31. friend std::ostream &operator << (std::ostream &stream, const EMSchur3D &ob) {
  32. stream << ob.Serialize() << "\n"; // End of doc
  33. return stream;
  34. }
  35. //friend std::ostream &operator<<(std::ostream &stream,
  36. // const EMSchur3D &ob);
  37. public:
  38. // ==================== LIFECYCLE =======================
  39. /**
  40. * @copybrief LemmaObject::New()
  41. * @copydetails LemmaObject::New()
  42. */
  43. static std::shared_ptr< EMSchur3D > NewSP() {
  44. return std::make_shared< EMSchur3D<Solver> >( ctor_key() );
  45. //return std::make_shared< EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > > >( ctor_key() ) ;
  46. }
  47. /** Default protected constructor, use New */
  48. explicit EMSchur3D ( const ctor_key& key ) : EMSchur3DBase( key ), CSolver( nullptr ) {
  49. }
  50. /** Locked DeDerializing constructor, use factory DeSerialize method*/
  51. EMSchur3D (const YAML::Node& node, const ctor_key& key): EMSchur3DBase(node, key), CSolver( nullptr ) {
  52. }
  53. /** Default protected destructor, use Delete */
  54. virtual ~EMSchur3D () {
  55. // TODO delete arrays
  56. }
  57. /**
  58. * Uses YAML to serialize this object.
  59. * @return a YAML::Node
  60. */
  61. YAML::Node Serialize() const {
  62. YAML::Node node = EMSchur3DBase::Serialize();
  63. //node["NumberOfLayers"] = NumberOfLayers;
  64. node.SetTag( this->GetName() );
  65. return node;
  66. }
  67. /**
  68. * Constructs an object from a YAML::Node.
  69. */
  70. static EMSchur3D* DeSerialize(const YAML::Node& node);
  71. // ==================== OPERATORS =======================
  72. // ==================== OPERATIONS =======================
  73. /** Solves a single source problem. This method is thread safe.
  74. * @param[in] Source is the source term for generating primary fields
  75. * @param[in] isource is the source index
  76. */
  77. void SolveSource( std::shared_ptr<DipoleSource> Source , const int& isource);
  78. /** Builds the solver for the C matrix */
  79. void BuildCDirectSolver( );
  80. // ==================== ACCESS =======================
  81. virtual std::string GetName() const {
  82. return this->CName;
  83. }
  84. // ==================== INQUIRY =======================
  85. protected:
  86. // ==================== LIFECYCLE =======================
  87. private:
  88. /** Copy constructor */
  89. EMSchur3D( const EMSchur3D& ) = delete;
  90. // ==================== DATA MEMBERS =========================
  91. /** The templated solver for C */
  92. Solver* CSolver;
  93. Eigen::SparseMatrix<Complex> Csym;
  94. static constexpr auto CName = "EMSchur3D";
  95. }; // ----- end of class EMSchur3D -----
  96. ////////////////////////////////////////////////////////////////////////////////////////
  97. // Implimentation and Specialisations //
  98. ////////////////////////////////////////////////////////////////////////////////////////
  99. //--------------------------------------------------------------------------------------
  100. // Class: EMSchur3D
  101. // Method: SolveSource
  102. //--------------------------------------------------------------------------------------
  103. template < class Solver >
  104. void EMSchur3D<Solver>::SolveSource ( std::shared_ptr<DipoleSource> Source, const int& isource ) {
  105. std::cout << "In vanilla SolveSource" << std::endl;
  106. // figure out which omega we are working with
  107. int iw = -1;
  108. for (int iiw=0; iiw<Omegas.size(); ++iiw) {
  109. if (Omegas[iiw] - Source->GetAngularFrequency(0) < 1e-3 ) {
  110. iw = iiw;
  111. }
  112. }
  113. if (iw == -1) {
  114. std::cerr << "FREQUENCY DOOM IN EMSchur3D::SolveSource \n";
  115. exit(EXIT_FAILURE);
  116. }
  117. ///////////////////////////////////
  118. // Set up primary fields
  119. // TODO, this is a little stupid as they all share the same points. We need to extend
  120. // EmEARTH to be able to input a grid so that points are not explicitly needed like
  121. // this. This requires some care as calcs are made on faces.
  122. // Alternatively, the bins function of ReceiverPoints could be extended quite easily.
  123. // This may be the way to do this.
  124. //Lemma::ReceiverPoints* dpoint = Lemma::ReceiverPoints::New();
  125. std::shared_ptr< FieldPoints > dpoint = FieldPoints::NewSP();
  126. FillPoints(dpoint);
  127. PrimaryField(Source, dpoint);
  128. std::cout << "Done with primary field" << std::endl;
  129. // Allocate a ton of memory
  130. VectorXcr Phi = VectorXcr::Zero(uns);
  131. VectorXcr ms(unx+uny+unz); // mu sigma
  132. // Vector potential (A) Vector and phi
  133. VectorXcr Se = VectorXcr::Zero(unx+uny+unz);
  134. //VectorXcr A = VectorXcr::Zero(unx+uny+unz);
  135. VectorXcr E = VectorXcr::Zero(unx+uny+unz);
  136. VectorXcr E0 = VectorXcr::Zero(unx+uny+unz);
  137. // Lets get cracking
  138. std::cout << "Filling source terms" << std::endl;
  139. FillSourceTerms(ms, Se, E0, dpoint, Omegas[iw]);
  140. std::cout << "Done source terms" << std::endl;
  141. /////////////////////////////////////////////////
  142. // LOG File
  143. std::string logfile (ResFile);
  144. logfile += to_string(isource) + std::string(".log");
  145. ofstream logio(logfile.c_str());
  146. std::cout << "just logging" << std::endl;
  147. // logio << *Source << std::endl;
  148. logio << *Grid << std::endl;
  149. logio << *LayModel << std::endl;
  150. std::cout << "dun logging" << std::endl;
  151. // solve for RHS
  152. int max_it(nx*ny*nz), iter_done(0);
  153. Real tol(3e-16), errorn(0);
  154. logio << "solving RHS for source " << isource << std::endl;
  155. // TODO, this is stupid, try and get rid of this copy!
  156. Eigen::SparseMatrix<Complex> Cc = Cvec[iw];
  157. jsw_timer timer;
  158. jsw_timer timer2;
  159. timer.begin();
  160. timer2.begin();
  161. /////////////////////////////////////////
  162. // Solve for RHS
  163. VectorXcr A = CSolver[iw].solve(Se);
  164. // // Solve Real system instead
  165. // The Real system is quasi-definite, though an LDLT decomposition exists, CHOLMOD doesn't find it.
  166. // An LU can be done on this, but compute performance is very similiar to the complex system, and diagonal pivoting
  167. // cannot be assumed to be best, hurting solve time.
  168. // /* EXPERIMENTAL */
  169. // VectorXr b2 = VectorXr::Zero(2*(unx+uny+unz));
  170. // b2.head(unx+uny+unz) = Se.real();
  171. // b2.tail(unx+uny+unz) = Se.imag();
  172. // VectorXr A2 = CReSolver[iw].solve(b2);
  173. // A.real() = A2.head( unx+uny+unz );
  174. // A.imag() = -A2.tail( unx+uny+unz ); // Due to decomp. negative!
  175. // /* END EXPERIMENTAL */
  176. VectorXcr ADiv = D*A; // ADiv == RHS == D C^I Se
  177. VectorXcr Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() - Se.array());
  178. logio << "|| Div(A) || = " << ADiv.norm()
  179. // << " in " << iter_done << " iterations"
  180. //<< " with error " << errorn << "\t"
  181. << "\tInital solution error "<< Error.norm() // Iteritive info
  182. << "\ttime " << timer.end() << std::endl;
  183. //VectorXcr ADivMAC = ADiv.array() * MAC.array().cast<Complex>();
  184. //logio << "|| Div(A) || on MAC grid " << ADivMAC.norm() << std::endl;
  185. /////////////////////
  186. // Solve for Phi
  187. logio << "Solving for Phi " << std::flush;
  188. timer.begin();
  189. tol = 1e-18;
  190. int success(2);
  191. success = implicitbicgstab(D, idx, ms, ADiv, Phi, CSolver[iw], max_it, tol, errorn, iter_done, logio);
  192. //Phi.array() *= MAC.array().cast<Complex>(); // remove phi from air regions
  193. /* Restart if necessary */
  194. int nrestart(1);
  195. // TODO send MAC to implicitbicgstab?
  196. while (success == 2 && nrestart < 18 && iter_done > 1) {
  197. success = implicitbicgstab(D, idx, ms, ADiv, Phi, CSolver[iw], max_it, tol, errorn, iter_done, logio);
  198. //Phi.array() *= MAC.array().cast<Complex>(); // remove phi from air regions
  199. nrestart += 1;
  200. }
  201. logio << "Implicit BiCGStab solution in " << iter_done << " iterations."
  202. << " with error " << std::setprecision(8) << std::scientific << errorn << std::endl;
  203. logio << "time "<< timer.end() << " [s]" << std::endl;
  204. E = ms.array()*(D.transpose()*Phi).array(); // Temp, field due to charge
  205. /////////////////////////////////////
  206. // Compute A
  207. /////////////////////////////////////
  208. logio << "Solving for A using phi" << std::endl;
  209. std::cout << "Solving for A" << std::endl;
  210. max_it = nx*ny*nz;
  211. tol = 5e-16;
  212. errorn = 0;
  213. iter_done = 0;
  214. timer.begin();
  215. A = CSolver[iw].solve( (Se-E).eval() ); // UmfPack requires eval?
  216. VectorXcr ADiv2 = D*A;
  217. logio << "|| Div(A) || = " << ADiv2.norm() ;
  218. //" in " << iter_done << " iterations"
  219. //<< " with error " << errorn << "\t";
  220. // Report error of solutions
  221. Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() + E.array() - Se.array());
  222. logio << "\tsolution error " << Error.norm()
  223. << std::fixed << std::setprecision(2) << "\ttime " << timer.end() << "\ttotal time " << timer2.end() << std::endl;
  224. logio.close();
  225. //////////////////////////////////////
  226. // Update Fields and report
  227. E.array() = Complex(0,-Omegas[iw])*A.array() - (D.transpose()*Phi).array(); // Secondary Field Only
  228. VectorXcr B = StaggeredGridCurl(A);
  229. WriteVTKResults( ResFile+ to_string(isource), A, Se, E0, E , Phi, ADiv, ADiv2, B);
  230. //dpoint->Delete();
  231. return ;
  232. } // ----- end of method EMSchur3D::SolveSource -----
  233. //--------------------------------------------------------------------------------------
  234. // Class: EMSchur3DBase
  235. // Method: BuildCDirectSolver
  236. //--------------------------------------------------------------------------------------
  237. template < class Solver >
  238. void EMSchur3D<Solver>::BuildCDirectSolver ( ) {
  239. CSolver = new Solver[Omegas.size()];
  240. for (int iw=0; iw<Omegas.size(); ++iw) {
  241. jsw_timer timer;
  242. timer.begin();
  243. /* Complex system */
  244. /*
  245. std::cout << "Generic solver pattern analyzing C_" << iw << ",";
  246. std::cout.flush();
  247. CSolver[iw].analyzePattern( Cvec[iw].selfadjointView< Eigen::Lower>() );
  248. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  249. // factorize
  250. timer.begin();
  251. std::cout << "Generic solver factorising C_" << iw << ", ";
  252. std::cout.flush();
  253. CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
  254. */
  255. std::cerr << "No solver Specified!" << iw << ",";
  256. exit(EXIT_FAILURE);
  257. //CSolver[iw].compute( Cvec[iw].selfadjointView< Eigen::Lower>() );
  258. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  259. }
  260. }
  261. #ifdef HAVE_SUPERLUMT
  262. template<>
  263. void EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::BuildCDirectSolver() {
  264. CSolver = new Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > [Omegas.size()];
  265. for (int iw=0; iw<Omegas.size(); ++iw) {
  266. jsw_timer timer;
  267. timer.begin();
  268. /* SuperLU */
  269. //CSolver[iw].options().DiagPivotThresh = 0.01;
  270. //CSolver[iw].options().SymmetricMode = YES;
  271. //CSolver[iw].options().ColPerm = MMD_AT_PLUS_A;
  272. //CSolver[iw].options().Trans = NOTRANS;
  273. //CSolver[iw].options().ConditionNumber = NO;
  274. //std::cout << "SuperLU options:\n";
  275. //std::cout << "\tPivot Threshold: " << CSolver[iw].options().DiagPivotThresh << std::endl;
  276. //std::cout << "\tSymmetric mode: " << CSolver[iw].options().SymmetricMode << std::endl;
  277. //std::cout << "\tEquilibrate: " << CSolver[iw].options().Equil << std::endl;
  278. //std::cout << "\tCol Permutation: " << CSolver[iw].options().ColPerm << std::endl;
  279. //std::cout << "\tTrans: " << CSolver[iw].options().Trans << std::endl;
  280. //std::cout << "\tCondition Number: " << CSolver[iw].options().ConditionNumber << std::endl;
  281. /* Complex system */
  282. std::cout << "SuperLU_MT pattern analyzing C_" << iw << ",";
  283. std::cout.flush();
  284. CSolver[iw].analyzePattern( Cvec[iw].selfadjointView< Eigen::Lower>() );
  285. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  286. // factorize
  287. timer.begin();
  288. std::cout << "SuperLU_MT factorising C_" << iw << ", ";
  289. std::cout.flush();
  290. CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
  291. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  292. }
  293. }
  294. #endif
  295. template<>
  296. void EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::COLAMDOrdering<int> > >::BuildCDirectSolver() {
  297. CSolver = new Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::COLAMDOrdering<int> > [Omegas.size()];
  298. for (int iw=0; iw<Omegas.size(); ++iw) {
  299. jsw_timer timer;
  300. timer.begin();
  301. CSolver[iw].isSymmetric(true);
  302. CSolver[iw].setPivotThreshold(0.0);
  303. /* Complex system */
  304. std::cout << "SparseLU pattern analyzing C_" << iw << ",";
  305. std::cout.flush();
  306. CSolver[iw].analyzePattern( Cvec[iw].selfadjointView< Eigen::Lower>() );
  307. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  308. // factorize
  309. timer.begin();
  310. std::cout << "SparseLU factorising C_" << iw << ", ";
  311. std::cout.flush();
  312. CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
  313. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  314. }
  315. }
  316. // template<>
  317. // void EMSchur3D< Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
  318. // CSolver = new Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > [Omegas.size()];
  319. // for (int iw=0; iw<Omegas.size(); ++iw) {
  320. // Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  321. // jsw_timer timer;
  322. // timer.begin();
  323. // /* Complex system */
  324. // std::cout << "CholmodSupernodalLLT pattern analyzing C_" << iw << ",";
  325. // std::cout.flush();
  326. // CSolver[iw].analyzePattern( Csym );
  327. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  328. // /* factorize */
  329. // timer.begin();
  330. // std::cout << "CholmodSupernodalLLT factorising C_" << iw << ", ";
  331. // std::cout.flush();
  332. // CSolver[iw].factorize( Csym );
  333. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  334. // }
  335. // }
  336. // template<>
  337. // void EMSchur3D< Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::NaturalOrdering<int> > > ::BuildCDirectSolver() {
  338. // CSolver = new Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::NaturalOrdering<int> > [Omegas.size()];
  339. // for (int iw=0; iw<Omegas.size(); ++iw) {
  340. // Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  341. // jsw_timer timer;
  342. // timer.begin();
  343. // /* Complex system */
  344. // std::cout << "CSymSimplicialLLT<NaturalOrdering> pattern analyzing C_" << iw << ",";
  345. // std::cout.flush();
  346. // CSolver[iw].analyzePattern( Csym );
  347. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  348. // /* factorize */
  349. // timer.begin();
  350. // std::cout << "CSymSimplicialLLT<NaturalOrdering> factorising C_" << iw << ", ";
  351. // std::cout.flush();
  352. // CSolver[iw].factorize( Csym );
  353. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  354. // }
  355. // }
  356. //
  357. // template<>
  358. // void EMSchur3D< Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > > ::BuildCDirectSolver() {
  359. // CSolver = new Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > [Omegas.size()];
  360. // for (int iw=0; iw<Omegas.size(); ++iw) {
  361. // //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  362. // jsw_timer timer;
  363. // timer.begin();
  364. // /* Complex system */
  365. // std::cout << "CSymSimplicialLLT<AMDOrdering> pattern analyzing C_" << iw << ",";
  366. // std::cout.flush();
  367. // CSolver[iw].analyzePattern( Cvec[iw] );
  368. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  369. // /* factorize */
  370. // timer.begin();
  371. // std::cout << "CSymSimplicialLLT<AMDOrdering> factorising C_" << iw << ", ";
  372. // std::cout.flush();
  373. // CSolver[iw].factorize( Cvec[iw] );
  374. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  375. // }
  376. // }
  377. //
  378. // template<>
  379. // void EMSchur3D< Eigen::CSymSimplicialLDLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > > ::BuildCDirectSolver() {
  380. // CSolver = new Eigen::CSymSimplicialLDLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > [Omegas.size()];
  381. // for (int iw=0; iw<Omegas.size(); ++iw) {
  382. // Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  383. // jsw_timer timer;
  384. // timer.begin();
  385. // /* Complex system */
  386. // std::cout << "CSymSimplicialLDLT<AMDOrdering> pattern analyzing C_" << iw << ",";
  387. // std::cout.flush();
  388. // CSolver[iw].analyzePattern( Csym );
  389. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  390. // /* factorize */
  391. // timer.begin();
  392. // std::cout << "CSymSimplicialLDLT<AMDOrdering> factorising C_" << iw << ", ";
  393. // std::cout.flush();
  394. // CSolver[iw].factorize( Csym );
  395. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  396. // }
  397. // }
  398. template<>
  399. void EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::IncompleteLUT<Complex> > > ::BuildCDirectSolver() {
  400. CSolver = new Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::IncompleteLUT<Complex> > [Omegas.size()];
  401. for (int iw=0; iw<Omegas.size(); ++iw) {
  402. Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  403. CSolver[iw].preconditioner().setDroptol(1e-12);
  404. CSolver[iw].preconditioner().setFillfactor(1e2);
  405. jsw_timer timer;
  406. timer.begin();
  407. /* Complex system */
  408. std::cout << "BiCGSTAB(ILU) pattern analyzing C_" << iw << ",";
  409. std::cout.flush();
  410. CSolver[iw].analyzePattern( Csym );
  411. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  412. /* factorize */
  413. timer.begin();
  414. std::cout << "BiCGSTAB(ILU) factorising C_" << iw << ", ";
  415. std::cout.flush();
  416. CSolver[iw].factorize( Csym );
  417. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  418. }
  419. }
  420. template<>
  421. void EMSchur3D< Eigen::BiCGSTAB< Eigen::SparseMatrix<Complex, Eigen::RowMajor> > > ::BuildCDirectSolver() {
  422. CSolver = new Eigen::BiCGSTAB< Eigen::SparseMatrix<Complex, Eigen::RowMajor> > [Omegas.size()];
  423. for (int iw=0; iw<Omegas.size(); ++iw) {
  424. Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  425. jsw_timer timer;
  426. timer.begin();
  427. /* Complex system */
  428. std::cout << "BiCGSTAB pattern analyzing C_" << iw << ",";
  429. std::cout.flush();
  430. CSolver[iw].analyzePattern( Csym );
  431. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  432. // factorize
  433. timer.begin();
  434. std::cout << "BiCGSTAB factorising C_" << iw << ", ";
  435. std::cout.flush();
  436. CSolver[iw].factorize( Csym );
  437. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  438. }
  439. }
  440. // template<>
  441. // void EMSchur3D< Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
  442. // CSolver = new Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > [Omegas.size()];
  443. // for (int iw=0; iw<Omegas.size(); ++iw) {
  444. // //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  445. // jsw_timer timer;
  446. // timer.begin();
  447. // /* Complex system */
  448. // std::cout << "ConjugateGradient pattern analyzing C_" << iw << ",";
  449. // std::cout.flush();
  450. // CSolver[iw].analyzePattern( Cvec[iw] );
  451. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  452. // // factorize
  453. // timer.begin();
  454. // std::cout << "ConjugateGradient factorising C_" << iw << ", ";
  455. // std::cout.flush();
  456. // CSolver[iw].factorize( Cvec[iw] );
  457. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  458. // }
  459. // }
  460. // template<>
  461. // void EMSchur3D< Eigen::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
  462. // CSolver = new Eigen::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > [Omegas.size()];
  463. // //MPI_Init(NULL, NULL);
  464. // for (int iw=0; iw<Omegas.size(); ++iw) {
  465. // //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  466. // jsw_timer timer;
  467. // timer.begin();
  468. // /* Complex system */
  469. // std::cout << "PaStiX LLT pattern analyzing C_" << iw << ",";
  470. // std::cout.flush();
  471. // CSolver[iw].analyzePattern( Cvec[iw] );
  472. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  473. // // factorize
  474. // timer.begin();
  475. // std::cout << "PaStiX LLT factorising C_" << iw << ", ";
  476. // std::cout.flush();
  477. // CSolver[iw].factorize( Cvec[iw] );
  478. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  479. // }
  480. // }
  481. //
  482. // template<>
  483. // void EMSchur3D< Eigen::PastixLDLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
  484. // CSolver = new Eigen::PastixLDLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > [Omegas.size()];
  485. // //MPI_Init(NULL, NULL);
  486. // for (int iw=0; iw<Omegas.size(); ++iw) {
  487. // //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  488. // jsw_timer timer;
  489. // timer.begin();
  490. // /* Complex system */
  491. // std::cout << "PaStiX LDLT pattern analyzing C_" << iw << ",";
  492. // std::cout.flush();
  493. // CSolver[iw].analyzePattern( Cvec[iw] );
  494. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  495. // // factorize
  496. // timer.begin();
  497. // std::cout << "PaStiX LDLT factorising C_" << iw << ", ";
  498. // std::cout.flush();
  499. // CSolver[iw].factorize( Cvec[iw] );
  500. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  501. // std::cout << "INFO " << CSolver[iw].info( ) << std::endl;
  502. // }
  503. // }
  504. //
  505. // template<>
  506. // void EMSchur3D< Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, true > > ::BuildCDirectSolver() {
  507. // CSolver = new Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, true > [Omegas.size()];
  508. // //MPI_Init(NULL, NULL);
  509. // for (int iw=0; iw<Omegas.size(); ++iw) {
  510. // Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  511. // jsw_timer timer;
  512. // timer.begin();
  513. // /* Complex system */
  514. // std::cout << "PaStiX LU pattern analyzing C_" << iw << ",";
  515. // std::cout.flush();
  516. // CSolver[iw].compute( Csym );
  517. // std::cout << "PaStiX LU Done C_" << iw << std::endl;;
  518. // // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  519. // // // factorize
  520. // // timer.begin();
  521. // // std::cout << "PaStiX LU factorising C_" << iw << ", ";
  522. // // std::cout.flush();
  523. // // CSolver[iw].factorize( Csym );
  524. // // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  525. // }
  526. // }
  527. } // ----- end of Lemma name -----
  528. #endif // ----- #ifndef EMSCHUR3D_INC -----