From ed895f48b73fcf437fa972d59222984f348643a8 Mon Sep 17 00:00:00 2001 From: Dylan Copeland Date: Wed, 10 Jul 2024 21:47:08 -0700 Subject: [PATCH] Introduced smart pointers in some functions, removed some unnecessary functions with pointer arguments, and fixed a memory leak. --- .../dg_advection_local_rom_matrix_interp.cpp | 8 +- lib/algo/DMD.cpp | 18 ++-- lib/algo/DMD.h | 2 +- lib/algo/DMDc.cpp | 22 ++--- lib/algo/DMDc.h | 2 +- lib/algo/greedy/GreedySampler.cpp | 1 - .../manifold_interp/VectorInterpolator.cpp | 7 +- lib/algo/manifold_interp/VectorInterpolator.h | 3 +- lib/linalg/BasisGenerator.cpp | 4 +- lib/linalg/Matrix.h | 24 ----- lib/linalg/Vector.cpp | 59 +------------ lib/linalg/Vector.h | 88 ++----------------- lib/linalg/svd/IncrementalSVD.cpp | 4 +- lib/linalg/svd/IncrementalSVDBrand.cpp | 9 +- 14 files changed, 50 insertions(+), 201 deletions(-) mode change 100755 => 100644 lib/algo/DMD.cpp diff --git a/examples/prom/dg_advection_local_rom_matrix_interp.cpp b/examples/prom/dg_advection_local_rom_matrix_interp.cpp index abbe6621c..4958cb49d 100644 --- a/examples/prom/dg_advection_local_rom_matrix_interp.cpp +++ b/examples/prom/dg_advection_local_rom_matrix_interp.cpp @@ -697,7 +697,7 @@ int main(int argc, char *argv[]) CAROM::Matrix *M_hat_carom, *K_hat_carom; DenseMatrix *M_hat, *K_hat; - CAROM::Vector *b_hat_carom, *u_init_hat_carom; + std::shared_ptr b_hat_carom, u_init_hat_carom; Vector *b_hat, *u_init_hat; Vector u_init(*U); @@ -770,11 +770,11 @@ int main(int argc, char *argv[]) Vector b_vec = *B; CAROM::Vector b_carom(b_vec.GetData(), b_vec.Size(), true); - b_hat_carom = spatialbasis->transposeMult(&b_carom); + b_hat_carom.reset(spatialbasis->transposeMult(&b_carom)); if (interp_prep) b_hat_carom->write("b_hat_" + std::to_string(f_factor)); b_hat = new Vector(b_hat_carom->getData(), b_hat_carom->dim()); - u_init_hat_carom = new CAROM::Vector(numColumnRB, false); + u_init_hat_carom.reset(new CAROM::Vector(numColumnRB, false)); ComputeCtAB_vec(K_mat, *U, *spatialbasis, *u_init_hat_carom); if (interp_prep) u_init_hat_carom->write("u_init_hat_" + std::to_string( f_factor)); @@ -1012,8 +1012,6 @@ int main(int argc, char *argv[]) delete K_hat_carom; delete M_hat; delete K_hat; - delete b_hat_carom; - delete u_init_hat_carom; delete b_hat; delete u_init_hat; delete u_in; diff --git a/lib/algo/DMD.cpp b/lib/algo/DMD.cpp old mode 100755 new mode 100644 index a66f3b65b..2ea11e77e --- a/lib/algo/DMD.cpp +++ b/lib/algo/DMD.cpp @@ -144,7 +144,6 @@ DMD::~DMD() delete d_phi_imaginary; delete d_phi_real_squared_inverse; delete d_phi_imaginary_squared_inverse; - delete d_projected_init_real; delete d_projected_init_imaginary; } @@ -607,15 +606,15 @@ DMD::projectInitialCondition(const Vector* init, double t_offset) Vector* rhs_real = d_phi_real->transposeMult(init); Vector* rhs_imaginary = d_phi_imaginary->transposeMult(init); - Vector* d_projected_init_real_1 = d_phi_real_squared_inverse->mult(rhs_real); + Vector* d_projected_init_real_1 = d_phi_real_squared_inverse->mult(*rhs_real); Vector* d_projected_init_real_2 = d_phi_imaginary_squared_inverse->mult( - rhs_imaginary); - d_projected_init_real = d_projected_init_real_1->plus(d_projected_init_real_2); + *rhs_imaginary); + d_projected_init_real = d_projected_init_real_1->plus(*d_projected_init_real_2); Vector* d_projected_init_imaginary_1 = d_phi_real_squared_inverse->mult( - rhs_imaginary); + *rhs_imaginary); Vector* d_projected_init_imaginary_2 = d_phi_imaginary_squared_inverse->mult( - rhs_real); + *rhs_real); d_projected_init_imaginary = d_projected_init_imaginary_2->minus( d_projected_init_imaginary_1); @@ -655,9 +654,10 @@ DMD::predict(double t, int deg) Matrix* d_phi_mult_eigs_imaginary = d_phi_pair.second; Vector* d_predicted_state_real_1 = d_phi_mult_eigs_real->mult( - d_projected_init_real); + *d_projected_init_real); + Vector* d_predicted_state_real_2 = d_phi_mult_eigs_imaginary->mult( - d_projected_init_imaginary); + *d_projected_init_imaginary); Vector* d_predicted_state_real = d_predicted_state_real_1->minus( d_predicted_state_real_2); addOffset(d_predicted_state_real, t, deg); @@ -820,7 +820,7 @@ DMD::load(std::string base_file_name) d_phi_imaginary_squared_inverse->read(full_file_name); full_file_name = base_file_name + "_projected_init_real"; - d_projected_init_real = new Vector(); + d_projected_init_real.reset(new Vector()); d_projected_init_real->read(full_file_name); full_file_name = base_file_name + "_projected_init_imaginary"; diff --git a/lib/algo/DMD.h b/lib/algo/DMD.h index ba3750251..d5cd7b5f9 100644 --- a/lib/algo/DMD.h +++ b/lib/algo/DMD.h @@ -435,7 +435,7 @@ class DMD /** * @brief The real part of the projected initial condition. */ - Vector* d_projected_init_real = NULL; + std::shared_ptr d_projected_init_real; /** * @brief The imaginary part of the projected initial condition. diff --git a/lib/algo/DMDc.cpp b/lib/algo/DMDc.cpp index 7e4591767..d92da293d 100644 --- a/lib/algo/DMDc.cpp +++ b/lib/algo/DMDc.cpp @@ -150,7 +150,6 @@ DMDc::~DMDc() delete d_phi_imaginary; delete d_phi_real_squared_inverse; delete d_phi_imaginary_squared_inverse; - delete d_projected_init_real; delete d_projected_init_imaginary; } @@ -659,15 +658,15 @@ DMDc::project(const Vector* init, const Matrix* controls, double t_offset) Vector* init_real = d_phi_real->transposeMult(init); Vector* init_imaginary = d_phi_imaginary->transposeMult(init); - Vector* d_projected_init_real_1 = d_phi_real_squared_inverse->mult(init_real); + Vector* d_projected_init_real_1 = d_phi_real_squared_inverse->mult(*init_real); Vector* d_projected_init_real_2 = d_phi_imaginary_squared_inverse->mult( - init_imaginary); - d_projected_init_real = d_projected_init_real_1->plus(d_projected_init_real_2); + *init_imaginary); + d_projected_init_real = d_projected_init_real_1->plus(*d_projected_init_real_2); Vector* d_projected_init_imaginary_1 = d_phi_real_squared_inverse->mult( - init_imaginary); + *init_imaginary); Vector* d_projected_init_imaginary_2 = d_phi_imaginary_squared_inverse->mult( - init_real); + *init_real); d_projected_init_imaginary = d_projected_init_imaginary_2->minus( d_projected_init_imaginary_1); @@ -730,9 +729,10 @@ DMDc::predict(double t) Matrix* d_phi_mult_eigs_imaginary = d_phi_pair.second; Vector* d_predicted_state_real_1 = d_phi_mult_eigs_real->mult( - d_projected_init_real); + *d_projected_init_real); + Vector* d_predicted_state_real_2 = d_phi_mult_eigs_imaginary->mult( - d_projected_init_imaginary); + *d_projected_init_imaginary); Vector* d_predicted_state_real = d_predicted_state_real_1->minus( d_predicted_state_real_2); addOffset(d_predicted_state_real); @@ -754,9 +754,9 @@ DMDc::predict(double t) d_projected_controls_real->getColumn(k, *f_control_real); d_projected_controls_imaginary->getColumn(k, *f_control_imaginary); d_predicted_state_real_1 = d_phi_mult_eigs_real->mult( - f_control_real); + *f_control_real); d_predicted_state_real_2 = d_phi_mult_eigs_imaginary->mult( - f_control_imaginary); + *f_control_imaginary); *d_predicted_state_real += *d_predicted_state_real_1; *d_predicted_state_real -= *d_predicted_state_real_2; @@ -921,7 +921,7 @@ DMDc::load(std::string base_file_name) d_phi_imaginary_squared_inverse->read(full_file_name); full_file_name = base_file_name + "_projected_init_real"; - d_projected_init_real = new Vector(); + d_projected_init_real.reset(new Vector()); d_projected_init_real->read(full_file_name); full_file_name = base_file_name + "_projected_init_imaginary"; diff --git a/lib/algo/DMDc.h b/lib/algo/DMDc.h index f298c4cae..7a8c4e642 100644 --- a/lib/algo/DMDc.h +++ b/lib/algo/DMDc.h @@ -416,7 +416,7 @@ class DMDc /** * @brief The real part of the projected initial condition. */ - Vector* d_projected_init_real = NULL; + std::shared_ptr d_projected_init_real; /** * @brief The imaginary part of the projected initial condition. diff --git a/lib/algo/greedy/GreedySampler.cpp b/lib/algo/greedy/GreedySampler.cpp index ec41da39e..d3b8a3e41 100644 --- a/lib/algo/greedy/GreedySampler.cpp +++ b/lib/algo/greedy/GreedySampler.cpp @@ -1239,7 +1239,6 @@ GreedySampler::generateRandomPoints(int num_points) std::shared_ptr GreedySampler::getNearestROM(Vector point) { - CAROM_VERIFY(point.dim() == d_parameter_points[0].dim()); double closest_dist_to_points = INT_MAX; diff --git a/lib/algo/manifold_interp/VectorInterpolator.cpp b/lib/algo/manifold_interp/VectorInterpolator.cpp index c11f4d8dd..7b09dcd62 100644 --- a/lib/algo/manifold_interp/VectorInterpolator.cpp +++ b/lib/algo/manifold_interp/VectorInterpolator.cpp @@ -104,7 +104,7 @@ Vector* VectorInterpolator::obtainLogInterpolatedVector( return obtainInterpolatedVector(d_gammas, d_lambda_T, d_interp_method, rbf); } -Vector* VectorInterpolator::interpolate(Vector* point) +std::shared_ptr VectorInterpolator::interpolate(Vector* point) { if (d_gammas.size() == 0) { @@ -139,8 +139,9 @@ Vector* VectorInterpolator::interpolate(Vector* point) Vector* log_interpolated_vector = obtainLogInterpolatedVector(rbf); // The exp mapping is X + the interpolated gamma - Vector* interpolated_vector = d_rotated_reduced_vectors[d_ref_point]->plus( - log_interpolated_vector); + std::shared_ptr interpolated_vector = + d_rotated_reduced_vectors[d_ref_point]->plus( + *log_interpolated_vector); delete log_interpolated_vector; return interpolated_vector; } diff --git a/lib/algo/manifold_interp/VectorInterpolator.h b/lib/algo/manifold_interp/VectorInterpolator.h index a1fdf5fe3..fb213ff5c 100644 --- a/lib/algo/manifold_interp/VectorInterpolator.h +++ b/lib/algo/manifold_interp/VectorInterpolator.h @@ -14,6 +14,7 @@ #define included_VectorInterpolator_h #include "Interpolator.h" +#include #include #include @@ -68,7 +69,7 @@ class VectorInterpolator : public Interpolator * * @param[in] point The unsampled parameter point. */ - Vector* interpolate(Vector* point); + std::shared_ptr interpolate(Vector* point); private: diff --git a/lib/linalg/BasisGenerator.cpp b/lib/linalg/BasisGenerator.cpp index 9e2cec4ae..fa5cda5f6 100644 --- a/lib/linalg/BasisGenerator.cpp +++ b/lib/linalg/BasisGenerator.cpp @@ -227,7 +227,7 @@ BasisGenerator::computeNextSampleTime( Vector* l = basis->transposeMult(u_vec); // basisl = basis * l - Vector* basisl = basis->mult(l); + Vector* basisl = basis->mult(*l); // Compute u - basisl. Vector* eta = u_vec.minus(basisl); @@ -240,7 +240,7 @@ BasisGenerator::computeNextSampleTime( l = basis->transposeMult(rhs_vec); // basisl = basis * l - basisl = basis->mult(l); + basisl = basis->mult(*l); // Compute rhs - basisl. Vector* eta_dot = rhs_vec.minus(basisl); diff --git a/lib/linalg/Matrix.h b/lib/linalg/Matrix.h index 68414b4dc..7bcb092b8 100644 --- a/lib/linalg/Matrix.h +++ b/lib/linalg/Matrix.h @@ -380,30 +380,6 @@ class Matrix return result; } - /** - * @brief Multiplies this Matrix with other and returns the product, - * pointer version. - * - * Supports multiplication of an undistributed Matrix and Vector - * returning an undistributed Vector, and multiplication of a distributed - * Matrix and an undistributed Vector returning a distributed Vector. - * - * @pre other != 0 - * @pre !other->distributed() - * @pre numColumns() == other->dim() - * - * @param[in] other The Vector to multiply with this. - * - * @return The product Vector. - */ - Vector* - mult( - const Vector* other) const - { - CAROM_VERIFY(other != 0); - return mult(*other); - } - /** * @brief Multiplies this Matrix with other and fills result with the * answer. diff --git a/lib/linalg/Vector.cpp b/lib/linalg/Vector.cpp index 490ee31ca..7cb11ff6b 100644 --- a/lib/linalg/Vector.cpp +++ b/lib/linalg/Vector.cpp @@ -164,21 +164,6 @@ Vector::transform(Vector& result, transformer(d_dim, result.d_vec); } -void -Vector::transform(Vector*& result, - std::function transformer) const { - // If the result has not been allocated then do so. Otherwise size it - // correctly. - if (result == 0) { - result = new Vector(d_dim, d_distributed); - } - else { - result->setSize(d_dim); - result->d_distributed = d_distributed; - } - transformer(d_dim, result->d_vec); -} - Vector& Vector::transform( std::function @@ -200,24 +185,6 @@ Vector::transform(Vector& result, delete origVector; } -void -Vector::transform(Vector*& result, - std::function - transformer) const { - // If the result has not been allocated then do so. Otherwise size it - // correctly. - Vector* origVector = new Vector(*this); - if (result == 0) { - result = new Vector(d_dim, d_distributed); - } - else { - result->setSize(d_dim); - result->d_distributed = d_distributed; - } - transformer(d_dim, origVector->d_vec, result->d_vec); - delete origVector; -} - double Vector::inner_product( const Vector& other) const @@ -248,7 +215,7 @@ Vector::norm() const double Vector::norm2() const { - double norm2 = inner_product(this); + double norm2 = inner_product(*this); return norm2; } @@ -262,30 +229,6 @@ Vector::normalize() return Norm; } -void -Vector::plus( - const Vector& other, - Vector*& result) const -{ - CAROM_ASSERT(result == 0 || result->distributed() == distributed()); - CAROM_ASSERT(distributed() == other.distributed()); - CAROM_VERIFY(dim() == other.dim()); - - // If the result has not been allocated then do so. Otherwise size it - // correctly. - if (result == 0) { - result = new Vector(d_dim, d_distributed); - } - else { - result->setSize(d_dim); - } - - // Do the addition. - for (int i = 0; i < d_dim; ++i) { - result->d_vec[i] = d_vec[i] + other.d_vec[i]; - } -} - void Vector::plus( const Vector& other, diff --git a/lib/linalg/Vector.h b/lib/linalg/Vector.h index 0f7624f65..72d3b6c94 100644 --- a/lib/linalg/Vector.h +++ b/lib/linalg/Vector.h @@ -16,6 +16,7 @@ #define included_Vector_h #include "utils/Utilities.h" +#include #include #include @@ -204,21 +205,6 @@ class Vector std::function transformer) const; - /** - * @brief Transform a vector using a supplied function and store the - * results in another vector. - * - * @param[out] result A vector which will store the transformed result. - * - * @param[in] transformer A transformer function which takes in as input a - * size and transforms the origVector and stores the result in - * resultVector. - */ - void - transform(Vector*& result, - std::function - transformer) const; - /** * @brief Sets the length of the vector and reallocates storage if * needed. All values are initialized to zero. @@ -228,8 +214,7 @@ class Vector * the Vector on this processor. */ void - setSize( - int dim) + setSize(int dim) { if (dim > d_alloc_size) { if (!d_owns_data) { @@ -284,27 +269,6 @@ class Vector inner_product( const Vector& other) const; - /** - * @brief Inner product, pointer version. - * - * For distributed Vectors this is a parallel operation. - * - * @pre other != 0 - * @pre dim() == other->dim() - * @pre distributed() == other->distributed() - * - * @param[in] other The Vector to form the inner product with this. - * - * @return The inner product of this and other. - */ - double - inner_product( - const Vector* other) const - { - CAROM_VERIFY(other != 0); - return inner_product(*other); - } - /** * @brief Form the norm of this. * @@ -345,52 +309,14 @@ class Vector * * @return this + other */ - Vector* - plus( - const Vector& other) const - { - Vector* result = 0; - plus(other, result); - return result; - } - - /** - * @brief Adds other and this and returns the result, pointer version. - * - * @pre other != 0 - * @pre distributed() == other->distributed() - * @pre dim() == other->dim() - * - * @param[in] other The other summand. - * - * @return this + other - */ - Vector* - plus( - const Vector* other) const + std::unique_ptr + plus(const Vector& other) const { - CAROM_VERIFY(other != 0); - return plus(*other); + Vector *result = new Vector(d_dim, d_distributed); + plus(other, *result); + return std::unique_ptr(result); } - /** - * @brief Adds other and this and fills result with the answer. - * - * Result will be allocated if unallocated or resized appropriately if - * already allocated. - * - * @pre result == 0 || result->distributed() == distributed() - * @pre distributed() == other.distributed() - * @pre dim() == other.dim() - * - * @param[in] other The other summand. - * @param[out] result this + other - */ - void - plus( - const Vector& other, - Vector*& result) const; - /** * @brief Adds other and this and fills result with the answer. * diff --git a/lib/linalg/svd/IncrementalSVD.cpp b/lib/linalg/svd/IncrementalSVD.cpp index 5c5e912c4..556a411aa 100644 --- a/lib/linalg/svd/IncrementalSVD.cpp +++ b/lib/linalg/svd/IncrementalSVD.cpp @@ -288,13 +288,13 @@ IncrementalSVD::buildIncrementalSVD( Vector* l = d_basis->transposeMult(u_vec); // basisl = basis * l - Vector* basisl = d_basis->mult(l); + Vector* basisl = d_basis->mult(*l); // Computing as k = sqrt(u.u - 2.0*l.l + basisl.basisl) // results in catastrophic cancellation, and must be avoided. // Instead we compute as k = sqrt((u-basisl).(u-basisl)). Vector* e_proj = u_vec.minus(basisl); - double k = e_proj->inner_product(e_proj); + double k = e_proj->inner_product(*e_proj); delete e_proj; if (k <= 0) { diff --git a/lib/linalg/svd/IncrementalSVDBrand.cpp b/lib/linalg/svd/IncrementalSVDBrand.cpp index b4648ab58..de32c67c1 100644 --- a/lib/linalg/svd/IncrementalSVDBrand.cpp +++ b/lib/linalg/svd/IncrementalSVDBrand.cpp @@ -148,8 +148,13 @@ IncrementalSVDBrand::buildIncrementalSVD( // (accurate down to the machine precision) Vector u_vec(u, d_dim, true); Vector e_proj(u, d_dim, true); - e_proj -= *(d_U->mult(d_U->transposeMult(e_proj))); // Gram-Schmidt - e_proj -= *(d_U->mult(d_U->transposeMult(e_proj))); // Re-orthogonalization + + Vector *tmp = d_U->transposeMult(e_proj); + + e_proj -= *(d_U->mult(*tmp)); // Gram-Schmidt + e_proj -= *(d_U->mult(*tmp)); // Re-orthogonalization + + delete tmp; double k = e_proj.inner_product(e_proj); if (k <= 0) {