diff --git a/MathLib/LinAlg/BLAS.cpp b/MathLib/LinAlg/BLAS.cpp index 51f2add2370cb4fd1169c593476234958b4013f8..c72cb1684840ee6444aec83bb7dc260883ab583a 100644 --- a/MathLib/LinAlg/BLAS.cpp +++ b/MathLib/LinAlg/BLAS.cpp @@ -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 diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.cpp b/MathLib/LinAlg/PETSc/PETScMatrix.cpp index a0c56f6b21bc14049aaabd7fcb22847106bb6286..2695a1b755a75aca736e6fae3449d82e5deaa351 100644 --- a/MathLib/LinAlg/PETSc/PETScMatrix.cpp +++ b/MathLib/LinAlg/PETSc/PETScMatrix.cpp @@ -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) diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.h b/MathLib/LinAlg/PETSc/PETScMatrix.h index d063d58503758c19231e84da34e6b91c2ea358ee..1b696bf9e0088cc14a2cbece1385c2745950239d 100644 --- a/MathLib/LinAlg/PETSc/PETScMatrix.h +++ b/MathLib/LinAlg/PETSc/PETScMatrix.h @@ -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); }; /*! diff --git a/MathLib/LinAlg/PETSc/PETScVector.cpp b/MathLib/LinAlg/PETSc/PETScVector.cpp index 0ee62a9c8b9a9c81eceb35558a238645256e0078..b56c4031a0a2d6283683d42e849e0244b9baa6bc 100644 --- a/MathLib/LinAlg/PETSc/PETScVector.cpp +++ b/MathLib/LinAlg/PETSc/PETScVector.cpp @@ -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) diff --git a/MathLib/LinAlg/PETSc/PETScVector.h b/MathLib/LinAlg/PETSc/PETScVector.h index ce15cbe8f4614dd990a96d24e8191f1644bfe091..f2abd9aa187fb70898c55a410c8f5ca1da7e52e5 100644 --- a/MathLib/LinAlg/PETSc/PETScVector.h +++ b/MathLib/LinAlg/PETSc/PETScVector.h @@ -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;