More options for eigen solvers

diff --git a/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx b/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx
index 9bb804970ab82f301c7f3bcc79a8b6603624f4d5..e412439b6681b9c4be5c5cfbdd8e2dcecfed72d5 100755 (executable)
--- 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 91412cdf858b17cf3edddbecdb8824ba7b42cfb6..30afd8ea9031ded0eec62399f404c32cac7c12c8 100755 (executable)
--- 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<nsv)
                        throw CdmathException("SparseMatrixPetsc::getConditionNumber could not find the smallest singular value");
            sigma_min=valS[nsv-1];
@@ -935,7 +934,7 @@ SparseMatrixPetsc::getConditionNumber(bool isSingular, double tol) const
            
            /*** Largest singular value ***/
            nsv=1;
-               nconv=computeSVD(nsv, &valS, &vecS, SVD_LARGEST, tol);
+               nconv=computeSVD(nsv, &valS, &vecS, SVD_LARGEST, tol,SVDCYCLIC);
                if(nconv<nsv)
                        throw CdmathException("SparseMatrixPetsc::getConditionNumber could not find the largest singular value");
            sigma_max=valS[nsv-1];
@@ -946,13 +945,13 @@ SparseMatrixPetsc::getConditionNumber(bool isSingular, double tol) const
 }
 
 std::vector< double > 
-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);