From ab068405f4bf7c07c733e73a6c0cc3625cec1872 Mon Sep 17 00:00:00 2001 From: michael Date: Fri, 30 Apr 2021 23:03:19 +0200 Subject: [PATCH] Added computation of singular vectors --- CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx | 7 +- CDMATH/linearsolver/src/SparseMatrixPetsc.cxx | 74 ++++++++++++++++--- .../tests/cdmath/SparseMatrixPetscTests.cxx | 5 ++ 3 files changed, 71 insertions(+), 15 deletions(-) diff --git a/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx b/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx index d06cd26..37449c4 100755 --- a/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx +++ b/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx @@ -109,10 +109,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< 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< 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; double getConditionNumber(bool isSingular=false, double tol=1e-6) const; bool isSymmetric(double tol=1.e-6) const ; @@ -126,7 +127,7 @@ 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, SVDWhich which=SVD_SMALLEST , double tol=1e-6) const; + int computeSVD (int nsv, double ** valS, double ***vecS, SVDWhich which=SVD_SMALLEST , double tol=1e-6) 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 7041acc..58ab952 100755 --- a/CDMATH/linearsolver/src/SparseMatrixPetsc.cxx +++ b/CDMATH/linearsolver/src/SparseMatrixPetsc.cxx @@ -503,8 +503,8 @@ SparseMatrixPetsc::computeSpectrum(int nev, double ** valP, double ***vecP, EPSW MatAssemblyBegin(_mat,MAT_FINAL_ASSEMBLY); MatAssemblyEnd(_mat,MAT_FINAL_ASSEMBLY); - MatCreateVecs(_mat,NULL,&xr); - MatCreateVecs(_mat,NULL,&xi); + MatCreateVecs(_mat,&xr,NULL); + MatCreateVecs(_mat,&xi,NULL); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create the eigensolver and set various options @@ -616,19 +616,23 @@ SparseMatrixPetsc::computeSpectrum(int nev, double ** valP, double ***vecP, EPSW } int -SparseMatrixPetsc::computeSVD(int nsv, double ** valS, SVDWhich which, double tol) const +SparseMatrixPetsc::computeSVD(int nsv, double ** valS, double ***vecS, SVDWhich which, double tol) const { SVD svd; /* Singular value decomposition solver context */ SVDType type; PetscReal error; PetscScalar sigma; + Vec u,v; PetscInt i,maxit,its, nconv; SlepcInitialize(0, (char ***)"", PETSC_NULL, PETSC_NULL); MatAssemblyBegin(_mat,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(_mat,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd( _mat,MAT_FINAL_ASSEMBLY); + MatCreateVecs(_mat,&u,NULL); + MatCreateVecs(_mat,NULL,&v); + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create the singular value solver and set various options - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ @@ -677,7 +681,9 @@ SparseMatrixPetsc::computeSVD(int nsv, double ** valS, SVDWhich which, double to PetscPrintf(PETSC_COMM_WORLD," Number of converged singular values: %D\n\n",nconv); *valS=new double[nconv]; - + *vecS=new double * [2*nconv]; + double * myVecS; + if (nconv>0) { /* Display eigenvalues and relative errors @@ -690,19 +696,28 @@ SparseMatrixPetsc::computeSVD(int nsv, double ** valS, SVDWhich which, double to /* Get converged singular values: i-th eigenvalue is stored in valS */ - SVDGetSingularTriplet(svd,i,*valS+i,NULL,NULL); + SVDGetSingularTriplet(svd,i,*valS+i, u, v); /* Compute the relative error associated to each singular value */ SVDComputeError( svd, i, SVD_ERROR_RELATIVE, &error ); PetscPrintf(PETSC_COMM_WORLD," %12f %12g\n",(double)*(*valS+i),(double)error); + + VecGetArray(u,&myVecS); + *(*vecS + i)=new double [_numberOfRows]; + memcpy(*(*vecS + i),myVecS,_numberOfRows*sizeof(double)) ; + VecGetArray(v,&myVecS); + *(*vecS + i + nconv)=new double [_numberOfColumns]; + memcpy(*(*vecS + i + nconv),myVecS,_numberOfColumns*sizeof(double)) ; } PetscPrintf(PETSC_COMM_WORLD,"\n"); /* Free work space */ SVDDestroy(&svd); + VecDestroy(&u); + VecDestroy(&v); SlepcFinalize(); return nconv; @@ -713,6 +728,8 @@ SparseMatrixPetsc::computeSVD(int nsv, double ** valS, SVDWhich which, double to Free work space */ SVDDestroy(&svd); + VecDestroy(&u); + VecDestroy(&v); SlepcFinalize(); throw CdmathException("SparseMatrixPetsc::computeSVD : singular value decomposition failed"); @@ -834,6 +851,7 @@ SparseMatrixPetsc::getConditionNumber(bool isSingular, double tol) const { int nsv, nconv; double * valS; + double ** vecS; double sigma_max, sigma_min; /*** Lowest singular value, first check if matrix is singular ***/ @@ -842,19 +860,19 @@ SparseMatrixPetsc::getConditionNumber(bool isSingular, double tol) const else nsv=1 ; - nconv=computeSVD(nsv, &valS, SVD_SMALLEST, tol); + nconv=computeSVD(nsv, &valS, &vecS, SVD_SMALLEST, tol); if(nconv result(nconv); for (int i=0;i +SparseMatrixPetsc::getSingularVectors(int nsv, SVDWhich which, double tol) const +{ + int nconv; + double * valS; + double **vecS; + + nconv=computeSVD(nsv, &valS, &vecS, which, tol); + + std::vector< Vector > result(2*nconv); + + for (int i=0;i vp = A.getEigenvalues(2); CPPUNIT_ASSERT_DOUBLES_EQUAL( 0, vp[0],1.e-5); CPPUNIT_ASSERT_DOUBLES_EQUAL( 2, vp[1],1.e-5); @@ -146,6 +147,10 @@ SparseMatrixPetscTests::testClassSparseMatrixPetsc( void ) CPPUNIT_ASSERT_DOUBLES_EQUAL( 0, sigma[0],1.e-5); CPPUNIT_ASSERT_DOUBLES_EQUAL( 2, sigma[1],1.e-5); + std::vector< Vector > Vs = A.getSingularVectors(2); + CPPUNIT_ASSERT_DOUBLES_EQUAL( 0, abs(Vs[0][0]) - abs(Vs[0][1]),1.e-5); + CPPUNIT_ASSERT_DOUBLES_EQUAL( 0, abs(Vs[1][0]) - abs(Vs[1][1]),1.e-5); + /* Condition number of a symmetric non singular matrix */ A.setValue(0,0, 0.); A.setValue(0,1, 1.); -- 2.39.2