From 7346cde0839867c658ebc6fee1187c96c5aba967 Mon Sep 17 00:00:00 2001 From: michael Date: Fri, 10 Dec 2021 14:13:20 +0100 Subject: [PATCH] More options for eigen solvers --- CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx | 14 ++++----- CDMATH/linearsolver/src/SparseMatrixPetsc.cxx | 31 +++++++++---------- 2 files changed, 22 insertions(+), 23 deletions(-) diff --git a/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx b/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx index 9bb8049..e412439 100755 --- a/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx +++ b/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx @@ -114,11 +114,11 @@ public: void diagonalShift(double lambda); void zeroEntries();//sets the matrix coefficients to zero - std::vector< double > getEigenvalues( int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const; - std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const; - MEDCoupling::DataArrayDouble * getEigenvectorsDataArrayDouble(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const; - std::vector< double > getSingularValues( int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const; - std::vector< Vector > getSingularVectors(int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const; + std::vector< double > getEigenvalues( int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6, EPSType type = EPSKRYLOVSCHUR) const; + std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6, EPSType type = EPSKRYLOVSCHUR) const; + MEDCoupling::DataArrayDouble * getEigenvectorsDataArrayDouble(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6, EPSType type = EPSKRYLOVSCHUR) const; + std::vector< double > getSingularValues( int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6, SVDType type = SVDCYCLIC) const; + std::vector< Vector > getSingularVectors(int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6, SVDType type = SVDCYCLIC) const; double getConditionNumber(bool isSingular=false, double tol=1e-6) const; bool isSymmetric(double tol=1.e-6) const ; @@ -131,8 +131,8 @@ private: int _numberOfNonZeros ;//The maximum number of nonzeros coefficients per line (or an upper bound) - int computeSpectrum(int nev, double ** valP, double ***vecP, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const; - int computeSVD (int nsv, double ** valS, double ***vecS, SVDWhich which=SVD_SMALLEST , double tol=1e-6) const; + int computeSpectrum(int nev, double ** valP, double ***vecP, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6, EPSType type = EPSKRYLOVSCHUR) const; + int computeSVD (int nsv, double ** valS, double ***vecS, SVDWhich which=SVD_SMALLEST , double tol=1e-6, SVDType type = SVDCYCLIC) const; Vector vecToVector(const Vec& vec) const ; Vec vectorToVec( const Vector& myVector ) const ; diff --git a/CDMATH/linearsolver/src/SparseMatrixPetsc.cxx b/CDMATH/linearsolver/src/SparseMatrixPetsc.cxx index 91412cd..30afd8e 100755 --- a/CDMATH/linearsolver/src/SparseMatrixPetsc.cxx +++ b/CDMATH/linearsolver/src/SparseMatrixPetsc.cxx @@ -536,10 +536,9 @@ bool SparseMatrixPetsc::isSymmetric(double tol) const } int -SparseMatrixPetsc::computeSpectrum(int nev, double ** valP, double ***vecP, EPSWhich which, double tol) const +SparseMatrixPetsc::computeSpectrum(int nev, double ** valP, double ***vecP, EPSWhich which, double tol, EPSType type) const { EPS eps; /* eigenproblem solver context */ - EPSType type; PetscReal error; PetscScalar kr,ki; Vec xr,xi; @@ -572,7 +571,7 @@ SparseMatrixPetsc::computeSpectrum(int nev, double ** valP, double ***vecP, EPSW EPSSetWhichEigenpairs(eps,which); EPSSetDimensions(eps,nev,PETSC_DEFAULT,PETSC_DEFAULT); EPSSetTolerances(eps,tol,PETSC_DEFAULT); - + EPSSetType(eps,type); /* Set solver parameters at runtime */ @@ -673,10 +672,9 @@ SparseMatrixPetsc::computeSpectrum(int nev, double ** valP, double ***vecP, EPSW } int -SparseMatrixPetsc::computeSVD(int nsv, double ** valS, double ***vecS, SVDWhich which, double tol) const +SparseMatrixPetsc::computeSVD(int nsv, double ** valS, double ***vecS, SVDWhich which, double tol, SVDType type) const { SVD svd; /* Singular value decomposition solver context */ - SVDType type; PetscReal error; PetscScalar sigma; Vec u,v; @@ -709,6 +707,7 @@ SparseMatrixPetsc::computeSVD(int nsv, double ** valS, double ***vecS, SVDWhich SVDSetWhichSingularTriplets(svd,which); SVDSetDimensions(svd,nsv,PETSC_DEFAULT,PETSC_DEFAULT); + SVDSetType(svd, type); SVDSetTolerances(svd,tol,PETSC_DEFAULT); /* @@ -804,13 +803,13 @@ SparseMatrixPetsc::computeSVD(int nsv, double ** valS, double ***vecS, SVDWhich } std::vector< double > -SparseMatrixPetsc::getEigenvalues(int nev, EPSWhich which, double tol) const +SparseMatrixPetsc::getEigenvalues(int nev, EPSWhich which, double tol, EPSType type) const { int nconv; double * valP; double **vecP; - nconv=computeSpectrum(nev, &valP, &vecP, which, tol); + nconv=computeSpectrum(nev, &valP, &vecP, which, tol,type); std::vector< double > result(nconv); @@ -826,13 +825,13 @@ SparseMatrixPetsc::getEigenvalues(int nev, EPSWhich which, double tol) const } std::vector< Vector > -SparseMatrixPetsc::getEigenvectors(int nev, EPSWhich which, double tol) const +SparseMatrixPetsc::getEigenvectors(int nev, EPSWhich which, double tol, EPSType type) const { int nconv; double * valP; double **vecP; - nconv=computeSpectrum(nev, &valP, &vecP, which, tol); + nconv=computeSpectrum(nev, &valP, &vecP, which, tol,type); std::vector< Vector > result(nconv); @@ -853,7 +852,7 @@ SparseMatrixPetsc::getEigenvectors(int nev, EPSWhich which, double tol) const } MEDCoupling::DataArrayDouble * -SparseMatrixPetsc::getEigenvectorsDataArrayDouble(int nev, EPSWhich which, double tol) const +SparseMatrixPetsc::getEigenvectorsDataArrayDouble(int nev, EPSWhich which, double tol, EPSType type) const { int nconv; double * valP; @@ -927,7 +926,7 @@ SparseMatrixPetsc::getConditionNumber(bool isSingular, double tol) const else nsv=1 ; - nconv=computeSVD(nsv, &valS, &vecS, SVD_SMALLEST, tol); + nconv=computeSVD(nsv, &valS, &vecS, SVD_SMALLEST, tol,SVDCYCLIC); if(nconv -SparseMatrixPetsc::getSingularValues(int nsv, SVDWhich which, double tol) const +SparseMatrixPetsc::getSingularValues(int nsv, SVDWhich which, double tol, SVDType type) const { int nconv; double * valS; double **vecS; - nconv=computeSVD(nsv, &valS, &vecS, which, tol); + nconv=computeSVD(nsv, &valS, &vecS, which, tol, type); std::vector< double > result(nconv); @@ -965,13 +964,13 @@ SparseMatrixPetsc::getSingularValues(int nsv, SVDWhich which, double tol) const } std::vector< Vector > -SparseMatrixPetsc::getSingularVectors(int nsv, SVDWhich which, double tol) const +SparseMatrixPetsc::getSingularVectors(int nsv, SVDWhich which, double tol, SVDType type) const { int nconv; double * valS; double **vecS; - nconv=computeSVD(nsv, &valS, &vecS, which, tol); + nconv=computeSVD(nsv, &valS, &vecS, which, tol, type); std::vector< Vector > result(2*nconv); -- 2.39.2