From d24ee6bdea4f04ca85a18b6860faa1affaaeb9d6 Mon Sep 17 00:00:00 2001 From: michael Date: Thu, 13 Jan 2022 17:53:27 +0100 Subject: [PATCH] Added functions to visualise and save pictures of matrices eigenvalues and non zero structure --- CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx | 11 +-- CDMATH/linearsolver/src/SparseMatrixPetsc.cxx | 76 ++++++++++++++++--- 2 files changed, 71 insertions(+), 16 deletions(-) diff --git a/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx b/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx index e412439..2296477 100755 --- a/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx +++ b/CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx @@ -101,7 +101,7 @@ public: friend SparseMatrixPetsc operator*(const SparseMatrixPetsc& M, const SparseMatrixPetsc& N) ; - void viewMatrix() const ; + void viewMatrix(bool useXWindow=false, double pause_lenght=0, std::string matrixName="") const ; //Save matrix coefficients into a file in ascii or binary mode void saveMatrix(std::string filename, bool binaryMode=false) const ; double getMatrixCoeff(int i, int j) const; @@ -114,10 +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, EPSType type = EPSKRYLOVSCHUR) const; + std::vector< double > getEigenvalues( int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6, EPSType type = EPSKRYLOVSCHUR, bool viewEigenvaluesInXWindows=false, double pause_lenght=0, std::string matrixName="") const; std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6, EPSType type = EPSKRYLOVSCHUR) const; + void plotEigenvalues(std::string matrixName="", int nev=-1, double pause_lenght=0, double tol=1e-6, EPSWhich which=EPS_SMALLEST_MAGNITUDE, 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< double > getSingularValues( int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6, SVDType type = SVDCYCLIC, bool viewSingularValuesInXWindows=false, double pause_lenght=0, std::string matrixName="") 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; @@ -131,8 +132,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, EPSType type = EPSKRYLOVSCHUR) const; - int computeSVD (int nsv, double ** valS, double ***vecS, SVDWhich which=SVD_SMALLEST , double tol=1e-6, SVDType type = SVDCYCLIC) const; + int computeSpectrum(int nev, double ** valP, double ***vecP, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6, EPSType type = EPSKRYLOVSCHUR, bool viewEigenvaluesInXWindows=false, double pause_lenght=0, std::string matrixName="") const; + int computeSVD (int nsv, double ** valS, double ***vecS, SVDWhich which=SVD_SMALLEST , double tol=1e-6, SVDType type = SVDCYCLIC, bool viewSingularValuesInXWindows=false, double pause_lenght=0, std::string matrixName="") 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 30afd8e..e663b7c 100755 --- a/CDMATH/linearsolver/src/SparseMatrixPetsc.cxx +++ b/CDMATH/linearsolver/src/SparseMatrixPetsc.cxx @@ -406,12 +406,23 @@ SparseMatrixPetsc::operator/= (double value) } void -SparseMatrixPetsc::viewMatrix() const +SparseMatrixPetsc::viewMatrix(bool useXWindow, double pause_lenght, std::string matrixName) const { MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY); - MatView(_mat,PETSC_VIEWER_STDOUT_SELF); + if(!useXWindow) + MatView(_mat,PETSC_VIEWER_STDOUT_WORLD); + else + { + PetscViewer viewer = PETSC_VIEWER_DRAW_WORLD; + PetscDraw draw; + PetscViewerDrawGetDraw(viewer, 0, &draw); + PetscDrawSetPause(draw, pause_lenght); // Wait for user + std::string filename="MatrixNonZeroPlot_"+matrixName+".ppm"; + PetscDrawSetSave( draw, filename.c_str()); + MatView(_mat, viewer); + } } void @@ -536,7 +547,7 @@ bool SparseMatrixPetsc::isSymmetric(double tol) const } int -SparseMatrixPetsc::computeSpectrum(int nev, double ** valP, double ***vecP, EPSWhich which, double tol, EPSType type) const +SparseMatrixPetsc::computeSpectrum(int nev, double ** valP, double ***vecP, EPSWhich which, double tol, EPSType type, bool viewEigenvaluesInXWindows, double pause_lenght, std::string matrixName) const { EPS eps; /* eigenproblem solver context */ PetscReal error; @@ -582,8 +593,20 @@ SparseMatrixPetsc::computeSpectrum(int nev, double ** valP, double ***vecP, EPSW - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ EPSSolve(eps); + + if(viewEigenvaluesInXWindows) + { + PetscViewer viewer = PETSC_VIEWER_DRAW_WORLD; + PetscDraw draw; + PetscViewerDrawGetDraw(viewer, 0, &draw); + PetscDrawSetPause(draw, pause_lenght); // time duration of the display. if pause_lenght = -1 then wait for user to press a key + std::string filename="MatrixEigenvaluePlot_"+matrixName+".ppm"; + PetscDrawSetSave( draw, filename.c_str()); + EPSValuesView(eps, viewer); + } + /* - Optional: Get some information from the solver and display it + Optional: Get some information from the solver and display it on the screen */ EPSGetIterationNumber(eps,&its); PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its); @@ -672,7 +695,7 @@ SparseMatrixPetsc::computeSpectrum(int nev, double ** valP, double ***vecP, EPSW } int -SparseMatrixPetsc::computeSVD(int nsv, double ** valS, double ***vecS, SVDWhich which, double tol, SVDType type) const +SparseMatrixPetsc::computeSVD(int nsv, double ** valS, double ***vecS, SVDWhich which, double tol, SVDType type, bool viewSingularValuesInXWindows, double pause_lenght, std::string matrixName) const { SVD svd; /* Singular value decomposition solver context */ PetscReal error; @@ -720,6 +743,18 @@ SparseMatrixPetsc::computeSVD(int nsv, double ** valS, double ***vecS, SVDWhich - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ SVDSolve(svd); + + if(viewSingularValuesInXWindows) + { + PetscViewer viewer = PETSC_VIEWER_DRAW_WORLD; + PetscDraw draw; + PetscViewerDrawGetDraw(viewer, 0, &draw); + PetscDrawSetPause(draw, pause_lenght); // time duration of the display. If pause_lenght = -1 then wait for user to press a key + std::string filename="MatrixSingularValuePlot_"+matrixName+".ppm"; + PetscDrawSetSave( draw, filename.c_str()); + SVDValuesView(svd, viewer); + } + /* Optional: Get some information from the solver and display it */ @@ -803,13 +838,13 @@ SparseMatrixPetsc::computeSVD(int nsv, double ** valS, double ***vecS, SVDWhich } std::vector< double > -SparseMatrixPetsc::getEigenvalues(int nev, EPSWhich which, double tol, EPSType type) const +SparseMatrixPetsc::getEigenvalues(int nev, EPSWhich which, double tol, EPSType type, bool viewEigenvaluesInXWindows, double pause_lenght, std::string matrixName) const { int nconv; double * valP; double **vecP; - nconv=computeSpectrum(nev, &valP, &vecP, which, tol,type); + nconv=computeSpectrum(nev, &valP, &vecP, which, tol,type, viewEigenvaluesInXWindows, pause_lenght, matrixName); std::vector< double > result(nconv); @@ -824,6 +859,23 @@ SparseMatrixPetsc::getEigenvalues(int nev, EPSWhich which, double tol, EPSType t return result; } +void +SparseMatrixPetsc::plotEigenvalues(std::string matrixName, int nev, double pause_lenght, double tol, EPSWhich which, EPSType type) const +{ + double * valP; + double **vecP; + + if(nev <=0) + computeSpectrum(_numberOfRows, &valP, &vecP, which, tol,type, true, pause_lenght, matrixName); + else + computeSpectrum(nev, &valP, &vecP, which, tol,type, true, pause_lenght, matrixName); + + delete[] valP; + for (int i=0;i<_numberOfRows;i++) + delete[] vecP[i]; + delete[] vecP; +} + std::vector< Vector > SparseMatrixPetsc::getEigenvectors(int nev, EPSWhich which, double tol, EPSType type) const { @@ -926,7 +978,7 @@ SparseMatrixPetsc::getConditionNumber(bool isSingular, double tol) const else nsv=1 ; - nconv=computeSVD(nsv, &valS, &vecS, SVD_SMALLEST, tol,SVDCYCLIC); + nconv=computeSVD(nsv, &valS, &vecS, SVD_SMALLEST, tol,SVDCROSS); if(nconv -SparseMatrixPetsc::getSingularValues(int nsv, SVDWhich which, double tol, SVDType type) const +SparseMatrixPetsc::getSingularValues(int nsv, SVDWhich which, double tol, SVDType type, bool viewSingularValuesInXWindows, double pause_lenght, std::string matrixName) const { int nconv; double * valS; double **vecS; - nconv=computeSVD(nsv, &valS, &vecS, which, tol, type); + nconv=computeSVD(nsv, &valS, &vecS, which, tol, type, viewSingularValuesInXWindows, pause_lenght, matrixName); std::vector< double > result(nconv); @@ -994,6 +1046,8 @@ SparseMatrixPetsc::getSingularVectors(int nsv, SVDWhich which, double tol, SVDTy return result; } + + void SparseMatrixPetsc::leftDiagonalScale(Vector v) { -- 2.39.2