diff --git a/src/Algorithm/LinearSolvers/IpKLUSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpKLUSolverInterface.cpp index 8e45bfa3..45e764ff 100644 --- a/src/Algorithm/LinearSolvers/IpKLUSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpKLUSolverInterface.cpp @@ -16,7 +16,7 @@ namespace Ipopt static const Index dbg_verbosity = 0; #endif -KLUSolverInterface::KLUSolverInterface() : _val(NULL), _Numeric(NULL) +KLUSolverInterface::KLUSolverInterface() : _val(NULL), _Symbolic(NULL), _Numeric(NULL) { DBG_START_METH("KLUSolverInterface::KLUSolverInterface()", dbg_verbosity); _rcond_val = 1e-128; diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index 81f832a8..1443d033 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -7,8 +7,10 @@ #include "IpReSolveSolverInterface.hpp" #include "IpoptConfig.h" +#include #include -#include +#include +#include namespace Ipopt { @@ -16,66 +18,80 @@ namespace Ipopt static const Index dbg_verbosity = 0; #endif -ReSolveSolverInterface::ReSolveSolverInterface() : val_(NULL) -{ - DBG_START_METH("ReSolveSolverInterface::ReSolveSolverInterface()", dbg_verbosity); - rcond_val_ = 1e-128; - factor_by_t_ = 2; - initialized_resolve_ = false; - initialized_ = false; - -#if RESOLVE_WITH_GPU - printf("Resolve with GPU.\n"); -# if RESOLVE_WITH_CUDA - printf("Resolve with CUDA.\n"); - NVMLHelper::getAvailableGPUMemory(); -# else - printf("Resolve with HIP.\n"); +static_assert( + std::is_same::value, + "Ipopt and ReSolve must use the same floating-point type."); + +static_assert( + std::is_same::value, + "Ipopt and ReSolve must use the same index type."); + +ReSolveSolverInterface::ReSolveSolverInterface() + : nonzeros_(0), + initialized_(false), + initialized_resolve_(false), + ndim_(0), + val_(NULL), + numneg_(0), + re_factorize_(false), + factorize_(false), + method_(resolve_klu), + n_iteration_(0), + k_(1), + pivot_tol_(0.1), + ordering_(1), + halt_if_singular_(false), + rcond_val_(1e-128), + use_rcond_(false), + resolve_KLU_(NULL), + workspace_CPU_(NULL), + A_(NULL), + vec_rhs_(NULL), + vec_x_(NULL), +#ifdef RESOLVE_USE_GPU + workspace_GPU_(NULL), + resolve_Rf_(NULL), +# ifdef RESOLVE_USE_CUDA + resolve_GLU_(NULL), # endif -#else - printf("Resolve with CPU. CUDA or HIP Unavailable.\n"); + GS_(NULL), + resolve_preconditioner_(NULL), + resolve_FGMRES_(NULL), #endif + matrix_handler_(NULL), + vector_handler_(NULL) +{ + DBG_START_METH("ReSolveSolverInterface::ReSolveSolverInterface()", dbg_verbosity); } ReSolveSolverInterface::~ReSolveSolverInterface() { + DBG_START_METH("ReSolveSolverInterface::~ReSolveSolverInterface()", dbg_verbosity); +#ifdef RESOLVE_USE_GPU + delete resolve_FGMRES_; + delete resolve_preconditioner_; + delete GS_; -// printf("Begin of Destructor\n"); -// NVMLHelper::getAvailableGPUMemory(); +# ifdef RESOLVE_USE_CUDA + delete resolve_GLU_; +# endif + delete resolve_Rf_; +#endif - DBG_START_METH("ReSolveSolverInterface::~ReSolveSolverInterface()", dbg_verbosity); - delete[] val_; + delete resolve_KLU_; delete matrix_handler_; delete vector_handler_; - delete resolve_KLU_; -#if RESOLVE_WITH_GPU +#ifdef RESOLVE_USE_GPU delete workspace_GPU_; - - if (method_ == resolve_rf || method_ == resolve_rf_fgmres) - { - delete resolve_Rf_; - } -# if RESOLVE_WITH_CUDA - else if (method_ == resolve_glu) - { - delete resolve_GLU_; - } -# endif - - if (method_ == resolve_rf_fgmres) - { - delete GS_; - delete resolve_FGMRES_; - } #endif + delete workspace_CPU_; delete vec_rhs_; delete vec_x_; delete A_; - -// printf("End of Destructor\n"); + delete[] val_; } void ReSolveSolverInterface::RegisterOptions(SmartPtr roptions) @@ -86,8 +102,8 @@ void ReSolveSolverInterface::RegisterOptions(SmartPtr roption options.push_back(resolve_klu); descrs.push_back("Use KLU"); -#if RESOLVE_WITH_GPU -# if RESOLVE_WITH_CUDA +#ifdef RESOLVE_USE_GPU +# ifdef RESOLVE_USE_CUDA options.push_back(resolve_glu); descrs.push_back("Use GLU"); # endif @@ -120,28 +136,17 @@ void ReSolveSolverInterface::RegisterOptions(SmartPtr roption "ordering if P and Q are NULL), or 3 for the user order function.", false); - roptions->AddIntegerOption("resolve_btf", // - "Use BTF", // - 1, // - "if nonzero, then BTF is used to permute the input matrix into block upper triangular form.", // - false); - - roptions->AddIntegerOption("resolve_scale", // - "Whether or not the matrix should be scaled", // - 2, // - "If scale < 0, then no scaling is performed and the input matrix is not checked for errors. If scale >= 0, the input matrix is check for errors. If scale=0, then no scaling is performed. If scale=1, then each row of A is " - "divided by the sum of the absolute values in that row. If scale=2, then each row of A is divided by the maximum absolute value in that row. Default: 2.", // - false); - roptions->AddBoolOption("resolve_halt_if_singular", // "how to handle a singular matrix", // false, // Default is False in ReSolve "FALSE: keep going, TRUE: stop quickly.", false); - roptions->AddIntegerOption("resolve_n_skip_refactoring", // + roptions->AddLowerBoundedIntegerOption("resolve_n_skip_refactoring", // "How many iterations to skip refactoring", // 1, // - "Integer, Start Refactoring after k-th iteration", // + 1, // + "Number of initial KLU iterations before switching to the selected refactorization method. " // + "At least one KLU iteration is required to construct the factors and permutations used during setup.", // false); roptions->AddBoolOption("resolve_use_rcond", // @@ -157,35 +162,49 @@ void ReSolveSolverInterface::RegisterOptions(SmartPtr roption bool ReSolveSolverInterface::InitializeImpl(const OptionsList& options, const std::string& prefix) { + options.GetNumericValue("resolve_tol", pivot_tol_, prefix); - // printf("ReSolveSolverInterface::InitializeImpl is Called\n"); - - Number tol; - options.GetNumericValue("resolve_tol", tol, prefix); - - Index order_method; - options.GetIntegerValue("resolve_ordering", order_method, prefix); - - Index btf; - options.GetIntegerValue("resolve_btf", btf, prefix); - - Index scale; - options.GetIntegerValue("resolve_scale", scale, prefix); + options.GetIntegerValue("resolve_ordering", ordering_, prefix); Index n_skip_refactoring; options.GetIntegerValue("resolve_n_skip_refactoring", n_skip_refactoring, prefix); + + // ReSolve's GPU refactorization methods require an initial KLU solve + // to construct the factors and permutations used during setup. k_ = n_skip_refactoring; - bool halt_if_singular; - options.GetBoolValue("resolve_halt_if_singular", halt_if_singular, prefix); + options.GetBoolValue("resolve_halt_if_singular", halt_if_singular_, prefix); - std::string method; - options.GetStringValue("resolve_method", method, prefix); - method_ = method; + options.GetStringValue("resolve_method", method_, prefix); options.GetNumericValue("resolve_rcond_val", rcond_val_, prefix); options.GetBoolValue("resolve_use_rcond", use_rcond_, prefix); + bool method_available = (method_ == resolve_klu); + +#ifdef RESOLVE_USE_GPU + method_available = + method_available + || method_ == resolve_rf + || method_ == resolve_rf_fgmres; + +# ifdef RESOLVE_USE_CUDA + method_available = + method_available + || method_ == resolve_glu; +# endif +#endif + + if (!method_available) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve method '%s' is not available for this build.\n", + method_.c_str()); + return false; + } + return true; } @@ -195,11 +214,13 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no DBG_START_METH("ReSolveSolverInterface::InitializeStructure", dbg_verbosity); ESymSolverStatus retval = SYMSOLVER_SUCCESS; - printf("ReSolveSolverInterface::InitializeStructure Called: dim: %d, nonzeros %d\n", dim, nonzeros); if (!initialized_resolve_) { resolve_KLU_ = new ReSolve::LinSolverDirectKLU(); + resolve_KLU_->setPivotThreshold(pivot_tol_); + resolve_KLU_->setOrdering(static_cast(ordering_)); + resolve_KLU_->setHaltIfSingular(halt_if_singular_); if (method_ == resolve_klu) { @@ -208,67 +229,130 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no vector_handler_ = new ReSolve::VectorHandler(workspace_CPU_); } -#if RESOLVE_WITH_GPU - workspace_GPU_ = new workspace_type(); - workspace_GPU_->initializeHandles(); - - matrix_handler_ = new ReSolve::MatrixHandler(workspace_GPU_); - vector_handler_ = new ReSolve::VectorHandler(workspace_GPU_); +#ifdef RESOLVE_USE_GPU + else + { + workspace_GPU_ = new workspace_type(); + workspace_GPU_->initializeHandles(); + matrix_handler_ = new ReSolve::MatrixHandler(workspace_GPU_); + vector_handler_ = new ReSolve::VectorHandler(workspace_GPU_); - if (method_ == resolve_rf || method_ == resolve_rf_fgmres) - { - resolve_Rf_ = new rf_solver(workspace_GPU_); - } -# if RESOLVE_WITH_CUDA - else if (method_ == resolve_glu) - { - resolve_GLU_ = new ReSolve::LinSolverDirectCuSolverGLU(workspace_GPU_); - } + if (method_ == resolve_rf || method_ == resolve_rf_fgmres) + { + resolve_Rf_ = new rf_solver(workspace_GPU_); + } +# ifdef RESOLVE_USE_CUDA + else if (method_ == resolve_glu) + { + resolve_GLU_ = + new ReSolve::LinSolverDirectCuSolverGLU(workspace_GPU_); + } # endif - - if (method_ == resolve_rf_fgmres) - { - GS_ = new ReSolve::GramSchmidt(vector_handler_, ReSolve::GramSchmidt::CGS2); - resolve_FGMRES_ = new ReSolve::LinSolverIterativeFGMRES(matrix_handler_, vector_handler_, GS_); + + if (method_ == resolve_rf_fgmres) + { + GS_ = new ReSolve::GramSchmidt( + vector_handler_, + ReSolve::GramSchmidt::CGS2); + + resolve_FGMRES_ = new ReSolve::LinSolverIterativeFGMRES( + matrix_handler_, + vector_handler_, + GS_); + + resolve_preconditioner_ = + new ReSolve::PreconditionerLU(resolve_Rf_); + + if (resolve_FGMRES_->setPreconditioner( + resolve_preconditioner_) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to attach the ReSolve RF preconditioner " + "to FGMRES.\n"); + return SYMSOLVER_FATAL_ERROR; + } + } } #endif + initialized_resolve_ = true; } - initialized_resolve_ = true; - // Store size for later use ndim_ = dim; nonzeros_ = nonzeros; - printf("Using Refactorization after %d iterations\n\n", k_); if (!initialized_) { A_ = new ReSolve::matrix::Csr(dim, dim, nonzeros); - if (val_ != NULL) + delete[] val_; + val_ = new Number[static_cast(nonzeros)]; + + // ReSolve borrows the matrix storage and does not take ownership. + // Ipopt owns ia/ja, while this interface owns val_. + if( A_->setDataPointers(const_cast(ia), const_cast(ja), val_, ReSolve::memory::HOST) != 0) { - delete[] val_; + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to attach Ipopt matrix storage to ReSolve.\n"); + return SYMSOLVER_FATAL_ERROR; } - val_ = new Number[nonzeros]; - - A_->setDataPointers(const_cast(ia), const_cast(ja), val_, ReSolve::memory::HOST); - resolve_KLU_->setup(A_); - - vec_rhs_ = new ReSolve::vector::Vector(A_->getNumRows()); - vec_x_ = new ReSolve::vector::Vector(A_->getNumRows()); +#ifdef RESOLVE_USE_GPU + if( method_ != resolve_klu + && A_->allocateMatrixData(ReSolve::memory::DEVICE) != 0 ) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to allocate ReSolve matrix device storage.\n"); + return SYMSOLVER_FATAL_ERROR; + } +#endif - vec_x_->allocate(ReSolve::memory::HOST); // for KLU - // vec_x_->allocate(ReSolve::memory::DEVICE); + if( resolve_KLU_->setup(A_) != 0 ) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to set up the ReSolve KLU solver.\n"); + return SYMSOLVER_FATAL_ERROR; + } + + vec_rhs_ = new ReSolve::vector::Vector(A_->getNumRows()); + vec_x_ = new ReSolve::vector::Vector(A_->getNumRows()); + + if( vec_rhs_->allocate(ReSolve::memory::HOST) != 0 + || vec_x_->allocate(ReSolve::memory::HOST) != 0 ) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to allocate ReSolve host vectors.\n"); + return SYMSOLVER_FATAL_ERROR; + } + +#ifdef RESOLVE_USE_GPU + if( method_ != resolve_klu ) + { + if( vec_rhs_->allocate(ReSolve::memory::DEVICE) != 0 + || vec_x_->allocate(ReSolve::memory::DEVICE) != 0 ) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to allocate ReSolve device vectors.\n"); + return SYMSOLVER_FATAL_ERROR; + } + } +#endif } - factorize_ = true; n_iteration_ = 0; initialized_ = true; - pivtol_changed_ = false; - -// printf("After Initialize\n"); -// NVMLHelper::getAvailableGPUMemory(); return retval; } @@ -277,27 +361,51 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index { DBG_START_METH("ReSolveSolverInterface::MultiSolve", dbg_verbosity); + std::vector solution_vals; + + if (nrhs > 1) + { + solution_vals.resize( + static_cast(nrhs) * static_cast(ndim_)); + } + int status; - int status_refactor = 0; bool full_factor_done = false; - // Get Data from CPU and update the A Matrix - A_->copyDataFrom(A_->getRowData(ReSolve::memory::HOST), A_->getColData(ReSolve::memory::HOST), A_->getValues(ReSolve::memory::HOST), ReSolve::memory::HOST, ReSolve::memory::DEVICE); - - // FACTORIZE + (void)ia; + (void)ja; + (void)check_NegEVals; + (void)numberOfNegEVals; - // Every factor_by_t_ iteration do a Factorization!!! - if (n_iteration_ % factor_by_t_ == 0) + // Ipopt updates val_ directly, so mark the host matrix current before + // synchronizing updated values to the GPU. + if (new_matrix) { - // factorize_ = true; - // re_factorize_ = true; + if( A_->setUpdated(ReSolve::memory::HOST) != 0 ) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to mark the ReSolve host matrix as updated.\n"); + return SYMSOLVER_FATAL_ERROR; + } + +#ifdef RESOLVE_USE_GPU + if( method_ != resolve_klu + && A_->syncData(ReSolve::memory::DEVICE) != 0 ) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to synchronize the ReSolve matrix " + "to the device.\n"); + return SYMSOLVER_FATAL_ERROR; + } +#endif } - if (factorize_ && (new_matrix || re_factorize_)) - { - // printf("Iteration: %d: Performing KLU Factorization\n", n_iteration_); - + if( n_iteration_ == 0){ // Symbolic Factorization if (HaveIpData()) { @@ -306,15 +414,21 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index status = resolve_KLU_->analyze(); if (status != 0) { - printf("Symbolic_ factorization crashed with Common_.status = %d \n", status); - printf("%s:0 Singular\n", __func__); + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve KLU symbolic analysis failed with status %d.\n", + status); return SYMSOLVER_SINGULAR; } if (HaveIpData()) { IpData().TimingStats().LinearSystemSymbolicFactorization().End(); } + } + if (factorize_ && (new_matrix || re_factorize_)) + { // perform the factorization if (HaveIpData()) { @@ -322,35 +436,47 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } // First Factorization is always done by KLU - std::cout << "%" << n_iteration_ << "%" << "FULL FACTORIZATIOM" << std::endl; status = resolve_KLU_->factorize(); full_factor_done = true; if (status != 0) { - DBG_PRINT((1, "FACTORIZATION FAILED!\n")); - printf("%s: Singular\n", __func__); + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve KLU factorization failed with status %d.\n", + status); return SYMSOLVER_SINGULAR; // Matrix singular or error occurred } - // printf("Iteration: %d: Done KLU Factorization\n", n_iteration_); // GLU can be setup as early as possible if (n_iteration_ == k_ - 1) { -#if RESOLVE_WITH_CUDA +#ifdef RESOLVE_USE_CUDA if (method_ == resolve_glu) { - printf("Iteration: %d: Setting Up GLU\n", n_iteration_); - ReSolve::matrix::Sparse* L = resolve_KLU_->getLFactor(); ReSolve::matrix::Sparse* U = resolve_KLU_->getUFactor(); - if (L == nullptr) - { - printf("ERROR"); - } ReSolve::index_type* P = resolve_KLU_->getPOrdering(); ReSolve::index_type* Q = resolve_KLU_->getQOrdering(); - resolve_GLU_->setup(A_, L, U, P, Q); + if (L == nullptr || U == nullptr || P == nullptr || Q == nullptr) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to obtain KLU factors or permutations for ReSolve CUDA GLU setup.\n"); + return SYMSOLVER_FATAL_ERROR; + } + status = resolve_GLU_->setup(A_, L, U, P, Q); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve CUDA GLU setup failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } } #endif } @@ -364,22 +490,6 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // Stop doing factorization after iteration _k factorize_ = (n_iteration_ >= (k_ - 1)) ? false : true; re_factorize_ = false; - // printf("Iteration: %d: Ending Factorization Section\n", n_iteration_); - } - - if (pivtol_changed_) - { - DBG_PRINT((1, "Pivot tolerance has changed.\n")); - pivtol_changed_ = false; - // If the pivot tolerance has been changed but the matrix is not - // new, we have to request the values for the matrix again to do - // the factorization again. - if (!new_matrix) - { - DBG_PRINT((1, "Ask caller to call again.\n")); - factorize_ = true; - return SYMSOLVER_CALL_AGAIN; - } } // REFACTORIZE @@ -393,42 +503,61 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index IpData().TimingStats().LinearSystemFactorization().Start(); } -# if RESOLVE_WITH_CUDA +# ifdef RESOLVE_USE_CUDA // Actual Refactorize if (method_ == resolve_glu) { status = resolve_GLU_->refactorize(); if (status != 0) { - std::cout << "CUSOLVER GLU refactorization status: " << status << std::endl; + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve CUDA GLU refactorization failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; } } else if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - status_refactor = resolve_Rf_->refactorize(); - if (status != 0) + int status_refactor = resolve_Rf_->refactorize(); + if (status_refactor != 0) { - std::cout << "CUSOLVER RF refactorization status: " << status_refactor << std::endl; + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve CUDA RF refactorization failed with status %d.\n", + status_refactor); + return SYMSOLVER_FATAL_ERROR; } } -# else +# elif defined(RESOLVE_USE_HIP) if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - int status = resolve_Rf_->refactorize(); + status = resolve_Rf_->refactorize(); if (status != 0) { - std::cout << "ROCSOLVER RF refactorization status: " << status << std::endl; + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve HIP RF refactorization failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; } } #endif if (method_ == resolve_klu) { - std::cout << "%" << n_iteration_ << "%" << "RE-FACTORIZATIOM" << std::endl; status = resolve_KLU_->refactorize(); if (status != 0) { - std::cout << "KLU refactorization status: " << status << std::endl; + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve KLU refactorization failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; } } re_factorize_ = false; @@ -452,12 +581,14 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index if (use_rcond_) { Number rcond_val = resolve_KLU_->getMatrixConditionNumber(); - printf("RCond: %12.8e\n", rcond_val); if (rcond_val < rcond_val_) { if (full_factor_done) { - printf("%s:1 Singular\n", __func__); + Jnlst().Printf( + J_DETAILED, + J_LINEAR_ALGEBRA, + "ReSolve KLU reciprocal condition estimate is below the configured threshold.\n"); return SYMSOLVER_SINGULAR; } else @@ -466,70 +597,163 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // and do full factorization factorize_ = true; re_factorize_ = true; - printf("%s:1 Need to do full factorization again.\n", __func__); - DBG_PRINT((1, "Ask caller to call again.\n")) + Jnlst().Printf( + J_DETAILED, + J_LINEAR_ALGEBRA, + "ReSolve KLU reciprocal condition estimate is below the configured threshold; " + "requesting a full factorization.\n"); return SYMSOLVER_CALL_AGAIN; } } } - // Copy rhs_vals to vec_rhs - vec_rhs_->copyDataFrom(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::HOST); - status = resolve_KLU_->solve(vec_rhs_, vec_x_); - if (status != 0) + for (Index irhs = 0; irhs < nrhs; ++irhs) { - std::cout << "KLU solve status: " << status << std::endl; + Number* rhs = rhs_vals + irhs * ndim_; + // Copy the current RHS to ReSolve + if (vec_rhs_->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to copy rhs data to ReSolve host storage.\n"); + return SYMSOLVER_FATAL_ERROR; + } + status = resolve_KLU_->solve(vec_rhs_, vec_x_); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve KLU solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } + if (nrhs > 1) + { + Number* solution = solution_vals.data() + irhs * ndim_; + if (vec_x_->copyToExternal(solution, ReSolve::memory::HOST, ReSolve::memory::HOST) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to copy the ReSolve solution to temporary storage.\n"); + return SYMSOLVER_FATAL_ERROR; + } + } } } // Setup RF here if (n_iteration_ == (k_ - 1)) { -# if RESOLVE_WITH_CUDA +# ifdef RESOLVE_USE_CUDA if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - printf("Iteration: %d: Setting up %s\n", n_iteration_, method_.c_str()); - ReSolve::matrix::Csc* L_csc = (ReSolve::matrix::Csc*)resolve_KLU_->getLFactor(); - ReSolve::matrix::Csc* U_csc = (ReSolve::matrix::Csc*)resolve_KLU_->getUFactor(); - ReSolve::matrix::Csr* L = new ReSolve::matrix::Csr(L_csc->getNumRows(), L_csc->getNumColumns(), L_csc->getNnz()); - ReSolve::matrix::Csr* U = new ReSolve::matrix::Csr(U_csc->getNumRows(), U_csc->getNumColumns(), U_csc->getNnz()); - matrix_handler_->csc2csr(L_csc, L, ReSolve::memory::DEVICE); - matrix_handler_->csc2csr(U_csc, U, ReSolve::memory::DEVICE); - if (L == nullptr) - { - printf("ERROR"); - } + ReSolve::matrix::Sparse* L = resolve_KLU_->getLFactor(); + ReSolve::matrix::Sparse* U = resolve_KLU_->getUFactor(); ReSolve::index_type* P = resolve_KLU_->getPOrdering(); ReSolve::index_type* Q = resolve_KLU_->getQOrdering(); - resolve_Rf_->setup(A_, L, U, P, Q); - delete L; - delete U; + if (L == nullptr || U == nullptr || P == nullptr || Q == nullptr) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to obtain KLU factors or permutations for ReSolve CUDA RF setup.\n"); + return SYMSOLVER_FATAL_ERROR; + } + + status = resolve_Rf_->setup(A_, L, U, P, Q); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve CUDA RF setup failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } + + status = resolve_Rf_->refactorize(); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve CUDA RF initial refactorization failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } } - - if (method_ == resolve_rf_fgmres) + + if (method_ == resolve_rf_fgmres && resolve_FGMRES_->setup(A_) != 0) { - resolve_FGMRES_->setup(A_); - resolve_FGMRES_->setupPreconditioner("CuSolverRf", resolve_Rf_); + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to set up ReSolve FGMRES.\n"); + return SYMSOLVER_FATAL_ERROR; } -# else - +# elif defined(RESOLVE_USE_HIP) if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - ReSolve::matrix::Csc* L = (ReSolve::matrix::Csc*)resolve_KLU_->getLFactor(); - ReSolve::matrix::Csc* U = (ReSolve::matrix::Csc*)resolve_KLU_->getUFactor(); + ReSolve::matrix::Sparse* L = resolve_KLU_->getLFactor(); + ReSolve::matrix::Sparse* U = resolve_KLU_->getUFactor(); ReSolve::index_type* P = resolve_KLU_->getPOrdering(); ReSolve::index_type* Q = resolve_KLU_->getQOrdering(); - vec_rhs_->copyDataFrom(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - resolve_Rf_->setup(A_, L, U, P, Q, vec_rhs_); + if (L == nullptr || U == nullptr || P == nullptr || Q == nullptr) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to obtain KLU factors or permutations for ReSolve HIP RF setup.\n"); + return SYMSOLVER_FATAL_ERROR; + } + + if (vec_rhs_->copyFromExternal( + rhs_vals, + ReSolve::memory::HOST, + ReSolve::memory::DEVICE) + != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to copy rhs data to ReSolve device storage.\n"); + return SYMSOLVER_FATAL_ERROR; + } + + status = resolve_Rf_->setup(A_, L, U, P, Q, vec_rhs_); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve HIP RF setup failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } + + status = resolve_Rf_->refactorize(); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve HIP RF initial refactorization failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } } - if (method_ == resolve_rf_fgmres) + if (method_ == resolve_rf_fgmres && resolve_FGMRES_->setup(A_) != 0) { - std::cout << "about to set FGMRES" << std::endl; - GS_->setup(A_->getNumRows(), resolve_FGMRES_->getRestart()); - resolve_FGMRES_->setup(A_); - resolve_FGMRES_->setupPreconditioner("LU", resolve_Rf_); + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to set up ReSolve FGMRES.\n"); + return SYMSOLVER_FATAL_ERROR; } # endif } @@ -537,75 +761,190 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index else // After k iteration only solve { // Every 10 iteration setup again with full factorization? Already factorization done. -# if RESOLVE_WITH_CUDA - +# ifdef RESOLVE_USE_CUDA if (method_ == resolve_glu || method_ == resolve_rf || method_ == resolve_rf_fgmres) { - // Copy rhs_vals to vec_rhs cuda - vec_rhs_->copyDataFrom(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - - if (method_ == resolve_glu) - { - status = resolve_GLU_->solve(vec_rhs_, vec_x_); - if (status != 0) + for (Index irhs = 0; irhs < nrhs; ++irhs) + { + Number* rhs = rhs_vals + irhs * ndim_; + // Copy the current RHS to ReSolve device storage + if (vec_rhs_->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE) != 0) { - std::cout << "GLU solve status: " << status << std::endl; + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to copy rhs data to ReSolve device storage.\n"); + return SYMSOLVER_FATAL_ERROR; } - } - else if (method_ == resolve_rf) - { - status = resolve_Rf_->solve(vec_rhs_, vec_x_); - if (status != 0) + + if (method_ == resolve_glu) { - std::cout << "RF solve status: " << status << std::endl; + status = resolve_GLU_->solve(vec_rhs_, vec_x_); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve CUDA GLU solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } } - } - else if (method_ == resolve_rf_fgmres) - { - int status = resolve_Rf_->solve(vec_rhs_, vec_x_); - if (status != 0) + else if (method_ == resolve_rf) { - std::cout << "RF solve status: " << status << std::endl; + status = resolve_Rf_->solve(vec_rhs_, vec_x_); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve CUDA RF solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } } + else if (method_ == resolve_rf_fgmres) + { + status = resolve_Rf_->solve(vec_rhs_, vec_x_); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve CUDA RF initial solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } - resolve_FGMRES_->resetMatrix(A_); - status = resolve_FGMRES_->solve(vec_rhs_, vec_x_); - if (status != 0) + status = resolve_FGMRES_->resetMatrix(A_); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to reset the ReSolve CUDA FGMRES matrix " + "with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } + + status = resolve_FGMRES_->solve(vec_rhs_, vec_x_); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve CUDA FGMRES solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } + } + + // GPU solvers leave the current solution on the device; synchronize it + // to host memory before copying the solution back to Ipopt. + if (vec_x_->syncData(ReSolve::memory::HOST) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to synchronize the ReSolve solution to the host.\n"); + return SYMSOLVER_FATAL_ERROR; + } + if (nrhs > 1) { - std::cout << "RF_FGMRES solve status: " << status << std::endl; + Number* solution = solution_vals.data() + irhs * ndim_; + if (vec_x_->copyToExternal(solution, ReSolve::memory::HOST, ReSolve::memory::HOST) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to copy the ReSolve solution to temporary storage.\n"); + return SYMSOLVER_FATAL_ERROR; + } } } - - // Copy vec_x cuda to vec_x in cpu - vec_x_->copyDataFrom(vec_x_->getData(ReSolve::memory::DEVICE), ReSolve::memory::DEVICE, ReSolve::memory::HOST); matrix_handler_->setValuesChanged(true, ReSolve::memory::DEVICE); } -# else +# elif defined(RESOLVE_USE_HIP) if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - // Copy rhs_vals to vec_rhs cuda - vec_rhs_->copyDataFrom(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - - int status = resolve_Rf_->solve(vec_rhs_, vec_x_); - if (status != 0) - { - std::cout << "RF solve status: " << status << std::endl; - } - - if (method_ == resolve_rf_fgmres) + for (Index irhs = 0; irhs < nrhs; ++irhs) { + Number* rhs = rhs_vals + irhs * ndim_; + // Copy the current RHS to ReSolve device storage + if (vec_rhs_->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to copy rhs data to ReSolve device storage.\n"); + return SYMSOLVER_FATAL_ERROR; + } - resolve_FGMRES_->resetMatrix(A_); - vec_rhs_->copyDataFrom(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - status = resolve_FGMRES_->solve(vec_rhs_, vec_x_); + status = resolve_Rf_->solve(vec_rhs_, vec_x_); if (status != 0) { - std::cout << "RF_FGMRES solve status: " << status << std::endl; + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve HIP RF solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } + + if (method_ == resolve_rf_fgmres) + { + status = resolve_FGMRES_->resetMatrix(A_); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to reset the ReSolve HIP FGMRES matrix " + "with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } + + status = resolve_FGMRES_->solve(vec_rhs_, vec_x_); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve HIP FGMRES solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } } - } - // Copy vec_x cuda to vec_x in cpu - vec_x_->copyDataFrom(vec_x_->getData(ReSolve::memory::DEVICE), ReSolve::memory::DEVICE, ReSolve::memory::HOST); + // GPU solvers leave the current solution on the device; synchronize it + // to host memory before copying the solution back to Ipopt. + if (vec_x_->syncData(ReSolve::memory::HOST) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to synchronize the ReSolve solution to the host.\n"); + return SYMSOLVER_FATAL_ERROR; + } + if (nrhs > 1) + { + Number* solution = solution_vals.data() + irhs * ndim_; + + if (vec_x_->copyToExternal( + solution, + ReSolve::memory::HOST, + ReSolve::memory::HOST) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to copy the ReSolve solution to temporary storage.\n"); + return SYMSOLVER_FATAL_ERROR; + } + } + } } #endif @@ -615,12 +954,14 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index if (use_rcond_) { Number rcond_val = resolve_KLU_->getMatrixConditionNumber(); - printf("RCond: %12.8e\n", rcond_val); if (rcond_val < rcond_val_) { if (full_factor_done) { - printf("%s:2 Singular\n", __func__); + Jnlst().Printf( + J_DETAILED, + J_LINEAR_ALGEBRA, + "ReSolve KLU reciprocal condition estimate is below the configured threshold.\n"); return SYMSOLVER_SINGULAR; } else @@ -629,25 +970,76 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // and do full factorization factorize_ = true; re_factorize_ = true; - printf("Need to do full factorization again.\n"); - DBG_PRINT((1, "Ask caller to call again.\n")) + Jnlst().Printf( + J_DETAILED, + J_LINEAR_ALGEBRA, + "ReSolve KLU reciprocal condition estimate is below the configured threshold; " + "requesting a full factorization.\n"); return SYMSOLVER_CALL_AGAIN; } } } - // Copy rhs_vals to vec_rhs cuda - vec_rhs_->copyDataFrom(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::HOST); - int status = resolve_KLU_->solve(vec_rhs_, vec_x_); - if (status != 0) + for (Index irhs = 0; irhs < nrhs; ++irhs) { - std::cout << "KLU solve status: " << status << std::endl; + Number* rhs = rhs_vals + irhs * ndim_; + // Copy the current RHS to ReSolve host storage + if (vec_rhs_->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to copy rhs data to ReSolve host storage.\n"); + return SYMSOLVER_FATAL_ERROR; + } + status = resolve_KLU_->solve(vec_rhs_, vec_x_); + if (status != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve KLU solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; + } + if (nrhs > 1) + { + Number* solution = solution_vals.data() + irhs * ndim_; + + if (vec_x_->copyToExternal( + solution, + ReSolve::memory::HOST, + ReSolve::memory::HOST) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to copy the ReSolve solution to temporary storage.\n"); + return SYMSOLVER_FATAL_ERROR; + } + } } } } - // copy vec_x to rhs_vals - memcpy(rhs_vals, vec_x_->getData(ReSolve::memory::HOST), (ndim_) * sizeof(ReSolve::real_type)); + if (nrhs == 1) + { + if (vec_x_->copyToExternal(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::HOST) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to copy the Resolve solution to Ipopt.\n"); + return SYMSOLVER_FATAL_ERROR; + } + } + else + { + std::copy( + solution_vals.begin(), + solution_vals.end(), + rhs_vals); + } if (HaveIpData()) { @@ -674,7 +1066,8 @@ Index ReSolveSolverInterface::NumberOfNegEVals() const bool ReSolveSolverInterface::IncreaseQuality() { - return true; + // Dynamic quality increases are not currently supported by this interface. + return false; } } // namespace Ipopt diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp index ebcafa08..efbc91ef 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp @@ -9,29 +9,29 @@ #include #include -#include #include "IpoptConfig.h" +#include -#if RESOLVE_WITH_CUDA -#include "NVMLHelper.hpp" +#if defined(RESOLVE_USE_GPU) && !defined(RESOLVE_USE_CUDA) && !defined(RESOLVE_USE_HIP) +# error "ReSolve GPU support requires either CUDA or HIP." #endif -#include -#include #include #include #include #include #include +#include +#include #include #include -#if RESOLVE_WITH_GPU +#ifdef RESOLVE_USE_GPU #include -#if RESOLVE_WITH_CUDA +#ifdef RESOLVE_USE_CUDA #include #include using workspace_type = ReSolve::LinAlgWorkspaceCUDA; @@ -59,8 +59,8 @@ namespace Ipopt static const std::string resolve_klu = "klu"; -#if RESOLVE_WITH_GPU -#if RESOLVE_WITH_CUDA +#ifdef RESOLVE_USE_GPU +#ifdef RESOLVE_USE_CUDA static const std::string resolve_glu = "glu"; # endif static const std::string resolve_rf = "rf"; @@ -138,18 +138,19 @@ class ReSolveSolverInterface : public SparseSymLinearSolverInterface Index ndim_; ///< Number of dimensions Number* val_; ///< Storage for variables Index numneg_; ///< Number of negative pivots in last factorization - bool pivtol_changed_; ///< indicates if pivtol has been changed bool re_factorize_; bool factorize_; std::string method_; int n_iteration_; int k_; + Number pivot_tol_; + Index ordering_; + bool halt_if_singular_; + Number rcond_val_; bool use_rcond_; - int factor_by_t_; - ReSolve::LinSolverDirectKLU* resolve_KLU_; ReSolve::LinAlgWorkspaceCpu* workspace_CPU_; @@ -158,13 +159,14 @@ class ReSolveSolverInterface : public SparseSymLinearSolverInterface ReSolve::vector::Vector* vec_x_; -#if RESOLVE_WITH_GPU +#ifdef RESOLVE_USE_GPU workspace_type* workspace_GPU_; rf_solver* resolve_Rf_; -# if RESOLVE_WITH_CUDA +# ifdef RESOLVE_USE_CUDA ReSolve::LinSolverDirectCuSolverGLU* resolve_GLU_; # endif ReSolve::GramSchmidt* GS_; + ReSolve::PreconditionerLU* resolve_preconditioner_; ReSolve::LinSolverIterativeFGMRES* resolve_FGMRES_; #endif diff --git a/src/Algorithm/LinearSolvers/NVMLHelper.hpp b/src/Algorithm/LinearSolvers/NVMLHelper.hpp deleted file mode 100644 index bc10329e..00000000 --- a/src/Algorithm/LinearSolvers/NVMLHelper.hpp +++ /dev/null @@ -1,70 +0,0 @@ -#pragma once - -#include -#include - -class NVMLHelper -{ -public: - static int getAvailableGPUMemory() - { - nvmlReturn_t result; - unsigned int device_count, i; - nvmlDevice_t device; - - // Initialize NVML - result = nvmlInit(); - if (result != NVML_SUCCESS) - { - printf("Failed to initialize NVML: %s\n", nvmlErrorString(result)); - return 1; - } - - // Get the number of GPUs - result = nvmlDeviceGetCount(&device_count); - if (result != NVML_SUCCESS) - { - printf("Failed to get device count: %s\n", nvmlErrorString(result)); - nvmlShutdown(); - return 1; - } - - // Loop through each GPU and get its memory information - for (i = 0; i < device_count; i++) - { - result = nvmlDeviceGetHandleByIndex(i, &device); - if (result != NVML_SUCCESS) - { - printf("Failed to get handle for device %d: %s\n", i, nvmlErrorString(result)); - nvmlShutdown(); - return 1; - } - - // Get memory information - nvmlMemory_t memory_info; - result = nvmlDeviceGetMemoryInfo(device, &memory_info); - if (result != NVML_SUCCESS) - { - printf("Failed to get memory info for device %d: %s\n", i, nvmlErrorString(result)); - nvmlShutdown(); - return 1; - } - - // Print total, free, and used memory - printf("GPU %d:\n", i); - printf(" Total memory: %ld\n", memory_info.total / (1024 * 1024)); - printf(" Free memory: %ld\n", memory_info.free / (1024 * 1024)); - printf(" Used memory: %ld\n", memory_info.used / (1024 * 1024)); - } - - // Shutdown NVML - result = nvmlShutdown(); - if (result != NVML_SUCCESS) - { - printf("Failed to shutdown NVML: %s\n", nvmlErrorString(result)); - return 1; - } - - return 0; - } -}; \ No newline at end of file diff --git a/src/Makefile.am b/src/Makefile.am index c0a8df20..1dcc4ff7 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -214,7 +214,7 @@ if COIN_HAS_KLU libipopt_la_SOURCES += Algorithm/LinearSolvers/IpKLUSolverInterface.cpp endif -if COIN_HAS_KLU +if COIN_HAS_RESOLVE libipopt_la_SOURCES += Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp endif diff --git a/src/Makefile.in b/src/Makefile.in index e0852917..88578501 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -110,10 +110,11 @@ host_triplet = @host@ @HAVE_MA28_TRUE@@IPOPT_INT64_FALSE@@IPOPT_SINGLE_TRUE@am__append_4 = Algorithm/LinearSolvers/IpMa28sPartition.F @HAVE_MA28_TRUE@@IPOPT_INT64_FALSE@@IPOPT_SINGLE_FALSE@am__append_5 = Algorithm/LinearSolvers/IpMa28Partition.F @COIN_HAS_KLU_TRUE@am__append_6 = Algorithm/LinearSolvers/IpKLUSolverInterface.cpp -@HAVE_WSMP_TRUE@am__append_7 = Algorithm/LinearSolvers/IpWsmpSolverInterface.cpp Algorithm/LinearSolvers/IpIterativeWsmpSolverInterface.cpp -@COIN_HAS_MUMPS_TRUE@am__append_8 = Algorithm/LinearSolvers/IpMumpsSolverInterface.cpp -@COIN_HAS_SPRAL_TRUE@am__append_9 = Algorithm/LinearSolvers/IpSpralSolverInterface.cpp -@BUILD_INEXACT_TRUE@am__append_10 = \ +@COIN_HAS_RESOLVE_TRUE@am__append_7 = Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +@HAVE_WSMP_TRUE@am__append_8 = Algorithm/LinearSolvers/IpWsmpSolverInterface.cpp Algorithm/LinearSolvers/IpIterativeWsmpSolverInterface.cpp +@COIN_HAS_MUMPS_TRUE@am__append_9 = Algorithm/LinearSolvers/IpMumpsSolverInterface.cpp +@COIN_HAS_SPRAL_TRUE@am__append_10 = Algorithm/LinearSolvers/IpSpralSolverInterface.cpp +@BUILD_INEXACT_TRUE@am__append_11 = \ @BUILD_INEXACT_TRUE@ Algorithm/Inexact/IpInexactAlgBuilder.cpp \ @BUILD_INEXACT_TRUE@ Algorithm/Inexact/IpInexactCq.cpp \ @BUILD_INEXACT_TRUE@ Algorithm/Inexact/IpInexactData.cpp \ @@ -129,8 +130,7 @@ host_triplet = @host@ @BUILD_INEXACT_TRUE@ Algorithm/Inexact/IpIterativePardisoSolverInterface.cpp \ @BUILD_INEXACT_TRUE@ Algorithm/Inexact/IpIterativeSolverTerminationTester.cpp -@BUILD_JAVA_TRUE@am__append_11 = Interfaces/IpStdJInterface.cpp org_coinor_Ipopt.h -@COIN_HAS_RESOLVE_TRUE@am__append_12 = Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +@BUILD_JAVA_TRUE@am__append_12 = Interfaces/IpStdJInterface.cpp org_coinor_Ipopt.h subdir = src ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 am__aclocal_m4_deps = $(top_srcdir)/configure.ac @@ -187,11 +187,12 @@ am__dirstamp = $(am__leading_dot)dirstamp @HAVE_MA28_TRUE@@IPOPT_INT64_FALSE@@IPOPT_SINGLE_TRUE@am__objects_4 = Algorithm/LinearSolvers/IpMa28sPartition.lo @HAVE_MA28_TRUE@@IPOPT_INT64_FALSE@@IPOPT_SINGLE_FALSE@am__objects_5 = Algorithm/LinearSolvers/IpMa28Partition.lo @COIN_HAS_KLU_TRUE@am__objects_6 = Algorithm/LinearSolvers/IpKLUSolverInterface.lo -@HAVE_WSMP_TRUE@am__objects_7 = Algorithm/LinearSolvers/IpWsmpSolverInterface.lo \ +@COIN_HAS_RESOLVE_TRUE@am__objects_7 = Algorithm/LinearSolvers/IpReSolveSolverInterface.lo +@HAVE_WSMP_TRUE@am__objects_8 = Algorithm/LinearSolvers/IpWsmpSolverInterface.lo \ @HAVE_WSMP_TRUE@ Algorithm/LinearSolvers/IpIterativeWsmpSolverInterface.lo -@COIN_HAS_MUMPS_TRUE@am__objects_8 = Algorithm/LinearSolvers/IpMumpsSolverInterface.lo -@COIN_HAS_SPRAL_TRUE@am__objects_9 = Algorithm/LinearSolvers/IpSpralSolverInterface.lo -@BUILD_INEXACT_TRUE@am__objects_10 = \ +@COIN_HAS_MUMPS_TRUE@am__objects_9 = Algorithm/LinearSolvers/IpMumpsSolverInterface.lo +@COIN_HAS_SPRAL_TRUE@am__objects_10 = Algorithm/LinearSolvers/IpSpralSolverInterface.lo +@BUILD_INEXACT_TRUE@am__objects_11 = \ @BUILD_INEXACT_TRUE@ Algorithm/Inexact/IpInexactAlgBuilder.lo \ @BUILD_INEXACT_TRUE@ Algorithm/Inexact/IpInexactCq.lo \ @BUILD_INEXACT_TRUE@ Algorithm/Inexact/IpInexactData.lo \ @@ -206,8 +207,7 @@ am__dirstamp = $(am__leading_dot)dirstamp @BUILD_INEXACT_TRUE@ Algorithm/Inexact/IpInexactTSymScalingMethod.lo \ @BUILD_INEXACT_TRUE@ Algorithm/Inexact/IpIterativePardisoSolverInterface.lo \ @BUILD_INEXACT_TRUE@ Algorithm/Inexact/IpIterativeSolverTerminationTester.lo -@BUILD_JAVA_TRUE@am__objects_11 = Interfaces/IpStdJInterface.lo -@COIN_HAS_RESOLVE_TRUE@am__objects_12 = Algorithm/LinearSolvers/IpReSolveSolverInterface.lo +@BUILD_JAVA_TRUE@am__objects_12 = Interfaces/IpStdJInterface.lo am_libipopt_la_OBJECTS = Common/IpDebug.lo Common/IpJournalist.lo \ Common/IpObserver.lo Common/IpOptionsList.lo \ Common/IpRegOptions.lo Common/IpTaggedObject.lo \ @@ -285,7 +285,8 @@ am_libipopt_la_OBJECTS = Common/IpDebug.lo Common/IpJournalist.lo \ Interfaces/IpTNLPReducer.lo $(am__objects_1) $(am__objects_2) \ $(am__objects_3) $(am__objects_4) $(am__objects_5) \ $(am__objects_6) $(am__objects_7) $(am__objects_8) \ - $(am__objects_9) $(am__objects_10) $(am__objects_11) $(am__objects_12) + $(am__objects_9) $(am__objects_10) $(am__objects_11) \ + $(am__objects_12) libipopt_la_OBJECTS = $(am_libipopt_la_OBJECTS) AM_V_lt = $(am__v_lt_@AM_V@) am__v_lt_ = $(am__v_lt_@AM_DEFAULT_V@) @@ -367,7 +368,6 @@ am__depfiles_remade = Algorithm/$(DEPDIR)/IpAdaptiveMuUpdate.Plo \ Algorithm/Inexact/$(DEPDIR)/IpIterativeSolverTerminationTester.Plo \ Algorithm/LinearSolvers/$(DEPDIR)/IpIterativeWsmpSolverInterface.Plo \ Algorithm/LinearSolvers/$(DEPDIR)/IpKLUSolverInterface.Plo \ - Algorithm/LinearSolvers/$(DEPDIR)/IpReSolveSolverInterface.Plo \ Algorithm/LinearSolvers/$(DEPDIR)/IpLinearSolvers.Plo \ Algorithm/LinearSolvers/$(DEPDIR)/IpLinearSolversRegOp.Plo \ Algorithm/LinearSolvers/$(DEPDIR)/IpMa27TSolverInterface.Plo \ @@ -380,6 +380,7 @@ am__depfiles_remade = Algorithm/$(DEPDIR)/IpAdaptiveMuUpdate.Plo \ Algorithm/LinearSolvers/$(DEPDIR)/IpMumpsSolverInterface.Plo \ Algorithm/LinearSolvers/$(DEPDIR)/IpPardisoMKLSolverInterface.Plo \ Algorithm/LinearSolvers/$(DEPDIR)/IpPardisoSolverInterface.Plo \ + Algorithm/LinearSolvers/$(DEPDIR)/IpReSolveSolverInterface.Plo \ Algorithm/LinearSolvers/$(DEPDIR)/IpSlackBasedTSymScalingMethod.Plo \ Algorithm/LinearSolvers/$(DEPDIR)/IpSpralSolverInterface.Plo \ Algorithm/LinearSolvers/$(DEPDIR)/IpTSymDependencyDetector.Plo \ @@ -864,7 +865,8 @@ libipopt_la_SOURCES = Common/IpDebug.cpp Common/IpJournalist.cpp \ Interfaces/IpTNLPReducer.cpp $(am__append_1) $(am__append_2) \ $(am__append_3) $(am__append_4) $(am__append_5) \ $(am__append_6) $(am__append_7) $(am__append_8) \ - $(am__append_9) $(am__append_10) $(am__append_11) $(am__append_12) + $(am__append_9) $(am__append_10) $(am__append_11) \ + $(am__append_12) AM_CPPFLAGS = \ -I$(srcdir)/Common \ -I$(srcdir)/LinAlg \ @@ -1396,7 +1398,6 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/Inexact/$(DEPDIR)/IpIterativeSolverTerminationTester.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpIterativeWsmpSolverInterface.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpKLUSolverInterface.Plo@am__quote@ # am--include-marker -@AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpReSolveSolverInterface.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpLinearSolvers.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpLinearSolversRegOp.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpMa27TSolverInterface.Plo@am__quote@ # am--include-marker @@ -1409,6 +1410,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpMumpsSolverInterface.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpPardisoMKLSolverInterface.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpPardisoSolverInterface.Plo@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpReSolveSolverInterface.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpSlackBasedTSymScalingMethod.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpSpralSolverInterface.Plo@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@Algorithm/LinearSolvers/$(DEPDIR)/IpTSymDependencyDetector.Plo@am__quote@ # am--include-marker @@ -1813,7 +1815,6 @@ distclean: distclean-recursive -rm -f Algorithm/Inexact/$(DEPDIR)/IpIterativeSolverTerminationTester.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpIterativeWsmpSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpKLUSolverInterface.Plo - -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpReSolveSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpLinearSolvers.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpLinearSolversRegOp.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpMa27TSolverInterface.Plo @@ -1826,6 +1827,7 @@ distclean: distclean-recursive -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpMumpsSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpPardisoMKLSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpPardisoSolverInterface.Plo + -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpReSolveSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpSlackBasedTSymScalingMethod.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpSpralSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpTSymDependencyDetector.Plo @@ -1989,7 +1991,6 @@ maintainer-clean: maintainer-clean-recursive -rm -f Algorithm/Inexact/$(DEPDIR)/IpIterativeSolverTerminationTester.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpIterativeWsmpSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpKLUSolverInterface.Plo - -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpReSolveSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpLinearSolvers.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpLinearSolversRegOp.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpMa27TSolverInterface.Plo @@ -2002,6 +2003,7 @@ maintainer-clean: maintainer-clean-recursive -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpMumpsSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpPardisoMKLSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpPardisoSolverInterface.Plo + -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpReSolveSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpSlackBasedTSymScalingMethod.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpSpralSolverInterface.Plo -rm -f Algorithm/LinearSolvers/$(DEPDIR)/IpTSymDependencyDetector.Plo diff --git a/test/Makefile.am b/test/Makefile.am index 58d59f48..cd1342dd 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -10,6 +10,10 @@ noinst_PROGRAMS = hs071_cpp hs071_c emptynlp getcurr +if COIN_HAS_RESOLVE +noinst_PROGRAMS += resolve_multirhs +endif + if COIN_HAS_F77 noinst_PROGRAMS += hs071_f endif @@ -30,6 +34,12 @@ emptynlp_LDADD = ../src/libipopt.la nodist_getcurr_SOURCES = getcurr.cpp getcurr_LDADD = ../src/libipopt.la +if COIN_HAS_RESOLVE +resolve_multirhs_SOURCES = resolve_multirhs.cpp +resolve_multirhs_CPPFLAGS = $(AM_CPPFLAGS) $(IPOPTLIB_CFLAGS_NOPC) +resolve_multirhs_LDADD = ../src/libipopt.la +endif + if !IPOPT_SINGLE nodist_hs071_f_SOURCES = hs071_f.f else @@ -65,6 +75,7 @@ AM_CPPFLAGS = \ -I$(srcdir)/../src/LinAlg \ -I$(srcdir)/../src/LinAlg/TMatrices \ -I$(srcdir)/../src/Algorithm \ + -I$(srcdir)/../src/Algorithm/LinearSolvers \ -I$(srcdir)/../src/Interfaces \ -I$(srcdir)/../contrib/sIPOPT/src \ -I$(srcdir)/../contrib/sIPOPT/examples/parametric_cpp \ @@ -83,6 +94,10 @@ endif chmod u+x ./run_unitTests ./run_unitTests +if COIN_HAS_RESOLVE + ./resolve_multirhs +endif + unitTest: test .PHONY: test unitTest diff --git a/test/Makefile.in b/test/Makefile.in index 313927e5..ee54e2d7 100644 --- a/test/Makefile.in +++ b/test/Makefile.in @@ -100,9 +100,10 @@ build_triplet = @build@ host_triplet = @host@ noinst_PROGRAMS = hs071_cpp$(EXEEXT) hs071_c$(EXEEXT) \ emptynlp$(EXEEXT) getcurr$(EXEEXT) $(am__EXEEXT_1) \ - $(am__EXEEXT_2) -@COIN_HAS_F77_TRUE@am__append_1 = hs071_f -@BUILD_SIPOPT_TRUE@am__append_2 = parametric_cpp redhess_cpp + $(am__EXEEXT_2) $(am__EXEEXT_3) +@COIN_HAS_RESOLVE_TRUE@am__append_1 = resolve_multirhs +@COIN_HAS_F77_TRUE@am__append_2 = hs071_f +@BUILD_SIPOPT_TRUE@am__append_3 = parametric_cpp redhess_cpp subdir = test ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 am__aclocal_m4_deps = $(top_srcdir)/configure.ac @@ -116,8 +117,9 @@ CONFIG_CLEAN_FILES = run_unitTests hs071_main.cpp hs071_nlp.cpp \ hs071_nlp.hpp hs071_c.c parametric_driver.cpp \ parametricTNLP.cpp MySensTNLP.cpp redhess_cpp.cpp CONFIG_CLEAN_VPATH_FILES = -@COIN_HAS_F77_TRUE@am__EXEEXT_1 = hs071_f$(EXEEXT) -@BUILD_SIPOPT_TRUE@am__EXEEXT_2 = parametric_cpp$(EXEEXT) \ +@COIN_HAS_RESOLVE_TRUE@am__EXEEXT_1 = resolve_multirhs$(EXEEXT) +@COIN_HAS_F77_TRUE@am__EXEEXT_2 = hs071_f$(EXEEXT) +@BUILD_SIPOPT_TRUE@am__EXEEXT_3 = parametric_cpp$(EXEEXT) \ @BUILD_SIPOPT_TRUE@ redhess_cpp$(EXEEXT) PROGRAMS = $(noinst_PROGRAMS) nodist_emptynlp_OBJECTS = emptynlp.$(OBJEXT) @@ -149,6 +151,10 @@ nodist_redhess_cpp_OBJECTS = MySensTNLP.$(OBJEXT) \ redhess_cpp.$(OBJEXT) redhess_cpp_OBJECTS = $(nodist_redhess_cpp_OBJECTS) redhess_cpp_DEPENDENCIES = ../contrib/sIPOPT/src/libsipopt.la +@COIN_HAS_RESOLVE_TRUE@am_resolve_multirhs_OBJECTS = resolve_multirhs-resolve_multirhs.$(OBJEXT) +resolve_multirhs_OBJECTS = $(am_resolve_multirhs_OBJECTS) +@COIN_HAS_RESOLVE_TRUE@resolve_multirhs_DEPENDENCIES = \ +@COIN_HAS_RESOLVE_TRUE@ ../src/libipopt.la AM_V_P = $(am__v_P_@AM_V@) am__v_P_ = $(am__v_P_@AM_DEFAULT_V@) am__v_P_0 = false @@ -168,7 +174,8 @@ am__depfiles_remade = ./$(DEPDIR)/MySensTNLP.Po \ ./$(DEPDIR)/emptynlp.Po ./$(DEPDIR)/getcurr.Po \ ./$(DEPDIR)/hs071_c.Po ./$(DEPDIR)/hs071_main.Po \ ./$(DEPDIR)/hs071_nlp.Po ./$(DEPDIR)/parametricTNLP.Po \ - ./$(DEPDIR)/parametric_driver.Po ./$(DEPDIR)/redhess_cpp.Po + ./$(DEPDIR)/parametric_driver.Po ./$(DEPDIR)/redhess_cpp.Po \ + ./$(DEPDIR)/resolve_multirhs-resolve_multirhs.Po am__mv = mv -f COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \ $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) @@ -224,7 +231,7 @@ am__v_F77LD_1 = SOURCES = $(nodist_emptynlp_SOURCES) $(nodist_getcurr_SOURCES) \ $(nodist_hs071_c_SOURCES) $(nodist_hs071_cpp_SOURCES) \ $(nodist_hs071_f_SOURCES) $(nodist_parametric_cpp_SOURCES) \ - $(nodist_redhess_cpp_SOURCES) + $(nodist_redhess_cpp_SOURCES) $(resolve_multirhs_SOURCES) am__can_run_installinfo = \ case $$AM_UPDATE_INFO_DIR in \ n|no|NO) false;; \ @@ -430,6 +437,9 @@ nodist_emptynlp_SOURCES = emptynlp.cpp emptynlp_LDADD = ../src/libipopt.la nodist_getcurr_SOURCES = getcurr.cpp getcurr_LDADD = ../src/libipopt.la +@COIN_HAS_RESOLVE_TRUE@resolve_multirhs_SOURCES = resolve_multirhs.cpp +@COIN_HAS_RESOLVE_TRUE@resolve_multirhs_CPPFLAGS = $(AM_CPPFLAGS) $(IPOPTLIB_CFLAGS_NOPC) +@COIN_HAS_RESOLVE_TRUE@resolve_multirhs_LDADD = ../src/libipopt.la @IPOPT_SINGLE_FALSE@nodist_hs071_f_SOURCES = hs071_f.f @IPOPT_SINGLE_TRUE@nodist_hs071_f_SOURCES = hs071_fs.f hs071_f_LDADD = ../src/libipopt.la $(CXXLIBS) @@ -447,6 +457,7 @@ AM_CPPFLAGS = \ -I$(srcdir)/../src/LinAlg \ -I$(srcdir)/../src/LinAlg/TMatrices \ -I$(srcdir)/../src/Algorithm \ + -I$(srcdir)/../src/Algorithm/LinearSolvers \ -I$(srcdir)/../src/Interfaces \ -I$(srcdir)/../contrib/sIPOPT/src \ -I$(srcdir)/../contrib/sIPOPT/examples/parametric_cpp \ @@ -526,6 +537,10 @@ redhess_cpp$(EXEEXT): $(redhess_cpp_OBJECTS) $(redhess_cpp_DEPENDENCIES) $(EXTRA @rm -f redhess_cpp$(EXEEXT) $(AM_V_CXXLD)$(CXXLINK) $(redhess_cpp_OBJECTS) $(redhess_cpp_LDADD) $(LIBS) +resolve_multirhs$(EXEEXT): $(resolve_multirhs_OBJECTS) $(resolve_multirhs_DEPENDENCIES) $(EXTRA_resolve_multirhs_DEPENDENCIES) + @rm -f resolve_multirhs$(EXEEXT) + $(AM_V_CXXLD)$(CXXLINK) $(resolve_multirhs_OBJECTS) $(resolve_multirhs_LDADD) $(LIBS) + mostlyclean-compile: -rm -f *.$(OBJEXT) @@ -541,6 +556,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/parametricTNLP.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/parametric_driver.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/redhess_cpp.Po@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/resolve_multirhs-resolve_multirhs.Po@am__quote@ # am--include-marker $(am__depfiles_remade): @$(MKDIR_P) $(@D) @@ -596,6 +612,20 @@ am--depfiles: $(am__depfiles_remade) @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LTCXXCOMPILE) -c -o $@ $< +resolve_multirhs-resolve_multirhs.o: resolve_multirhs.cpp +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(resolve_multirhs_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT resolve_multirhs-resolve_multirhs.o -MD -MP -MF $(DEPDIR)/resolve_multirhs-resolve_multirhs.Tpo -c -o resolve_multirhs-resolve_multirhs.o `test -f 'resolve_multirhs.cpp' || echo '$(srcdir)/'`resolve_multirhs.cpp +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/resolve_multirhs-resolve_multirhs.Tpo $(DEPDIR)/resolve_multirhs-resolve_multirhs.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='resolve_multirhs.cpp' object='resolve_multirhs-resolve_multirhs.o' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(resolve_multirhs_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o resolve_multirhs-resolve_multirhs.o `test -f 'resolve_multirhs.cpp' || echo '$(srcdir)/'`resolve_multirhs.cpp + +resolve_multirhs-resolve_multirhs.obj: resolve_multirhs.cpp +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(resolve_multirhs_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT resolve_multirhs-resolve_multirhs.obj -MD -MP -MF $(DEPDIR)/resolve_multirhs-resolve_multirhs.Tpo -c -o resolve_multirhs-resolve_multirhs.obj `if test -f 'resolve_multirhs.cpp'; then $(CYGPATH_W) 'resolve_multirhs.cpp'; else $(CYGPATH_W) '$(srcdir)/resolve_multirhs.cpp'; fi` +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/resolve_multirhs-resolve_multirhs.Tpo $(DEPDIR)/resolve_multirhs-resolve_multirhs.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='resolve_multirhs.cpp' object='resolve_multirhs-resolve_multirhs.obj' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(resolve_multirhs_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o resolve_multirhs-resolve_multirhs.obj `if test -f 'resolve_multirhs.cpp'; then $(CYGPATH_W) 'resolve_multirhs.cpp'; else $(CYGPATH_W) '$(srcdir)/resolve_multirhs.cpp'; fi` + .f.o: $(AM_V_F77)$(F77COMPILE) -c -o $@ $< @@ -713,6 +743,7 @@ distclean: distclean-am -rm -f ./$(DEPDIR)/parametricTNLP.Po -rm -f ./$(DEPDIR)/parametric_driver.Po -rm -f ./$(DEPDIR)/redhess_cpp.Po + -rm -f ./$(DEPDIR)/resolve_multirhs-resolve_multirhs.Po -rm -f Makefile distclean-am: clean-am distclean-compile distclean-generic \ distclean-tags @@ -767,6 +798,7 @@ maintainer-clean: maintainer-clean-am -rm -f ./$(DEPDIR)/parametricTNLP.Po -rm -f ./$(DEPDIR)/parametric_driver.Po -rm -f ./$(DEPDIR)/redhess_cpp.Po + -rm -f ./$(DEPDIR)/resolve_multirhs-resolve_multirhs.Po -rm -f Makefile maintainer-clean-am: distclean-am maintainer-clean-generic @@ -817,6 +849,8 @@ hs071_f.f hs071_fs.f: chmod u+x ./run_unitTests ./run_unitTests +@COIN_HAS_RESOLVE_TRUE@ ./resolve_multirhs + unitTest: test .PHONY: test unitTest diff --git a/test/resolve_multirhs.cpp b/test/resolve_multirhs.cpp new file mode 100644 index 00000000..5d6e324b --- /dev/null +++ b/test/resolve_multirhs.cpp @@ -0,0 +1,285 @@ +#include "IpOptionsList.hpp" +#include "IpReSolveSolverInterface.hpp" + +#include +#include +#include +#include + +using namespace Ipopt; + +namespace +{ + +bool check_solution( + const std::string& name, + const Number* actual, + const Number* expected, + Index n, + Number tolerance = 1e-8) +{ + for( Index i = 0; i < n; ++i ) + { + if( std::abs(actual[i] - expected[i]) > tolerance ) + { + std::cerr << "FAIL: " << name + << " at index " << i + << ": expected " << expected[i] + << ", got " << actual[i] << '\n'; + return false; + } + } + + return true; +} + +bool solve_succeeded( + const std::string& name, + ESymSolverStatus status) +{ + if( status != SYMSOLVER_SUCCESS ) + { + std::cerr << "FAIL: " << name + << " returned solver status " << status << '\n'; + return false; + } + + return true; +} + +bool configure_solver( + ReSolveSolverInterface& solver, + const std::string& method, + const char* label) +{ + OptionsList options; + + if( !options.SetNumericValue("resolve_tol", 1e-3) + || !options.SetIntegerValue("resolve_ordering", 0) + || !options.SetIntegerValue("resolve_n_skip_refactoring", 1) + || !options.SetStringValue("resolve_halt_if_singular", "no") + || !options.SetStringValue("resolve_method", method) + || !options.SetNumericValue("resolve_rcond_val", 1e-128) + || !options.SetStringValue("resolve_use_rcond", "no") ) + { + std::cerr << "FAIL: " << label << " options could not be set\n"; + return false; + } + + if( !solver.InitializeImpl(options, "") ) + { + std::cerr << "FAIL: " << label << " InitializeImpl failed\n"; + return false; + } + + return true; +} + +bool run_method( + const std::string& method, + const char* label) +{ + constexpr Index dim = 3; + constexpr Index nnz = 7; + + // Full CSR, zero-based indexing: + // + // [ * * 0 ] + // [ * * * ] + // [ 0 * * ] + const Index ia[dim + 1] = {0, 2, 5, 7}; + const Index ja[nnz] = {0, 1, 0, 1, 2, 1, 2}; + + const Number A1[nnz] = { + 4.0, 1.0, + 1.0, 3.0, 1.0, + 1.0, 2.0 + }; + + const Number A2[nnz] = { + 5.0, 1.0, + 1.0, 4.0, 1.0, + 1.0, 3.0 + }; + + const Number x1[dim] = {1.0, 2.0, 3.0}; + const Number x2[dim] = {-1.0, 0.0, 2.0}; + + const Number b1_A1[dim] = {6.0, 10.0, 8.0}; + const Number b2_A1[dim] = {-4.0, 1.0, 4.0}; + + const Number b1_A2[dim] = {7.0, 12.0, 11.0}; + const Number b2_A2[dim] = {-5.0, 1.0, 6.0}; + + ReSolveSolverInterface solver; + + if( !configure_solver(solver, method, label) ) + { + return false; + } + + if( solver.MatrixFormat() + != SparseSymLinearSolverInterface::CSR_Full_Format_0_Offset ) + { + std::cerr << "FAIL: " << label + << " returned an unexpected matrix format\n"; + return false; + } + + if( !solve_succeeded( + std::string(label) + " InitializeStructure", + solver.InitializeStructure(dim, nnz, ia, ja)) ) + { + return false; + } + + Number* values = solver.GetValuesArrayPtr(); + if( values == nullptr ) + { + std::cerr << "FAIL: " << label + << " GetValuesArrayPtr returned null\n"; + return false; + } + + // nrhs = 1 baseline. + std::copy(A1, A1 + nnz, values); + + Number rhs_single[dim]; + std::copy(b1_A1, b1_A1 + dim, rhs_single); + + if( !solve_succeeded( + std::string(label) + " nrhs=1", + solver.MultiSolve(true, ia, ja, 1, rhs_single, false, 0)) ) + { + return false; + } + + if( !check_solution( + std::string(label) + " nrhs=1", + rhs_single, + x1, + dim) ) + { + return false; + } + + // nrhs = 2 using the existing factorization. + Number rhs_multiple[2 * dim] = { + b1_A1[0], b1_A1[1], b1_A1[2], + b2_A1[0], b2_A1[1], b2_A1[2] + }; + + if( !solve_succeeded( + std::string(label) + " nrhs=2", + solver.MultiSolve(false, ia, ja, 2, rhs_multiple, false, 0)) ) + { + return false; + } + + if( !check_solution( + std::string(label) + " nrhs=2 rhs 1", + rhs_multiple, + x1, + dim) + || !check_solution( + std::string(label) + " nrhs=2 rhs 2", + rhs_multiple + dim, + x2, + dim) ) + { + return false; + } + + // Change numerical values while preserving the sparsity structure. + std::copy(A2, A2 + nnz, values); + + Number rhs_updated[2 * dim] = { + b1_A2[0], b1_A2[1], b1_A2[2], + b2_A2[0], b2_A2[1], b2_A2[2] + }; + + if( !solve_succeeded( + std::string(label) + " matrix update", + solver.MultiSolve(true, ia, ja, 2, rhs_updated, false, 0)) ) + { + return false; + } + + if( !check_solution( + std::string(label) + " matrix update rhs 1", + rhs_updated, + x1, + dim) + || !check_solution( + std::string(label) + " matrix update rhs 2", + rhs_updated + dim, + x2, + dim) ) + { + return false; + } + + // Reuse the updated factorization. + Number rhs_repeat[dim]; + std::copy(b1_A2, b1_A2 + dim, rhs_repeat); + + if( !solve_succeeded( + std::string(label) + " repeated solve", + solver.MultiSolve(false, ia, ja, 1, rhs_repeat, false, 0)) ) + { + return false; + } + + if( !check_solution( + std::string(label) + " repeated solve", + rhs_repeat, + x1, + dim) ) + { + return false; + } + + if( solver.IncreaseQuality() ) + { + std::cerr << "FAIL: " << label + << " IncreaseQuality should return false\n"; + return false; + } + + return true; +} + +} // namespace + +int main() +{ + std::cout << "Testing ReSolve Solver Interface...\n"; + + if( !run_method(resolve_klu, "KLU") ) + { + return 1; + } + +#ifdef RESOLVE_USE_GPU +# ifdef RESOLVE_USE_CUDA + if( !run_method(resolve_glu, "GLU") ) + { + return 1; + } +# endif + + if( !run_method(resolve_rf, "RF") ) + { + return 1; + } + + if( !run_method(resolve_rf_fgmres, "RF-FGMRES") ) + { + return 1; + } +#endif + + std::cout << " Test passed!\n"; + + return 0; +}