From 999753cd7693619d02205c8616187eac5eb7e6d6 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Wed, 24 Jul 2013 10:11:06 +0000
Subject: [PATCH] linearSystemPETSc : implement linearSystemPETScBlockDouble as
a template spetialization of linearSystemPETSc to avoid code duplication
---
Solver/linearSystemPETSc.cpp | 280 +++++------------------------------
Solver/linearSystemPETSc.h | 97 +++++-------
Solver/linearSystemPETSc.hpp | 66 ++++++---
wrappers/gmshpy/gmshSolver.i | 1 +
4 files changed, 115 insertions(+), 329 deletions(-)
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index 8f27fba74a..2516b5de5d 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -18,28 +18,18 @@ template class linearSystemPETSc<double>;
template class linearSystemPETSc<std::complex<double> >;
#endif
-void linearSystemPETScBlockDouble::_kspCreate()
+template<>
+int linearSystemPETSc<fullMatrix<double> >::_getBlockSizeFromParameters() const
{
- KSPCreate(_sequential ? PETSC_COMM_SELF : PETSC_COMM_WORLD, &_ksp);
- if (this->_parameters.count("petscPrefix"))
- KSPAppendOptionsPrefix(_ksp, this->_parameters["petscPrefix"].c_str());
- KSPSetFromOptions(_ksp);
- _kspAllocated = true;
+ if ( _parameters.find("blockSize") == _parameters.end())
+ Msg::Error ("'blockSize' parameters must be set for linearSystemPETScBlock");
+ int blockSize = strtol(_parameters.find("blockSize")->second.c_str(), NULL, 10);
+ return blockSize;
}
-linearSystemPETScBlockDouble::~linearSystemPETScBlockDouble()
-{
- if (_isAllocated) {
- MatDestroy(&_a);
- VecDestroy(&_b);
- VecDestroy(&_x);
- }
- if (_kspAllocated) {
- KSPDestroy(&_ksp);
- }
-}
-void linearSystemPETScBlockDouble::addToMatrix(int row, int col,
+template<>
+void linearSystemPETSc<fullMatrix<double> >::addToMatrix(int row, int col,
const fullMatrix<double> &val)
{
if (!_entriesPreAllocated)
@@ -73,30 +63,30 @@ void linearSystemPETScBlockDouble::addToMatrix(int row, int col,
modval(jj, ii) = buff;
}
#endif
- _matrixModified=true;
+ _valuesNotAssembled = true;
}
-void linearSystemPETScBlockDouble::addToRightHandSide(int row,
+template<>
+void linearSystemPETSc<fullMatrix<double> >::addToRightHandSide(int row,
const fullMatrix<double> &val)
{
- for (int ii = 0; ii < _blockSize; ii++) {
- PetscInt i = row * _blockSize + ii;
+ int blockSize;
+ _try(MatGetBlockSize(_a, &blockSize));
+ for (int ii = 0; ii < blockSize; ii++) {
+ PetscInt i = row * blockSize + ii;
PetscScalar v = val(ii, 0);
VecSetValues(_b, 1, &i, &v, ADD_VALUES);
}
}
-void linearSystemPETScBlockDouble::getFromMatrix(int row, int col,
- fullMatrix<double> &val ) const
-{
- Msg::Error("getFromMatrix not implemented for PETSc");
-}
-
-void linearSystemPETScBlockDouble::getFromRightHandSide(int row,
+template<>
+void linearSystemPETSc<fullMatrix<double > >::getFromRightHandSide(int row,
fullMatrix<double> &val) const
{
- for (int i = 0; i < _blockSize; i++) {
- int ii = row*_blockSize +i;
+ int blockSize;
+ _try(MatGetBlockSize(_a, &blockSize));
+ for (int i = 0; i < blockSize; i++) {
+ int ii = row*blockSize +i;
#ifdef PETSC_USE_COMPLEX
PetscScalar s;
VecGetValues ( _b, 1, &ii, &s);
@@ -107,19 +97,25 @@ void linearSystemPETScBlockDouble::getFromRightHandSide(int row,
}
}
-void linearSystemPETScBlockDouble::addToSolution(int row, const fullMatrix<double> &val)
+template<>
+void linearSystemPETSc<fullMatrix<double> >::addToSolution(int row, const fullMatrix<double> &val)
{
- for (int ii = 0; ii < _blockSize; ii++) {
- PetscInt i = row * _blockSize + ii;
+ int blockSize;
+ _try(MatGetBlockSize(_a, &blockSize));
+ for (int ii = 0; ii < blockSize; ii++) {
+ PetscInt i = row * blockSize + ii;
PetscScalar v = val(ii, 0);
VecSetValues(_x, 1, &i, &v, ADD_VALUES);
}
}
-void linearSystemPETScBlockDouble::getFromSolution(int row, fullMatrix<double> &val) const
+template<>
+void linearSystemPETSc<fullMatrix<double> >::getFromSolution(int row, fullMatrix<double> &val) const
{
- for (int i = 0; i < _blockSize; i++) {
- int ii = row*_blockSize +i;
+ int blockSize;
+ _try(MatGetBlockSize(_a, &blockSize));
+ for (int i = 0; i < blockSize; i++) {
+ int ii = row*blockSize +i;
#ifdef PETSC_USE_COMPLEX
PetscScalar s;
VecGetValues ( _x, 1, &ii, &s);
@@ -130,213 +126,5 @@ void linearSystemPETScBlockDouble::getFromSolution(int row, fullMatrix<double> &
}
}
-
-void linearSystemPETScBlockDouble::allocate(int nbRows)
-{
- MPI_Comm comm = _sequential ? PETSC_COMM_SELF: PETSC_COMM_WORLD;
- if (this->_parameters.count("petscOptions"))
- PetscOptionsInsertString(this->_parameters["petscOptions"].c_str());
- _blockSize = strtol (_parameters["blockSize"].c_str(), NULL, 10);
- if (_blockSize == 0)
- Msg::Error ("'blockSize' parameters must be set for linearSystemPETScBlock");
- clear();
- MatCreate(comm, &_a);
- MatSetSizes(_a,nbRows * _blockSize, nbRows * _blockSize, PETSC_DETERMINE, PETSC_DETERMINE);
- if (Msg::GetCommSize() > 1 && !_sequential) {
- MatSetType(_a, MATMPIBAIJ);
- }
- else {
- MatSetType(_a, MATSEQBAIJ);
- }
- if (_parameters.count("petscPrefix"))
- MatAppendOptionsPrefix(_a, _parameters["petscPrefix"].c_str());
- MatSetFromOptions(_a);
- //since PETSc 3.3 GetOwnershipRange and MatGetSize() cannot be called before SetPreallocation
- _localSize = nbRows;
- _localRowStart = 0;
- _localRowEnd = nbRows;
- _globalSize = _localSize;
- #ifdef HAVE_MPI
- if (!_sequential) {
- _localRowStart = 0;
- if (Msg::GetCommRank() != 0) {
- MPI_Status status;
- MPI_Recv((void*)&_localRowStart, 1, MPI_INT, Msg::GetCommRank() - 1, 1, MPI_COMM_WORLD, &status);
- }
- _localRowEnd = _localRowStart + nbRows;
- if (Msg::GetCommRank() != Msg::GetCommSize() - 1) {
- MPI_Send((void*)&_localRowEnd, 1, MPI_INT, Msg::GetCommRank() + 1, 1, MPI_COMM_WORLD);
- }
- MPI_Allreduce((void*)&_localSize, (void*)&_globalSize, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
- }
- #endif
- // override the default options with the ones from the option
- // database (if any)
- VecCreate(comm, &_x);
- VecSetSizes(_x, nbRows * _blockSize, PETSC_DETERMINE);
- // override the default options with the ones from the option
- // database (if any)
- if (_parameters.count("petscPrefix"))
- VecAppendOptionsPrefix(_x, _parameters["petscPrefix"].c_str());
- VecSetFromOptions(_x);
- VecDuplicate(_x, &_b);
- _isAllocated = true;
-}
-
-bool linearSystemPETScBlockDouble::isAllocated() const
-{
- return _isAllocated;
-}
-
-void linearSystemPETScBlockDouble::clear()
-{
- if(_isAllocated){
- MatDestroy(&_a);
- VecDestroy(&_x);
- VecDestroy(&_b);
- }
- _isAllocated = false;
-}
-
-void linearSystemPETScBlockDouble::print()
-{
- _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
- _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
- _try(VecAssemblyBegin(_b));
- _try(VecAssemblyEnd(_b));
- if(Msg::GetCommRank()==0)
- printf("a :\n");
- MatView(_a, PETSC_VIEWER_STDOUT_WORLD);
- if(Msg::GetCommRank()==0)
- printf("b :\n");
- VecView(_b, PETSC_VIEWER_STDOUT_WORLD);
- if(Msg::GetCommRank()==0)
- printf("x :\n");
- VecView(_x, PETSC_VIEWER_STDOUT_WORLD);
-}
-
-
-int linearSystemPETScBlockDouble::systemSolve()
-{
- if (!_kspAllocated)
- _kspCreate();
- if (!_matrixModified || _parameters["matrix_reuse"] == "same_matrix")
- KSPSetOperators(_ksp, _a, _a, SAME_PRECONDITIONER);
- else if (_parameters["matrix_reuse"] == "same_sparsity")
- KSPSetOperators(_ksp, _a, _a, SAME_NONZERO_PATTERN);
- else
- KSPSetOperators(_ksp, _a, _a, DIFFERENT_NONZERO_PATTERN);
- MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY);
- _matrixModified=false;
- VecAssemblyBegin(_b);
- VecAssemblyEnd(_b);
- KSPSolve(_ksp, _b, _x);
- return 1;
-}
-
-void linearSystemPETScBlockDouble::insertInSparsityPattern (int i, int j)
-{
- i -= _localRowStart;
- if (i<0 || i>= _localSize) return;
- _sparsity.insertEntry (i,j);
-}
-
-void linearSystemPETScBlockDouble::preAllocateEntries()
-{
- if (_entriesPreAllocated) return;
- if (!_isAllocated) Msg::Fatal("system must be allocated first");
- if (_sparsity.getNbRows() == 0) {
- PetscInt prealloc = 300;
- PetscBool set;
- PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set);
- if (_blockSize == 0) {
- MatSeqAIJSetPreallocation(_a, prealloc, PETSC_NULL);
- } else {
- MatSeqBAIJSetPreallocation(_a, _blockSize, 5, PETSC_NULL);
- }
- } else {
- std::vector<int> nByRowDiag (_localSize), nByRowOffDiag (_localSize);
- for (int i = 0; i < _localSize; i++) {
- int n;
- const int *r = _sparsity.getRow(i, n);
- for (int j = 0; j < n; j++) {
- if (r[j] >= _localRowStart && r[j] < _localRowEnd)
- nByRowDiag[i] ++;
- else
- nByRowOffDiag[i] ++;
- }
- }
- if (_blockSize == 0) {
- MatSeqAIJSetPreallocation(_a, 0, &nByRowDiag[0]);
- MatMPIAIJSetPreallocation(_a, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]);
- } else {
- MatSeqBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0]);
- MatMPIBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]);
- }
- _sparsity.clear();
- }
- _entriesPreAllocated = true;
-}
-
-void linearSystemPETScBlockDouble::zeroMatrix()
-{
- if (_isAllocated && _entriesPreAllocated) {
- MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY);
- MatZeroEntries(_a);
- }
-}
-
-void linearSystemPETScBlockDouble::zeroRightHandSide()
-{
- if (_isAllocated) {
- VecAssemblyBegin(_b);
- VecAssemblyEnd(_b);
- VecZeroEntries(_b);
- }
-}
-
-void linearSystemPETScBlockDouble::zeroSolution()
-{
- if (_isAllocated) {
- VecAssemblyBegin(_x);
- VecAssemblyEnd(_x);
- VecZeroEntries(_x);
- }
-}
-
-linearSystemPETScBlockDouble::linearSystemPETScBlockDouble(bool sequential)
-{
- _entriesPreAllocated = false;
- _isAllocated = false;
- _kspAllocated = false;
- _matrixModified=true;
- _sequential = sequential;
-}
-
-double linearSystemPETScBlockDouble::normInfRightHandSide() const
-{
- PetscReal nor;
- VecAssemblyBegin(_b);
- VecAssemblyEnd(_b);
- VecNorm(_b, NORM_INFINITY, &nor);
- return nor;
-}
-
-void linearSystemPETScBlockDouble::printMatlab(const char *filename) const
-{
- _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
- _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
- _try(VecAssemblyBegin(_b));
- _try(VecAssemblyEnd(_b));
-
- PetscViewer viewer;
- PetscViewerASCIIOpen(PETSC_COMM_WORLD, filename, &viewer);
- PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
- MatView(_a, viewer);
- PetscViewerDestroy(&viewer);
- return;
-}
-
+template class linearSystemPETSc<fullMatrix<double> >;
#endif // HAVE_PETSC
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index ea8579a35b..5cd933d066 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -65,7 +65,6 @@ typedef struct _p_KSP* KSP;
template <class scalar>
class linearSystemPETSc : public linearSystem<scalar> {
protected:
- int _blockSize; // for block Matrix
bool _isAllocated, _kspAllocated, _entriesPreAllocated;
bool _matrixChangedSinceLastSolve;
bool _valuesNotAssembled; //cannot use MatAssembled since MatAssembled return false for an empty matrix
@@ -79,25 +78,27 @@ class linearSystemPETSc : public linearSystem<scalar> {
#ifndef SWIG
MPI_Comm _comm;
#endif
+ int _getBlockSizeFromParameters() const;
public:
- virtual ~linearSystemPETSc();
+ ~linearSystemPETSc();
void insertInSparsityPattern (int i, int j);
- virtual bool isAllocated() const { return _isAllocated; }
- virtual void preAllocateEntries();
- virtual void allocate(int nbRows);
+ bool isAllocated() const { return _isAllocated; }
+ void preAllocateEntries();
+ void allocate(int nbRows);
void print();
- virtual void clear();
- virtual void getFromMatrix(int row, int col, scalar &val) const;
- virtual void addToRightHandSide(int row, const scalar &val);
- virtual void getFromRightHandSide(int row, scalar &val) const;
- virtual double normInfRightHandSide() const;
- virtual void addToMatrix(int row, int col, const scalar &val);
- virtual void addToSolution(int row, const scalar &val);
- virtual void getFromSolution(int row, scalar &val) const;
- virtual void zeroMatrix();
- virtual void zeroRightHandSide();
- virtual void zeroSolution();
- virtual int systemSolve();
+ void clear();
+ void getFromMatrix(int row, int col, scalar &val) const;
+ void addToRightHandSide(int row, const scalar &val);
+ void getFromRightHandSide(int row, scalar &val) const;
+ double normInfRightHandSide() const;
+ void addToMatrix(int row, int col, const scalar &val);
+ void addToSolution(int row, const scalar &val);
+ void getFromSolution(int row, scalar &val) const;
+ void zeroMatrix();
+ void zeroRightHandSide();
+ void zeroSolution();
+ void printMatlab(const char *filename) const;
+ int systemSolve();
Mat &getMatrix(){ return _a; }
//std::vector<scalar> getData();
//std::vector<int> getRowPointers();
@@ -108,39 +109,6 @@ class linearSystemPETSc : public linearSystem<scalar> {
#endif
linearSystemPETSc();
};
-
-class linearSystemPETScBlockDouble : public linearSystem<fullMatrix<double> > {
- bool _entriesPreAllocated, _isAllocated, _kspAllocated, _matrixModified;
- sparsityPattern _sparsity;
- Mat _a;
- Vec _b, _x;
- KSP _ksp;
- int _blockSize, _localRowStart, _localRowEnd, _globalSize, _localSize;
- bool _sequential;
- public:
- void _kspCreate();
- void print();
- void printMatlab(const char *filename) const;
- virtual void addToMatrix(int row, int col, const fullMatrix<double> &val);
- virtual void addToRightHandSide(int row, const fullMatrix<double> &val);
- virtual void addToSolution(int row, const fullMatrix<double> &val);
- virtual void getFromMatrix(int row, int col, fullMatrix<double> &val ) const;
- virtual void getFromRightHandSide(int row, fullMatrix<double> &val) const;
- virtual void getFromSolution(int row, fullMatrix<double> &val) const;
- void allocate(int nbRows);
- bool isAllocated() const;
- int systemSolve();
- void preAllocateEntries();
- void clear();
- void zeroMatrix();
- void zeroRightHandSide();
- void zeroSolution();
- double normInfRightHandSide() const;
- void insertInSparsityPattern (int i, int j);
- linearSystemPETScBlockDouble(bool sequential = false);
- ~linearSystemPETScBlockDouble();
-};
-
#else
template <class scalar>
@@ -150,20 +118,21 @@ class linearSystemPETSc : public linearSystem<scalar> {
{
Msg::Error("PETSc is not available in this version of Gmsh");
}
- virtual bool isAllocated() const { return false; }
- virtual void allocate(int nbRows) {}
- virtual void clear(){}
- virtual void addToMatrix(int row, int col, const scalar &val) {}
- virtual void getFromMatrix(int row, int col, scalar &val) const {}
- virtual void addToRightHandSide(int row, const scalar &val) {}
- virtual void addToSolution(int row, const scalar &val) {}
- virtual void getFromRightHandSide(int row, scalar &val) const {}
- virtual void getFromSolution(int row, scalar &val) const {}
- virtual void zeroMatrix() {}
- virtual void zeroRightHandSide() {}
- virtual void zeroSolution() {}
- virtual int systemSolve() { return 0; }
- virtual double normInfRightHandSide() const{return 0;}
+ bool isAllocated() const { return false; }
+ void allocate(int nbRows) {}
+ void clear(){}
+ void addToMatrix(int row, int col, const scalar &val) {}
+ void getFromMatrix(int row, int col, scalar &val) const {}
+ void addToRightHandSide(int row, const scalar &val) {}
+ void addToSolution(int row, const scalar &val) {}
+ void getFromRightHandSide(int row, scalar &val) const {}
+ void getFromSolution(int row, scalar &val) const {}
+ void zeroMatrix() {}
+ void zeroRightHandSide() {}
+ void zeroSolution() {}
+ void printMatlab(const char *filename) const{};
+ int systemSolve() { return 0; }
+ double normInfRightHandSide() const{return 0;}
};
#endif
#endif
diff --git a/Solver/linearSystemPETSc.hpp b/Solver/linearSystemPETSc.hpp
index 81024fadf3..5f7098f2a7 100644
--- a/Solver/linearSystemPETSc.hpp
+++ b/Solver/linearSystemPETSc.hpp
@@ -9,6 +9,12 @@ static void _try(int ierr)
CHKERRABORT(PETSC_COMM_WORLD, ierr);
}
+template <class scalar>
+int linearSystemPETSc<scalar>::_getBlockSizeFromParameters() const
+{
+ return 1;
+}
+
template <class scalar>
void linearSystemPETSc<scalar>::_kspCreate()
{
@@ -34,10 +40,10 @@ linearSystemPETSc<scalar>::linearSystemPETSc(MPI_Comm com)
{
_comm = com;
_isAllocated = false;
- _blockSize = 0;
_kspAllocated = false;
_matrixChangedSinceLastSolve = true;
_valuesNotAssembled = false;
+ _entriesPreAllocated = false;
}
template <class scalar>
@@ -45,10 +51,10 @@ linearSystemPETSc<scalar>::linearSystemPETSc()
{
_comm = PETSC_COMM_WORLD;
_isAllocated = false;
- _blockSize = 0;
_kspAllocated = false;
_matrixChangedSinceLastSolve = true;
_valuesNotAssembled = false;
+ _entriesPreAllocated = false;
}
template <class scalar>
@@ -84,17 +90,15 @@ void linearSystemPETSc<scalar>::preAllocateEntries()
{
if (_entriesPreAllocated) return;
if (!_isAllocated) Msg::Fatal("system must be allocated first");
+ int blockSize = _getBlockSizeFromParameters();
+ std::vector<int> nByRowDiag (_localSize), nByRowOffDiag (_localSize);
if (_sparsity.getNbRows() == 0) {
PetscInt prealloc = 300;
PetscBool set;
PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set);
- if (_blockSize == 0) {
- _try(MatSeqAIJSetPreallocation(_a, prealloc, PETSC_NULL));
- } else {
- _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 5, PETSC_NULL));
- }
+ nByRowDiag.resize(0);
+ nByRowDiag.resize(_localSize, prealloc);
} else {
- std::vector<int> nByRowDiag (_localSize), nByRowOffDiag (_localSize);
for (int i = 0; i < _localSize; i++) {
int n;
const int *r = _sparsity.getRow(i, n);
@@ -105,24 +109,34 @@ void linearSystemPETSc<scalar>::preAllocateEntries()
nByRowOffDiag[i] ++;
}
}
- if (_blockSize == 0) {
- _try(MatSeqAIJSetPreallocation(_a, 0, &nByRowDiag[0]));
- _try(MatMPIAIJSetPreallocation(_a, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]));
- } else {
- _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0]));
- _try(MatMPIBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]));
- }
_sparsity.clear();
}
+ _try(MatXAIJSetPreallocation(_a, blockSize, &nByRowDiag[0], &nByRowOffDiag[0], NULL, NULL));
_entriesPreAllocated = true;
}
template <class scalar>
void linearSystemPETSc<scalar>::allocate(int nbRows)
{
+ #ifdef HAVE_MPI
+ PetscMPIInt commSize;
+ MPI_Comm_size(_comm,&commSize);
+ #endif
+ int blockSize = _getBlockSizeFromParameters();
clear();
_try(MatCreate(_comm, &_a));
- _try(MatSetSizes(_a, nbRows, nbRows, PETSC_DETERMINE, PETSC_DETERMINE));
+ _try(MatSetSizes(_a, blockSize * nbRows, blockSize * nbRows, PETSC_DETERMINE, PETSC_DETERMINE));
+ if (blockSize > 1) {
+ #ifdef HAVE_MPI
+ if (commSize > 1) {
+ MatSetType(_a, MATMPIBAIJ);
+ }
+ else
+ #endif
+ {
+ MatSetType(_a, MATSEQBAIJ);
+ }
+ }
// override the default options with the ones from the option
// database (if any)
if (this->_parameters.count("petscOptions"))
@@ -133,8 +147,6 @@ void linearSystemPETSc<scalar>::allocate(int nbRows)
//since PETSc 3.3 GetOwnershipRange and MatGetSize cannot be called before MatXXXSetPreallocation
_localSize = nbRows;
#ifdef HAVE_MPI
- PetscMPIInt commSize;
- MPI_Comm_size(_comm,&commSize);
if (commSize>1){
_localRowStart = 0;
if (Msg::GetCommRank() != 0) {
@@ -159,7 +171,7 @@ void linearSystemPETSc<scalar>::allocate(int nbRows)
#endif
// preallocation option must be set after other options
_try(VecCreate(_comm, &_x));
- _try(VecSetSizes(_x, nbRows, PETSC_DETERMINE));
+ _try(VecSetSizes(_x, blockSize * nbRows, PETSC_DETERMINE));
// override the default options with the ones from the option
// database (if any)
if (this->_parameters.count("petscPrefix"))
@@ -345,6 +357,22 @@ int linearSystemPETSc<scalar>::systemSolve()
return 1;
}
+template <class scalar>
+void linearSystemPETSc<scalar>::printMatlab(const char *filename) const
+{
+ _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
+ _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+ _try(VecAssemblyBegin(_b));
+ _try(VecAssemblyEnd(_b));
+
+ PetscViewer viewer;
+ PetscViewerASCIIOpen(PETSC_COMM_WORLD, filename, &viewer);
+ PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
+ MatView(_a, viewer);
+ PetscViewerDestroy(&viewer);
+ return;
+}
+
/*template <class scalar>
std::vector<scalar> linearSystemPETSc<scalar>::getData()
{
diff --git a/wrappers/gmshpy/gmshSolver.i b/wrappers/gmshpy/gmshSolver.i
index cff31d8144..20b8e64462 100644
--- a/wrappers/gmshpy/gmshSolver.i
+++ b/wrappers/gmshpy/gmshSolver.i
@@ -35,6 +35,7 @@
#if defined(HAVE_PETSC)
%include "linearSystemPETSc.h"
%template(linearSystemPETScDouble) linearSystemPETSc<double>;
+%template(linearSystemPETScBlockDouble) linearSystemPETSc<fullMatrix<double> >;
#endif
%include "eigenSolver.h"
#endif
--
GitLab