Skip to content
Snippets Groups Projects
Commit 2feceb74 authored by Norihiro Watanabe's avatar Norihiro Watanabe
Browse files

remove std::unique_ptr from PETSc vec mat

parent 44ec38a4
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MathLib/LinAlg/BLAS.cpp
+ 10
− 10
View file @ 2feceb74
......@@ -75,28 +75,28 @@ void copy(PETScVector const& x, PETScVector& y)
void scale(PETScVector& x, double const a)
{
VecScale(*x.getRawVector(), a);
VecScale(x.getRawVector(), a);
}
// y = a*y + X
void aypx(PETScVector& y, double const a, PETScVector const& x)
{
// TODO check sizes
VecAYPX(*y.getRawVector(), a, *x.getRawVector());
VecAYPX(y.getRawVector(), a, x.getRawVector());
}
// y = a*x + y
void axpy(PETScVector& y, double const a, PETScVector const& x)
{
// TODO check sizes
VecAXPY(*y.getRawVector(), a, *x.getRawVector());
VecAXPY(y.getRawVector(), a, x.getRawVector());
}
// y = a*x + b*y
void axpby(PETScVector& y, double const a, double const b, PETScVector const& x)
{
// TODO check sizes
VecAXPBY(*y.getRawVector(), a, b, *x.getRawVector());
VecAXPBY(y.getRawVector(), a, b, x.getRawVector());
}
// Explicit specialization
......@@ -105,7 +105,7 @@ template<>
void componentwiseDivide(PETScVector& w,
PETScVector const& x, PETScVector const& y)
{
VecPointwiseDivide(*w.getRawVector(), *x.getRawVector(), *y.getRawVector());
VecPointwiseDivide(w.getRawVector(), x.getRawVector(), y.getRawVector());
}
// Explicit specialization
......@@ -173,7 +173,7 @@ void matMult(PETScMatrix const& A, PETScVector const& x, PETScVector& y)
// TODO check sizes
assert(&x != &y);
if (!y.getRawVector()) y.shallowCopy(x);
MatMult(A.getRawMatrix(), *x.getRawVector(), *y.getRawVector());
MatMult(A.getRawMatrix(), x.getRawVector(), y.getRawVector());
}
// v3 = A*v1 + v2
......@@ -183,7 +183,7 @@ void matMultAdd(PETScMatrix const& A, PETScVector const& v1,
// TODO check sizes
assert(&v1 != &v3);
if (!v3.getRawVector()) v3.shallowCopy(v1);
MatMultAdd(A.getRawMatrix(), *v1.getRawVector(), *v2.getRawVector(), *v3.getRawVector());
MatMultAdd(A.getRawMatrix(), v1.getRawVector(), v2.getRawVector(), v3.getRawVector());
}
void finalizeAssembly(PETScMatrix& A)
......@@ -225,21 +225,21 @@ void scale(EigenVector& x, double const a)
void aypx(EigenVector& y, double const a, EigenVector const& x)
{
// TODO: does that break anything?
y.getRawVector() = a*y.getRawVector() + x.getRawVector();
y.getRawVector() = a * y.getRawVector() + x.getRawVector();
}
// y = a*x + y
void axpy(EigenVector& y, double const a, EigenVector const& x)
{
// TODO: does that break anything?
y.getRawVector() += a*x.getRawVector();
y.getRawVector() += a * x.getRawVector();
}
// y = a*x + y
void axpby(EigenVector& y, double const a, double const b, EigenVector const& x)
{
// TODO: does that break anything?
y.getRawVector() = a*x.getRawVector() + b*y.getRawVector();
y.getRawVector() = a * x.getRawVector() + b * y.getRawVector();
}
// Explicit specialization
......
MathLib/LinAlg/PETSc/PETScMatrix.cpp
+ 19
− 24
View file @ 2feceb74
......@@ -49,16 +49,14 @@ PETScMatrix::PETScMatrix (const PetscInt nrows, const PetscInt ncols, const PETS
}
PETScMatrix::PETScMatrix(const PETScMatrix &A)
: _A(new PETSc_Mat)
, _nrows(A._nrows)
: _nrows(A._nrows)
, _ncols(A._ncols)
, _n_loc_rows(A._n_loc_rows)
, _n_loc_cols(A._n_loc_cols)
, _start_rank(A._start_rank)
, _end_rank(A._end_rank)
{
_A.reset(new PETSc_Mat);
MatConvert(*A._A, MATSAME, MAT_INITIAL_MATRIX, _A.get());
MatConvert(A._A, MATSAME, MAT_INITIAL_MATRIX, &_A);
}
PETScMatrix&
......@@ -73,11 +71,10 @@ PETScMatrix::operator=(PETScMatrix const& A)
if (_A) {
// TODO this is the slowest option for copying
MatCopy(*A._A, *_A, DIFFERENT_NONZERO_PATTERN);
MatCopy(A._A, _A, DIFFERENT_NONZERO_PATTERN);
} else {
destroy();
_A.reset(new PETSc_Mat);
MatConvert(*A._A, MATSAME, MAT_INITIAL_MATRIX, _A.get());
MatConvert(A._A, MATSAME, MAT_INITIAL_MATRIX, &_A);
}
return *this;
......@@ -93,11 +90,11 @@ void PETScMatrix::setRowsColumnsZero(std::vector<PetscInt> const& row_pos)
// This avoids all reductions in the zero row routines
// and thus improves performance for very large process counts.
// See PETSc doc about MAT_NO_OFF_PROC_ZERO_ROWS.
MatSetOption(*_A, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
MatSetOption(_A, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
if(nrows>0)
MatZeroRows(*_A, nrows, &row_pos[0], one, PETSC_NULL, PETSC_NULL);
MatZeroRows(_A, nrows, &row_pos[0], one, PETSC_NULL, PETSC_NULL);
else
MatZeroRows(*_A, 0, PETSC_NULL, one, PETSC_NULL, PETSC_NULL);
MatZeroRows(_A, 0, PETSC_NULL, one, PETSC_NULL, PETSC_NULL);
}
void PETScMatrix::viewer(const std::string &file_name, const PetscViewerFormat vw_format)
......@@ -108,8 +105,8 @@ void PETScMatrix::viewer(const std::string &file_name, const PetscViewerFormat v
finalizeAssembly();
PetscObjectSetName((PetscObject)*_A,"Stiffness_matrix");
MatView(*_A,viewer);
PetscObjectSetName((PetscObject)_A,"Stiffness_matrix");
MatView(_A,viewer);
// This preprocessor is only for debugging, e.g. dump the matrix and exit the program.
//#define EXIT_TEST
......@@ -123,22 +120,20 @@ void PETScMatrix::viewer(const std::string &file_name, const PetscViewerFormat v
void PETScMatrix::create(const PetscInt d_nz, const PetscInt o_nz)
{
_A.reset(new PETSc_Mat);
MatCreate(PETSC_COMM_WORLD, &_A);
MatSetSizes(_A, _n_loc_rows, _n_loc_cols, _nrows, _ncols);
MatCreate(PETSC_COMM_WORLD, _A.get());
MatSetSizes(*_A, _n_loc_rows, _n_loc_cols, _nrows, _ncols);
MatSetFromOptions(_A);
MatSetFromOptions(*_A);
MatSetType(*_A, MATMPIAIJ);
MatSeqAIJSetPreallocation(*_A, d_nz, PETSC_NULL);
MatMPIAIJSetPreallocation(*_A, d_nz, PETSC_NULL, o_nz, PETSC_NULL);
// If pre-allocation does not work one can use MatSetUp(*_A), which is much
MatSetType(_A, MATMPIAIJ);
MatSeqAIJSetPreallocation(_A, d_nz, PETSC_NULL);
MatMPIAIJSetPreallocation(_A, d_nz, PETSC_NULL, o_nz, PETSC_NULL);
// If pre-allocation does not work one can use MatSetUp(_A), which is much
// slower.
MatGetOwnershipRange(*_A, &_start_rank, &_end_rank);
MatGetSize(*_A, &_nrows, &_ncols);
MatGetLocalSize(*_A, &_n_loc_rows, &_n_loc_cols);
MatGetOwnershipRange(_A, &_start_rank, &_end_rank);
MatGetSize(_A, &_nrows, &_ncols);
MatGetLocalSize(_A, &_n_loc_rows, &_n_loc_cols);
}
bool finalizeMatrixAssembly(PETScMatrix &mat, const MatAssemblyType asm_type)
......
MathLib/LinAlg/PETSc/PETScMatrix.h
+ 11
− 12
View file @ 2feceb74
......@@ -15,7 +15,6 @@
#ifndef PETSCMATRIX_H_
#define PETSCMATRIX_H_
#include <memory>
#include <string>
#include <vector>
......@@ -72,8 +71,8 @@ class PETScMatrix
*/
void finalizeAssembly(const MatAssemblyType asm_type = MAT_FINAL_ASSEMBLY)
{
MatAssemblyBegin(*_A, asm_type);
MatAssemblyEnd(*_A, asm_type);
MatAssemblyBegin(_A, asm_type);
MatAssemblyEnd(_A, asm_type);
}
/// Get the number of rows.
......@@ -116,7 +115,7 @@ class PETScMatrix
/// Get matrix reference.
PETSc_Mat &getRawMatrix()
{
return *_A;
return _A;
}
/*! Get a matrix reference.
......@@ -127,13 +126,13 @@ class PETScMatrix
*/
PETSc_Mat const& getRawMatrix() const
{
return *_A;
return _A;
}
/// Set all entries to zero.
void setZero()
{
MatZeroEntries(*_A);
MatZeroEntries(_A);
}
/*
......@@ -155,7 +154,7 @@ class PETScMatrix
*/
void multiply(const PETScVector &vec, PETScVector &vec_r)
{
MatMult(*_A, vec.getData(), vec_r.getData() );
MatMult(_A, vec.getData(), vec_r.getData() );
}
/*!
......@@ -166,7 +165,7 @@ class PETScMatrix
*/
void set(const PetscInt i, const PetscInt j, const PetscScalar value)
{
MatSetValue(*_A, i, j, value, INSERT_VALUES);
MatSetValue(_A, i, j, value, INSERT_VALUES);
}
/*!
......@@ -177,7 +176,7 @@ class PETScMatrix
*/
void add(const PetscInt i, const PetscInt j, const PetscScalar value)
{
MatSetValue(*_A, i, j, value, ADD_VALUES);
MatSetValue(_A, i, j, value, ADD_VALUES);
}
/// Add sub-matrix at positions given by \c indices.
......@@ -237,10 +236,10 @@ class PETScMatrix
const PetscViewerFormat vw_format = PETSC_VIEWER_ASCII_MATLAB );
private:
void destroy() { if (_A) MatDestroy(_A.get()); _A.reset(nullptr); }
void destroy() { if (_A) MatDestroy(&_A); _A = nullptr; }
/// PETSc matrix
std::unique_ptr<PETSc_Mat> _A;
PETSc_Mat _A = nullptr;
/// Number of the global rows
PetscInt _nrows;
......@@ -286,7 +285,7 @@ void PETScMatrix::add(std::vector<PetscInt> const& row_pos,
const PetscInt nrows = static_cast<PetscInt> (row_pos.size());
const PetscInt ncols = static_cast<PetscInt> (col_pos.size());
MatSetValues(*_A, nrows, &row_pos[0], ncols, &col_pos[0], &sub_mat(0,0), ADD_VALUES);
MatSetValues(_A, nrows, &row_pos[0], ncols, &col_pos[0], &sub_mat(0,0), ADD_VALUES);
};
/*!
......
MathLib/LinAlg/PETSc/PETScVector.cpp
+ 30
− 36
View file @ 2feceb74
......@@ -28,15 +28,13 @@ namespace MathLib
{
PETScVector::PETScVector(const PetscInt vec_size, const bool is_global_size)
{
_v.reset(new PETSc_Vec);
if( is_global_size ) {
VecCreate(PETSC_COMM_WORLD, _v.get());
VecSetSizes(*_v, PETSC_DECIDE, vec_size);
VecCreate(PETSC_COMM_WORLD, &_v);
VecSetSizes(_v, PETSC_DECIDE, vec_size);
} else {
// Fix size partitioning
// the size can be associated to specific memory allocation of a matrix
VecCreateMPI(PETSC_COMM_WORLD, vec_size, PETSC_DECIDE, _v.get());
VecCreateMPI(PETSC_COMM_WORLD, vec_size, PETSC_DECIDE, &_v);
}
config();
......@@ -48,20 +46,18 @@ PETScVector::PETScVector(const PetscInt vec_size,
: _size_ghosts{static_cast<PetscInt>(ghost_ids.size())}
, _has_ghost_id{true}
{
_v.reset(new PETSc_Vec);
PetscInt nghosts = static_cast<PetscInt>( ghost_ids.size() );
if ( is_global_size )
{
VecCreateGhost(PETSC_COMM_WORLD, PETSC_DECIDE, vec_size, nghosts,
ghost_ids.data(), _v.get());
ghost_ids.data(), &_v);
}
else
{
VecCreate(PETSC_COMM_WORLD, _v.get());
VecSetType(*_v, VECMPI);
VecSetSizes(*_v, vec_size, PETSC_DECIDE);
VecMPISetGhost(*_v, nghosts, ghost_ids.data());
VecCreate(PETSC_COMM_WORLD, &_v);
VecSetType(_v, VECMPI);
VecSetSizes(_v, vec_size, PETSC_DECIDE);
VecMPISetGhost(_v, nghosts, ghost_ids.data());
}
config();
......@@ -74,7 +70,7 @@ PETScVector::PETScVector(const PETScVector &existing_vec, const bool deep_copy)
// Copy values
if(deep_copy)
{
VecCopy(*existing_vec._v, *_v);
VecCopy(existing_vec._v, _v);
}
}
......@@ -91,20 +87,20 @@ PETScVector::PETScVector(PETScVector &&other)
void PETScVector::config()
{
VecSetFromOptions(*_v);
// VecSetUp(*_v); // for petsc ver.>3.3
VecGetOwnershipRange(*_v, &_start_rank, &_end_rank);
VecSetFromOptions(_v);
// VecSetUp(_v); // for petsc ver.>3.3
VecGetOwnershipRange(_v, &_start_rank, &_end_rank);
VecGetLocalSize(*_v, &_size_loc);
VecGetSize(*_v, &_size);
VecGetLocalSize(_v, &_size_loc);
VecGetSize(_v, &_size);
VecSetOption(*_v, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
VecSetOption(_v, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
}
void PETScVector::finalizeAssembly()
{
VecAssemblyBegin(*_v);
VecAssemblyEnd(*_v);
VecAssemblyBegin(_v);
VecAssemblyEnd(_v);
}
void PETScVector::gatherLocalVectors( PetscScalar local_array[],
......@@ -143,7 +139,7 @@ void PETScVector::getGlobalVector(PetscScalar u[])
#endif
PetscScalar *xp = nullptr;
VecGetArray(*_v, &xp);
VecGetArray(_v, &xp);
gatherLocalVectors(xp, u);
......@@ -151,7 +147,7 @@ void PETScVector::getGlobalVector(PetscScalar u[])
// for a communication load balance:
//MPI_Barrier(PETSC_COMM_WORLD);
VecRestoreArray(*_v, &xp);
VecRestoreArray(_v, &xp);
//TEST
#ifdef TEST_MEM_PETSC
......@@ -174,13 +170,13 @@ PetscScalar* PETScVector::getLocalVector() const
PetscScalar *loc_array;
if (_has_ghost_id)
{
VecGhostUpdateBegin(*_v, INSERT_VALUES, SCATTER_FORWARD);
VecGhostUpdateEnd(*_v, INSERT_VALUES, SCATTER_FORWARD);
VecGhostGetLocalForm(*_v, &_v_loc);
VecGhostUpdateBegin(_v, INSERT_VALUES, SCATTER_FORWARD);
VecGhostUpdateEnd(_v, INSERT_VALUES, SCATTER_FORWARD);
VecGhostGetLocalForm(_v, &_v_loc);
VecGetArray(_v_loc, &loc_array);
}
else
VecGetArray(*_v, &loc_array);
VecGetArray(_v, &loc_array);
return loc_array;
}
......@@ -189,10 +185,10 @@ void PETScVector::restoreArray(PetscScalar* array) const
if (_has_ghost_id)
{
VecRestoreArray(_v_loc, &array);
VecGhostRestoreLocalForm(*_v, &_v_loc);
VecGhostRestoreLocalForm(_v, &_v_loc);
}
else
VecRestoreArray(*_v, &array);
VecRestoreArray(_v, &array);
}
PetscScalar PETScVector::getNorm(MathLib::VecNormType nmtype) const
......@@ -214,7 +210,7 @@ PetscScalar PETScVector::getNorm(MathLib::VecNormType nmtype) const
}
PetscScalar norm = 0.;
VecNorm(*_v, petsc_norm, &norm);
VecNorm(_v, petsc_norm, &norm);
return norm;
}
......@@ -224,8 +220,8 @@ void PETScVector::viewer(const std::string &file_name, const PetscViewerFormat v
PetscViewerASCIIOpen(PETSC_COMM_WORLD, file_name.c_str(), &viewer);
PetscViewerPushFormat(viewer, vw_format);
PetscObjectSetName((PetscObject)*_v, file_name.c_str());
VecView(*_v, viewer);
PetscObjectSetName((PetscObject)_v, file_name.c_str());
VecView(_v, viewer);
#define nEXIT_TEST
#ifdef EXIT_TEST
......@@ -240,9 +236,7 @@ void PETScVector::shallowCopy(const PETScVector &v)
{
destroy();
_v.reset(new PETSc_Vec);
VecDuplicate(*v._v, _v.get());
VecDuplicate(v.getRawVector(), &_v);
_start_rank = v._start_rank;
_end_rank = v._end_rank;
......@@ -251,7 +245,7 @@ void PETScVector::shallowCopy(const PETScVector &v)
_size_ghosts = v._size_ghosts;
_has_ghost_id = v._has_ghost_id;
VecSetOption(*_v, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
VecSetOption(_v, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
}
void finalizeVectorAssembly(PETScVector &vec)
......
MathLib/LinAlg/PETSc/PETScVector.h
+ 17
− 18
View file @ 2feceb74
......@@ -17,7 +17,6 @@
#ifndef PETSCVECTOR_H_
#define PETSCVECTOR_H_
#include <memory>
#include <string>
#include <vector>
......@@ -128,7 +127,7 @@ class PETScVector
*/
void set(const PetscInt i, const PetscScalar value)
{
VecSetValue(*_v, i, value, INSERT_VALUES);
VecSetValue(_v, i, value, INSERT_VALUES);
}
/*!
......@@ -138,7 +137,7 @@ class PETScVector
*/
void add(const PetscInt i, const PetscScalar value)
{
VecSetValue(*_v, i, value, ADD_VALUES);
VecSetValue(_v, i, value, ADD_VALUES);
}
/*!
......@@ -150,7 +149,7 @@ class PETScVector
template<class T_SUBVEC> void add(const std::vector<PetscInt> &e_idxs,
const T_SUBVEC &sub_vec)
{
VecSetValues(*_v, e_idxs.size(), &e_idxs[0], &sub_vec[0], ADD_VALUES);
VecSetValues(_v, e_idxs.size(), &e_idxs[0], &sub_vec[0], ADD_VALUES);
}
/*!
......@@ -162,7 +161,7 @@ class PETScVector
template<class T_SUBVEC> void set(const std::vector<PetscInt> &e_idxs,
const T_SUBVEC &sub_vec)
{
VecSetValues(*_v, e_idxs.size(), &e_idxs[0], &sub_vec[0], INSERT_VALUES);
VecSetValues(_v, e_idxs.size(), &e_idxs[0], &sub_vec[0], INSERT_VALUES);
}
/*!
......@@ -174,14 +173,14 @@ class PETScVector
template<class T_SUBVEC> void get(const std::vector<PetscInt> &e_idxs,
T_SUBVEC &sub_vec)
{
VecGetValues(*_v, e_idxs.size(), &e_idxs[0], &sub_vec[0]);
VecGetValues(_v, e_idxs.size(), &e_idxs[0], &sub_vec[0]);
}
// TODO preliminary
double operator[] (PetscInt idx) const
{
double value;
VecGetValues(*_v, 1, &idx, &value);
VecGetValues(_v, 1, &idx, &value);
return value;
}
......@@ -202,7 +201,7 @@ class PETScVector
PetscScalar get(const PetscInt idx) const
{
PetscScalar x;
VecGetValues(*_v, 1, &idx, &x);
VecGetValues(_v, 1, &idx, &x);
return x;
}
......@@ -210,13 +209,13 @@ class PETScVector
/// Get PETsc vector. Use it only for test purpose
const PETSc_Vec &getData() const
{
return *_v;
return _v;
}
/// Initialize the vector with a constant value
void operator = (const PetscScalar val)
{
VecSet(*_v, val);
VecSet(_v, val);
}
// TODO preliminary
......@@ -226,7 +225,7 @@ class PETScVector
PETScVector& operator = (const PETScVector &v_in)
{
if (!_v) shallowCopy(v_in);
VecCopy(*v_in._v, *_v);
VecCopy(v_in._v, _v);
return *this;
}
......@@ -240,18 +239,18 @@ class PETScVector
void operator += (const PETScVector& v_in)
{
if (!_v) shallowCopy(v_in);
VecAXPY(*_v, 1.0, *v_in._v);
VecAXPY(_v, 1.0, v_in._v);
}
/// Overloaded operator: subtract
void operator -= (const PETScVector& v_in)
{
if (!_v) shallowCopy(v_in);
VecAXPY(*_v, -1.0, *v_in._v);
VecAXPY(_v, -1.0, v_in._v);
}
//! Exposes the underlying PETSc vector.
PETSc_Vec* getRawVector() { return _v.get(); }
PETSc_Vec getRawVector() { return _v; }
/*! Exposes the underlying PETSc vector.
*
......@@ -259,7 +258,7 @@ class PETScVector
* This method is dangerous insofar as you can do arbitrary things also
* with a const PETSc vector!
*/
const PETSc_Vec* getRawVector() const {return _v.get(); }
PETSc_Vec getRawVector() const {return _v; }
/*! View the global vector for test purpose. Do not use it for output a big vector.
......@@ -291,12 +290,12 @@ class PETScVector
void shallowCopy(const PETScVector &v);
private:
void destroy() { if (_v) VecDestroy(_v.get()); _v.reset(nullptr); }
void destroy() { if (_v) VecDestroy(&_v); _v = nullptr; }
std::unique_ptr<PETSc_Vec> _v;
PETSc_Vec _v = nullptr;
/// Local vector, which is only for the case that _v is created
/// with ghost entries.
mutable PETSc_Vec _v_loc;
mutable PETSc_Vec _v_loc = nullptr;
/// Starting index in a rank
PetscInt _start_rank;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment