3D EM based on Schur decomposition
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

EMSchur3D.h 26KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620
  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---\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::ColMajor> > > >( 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. // figure out which omega we are working with
  106. int iw = -1;
  107. for (int iiw=0; iiw<Omegas.size(); ++iiw) {
  108. if (Omegas[iiw] - Source->GetAngularFrequency(0) < 1e-3 ) {
  109. iw = iiw;
  110. }
  111. }
  112. if (iw == -1) {
  113. std::cerr << "FREQUENCY DOOM IN EMSchur3D::SolveSource \n";
  114. exit(EXIT_FAILURE);
  115. }
  116. ///////////////////////////////////
  117. // Set up primary fields
  118. // TODO, this is a little stupid as they all share the same points. We need to extend
  119. // EmEARTH to be able to input a grid so that points are not explicitly needed like
  120. // this. This requires some care as calcs are made on faces.
  121. // Alternatively, the bins function of ReceiverPoints could be extended quite easily.
  122. // This may be the way to do this.
  123. //Lemma::ReceiverPoints* dpoint = Lemma::ReceiverPoints::New();
  124. std::shared_ptr< FieldPoints > dpoint = FieldPoints::NewSP();
  125. FillPoints(dpoint);
  126. PrimaryField(Source, dpoint);
  127. // Allocate a ton of memory
  128. VectorXcr Phi = VectorXcr::Zero(uns);
  129. VectorXcr ms(unx+uny+unz); // mu sigma
  130. // Vector potential (A) Vector and phi
  131. VectorXcr Se = VectorXcr::Zero(unx+uny+unz);
  132. //VectorXcr A = VectorXcr::Zero(unx+uny+unz);
  133. VectorXcr E = VectorXcr::Zero(unx+uny+unz);
  134. VectorXcr E0 = VectorXcr::Zero(unx+uny+unz);
  135. // Lets get cracking
  136. FillSourceTerms(ms, Se, E0, dpoint, Omegas[iw]);
  137. /////////////////////////////////////////////////
  138. // LOG File
  139. std::string logfile (ResFile);
  140. logfile += to_string(isource) + std::string(".log");
  141. ofstream logio(logfile.c_str());
  142. logio << *Source << std::endl;
  143. logio << *Grid << std::endl;
  144. logio << *LayModel << std::endl;
  145. // solve for RHS
  146. int max_it(nx*ny*nz), iter_done(0);
  147. Real tol(3e-16), errorn(0);
  148. logio << "solving RHS for source " << isource << std::endl;
  149. // TODO, this is stupid, try and get rid of this copy!
  150. Eigen::SparseMatrix<Complex> Cc = Cvec[iw];
  151. jsw_timer timer;
  152. jsw_timer timer2;
  153. timer.begin();
  154. timer2.begin();
  155. /////////////////////////////////////////
  156. // Solve for RHS
  157. VectorXcr A = CSolver[iw].solve(Se);
  158. // // Solve Real system instead
  159. // The Real system is quasi-definite, though an LDLT decomposition exists, CHOLMOD doesn't find it.
  160. // An LU can be done on this, but compute performance is very similiar to the complex system, and diagonal pivoting
  161. // cannot be assumed to be best, hurting solve time.
  162. // /* EXPERIMENTAL */
  163. // VectorXr b2 = VectorXr::Zero(2*(unx+uny+unz));
  164. // b2.head(unx+uny+unz) = Se.real();
  165. // b2.tail(unx+uny+unz) = Se.imag();
  166. // VectorXr A2 = CReSolver[iw].solve(b2);
  167. // A.real() = A2.head( unx+uny+unz );
  168. // A.imag() = -A2.tail( unx+uny+unz ); // Due to decomp. negative!
  169. // /* END EXPERIMENTAL */
  170. VectorXcr ADiv = D*A; // ADiv == RHS == D C^I Se
  171. VectorXcr Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() - Se.array());
  172. logio << "|| Div(A) || = " << ADiv.norm()
  173. // << " in " << iter_done << " iterations"
  174. //<< " with error " << errorn << "\t"
  175. << "\tInital solution error "<< Error.norm() // Iteritive info
  176. << "\ttime " << timer.end() << std::endl;
  177. //VectorXcr ADivMAC = ADiv.array() * MAC.array().cast<Complex>();
  178. //logio << "|| Div(A) || on MAC grid " << ADivMAC.norm() << std::endl;
  179. /////////////////////
  180. // Solve for Phi
  181. logio << "Solving for Phi " << std::flush;
  182. timer.begin();
  183. tol = 1e-18;
  184. int success(2);
  185. success = implicitbicgstab(D, idx, ms, ADiv, Phi, CSolver[iw], max_it, tol, errorn, iter_done, logio);
  186. //Phi.array() *= MAC.array().cast<Complex>(); // remove phi from air regions
  187. /* Restart if necessary */
  188. int nrestart(1);
  189. // TODO send MAC to implicitbicgstab?
  190. while (success == 2 && nrestart < 18 && iter_done > 1) {
  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. nrestart += 1;
  194. }
  195. logio << "Implicit BiCGStab solution in " << iter_done << " iterations."
  196. << " with error " << std::setprecision(8) << std::scientific << errorn << std::endl;
  197. logio << "time "<< timer.end() << " [s]" << std::endl;
  198. E = ms.array()*(D.transpose()*Phi).array(); // Temp, field due to charge
  199. /////////////////////////////////////
  200. // Compute A
  201. /////////////////////////////////////
  202. logio << "Solving for A using phi" << std::endl;
  203. std::cout << "Solving for A" << std::endl;
  204. max_it = nx*ny*nz;
  205. tol = 5e-16;
  206. errorn = 0;
  207. iter_done = 0;
  208. timer.begin();
  209. A = CSolver[iw].solve( (Se-E).eval() ); // UmfPack requires eval?
  210. VectorXcr ADiv2 = D*A;
  211. logio << "|| Div(A) || = " << ADiv2.norm() ;
  212. //" in " << iter_done << " iterations"
  213. //<< " with error " << errorn << "\t";
  214. // Report error of solutions
  215. Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() + E.array() - Se.array());
  216. logio << "\tsolution error " << Error.norm()
  217. << std::fixed << std::setprecision(2) << "\ttime " << timer.end() << "\ttotal time " << timer2.end() << std::endl;
  218. logio.close();
  219. //////////////////////////////////////
  220. // Update Fields and report
  221. E.array() = Complex(0,-Omegas[iw])*A.array() - (D.transpose()*Phi).array(); // Secondary Field Only
  222. VectorXcr B = StaggeredGridCurl(A);
  223. WriteVTKResults( ResFile+ to_string(isource), A, Se, E0, E , Phi, ADiv, ADiv2, B);
  224. //dpoint->Delete();
  225. return ;
  226. } // ----- end of method EMSchur3D::SolveSource -----
  227. //--------------------------------------------------------------------------------------
  228. // Class: EMSchur3DBase
  229. // Method: BuildCDirectSolver
  230. //--------------------------------------------------------------------------------------
  231. template < class Solver >
  232. void EMSchur3D<Solver>::BuildCDirectSolver ( ) {
  233. CSolver = new Solver[Omegas.size()];
  234. for (int iw=0; iw<Omegas.size(); ++iw) {
  235. jsw_timer timer;
  236. timer.begin();
  237. /* Complex system */
  238. /*
  239. std::cout << "Generic solver pattern analyzing C_" << iw << ",";
  240. std::cout.flush();
  241. CSolver[iw].analyzePattern( Cvec[iw].selfadjointView< Eigen::Lower>() );
  242. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  243. // factorize
  244. timer.begin();
  245. std::cout << "Generic solver factorising C_" << iw << ", ";
  246. std::cout.flush();
  247. CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
  248. */
  249. std::cerr << "No solver Specified!" << iw << ",";
  250. exit(EXIT_FAILURE);
  251. //CSolver[iw].compute( Cvec[iw].selfadjointView< Eigen::Lower>() );
  252. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  253. }
  254. }
  255. #ifdef HAVE_SUPERLUMT
  256. template<>
  257. void EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::BuildCDirectSolver() {
  258. CSolver = new Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > [Omegas.size()];
  259. for (int iw=0; iw<Omegas.size(); ++iw) {
  260. jsw_timer timer;
  261. timer.begin();
  262. /* SuperLU */
  263. //CSolver[iw].options().DiagPivotThresh = 0.01;
  264. //CSolver[iw].options().SymmetricMode = YES;
  265. //CSolver[iw].options().ColPerm = MMD_AT_PLUS_A;
  266. //CSolver[iw].options().Trans = NOTRANS;
  267. //CSolver[iw].options().ConditionNumber = NO;
  268. //std::cout << "SuperLU options:\n";
  269. //std::cout << "\tPivot Threshold: " << CSolver[iw].options().DiagPivotThresh << std::endl;
  270. //std::cout << "\tSymmetric mode: " << CSolver[iw].options().SymmetricMode << std::endl;
  271. //std::cout << "\tEquilibrate: " << CSolver[iw].options().Equil << std::endl;
  272. //std::cout << "\tCol Permutation: " << CSolver[iw].options().ColPerm << std::endl;
  273. //std::cout << "\tTrans: " << CSolver[iw].options().Trans << std::endl;
  274. //std::cout << "\tCondition Number: " << CSolver[iw].options().ConditionNumber << std::endl;
  275. /* Complex system */
  276. std::cout << "SuperLU_MT pattern analyzing C_" << iw << ",";
  277. std::cout.flush();
  278. CSolver[iw].analyzePattern( Cvec[iw].selfadjointView< Eigen::Lower>() );
  279. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  280. // factorize
  281. timer.begin();
  282. std::cout << "SuperLU_MT factorising C_" << iw << ", ";
  283. std::cout.flush();
  284. CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
  285. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  286. }
  287. }
  288. #endif
  289. template<>
  290. void EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::COLAMDOrdering<int> > >::BuildCDirectSolver() {
  291. CSolver = new Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::COLAMDOrdering<int> > [Omegas.size()];
  292. for (int iw=0; iw<Omegas.size(); ++iw) {
  293. jsw_timer timer;
  294. timer.begin();
  295. CSolver[iw].isSymmetric(true);
  296. CSolver[iw].setPivotThreshold(0.0);
  297. /* Complex system */
  298. std::cout << "SparseLU pattern analyzing C_" << iw << ",";
  299. std::cout.flush();
  300. CSolver[iw].analyzePattern( Cvec[iw].selfadjointView< Eigen::Lower>() );
  301. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  302. // factorize
  303. timer.begin();
  304. std::cout << "SparseLU factorising C_" << iw << ", ";
  305. std::cout.flush();
  306. CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
  307. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  308. }
  309. }
  310. // template<>
  311. // void EMSchur3D< Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
  312. // CSolver = new Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > [Omegas.size()];
  313. // for (int iw=0; iw<Omegas.size(); ++iw) {
  314. // Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  315. // jsw_timer timer;
  316. // timer.begin();
  317. // /* Complex system */
  318. // std::cout << "CholmodSupernodalLLT pattern analyzing C_" << iw << ",";
  319. // std::cout.flush();
  320. // CSolver[iw].analyzePattern( Csym );
  321. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  322. // /* factorize */
  323. // timer.begin();
  324. // std::cout << "CholmodSupernodalLLT factorising C_" << iw << ", ";
  325. // std::cout.flush();
  326. // CSolver[iw].factorize( Csym );
  327. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  328. // }
  329. // }
  330. // template<>
  331. // void EMSchur3D< Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::NaturalOrdering<int> > > ::BuildCDirectSolver() {
  332. // CSolver = new Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::NaturalOrdering<int> > [Omegas.size()];
  333. // for (int iw=0; iw<Omegas.size(); ++iw) {
  334. // Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  335. // jsw_timer timer;
  336. // timer.begin();
  337. // /* Complex system */
  338. // std::cout << "CSymSimplicialLLT<NaturalOrdering> pattern analyzing C_" << iw << ",";
  339. // std::cout.flush();
  340. // CSolver[iw].analyzePattern( Csym );
  341. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  342. // /* factorize */
  343. // timer.begin();
  344. // std::cout << "CSymSimplicialLLT<NaturalOrdering> factorising C_" << iw << ", ";
  345. // std::cout.flush();
  346. // CSolver[iw].factorize( Csym );
  347. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  348. // }
  349. // }
  350. //
  351. // template<>
  352. // void EMSchur3D< Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > > ::BuildCDirectSolver() {
  353. // CSolver = new Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > [Omegas.size()];
  354. // for (int iw=0; iw<Omegas.size(); ++iw) {
  355. // //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  356. // jsw_timer timer;
  357. // timer.begin();
  358. // /* Complex system */
  359. // std::cout << "CSymSimplicialLLT<AMDOrdering> pattern analyzing C_" << iw << ",";
  360. // std::cout.flush();
  361. // CSolver[iw].analyzePattern( Cvec[iw] );
  362. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  363. // /* factorize */
  364. // timer.begin();
  365. // std::cout << "CSymSimplicialLLT<AMDOrdering> factorising C_" << iw << ", ";
  366. // std::cout.flush();
  367. // CSolver[iw].factorize( Cvec[iw] );
  368. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  369. // }
  370. // }
  371. //
  372. // template<>
  373. // void EMSchur3D< Eigen::CSymSimplicialLDLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > > ::BuildCDirectSolver() {
  374. // CSolver = new Eigen::CSymSimplicialLDLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > [Omegas.size()];
  375. // for (int iw=0; iw<Omegas.size(); ++iw) {
  376. // Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  377. // jsw_timer timer;
  378. // timer.begin();
  379. // /* Complex system */
  380. // std::cout << "CSymSimplicialLDLT<AMDOrdering> pattern analyzing C_" << iw << ",";
  381. // std::cout.flush();
  382. // CSolver[iw].analyzePattern( Csym );
  383. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  384. // /* factorize */
  385. // timer.begin();
  386. // std::cout << "CSymSimplicialLDLT<AMDOrdering> factorising C_" << iw << ", ";
  387. // std::cout.flush();
  388. // CSolver[iw].factorize( Csym );
  389. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  390. // }
  391. // }
  392. template<>
  393. void EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::IncompleteLUT<Complex> > > ::BuildCDirectSolver() {
  394. CSolver = new Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::IncompleteLUT<Complex> > [Omegas.size()];
  395. for (int iw=0; iw<Omegas.size(); ++iw) {
  396. Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  397. jsw_timer timer;
  398. timer.begin();
  399. /* Complex system */
  400. std::cout << "BiCGSTAB(ILU) pattern analyzing C_" << iw << ",";
  401. std::cout.flush();
  402. CSolver[iw].analyzePattern( Csym );
  403. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  404. /* factorize */
  405. timer.begin();
  406. std::cout << "BiCGSTAB(ILU) factorising C_" << iw << ", ";
  407. std::cout.flush();
  408. CSolver[iw].factorize( Csym );
  409. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  410. }
  411. }
  412. template<>
  413. void EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > > ::BuildCDirectSolver() {
  414. CSolver = new Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > [Omegas.size()];
  415. for (int iw=0; iw<Omegas.size(); ++iw) {
  416. Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  417. jsw_timer timer;
  418. timer.begin();
  419. /* Complex system */
  420. std::cout << "BiCGSTAB pattern analyzing C_" << iw << ",";
  421. std::cout.flush();
  422. CSolver[iw].analyzePattern( Csym );
  423. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  424. // factorize
  425. timer.begin();
  426. std::cout << "BiCGSTAB factorising C_" << iw << ", ";
  427. std::cout.flush();
  428. CSolver[iw].factorize( Csym );
  429. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  430. }
  431. }
  432. template<>
  433. void EMSchur3D< Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
  434. CSolver = new Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > [Omegas.size()];
  435. for (int iw=0; iw<Omegas.size(); ++iw) {
  436. //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  437. jsw_timer timer;
  438. timer.begin();
  439. /* Complex system */
  440. std::cout << "ConjugateGradient pattern analyzing C_" << iw << ",";
  441. std::cout.flush();
  442. CSolver[iw].analyzePattern( Cvec[iw] );
  443. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  444. // factorize
  445. timer.begin();
  446. std::cout << "ConjugateGradient factorising C_" << iw << ", ";
  447. std::cout.flush();
  448. CSolver[iw].factorize( Cvec[iw] );
  449. std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  450. }
  451. }
  452. // template<>
  453. // void EMSchur3D< Eigen::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
  454. // CSolver = new Eigen::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > [Omegas.size()];
  455. // //MPI_Init(NULL, NULL);
  456. // for (int iw=0; iw<Omegas.size(); ++iw) {
  457. // //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  458. // jsw_timer timer;
  459. // timer.begin();
  460. // /* Complex system */
  461. // std::cout << "PaStiX LLT pattern analyzing C_" << iw << ",";
  462. // std::cout.flush();
  463. // CSolver[iw].analyzePattern( Cvec[iw] );
  464. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  465. // // factorize
  466. // timer.begin();
  467. // std::cout << "PaStiX LLT factorising C_" << iw << ", ";
  468. // std::cout.flush();
  469. // CSolver[iw].factorize( Cvec[iw] );
  470. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  471. // }
  472. // }
  473. //
  474. // template<>
  475. // void EMSchur3D< Eigen::PastixLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
  476. // CSolver = new Eigen::PastixLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > [Omegas.size()];
  477. // //MPI_Init(NULL, NULL);
  478. // for (int iw=0; iw<Omegas.size(); ++iw) {
  479. // //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  480. // jsw_timer timer;
  481. // timer.begin();
  482. // /* Complex system */
  483. // std::cout << "PaStiX LDLT pattern analyzing C_" << iw << ",";
  484. // std::cout.flush();
  485. // CSolver[iw].analyzePattern( Cvec[iw] );
  486. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  487. // // factorize
  488. // timer.begin();
  489. // std::cout << "PaStiX LDLT factorising C_" << iw << ", ";
  490. // std::cout.flush();
  491. // CSolver[iw].factorize( Cvec[iw] );
  492. // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  493. // std::cout << "INFO " << CSolver[iw].info( ) << std::endl;
  494. // }
  495. // }
  496. //
  497. // template<>
  498. // void EMSchur3D< Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, true > > ::BuildCDirectSolver() {
  499. // CSolver = new Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, true > [Omegas.size()];
  500. // //MPI_Init(NULL, NULL);
  501. // for (int iw=0; iw<Omegas.size(); ++iw) {
  502. // Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
  503. // jsw_timer timer;
  504. // timer.begin();
  505. // /* Complex system */
  506. // std::cout << "PaStiX LU pattern analyzing C_" << iw << ",";
  507. // std::cout.flush();
  508. // CSolver[iw].compute( Csym );
  509. // std::cout << "PaStiX LU Done C_" << iw << std::endl;;
  510. // // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  511. // // // factorize
  512. // // timer.begin();
  513. // // std::cout << "PaStiX LU factorising C_" << iw << ", ";
  514. // // std::cout.flush();
  515. // // CSolver[iw].factorize( Csym );
  516. // // std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
  517. // }
  518. // }
  519. } // ----- end of Lemma name -----
  520. #endif // ----- #ifndef EMSCHUR3D_INC -----