From 2e85a02ad57af0e552ac2dd27b7e1080b19859db Mon Sep 17 00:00:00 2001
From: sunwen <sunwen@realcan.com>
Date: Fri, 2 Jun 2023 13:47:31 +0800
Subject: [PATCH] Add eigs function.

---
 .../solvingEquationSystem/TVAL/TVAL.cpp       | 40 +++++++++++++++++++
 .../solvingEquationSystem/TVAL/TVAL.h         |  2 +
 2 files changed, 42 insertions(+)

diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
index 60addf8..715ace1 100644
--- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
+++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
@@ -11,8 +11,48 @@
 #include <iostream>
 #include <vector>
 
+#include <Eigen/Sparse>
+#include <Spectra/GenEigsSolver.h>
+#include <Spectra/MatOp/SparseGenMatProd.h>
+#include <Eigen/Core>
+#include <Eigen/SparseCore>
 namespace Recon
 {
+
+    double eigs(Aurora::Sparse& aM)
+    {
+        double result = NAN;
+        size_t size = aM.getM();
+        if(size < aM.getN())
+        {
+            size = aM.getN();
+        }
+        Eigen::SparseMatrix<double> M(size, size);
+        std::vector<Eigen::Triplet<double>> triplets;
+        Aurora::Matrix rows = aM.getRowVector();
+        Aurora::Matrix columns = aM.getColVector();
+        Aurora::Matrix values = aM.getValVector();
+        for (int i = 0; i < rows.getDataSize(); ++i)
+        {
+            triplets.push_back(Eigen::Triplet<double>(rows[i], columns[i], values[i]));
+        }
+
+        M.setFromTriplets(triplets.begin(), triplets.end());
+        Spectra::SparseGenMatProd<double> op(M);
+        Spectra::GenEigsSolver<Spectra::SparseGenMatProd<double>> eigs(op, 1, 6);
+        eigs.init();
+        int nconv = eigs.compute(Spectra::SortRule::LargestMagn);
+        Eigen::VectorXcd evalues;
+        if(eigs.info() == Spectra::CompInfo::Successful)
+        {
+            evalues = eigs.eigenvalues();
+            std::complex<double> complex = evalues[0];
+            result = complex.real();
+        }
+
+        return result;
+    }
+
     void checkAndScale(Aurora::Sparse& M, Aurora::Matrix& b,size_t n){
         return;
         //TODO:暂时啥都没做
diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h
index a2fbf67..35a8630 100644
--- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h
+++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h
@@ -7,6 +7,8 @@ namespace Recon {
 
 void checkAndScale(Aurora::Sparse& M, Aurora::Matrix& b,size_t n);
 
+double eigs(Aurora::Sparse& aM);
+
 Aurora::Matrix callTval3(Aurora::Sparse &M, Aurora::Matrix &b,
                          const Aurora::Matrix &dims, int device,
                          struct TVALOptions &options);