From ed7f6c2c780b82ba2f039812f29b18386b445c8e Mon Sep 17 00:00:00 2001
From: Arthur <arthur.alb53@gmail.com>
Date: Wed, 1 Jan 2025 12:17:40 +0100
Subject: [PATCH] PowerIteration class creation

---
 examples/main.cpp            | 29 ++++++++++++++++++++++---
 src/demeter/solve/solver.cpp | 35 ++++++++++++++++++++++++++++++
 src/demeter/solve/solver.hpp | 42 ++++++++++++++++++++++++++++++++++++
 3 files changed, 103 insertions(+), 3 deletions(-)

diff --git a/examples/main.cpp b/examples/main.cpp
index 5e2f75e..c568c7e 100644
--- a/examples/main.cpp
+++ b/examples/main.cpp
@@ -25,9 +25,34 @@ int main() {
   
   ArrayXd weights(nteta/2);
   weights = gauss.getweights();
-  std::cout << weights;
+  std::cout << weights.sum();
   std::cout << "\n";
 
+  // Test power iteration solver
+  // 1st problem, vp : (5, 2)
+  Eigen::MatrixXd A1(2,2);
+  A1 << 4, 1,
+        1, 3;
+  // PowerIteration PowIte(A);
+  // PowIte.solve(); 
+
+  // 2nd problem, vp : (11.25, 7.34, 5.11, 4.30)
+  Eigen::MatrixXd B1(4,4);
+  B1 << 10, 2, 0, 0,
+         2, 8, 1, 0,
+         0, 1, 6, 1,
+         0, 0, 1, 4;
+
+  // 3rd problem
+  // vp : (11.34, 11.28, 8.56, 6.68, 4.84, 2.20)
+  Eigen::MatrixXd C1(6,6);
+  C1 << 12,  3, 0, 0, 0, 0, 
+         3, 11, 2, 0, 0, 0,
+         0,  2, 9, 1, 0, 0,
+         0,  0, 1, 7, 1, 0,
+         0,  0, 0, 1, 5, 1,
+         0,  0, 0, 0, 1, 3;
+
   // Test linear solvers
   Eigen::Matrix2f A, b;
   A << 2, -1, -1, 3; 
@@ -39,8 +64,6 @@ int main() {
   LinSolver.solve_PartialPivLu();
   std::cout << '\n';
 
-
-  
   // Define cross sections
   ArrayXd  sigma_t(ngroups);
   ArrayXXd sigma_s(ngroups,ngroups);
diff --git a/src/demeter/solve/solver.cpp b/src/demeter/solve/solver.cpp
index 6851465..af88720 100644
--- a/src/demeter/solve/solver.cpp
+++ b/src/demeter/solve/solver.cpp
@@ -18,4 +18,39 @@ void Solver::solve() {
   // Calcul des sources
   // Methode de la puissance iteree
 }
+
+void PowerIteration::solve() {
+ 
+  size_t size;
+  double lambda_k1;
+  double lambda_k2;
+  double error_eigvalue;
+  
+  // Initialization
+  lambda_k1 = 0;
+  lambda_k2 = 0;
+  size = A_.rows();
+  Eigen::ArrayXd x1 = Eigen::ArrayXd::Ones(size);
+  Eigen::ArrayXd x2 = Eigen::ArrayXd::Ones(size);
+
+  // Iterations
+  while (error_eigvalue >= tol_eigvalue_ ) {
+
+     // Vector update 
+     // x2 = A_ * x1;
+     // x2 = x2 / |x2|
+
+     // Value update 
+     //lambda_k1 = 1;     
+     //lambda_k2 = <x2,x1> / (<x1,x1>;
+     //x2 = x2 / |x2|
+     
+     // Convergence test
+     error_eigvalue = lambda_k2-lambda_k1;
+     std::cout << "Lambda: " << lambda_k2 << std::endl;
+     } 
+
+  std::cout << "Eigenvalue: " << lambda_k2 << std::endl;
+
+}
 }  // namespace Demeter
diff --git a/src/demeter/solve/solver.hpp b/src/demeter/solve/solver.hpp
index e063f87..8d1cf9f 100644
--- a/src/demeter/solve/solver.hpp
+++ b/src/demeter/solve/solver.hpp
@@ -35,4 +35,46 @@ class Solver {
   double tolerence_flux_ = 1e-5;
   double tolerence_fiss_ = 1e-4;
 };
+
+/**
+ * @brief Class that solve a problem of the
+ * type Ax=lb, with l a scalar.
+ */
+class PowerIteration {
+ public:
+ 
+  /**
+  * @brief Standard eigenvalue problem constructor.
+  */
+  PowerIteration(const Eigen::MatrixXd& A) 
+      : A_(A) {
+      std::cout << "Initialisation" << std::endl;
+  }
+
+  /**
+  * @brief Generalized eigenvalue problem constructor. 
+  */
+  PowerIteration(const Eigen::MatrixXd& A, const Eigen::MatrixXd& B) 
+      : A_(A), B_(B) {
+      std::cout << "Initialisation" << std::endl;
+  }
+
+  /**
+  * @brief Solve the eigenvalue problem. 
+  */
+  void solve();
+
+ private:
+  
+  double lambda_;
+  double tol_eigvector_ = 1e-5;
+  double tol_eigvalue_  = 1e-5;
+
+  Eigen::MatrixXd A_; 
+  Eigen::MatrixXd B_;
+  Eigen::ArrayXd  x_;
+};
+
+
+
 }  // namespace Demeter