Add eigs function.

This commit is contained in:
sunwen
2023-06-02 13:47:31 +08:00
parent b5947c7181
commit 2e85a02ad5
2 changed files with 42 additions and 0 deletions

View File

@@ -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:暂时啥都没做

View File

@@ -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);