From bfb1efe9aa916620669d1d583b4c497b41cc2389 Mon Sep 17 00:00:00 2001 From: Maksudul Alam Date: Thu, 15 May 2025 17:46:11 +0000 Subject: [PATCH 01/19] Debug Fix + GPU Code Fix --- .../IpReSolveSolverInterface.cpp | 40 ++- src/Algorithm/LinearSolvers/test-csr.h | 257 ++++++++++++++++++ 2 files changed, 288 insertions(+), 9 deletions(-) create mode 100644 src/Algorithm/LinearSolvers/test-csr.h diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index 81f832a8..f0d4b92a 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -203,12 +203,14 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no if (method_ == resolve_klu) { + printf("Resolve::KLU Setup\n"); workspace_CPU_ = new ReSolve::LinAlgWorkspaceCpu(); matrix_handler_ = new ReSolve::MatrixHandler(workspace_CPU_); vector_handler_ = new ReSolve::VectorHandler(workspace_CPU_); } #if RESOLVE_WITH_GPU + printf("Resolve::GPU Setup\n"); workspace_GPU_ = new workspace_type(); workspace_GPU_->initializeHandles(); @@ -218,17 +220,20 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - resolve_Rf_ = new rf_solver(workspace_GPU_); + printf("Resolve::RF Setup\n"); + resolve_Rf_ = new rf_solver(); } # if RESOLVE_WITH_CUDA else if (method_ == resolve_glu) { + printf("Resolve::GLU Setup\n"); resolve_GLU_ = new ReSolve::LinSolverDirectCuSolverGLU(workspace_GPU_); } # endif if (method_ == resolve_rf_fgmres) { + printf("Resolve::FGMRES Setup\n"); GS_ = new ReSolve::GramSchmidt(vector_handler_, ReSolve::GramSchmidt::CGS2); resolve_FGMRES_ = new ReSolve::LinSolverIterativeFGMRES(matrix_handler_, vector_handler_, GS_); } @@ -283,7 +288,14 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index 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); + + A_->setDataPointers(const_cast(ia), const_cast(ja), const_cast(A_->getValues(ReSolve::memory::HOST)), ReSolve::memory::HOST); + +#if RESOLVE_WITH_GPU + if (method_ == resolve_rf || method_ == resolve_rf_fgmres || method_ == resolve_glu) { + A_->syncData(ReSolve::memory::DEVICE); + } +#endif // FACTORIZE @@ -294,10 +306,8 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // re_factorize_ = true; } - if (factorize_ && (new_matrix || re_factorize_)) - { - // printf("Iteration: %d: Performing KLU Factorization\n", n_iteration_); - + if( n_iteration_ == 0){ + printf("First Iteration: %d: Performing KLU Factorization\n", n_iteration_); // Symbolic Factorization if (HaveIpData()) { @@ -314,6 +324,11 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index { IpData().TimingStats().LinearSystemSymbolicFactorization().End(); } + } + + if (factorize_ && (new_matrix || re_factorize_)) + { + // printf("Iteration: %d: Performing KLU Factorization\n", n_iteration_); // perform the factorization if (HaveIpData()) @@ -322,7 +337,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } // First Factorization is always done by KLU - std::cout << "%" << n_iteration_ << "%" << "FULL FACTORIZATIOM" << std::endl; + std::cout << "%" << n_iteration_ << "%" << "KLU FULL FACTORIZATION" << std::endl; status = resolve_KLU_->factorize(); full_factor_done = true; @@ -397,6 +412,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // Actual Refactorize if (method_ == resolve_glu) { + std::cout << "%" << n_iteration_ << "%" << "GLU->refactorize()" << std::endl; status = resolve_GLU_->refactorize(); if (status != 0) { @@ -405,6 +421,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } else if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { + std::cout << "%" << n_iteration_ << "%" << "RF->refactorize()" << std::endl; status_refactor = resolve_Rf_->refactorize(); if (status != 0) { @@ -414,6 +431,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index # else if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { + std::cout << "%" << n_iteration_ << "%" << "RF->refactorize()" << std::endl; int status = resolve_Rf_->refactorize(); if (status != 0) { @@ -424,7 +442,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index if (method_ == resolve_klu) { - std::cout << "%" << n_iteration_ << "%" << "RE-FACTORIZATIOM" << std::endl; + std::cout << "%" << n_iteration_ << "%" << "KLU->refactorize()" << std::endl; status = resolve_KLU_->refactorize(); if (status != 0) { @@ -476,6 +494,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // 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_); + printf("Solving using KLU!\n"); if (status != 0) { std::cout << "KLU solve status: " << status << std::endl; @@ -489,10 +508,13 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index 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()); + L_csc->syncData(ReSolve::memory::DEVICE); + U_csc->syncData(ReSolve::memory::DEVICE); matrix_handler_->csc2csr(L_csc, L, ReSolve::memory::DEVICE); matrix_handler_->csc2csr(U_csc, U, ReSolve::memory::DEVICE); if (L == nullptr) @@ -615,7 +637,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index if (use_rcond_) { Number rcond_val = resolve_KLU_->getMatrixConditionNumber(); - printf("RCond: %12.8e\n", rcond_val); + //printf("RCond: %12.8e\n", rcond_val); if (rcond_val < rcond_val_) { if (full_factor_done) diff --git a/src/Algorithm/LinearSolvers/test-csr.h b/src/Algorithm/LinearSolvers/test-csr.h new file mode 100644 index 00000000..0ddeacad --- /dev/null +++ b/src/Algorithm/LinearSolvers/test-csr.h @@ -0,0 +1,257 @@ +#ifndef __TEST_CSR_H__ +#define __TEST_CSR_H__ + +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +namespace debug +{ + +struct CSR +{ + vector val; + vector col_ind; + vector row_ptr; +}; + +class Debug +{ +public: + Debug() + { + } + + template static T* load_b_data(std::string prefix, int iter, int& n) + { + std::fstream file; + + std::stringstream ss; + ss << prefix << "-" << iter << ".txt"; + std::string filename = ss.str(); + + file.open(filename, std::ios::in | std::fstream::binary); + file.read((char*)&n, sizeof(n)); + T* d = new T[n]; + + T p; + for (int i = 0; i < n; i++) + { + file.read((char*)&p, sizeof(T)); + d[i] = p; + } + + file.close(); + return d; + } + + template static void save_b_data(std::string prefix, int iter, int n, const T* data) + { + std::fstream file; + std::stringstream ss; + ss << prefix << "-" << iter << ".bin"; + std::string filename = ss.str(); + + file.open(filename, std::ios::out | std::fstream::binary); + + file.write((char*)&n, sizeof(n)); + // std::cout.precision(std::numeric_limits::max_digits10 - 1); + file.write((char*)(data), n * sizeof(T)); + file.close(); + } + + template static T* load_data(std::string prefix, int iter, int& n) + { + std::fstream file; + + std::stringstream ss; + ss << prefix << "-" << iter << ".txt"; + std::string filename = ss.str(); + + file.open(filename, std::ios::in); + file >> n; + T* d = new T[n]; + + T p; + for (int i = 0; i < n; i++) + { + file >> p; + d[i] = p; + } + + file.close(); + return d; + } + + template static void save_data(std::string prefix, int iter, int n, const T* data) + { + std::fstream file; + std::stringstream ss; + ss << prefix << "-" << iter << ".txt"; + std::string filename = ss.str(); + + file.open(filename, std::ios::out); + file << n << std::endl; + + // std::cout.precision(std::numeric_limits::max_digits10 - 1); + + for (int i = 0; i < n; i++) + { + // file << std::scientific << data[i] << std::endl; + file << data[i] << std::endl; + } + // file.write((char *)(data), n * sizeof(T)); + file.close(); + } + + static void save_coo_as_mm(std::string prefix, int iter, int dim, int nnz, const int* row_ptr, const int* col_ind, double* val) + { + printf("%s\t%d\n", __FILE__, __LINE__); + int cont = 0; + + std::stringstream ss; + ss << prefix << "-matrix-" << iter << ".mtx"; + std::string filename = ss.str(); + std::cout << "Filename: " << filename << std::endl; + + std::ofstream mat_file(filename); + + mat_file << "%%MatrixMarket matrix coordinate real general" << std::endl; + mat_file << "% ID: " << iter << std::endl; + mat_file << dim << "\t" << dim << "\t" << nnz << std::endl; + + int nnz_c = 0; + + for (int i = 0; i < nnz; i++) + { + mat_file << row_ptr[i] << "\t" << col_ind[i] << "\t" << val[i] << std::endl; + } + mat_file.close(); + printf("%d, %d==%d\n", dim, nnz_c, nnz); + } + + static void save_csc_as_csv(std::string prefix, int iter, int dim, int nnz, const int* row_ptr, const int* col_ind, double* val) + { + printf("%s\t%d\n", __FILE__, __LINE__); + int cont = 0; + + std::stringstream ss; + ss << prefix << "-matrix-" << iter << ".csv"; + std::string filename = ss.str(); + std::cout << "Filename: " << filename << std::endl; + + std::ofstream mat_file(filename); + + int nnz_c = 0; + + for (int i = 1; i <= dim; i++) + { + int row_start = row_ptr[i - 1]; + int row_end = row_ptr[i]; + + nnz_c += row_end - row_start; + // printf("%d, %d, %d\n", row_start, row_end, row_end - row_start); + + for (int jj = 0; jj < dim; jj++) + { + bool found = false; + for (int j = row_start; j < row_end; j++) + { + if (jj == col_ind[j]) + { + found = true; + break; + } + } + + if (found) + { + mat_file << val[cont] << ","; + cont++; + } + else + { + mat_file << "0,"; + } + } + mat_file << std::endl; + } + mat_file.close(); + printf("%d, %d==%d\n", dim, nnz_c, nnz); + } + + static void save_csc_as_mm(std::string prefix, int iter, int dim, int nnz, const int* row_ptr, const int* col_ind, double* val) + { + printf("%s\t%d\n", __FILE__, __LINE__); + int cont = 0; + + std::stringstream ss; + ss << prefix << "-matrix-" << iter << ".mtx"; + std::string filename = ss.str(); + std::cout << "Filename: " << filename << std::endl; + + std::ofstream mat_file(filename); + + mat_file << "%%MatrixMarket matrix coordinate real general" << std::endl; + mat_file << "% ID: " << iter << std::endl; + mat_file << dim << "\t" << dim << "\t" << nnz << std::endl; + + int nnz_c = 0; + + for (int i = 1; i <= dim; i++) + { + int row_start = row_ptr[i - 1]; + int row_end = row_ptr[i]; + + nnz_c += row_end - row_start; + // printf("%d, %d, %d\n", row_start, row_end, row_end - row_start); + + for (int j = row_start; j < row_end; j++) + { + // printf("%d ", col_ind[j]); + mat_file << i << "\t" << col_ind[j] + 1 << "\t" << val[cont] << std::endl; + cont++; + } + } + mat_file.close(); + printf("%d, %d==%d\n", dim, nnz_c, nnz); + } + + static void save_vec_as_csv(std::string prefix, int iter, int dim, int nnz, const int* row_ptr, const int* col_ind, double* val) + { + } + + static void save_vec_as_mm(std::string prefix, int iter, int dim, double* val) + { + printf("%s\t%d\n", __FILE__, __LINE__); + int cont = 0; + + std::stringstream ss; + ss << prefix << "-rhs-" << iter << ".mtx"; + std::string filename = ss.str(); + std::cout << "Filename: " << filename << std::endl; + + std::ofstream vec_file(filename); + + vec_file << "%%MatrixMarket matrix array real general" << std::endl; + vec_file << "% ID: " << iter << std::endl; + vec_file << dim << "\t" << 1 << std::endl; + + int nnz_c = 0; + + for (int i = 0; i < dim; i++) + { + vec_file << val[i] << std::endl; + } + vec_file.close(); + } +}; + +} // namespace debug +#endif From 6e243317ae29ec8247d874ca97d0fa08cd4dab1d Mon Sep 17 00:00:00 2001 From: Maksudul Alam Date: Thu, 31 Jul 2025 16:14:20 +0000 Subject: [PATCH 02/19] Cleanup --- .../IpReSolveSolverInterface.cpp | 28 +- .../IpReSolveSolverInterface.hpp | 4 - src/Algorithm/LinearSolvers/NVMLHelper.hpp | 70 ----- src/Algorithm/LinearSolvers/test-csr.h | 257 ------------------ 4 files changed, 9 insertions(+), 350 deletions(-) delete mode 100644 src/Algorithm/LinearSolvers/NVMLHelper.hpp delete mode 100644 src/Algorithm/LinearSolvers/test-csr.h diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index f0d4b92a..0819dfd9 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -28,7 +28,6 @@ ReSolveSolverInterface::ReSolveSolverInterface() : val_(NULL) printf("Resolve with GPU.\n"); # if RESOLVE_WITH_CUDA printf("Resolve with CUDA.\n"); - NVMLHelper::getAvailableGPUMemory(); # else printf("Resolve with HIP.\n"); # endif @@ -39,10 +38,6 @@ ReSolveSolverInterface::ReSolveSolverInterface() : val_(NULL) ReSolveSolverInterface::~ReSolveSolverInterface() { - -// printf("Begin of Destructor\n"); -// NVMLHelper::getAvailableGPUMemory(); - DBG_START_METH("ReSolveSolverInterface::~ReSolveSolverInterface()", dbg_verbosity); delete[] val_; @@ -74,8 +69,6 @@ ReSolveSolverInterface::~ReSolveSolverInterface() delete vec_rhs_; delete vec_x_; delete A_; - -// printf("End of Destructor\n"); } void ReSolveSolverInterface::RegisterOptions(SmartPtr roptions) @@ -263,7 +256,6 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no vec_x_ = new ReSolve::vector::Vector(A_->getNumRows()); vec_x_->allocate(ReSolve::memory::HOST); // for KLU - // vec_x_->allocate(ReSolve::memory::DEVICE); } factorize_ = true; @@ -272,9 +264,6 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no initialized_ = true; pivtol_changed_ = false; -// printf("After Initialize\n"); -// NVMLHelper::getAvailableGPUMemory(); - return retval; } @@ -337,7 +326,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } // First Factorization is always done by KLU - std::cout << "%" << n_iteration_ << "%" << "KLU FULL FACTORIZATION" << std::endl; + // std::cout << "%" << n_iteration_ << "%" << "KLU FULL FACTORIZATION" << std::endl; status = resolve_KLU_->factorize(); full_factor_done = true; @@ -412,7 +401,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // Actual Refactorize if (method_ == resolve_glu) { - std::cout << "%" << n_iteration_ << "%" << "GLU->refactorize()" << std::endl; + //std::cout << "%" << n_iteration_ << "%" << "GLU->refactorize()" << std::endl; status = resolve_GLU_->refactorize(); if (status != 0) { @@ -421,7 +410,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } else if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - std::cout << "%" << n_iteration_ << "%" << "RF->refactorize()" << std::endl; + //std::cout << "%" << n_iteration_ << "%" << "RF->refactorize()" << std::endl; status_refactor = resolve_Rf_->refactorize(); if (status != 0) { @@ -431,7 +420,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index # else if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - std::cout << "%" << n_iteration_ << "%" << "RF->refactorize()" << std::endl; + //std::cout << "%" << n_iteration_ << "%" << "RF->refactorize()" << std::endl; int status = resolve_Rf_->refactorize(); if (status != 0) { @@ -442,7 +431,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index if (method_ == resolve_klu) { - std::cout << "%" << n_iteration_ << "%" << "KLU->refactorize()" << std::endl; + //std::cout << "%" << n_iteration_ << "%" << "KLU->refactorize()" << std::endl; status = resolve_KLU_->refactorize(); if (status != 0) { @@ -470,7 +459,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index if (use_rcond_) { Number rcond_val = resolve_KLU_->getMatrixConditionNumber(); - printf("RCond: %12.8e\n", rcond_val); + //printf("RCond: %12.8e\n", rcond_val); if (rcond_val < rcond_val_) { if (full_factor_done) @@ -494,7 +483,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // 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_); - printf("Solving using KLU!\n"); + //printf("Solving using KLU!\n"); if (status != 0) { std::cout << "KLU solve status: " << status << std::endl; @@ -507,7 +496,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index # if RESOLVE_WITH_CUDA if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - printf("Iteration: %d: Setting up %s\n", n_iteration_, method_.c_str()); + printf("CUDA: 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(); @@ -538,6 +527,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { + printf("HIP: Iteration: %d: Setting up %s\n", n_iteration_, method_.c_str()); ReSolve::matrix::Csc* L = (ReSolve::matrix::Csc*)resolve_KLU_->getLFactor(); ReSolve::matrix::Csc* U = (ReSolve::matrix::Csc*)resolve_KLU_->getUFactor(); ReSolve::index_type* P = resolve_KLU_->getPOrdering(); diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp index ebcafa08..92d10bce 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp @@ -13,10 +13,6 @@ #include "IpoptConfig.h" -#if RESOLVE_WITH_CUDA -#include "NVMLHelper.hpp" -#endif - #include #include #include 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/Algorithm/LinearSolvers/test-csr.h b/src/Algorithm/LinearSolvers/test-csr.h deleted file mode 100644 index 0ddeacad..00000000 --- a/src/Algorithm/LinearSolvers/test-csr.h +++ /dev/null @@ -1,257 +0,0 @@ -#ifndef __TEST_CSR_H__ -#define __TEST_CSR_H__ - -#include -#include -#include -#include -#include -#include -#include - -using namespace std; - -namespace debug -{ - -struct CSR -{ - vector val; - vector col_ind; - vector row_ptr; -}; - -class Debug -{ -public: - Debug() - { - } - - template static T* load_b_data(std::string prefix, int iter, int& n) - { - std::fstream file; - - std::stringstream ss; - ss << prefix << "-" << iter << ".txt"; - std::string filename = ss.str(); - - file.open(filename, std::ios::in | std::fstream::binary); - file.read((char*)&n, sizeof(n)); - T* d = new T[n]; - - T p; - for (int i = 0; i < n; i++) - { - file.read((char*)&p, sizeof(T)); - d[i] = p; - } - - file.close(); - return d; - } - - template static void save_b_data(std::string prefix, int iter, int n, const T* data) - { - std::fstream file; - std::stringstream ss; - ss << prefix << "-" << iter << ".bin"; - std::string filename = ss.str(); - - file.open(filename, std::ios::out | std::fstream::binary); - - file.write((char*)&n, sizeof(n)); - // std::cout.precision(std::numeric_limits::max_digits10 - 1); - file.write((char*)(data), n * sizeof(T)); - file.close(); - } - - template static T* load_data(std::string prefix, int iter, int& n) - { - std::fstream file; - - std::stringstream ss; - ss << prefix << "-" << iter << ".txt"; - std::string filename = ss.str(); - - file.open(filename, std::ios::in); - file >> n; - T* d = new T[n]; - - T p; - for (int i = 0; i < n; i++) - { - file >> p; - d[i] = p; - } - - file.close(); - return d; - } - - template static void save_data(std::string prefix, int iter, int n, const T* data) - { - std::fstream file; - std::stringstream ss; - ss << prefix << "-" << iter << ".txt"; - std::string filename = ss.str(); - - file.open(filename, std::ios::out); - file << n << std::endl; - - // std::cout.precision(std::numeric_limits::max_digits10 - 1); - - for (int i = 0; i < n; i++) - { - // file << std::scientific << data[i] << std::endl; - file << data[i] << std::endl; - } - // file.write((char *)(data), n * sizeof(T)); - file.close(); - } - - static void save_coo_as_mm(std::string prefix, int iter, int dim, int nnz, const int* row_ptr, const int* col_ind, double* val) - { - printf("%s\t%d\n", __FILE__, __LINE__); - int cont = 0; - - std::stringstream ss; - ss << prefix << "-matrix-" << iter << ".mtx"; - std::string filename = ss.str(); - std::cout << "Filename: " << filename << std::endl; - - std::ofstream mat_file(filename); - - mat_file << "%%MatrixMarket matrix coordinate real general" << std::endl; - mat_file << "% ID: " << iter << std::endl; - mat_file << dim << "\t" << dim << "\t" << nnz << std::endl; - - int nnz_c = 0; - - for (int i = 0; i < nnz; i++) - { - mat_file << row_ptr[i] << "\t" << col_ind[i] << "\t" << val[i] << std::endl; - } - mat_file.close(); - printf("%d, %d==%d\n", dim, nnz_c, nnz); - } - - static void save_csc_as_csv(std::string prefix, int iter, int dim, int nnz, const int* row_ptr, const int* col_ind, double* val) - { - printf("%s\t%d\n", __FILE__, __LINE__); - int cont = 0; - - std::stringstream ss; - ss << prefix << "-matrix-" << iter << ".csv"; - std::string filename = ss.str(); - std::cout << "Filename: " << filename << std::endl; - - std::ofstream mat_file(filename); - - int nnz_c = 0; - - for (int i = 1; i <= dim; i++) - { - int row_start = row_ptr[i - 1]; - int row_end = row_ptr[i]; - - nnz_c += row_end - row_start; - // printf("%d, %d, %d\n", row_start, row_end, row_end - row_start); - - for (int jj = 0; jj < dim; jj++) - { - bool found = false; - for (int j = row_start; j < row_end; j++) - { - if (jj == col_ind[j]) - { - found = true; - break; - } - } - - if (found) - { - mat_file << val[cont] << ","; - cont++; - } - else - { - mat_file << "0,"; - } - } - mat_file << std::endl; - } - mat_file.close(); - printf("%d, %d==%d\n", dim, nnz_c, nnz); - } - - static void save_csc_as_mm(std::string prefix, int iter, int dim, int nnz, const int* row_ptr, const int* col_ind, double* val) - { - printf("%s\t%d\n", __FILE__, __LINE__); - int cont = 0; - - std::stringstream ss; - ss << prefix << "-matrix-" << iter << ".mtx"; - std::string filename = ss.str(); - std::cout << "Filename: " << filename << std::endl; - - std::ofstream mat_file(filename); - - mat_file << "%%MatrixMarket matrix coordinate real general" << std::endl; - mat_file << "% ID: " << iter << std::endl; - mat_file << dim << "\t" << dim << "\t" << nnz << std::endl; - - int nnz_c = 0; - - for (int i = 1; i <= dim; i++) - { - int row_start = row_ptr[i - 1]; - int row_end = row_ptr[i]; - - nnz_c += row_end - row_start; - // printf("%d, %d, %d\n", row_start, row_end, row_end - row_start); - - for (int j = row_start; j < row_end; j++) - { - // printf("%d ", col_ind[j]); - mat_file << i << "\t" << col_ind[j] + 1 << "\t" << val[cont] << std::endl; - cont++; - } - } - mat_file.close(); - printf("%d, %d==%d\n", dim, nnz_c, nnz); - } - - static void save_vec_as_csv(std::string prefix, int iter, int dim, int nnz, const int* row_ptr, const int* col_ind, double* val) - { - } - - static void save_vec_as_mm(std::string prefix, int iter, int dim, double* val) - { - printf("%s\t%d\n", __FILE__, __LINE__); - int cont = 0; - - std::stringstream ss; - ss << prefix << "-rhs-" << iter << ".mtx"; - std::string filename = ss.str(); - std::cout << "Filename: " << filename << std::endl; - - std::ofstream vec_file(filename); - - vec_file << "%%MatrixMarket matrix array real general" << std::endl; - vec_file << "% ID: " << iter << std::endl; - vec_file << dim << "\t" << 1 << std::endl; - - int nnz_c = 0; - - for (int i = 0; i < dim; i++) - { - vec_file << val[i] << std::endl; - } - vec_file.close(); - } -}; - -} // namespace debug -#endif From 4706b2c5f2c59d6a617336d0be4db293953e400c Mon Sep 17 00:00:00 2001 From: Maksudul Alam Date: Thu, 31 Jul 2025 20:04:13 +0000 Subject: [PATCH 03/19] Added GramSchmidt include --- src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp index 92d10bce..1d01ef0d 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include From 85767a62ea7b1c8ec35746b50f988b32bb7d5d81 Mon Sep 17 00:00:00 2001 From: Maksudul Alam Date: Thu, 7 Aug 2025 19:17:12 +0000 Subject: [PATCH 04/19] Added parameter for rf_solver --- src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index 0819dfd9..dd750d38 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -214,7 +214,7 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { printf("Resolve::RF Setup\n"); - resolve_Rf_ = new rf_solver(); + resolve_Rf_ = new rf_solver(workspace_GPU_); } # if RESOLVE_WITH_CUDA else if (method_ == resolve_glu) @@ -281,9 +281,15 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index A_->setDataPointers(const_cast(ia), const_cast(ja), const_cast(A_->getValues(ReSolve::memory::HOST)), ReSolve::memory::HOST); #if RESOLVE_WITH_GPU +#if RESOLVE_WITH_CUDA if (method_ == resolve_rf || method_ == resolve_rf_fgmres || method_ == resolve_glu) { A_->syncData(ReSolve::memory::DEVICE); } +#else + if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { + A_->syncData(ReSolve::memory::DEVICE); + } +#endif #endif // FACTORIZE From 0edf5fbb1a4f44b4623176b05dfa671c2f59570a Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Thu, 16 Jul 2026 15:50:35 -0700 Subject: [PATCH 05/19] Update ReSolve interface compatibility --- .../IpReSolveSolverInterface.cpp | 372 +++++++++++++----- .../IpReSolveSolverInterface.hpp | 6 + 2 files changed, 282 insertions(+), 96 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index dd750d38..8f6be162 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -9,6 +9,7 @@ #include #include +#include namespace Ipopt { @@ -16,14 +17,52 @@ namespace Ipopt static const Index dbg_verbosity = 0; #endif -ReSolveSolverInterface::ReSolveSolverInterface() : val_(NULL) +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), + pivtol_changed_(false), + 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), + factor_by_t_(2), + resolve_KLU_(NULL), + workspace_CPU_(NULL), + A_(NULL), + vec_rhs_(NULL), + vec_x_(NULL), +#if RESOLVE_WITH_GPU + workspace_GPU_(NULL), + resolve_Rf_(NULL), +# if RESOLVE_WITH_CUDA + resolve_GLU_(NULL), +# endif + GS_(NULL), + resolve_preconditioner_(NULL), + resolve_FGMRES_(NULL), +#endif + matrix_handler_(NULL), + vector_handler_(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 @@ -39,36 +78,31 @@ ReSolveSolverInterface::ReSolveSolverInterface() : val_(NULL) ReSolveSolverInterface::~ReSolveSolverInterface() { DBG_START_METH("ReSolveSolverInterface::~ReSolveSolverInterface()", dbg_verbosity); - delete[] val_; +#if RESOLVE_WITH_GPU + delete resolve_FGMRES_; + delete resolve_preconditioner_; + delete GS_; + +# if RESOLVE_WITH_CUDA + delete resolve_GLU_; +# endif + delete resolve_Rf_; +#endif + + delete resolve_KLU_; delete matrix_handler_; delete vector_handler_; - delete resolve_KLU_; #if RESOLVE_WITH_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_; + delete[] val_; } void ReSolveSolverInterface::RegisterOptions(SmartPtr roptions) @@ -153,11 +187,9 @@ bool ReSolveSolverInterface::InitializeImpl(const OptionsList& options, const st // printf("ReSolveSolverInterface::InitializeImpl is Called\n"); - Number tol; - options.GetNumericValue("resolve_tol", tol, prefix); + options.GetNumericValue("resolve_tol", pivot_tol_, prefix); - Index order_method; - options.GetIntegerValue("resolve_ordering", order_method, prefix); + options.GetIntegerValue("resolve_ordering", ordering_, prefix); Index btf; options.GetIntegerValue("resolve_btf", btf, prefix); @@ -169,12 +201,9 @@ bool ReSolveSolverInterface::InitializeImpl(const OptionsList& options, const st options.GetIntegerValue("resolve_n_skip_refactoring", n_skip_refactoring, prefix); 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); @@ -193,6 +222,9 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no 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) { @@ -203,38 +235,62 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no } #if RESOLVE_WITH_GPU - printf("Resolve::GPU Setup\n"); - workspace_GPU_ = new workspace_type(); - workspace_GPU_->initializeHandles(); + else + { + printf("Resolve::GPU Setup\n"); - matrix_handler_ = new ReSolve::MatrixHandler(workspace_GPU_); - vector_handler_ = new ReSolve::VectorHandler(workspace_GPU_); + 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) - { - printf("Resolve::RF Setup\n"); - resolve_Rf_ = new rf_solver(workspace_GPU_); - } + if (method_ == resolve_rf || method_ == resolve_rf_fgmres) + { + printf("Resolve::RF Setup\n"); + resolve_Rf_ = new rf_solver(workspace_GPU_); + } # if RESOLVE_WITH_CUDA - else if (method_ == resolve_glu) - { - printf("Resolve::GLU Setup\n"); - resolve_GLU_ = new ReSolve::LinSolverDirectCuSolverGLU(workspace_GPU_); - } + else if (method_ == resolve_glu) + { + printf("Resolve::GLU Setup\n"); + resolve_GLU_ = + new ReSolve::LinSolverDirectCuSolverGLU(workspace_GPU_); + } # endif - - if (method_ == resolve_rf_fgmres) - { - printf("Resolve::FGMRES Setup\n"); - GS_ = new ReSolve::GramSchmidt(vector_handler_, ReSolve::GramSchmidt::CGS2); - resolve_FGMRES_ = new ReSolve::LinSolverIterativeFGMRES(matrix_handler_, vector_handler_, GS_); + + if (method_ == resolve_rf_fgmres) + { + printf("Resolve::FGMRES Setup\n"); + + 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; @@ -243,21 +299,66 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no if (!initialized_) { A_ = new ReSolve::matrix::Csr(dim, dim, nonzeros); - if (val_ != NULL) - { - delete[] val_; - } + delete[] val_; val_ = new Number[nonzeros]; - A_->setDataPointers(const_cast(ia), const_cast(ja), val_, ReSolve::memory::HOST); - resolve_KLU_->setup(A_); + if( A_->setDataPointers(const_cast(ia), const_cast(ja), val_, ReSolve::memory::HOST) != 0) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to attach Ipopt matrix storage to ReSolve.\n"); + return SYMSOLVER_FATAL_ERROR; + } +#if RESOLVE_WITH_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_rhs_ = new ReSolve::vector::Vector(A_->getNumRows()); - vec_x_ = new ReSolve::vector::Vector(A_->getNumRows()); + 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; + } - vec_x_->allocate(ReSolve::memory::HOST); // for KLU +#if RESOLVE_WITH_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; @@ -276,21 +377,33 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index bool full_factor_done = false; - // Get Data from CPU and update the A Matrix + (void)ia; + (void)ja; - A_->setDataPointers(const_cast(ia), const_cast(ja), const_cast(A_->getValues(ReSolve::memory::HOST)), ReSolve::memory::HOST); + if( new_matrix ) + { + 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; + } #if RESOLVE_WITH_GPU -#if RESOLVE_WITH_CUDA - if (method_ == resolve_rf || method_ == resolve_rf_fgmres || method_ == resolve_glu) { - A_->syncData(ReSolve::memory::DEVICE); - } -#else - if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - A_->syncData(ReSolve::memory::DEVICE); - } -#endif + 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 + } // FACTORIZE @@ -487,7 +600,14 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } // Copy rhs_vals to vec_rhs - vec_rhs_->copyDataFrom(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::HOST); + if (vec_rhs_->copyFromExternal(rhs_vals, 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_); //printf("Solving using KLU!\n"); if (status != 0) @@ -524,10 +644,13 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index delete U; } - 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 @@ -538,16 +661,24 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index ReSolve::matrix::Csc* U = (ReSolve::matrix::Csc*)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); + 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; + } resolve_Rf_->setup(A_, L, U, P, Q, vec_rhs_); } - 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 } @@ -560,8 +691,15 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index 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 (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; + } + if (method_ == resolve_glu) { status = resolve_GLU_->solve(vec_rhs_, vec_x_); @@ -595,14 +733,28 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } // Copy vec_x cuda to vec_x in cpu - vec_x_->copyDataFrom(vec_x_->getData(ReSolve::memory::DEVICE), ReSolve::memory::DEVICE, ReSolve::memory::HOST); + 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; + } matrix_handler_->setValuesChanged(true, ReSolve::memory::DEVICE); } # else 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); + 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; + } int status = resolve_Rf_->solve(vec_rhs_, vec_x_); if (status != 0) @@ -614,7 +766,14 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index { resolve_FGMRES_->resetMatrix(A_); - vec_rhs_->copyDataFrom(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::DEVICE); + 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_FGMRES_->solve(vec_rhs_, vec_x_); if (status != 0) { @@ -623,7 +782,14 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } // Copy vec_x cuda to vec_x in cpu - vec_x_->copyDataFrom(vec_x_->getData(ReSolve::memory::DEVICE), ReSolve::memory::DEVICE, ReSolve::memory::HOST); + 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; + } } #endif @@ -655,7 +821,14 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } // Copy rhs_vals to vec_rhs cuda - vec_rhs_->copyDataFrom(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::HOST); + if (vec_rhs_->copyFromExternal(rhs_vals, 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; + } int status = resolve_KLU_->solve(vec_rhs_, vec_x_); if (status != 0) { @@ -665,7 +838,14 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } // copy vec_x to rhs_vals - memcpy(rhs_vals, vec_x_->getData(ReSolve::memory::HOST), (ndim_) * sizeof(ReSolve::real_type)); + 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; + } if (HaveIpData()) { diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp index 1d01ef0d..16456d95 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp @@ -21,6 +21,7 @@ #include #include #include +#include #include #include @@ -142,6 +143,10 @@ class ReSolveSolverInterface : public SparseSymLinearSolverInterface int n_iteration_; int k_; + Number pivot_tol_; + Index ordering_; + bool halt_if_singular_; + Number rcond_val_; bool use_rcond_; @@ -162,6 +167,7 @@ class ReSolveSolverInterface : public SparseSymLinearSolverInterface ReSolve::LinSolverDirectCuSolverGLU* resolve_GLU_; # endif ReSolve::GramSchmidt* GS_; + ReSolve::PreconditionerLU* resolve_preconditioner_; ReSolve::LinSolverIterativeFGMRES* resolve_FGMRES_; #endif From 091269a02db03223b4068ab62d9a6fea2074fa60 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Thu, 16 Jul 2026 20:25:36 -0700 Subject: [PATCH 06/19] Update ReSolve GPU factor setup --- .../IpReSolveSolverInterface.cpp | 91 +++++++++++++------ 1 file changed, 64 insertions(+), 27 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index 8f6be162..d4ca0daf 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -529,21 +529,29 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } else if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - //std::cout << "%" << n_iteration_ << "%" << "RF->refactorize()" << std::endl; status_refactor = resolve_Rf_->refactorize(); - if (status != 0) + 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 if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - //std::cout << "%" << n_iteration_ << "%" << "RF->refactorize()" << std::endl; - 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 @@ -624,26 +632,32 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index { printf("CUDA: 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()); - L_csc->syncData(ReSolve::memory::DEVICE); - U_csc->syncData(ReSolve::memory::DEVICE); - 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; + } } - + if (method_ == resolve_rf_fgmres && resolve_FGMRES_->setup(A_) != 0) { Jnlst().Printf( @@ -657,19 +671,42 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { printf("HIP: Iteration: %d: Setting up %s\n", n_iteration_, method_.c_str()); - 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(); - if (vec_rhs_->copyFromExternal(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::DEVICE) != 0) + 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, - "Failed to copy rhs data to ReSolve device storage.\n"); + "ReSolve HIP RF setup failed with status %d.\n", + status); return SYMSOLVER_FATAL_ERROR; } - resolve_Rf_->setup(A_, L, U, P, Q, vec_rhs_); } if (method_ == resolve_rf_fgmres && resolve_FGMRES_->setup(A_) != 0) From 4a16e48823af778ab5e3ff0dd9f7478646b299a3 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Thu, 16 Jul 2026 20:59:02 -0700 Subject: [PATCH 07/19] Handle ReSolve RF and FGMRES failures --- .../IpReSolveSolverInterface.cpp | 64 +++++++++++++++---- 1 file changed, 51 insertions(+), 13 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index d4ca0daf..bbba9eaa 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -750,22 +750,48 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index status = resolve_Rf_->solve(vec_rhs_, vec_x_); if (status != 0) { - std::cout << "RF solve status: " << status << std::endl; + 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) { - int status = resolve_Rf_->solve(vec_rhs_, vec_x_); + 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; + } + + status = resolve_FGMRES_->resetMatrix(A_); if (status != 0) { - std::cout << "RF solve status: " << status << std::endl; + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to reset the ReSolve CUDA FGMRES matrix " + "with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; } - resolve_FGMRES_->resetMatrix(A_); status = resolve_FGMRES_->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 CUDA FGMRES solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; } } @@ -793,28 +819,40 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index return SYMSOLVER_FATAL_ERROR; } - int status = resolve_Rf_->solve(vec_rhs_, vec_x_); + status = resolve_Rf_->solve(vec_rhs_, vec_x_); if (status != 0) { - std::cout << "RF 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) + if (method_ == resolve_rf_fgmres) { - - resolve_FGMRES_->resetMatrix(A_); - if (vec_rhs_->copyFromExternal(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::DEVICE) != 0) + status = resolve_FGMRES_->resetMatrix(A_); + if (status != 0) { Jnlst().Printf( J_ERROR, J_LINEAR_ALGEBRA, - "Failed to copy rhs data to ReSolve device storage.\n"); + "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) { - std::cout << "RF_FGMRES solve status: " << status << std::endl; + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve HIP FGMRES solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; } } From 0d29312255d951444e286a08ea4c6ee4604628cd Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Thu, 16 Jul 2026 21:24:10 -0700 Subject: [PATCH 08/19] Use ReSolve backend configuration --- .../IpReSolveSolverInterface.cpp | 44 +++++++++---------- .../IpReSolveSolverInterface.hpp | 13 +++--- 2 files changed, 28 insertions(+), 29 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index bbba9eaa..daef8b4a 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -49,10 +49,10 @@ ReSolveSolverInterface::ReSolveSolverInterface() A_(NULL), vec_rhs_(NULL), vec_x_(NULL), -#if RESOLVE_WITH_GPU +#ifdef RESOLVE_USE_GPU workspace_GPU_(NULL), resolve_Rf_(NULL), -# if RESOLVE_WITH_CUDA +# ifdef RESOLVE_USE_CUDA resolve_GLU_(NULL), # endif GS_(NULL), @@ -63,9 +63,9 @@ ReSolveSolverInterface::ReSolveSolverInterface() vector_handler_(NULL) { DBG_START_METH("ReSolveSolverInterface::ReSolveSolverInterface()", dbg_verbosity); -#if RESOLVE_WITH_GPU +#ifdef RESOLVE_USE_GPU printf("Resolve with GPU.\n"); -# if RESOLVE_WITH_CUDA +# ifdef RESOLVE_USE_CUDA printf("Resolve with CUDA.\n"); # else printf("Resolve with HIP.\n"); @@ -78,12 +78,12 @@ ReSolveSolverInterface::ReSolveSolverInterface() ReSolveSolverInterface::~ReSolveSolverInterface() { DBG_START_METH("ReSolveSolverInterface::~ReSolveSolverInterface()", dbg_verbosity); -#if RESOLVE_WITH_GPU +#ifdef RESOLVE_USE_GPU delete resolve_FGMRES_; delete resolve_preconditioner_; delete GS_; -# if RESOLVE_WITH_CUDA +# ifdef RESOLVE_USE_CUDA delete resolve_GLU_; # endif delete resolve_Rf_; @@ -94,7 +94,7 @@ ReSolveSolverInterface::~ReSolveSolverInterface() delete matrix_handler_; delete vector_handler_; -#if RESOLVE_WITH_GPU +#ifdef RESOLVE_USE_GPU delete workspace_GPU_; #endif delete workspace_CPU_; @@ -113,8 +113,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 @@ -234,7 +234,7 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no vector_handler_ = new ReSolve::VectorHandler(workspace_CPU_); } -#if RESOLVE_WITH_GPU +#ifdef RESOLVE_USE_GPU else { printf("Resolve::GPU Setup\n"); @@ -250,7 +250,7 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no printf("Resolve::RF Setup\n"); resolve_Rf_ = new rf_solver(workspace_GPU_); } -# if RESOLVE_WITH_CUDA +# ifdef RESOLVE_USE_CUDA else if (method_ == resolve_glu) { printf("Resolve::GLU Setup\n"); @@ -310,7 +310,7 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no "Failed to attach Ipopt matrix storage to ReSolve.\n"); return SYMSOLVER_FATAL_ERROR; } -#if RESOLVE_WITH_GPU +#ifdef RESOLVE_USE_GPU if( method_ != resolve_klu && A_->allocateMatrixData(ReSolve::memory::DEVICE) != 0 ) { @@ -344,7 +344,7 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no return SYMSOLVER_FATAL_ERROR; } -#if RESOLVE_WITH_GPU +#ifdef RESOLVE_USE_GPU if( method_ != resolve_klu ) { if( vec_rhs_->allocate(ReSolve::memory::DEVICE) != 0 @@ -391,7 +391,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index return SYMSOLVER_FATAL_ERROR; } -#if RESOLVE_WITH_GPU +#ifdef RESOLVE_USE_GPU if( method_ != resolve_klu && A_->syncData(ReSolve::memory::DEVICE) != 0 ) { @@ -460,7 +460,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // 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_); @@ -516,7 +516,7 @@ 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) { @@ -540,7 +540,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index return SYMSOLVER_FATAL_ERROR; } } -# else +# elif defined(RESOLVE_USE_HIP) if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { status = resolve_Rf_->refactorize(); @@ -627,7 +627,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // 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("CUDA: Iteration: %d: Setting up %s\n", n_iteration_, method_.c_str()); @@ -666,8 +666,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index "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) { printf("HIP: Iteration: %d: Setting up %s\n", n_iteration_, method_.c_str()); @@ -723,8 +722,7 @@ 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 @@ -806,7 +804,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } 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 diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp index 16456d95..c1ef0e58 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp @@ -12,6 +12,7 @@ #include #include "IpoptConfig.h" +#include #include #include @@ -26,10 +27,10 @@ #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; @@ -57,8 +58,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"; @@ -160,10 +161,10 @@ 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_; From 00a3146bf4083ba70a283321d0066e7ef16785f6 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 20 Jul 2026 14:03:54 -0700 Subject: [PATCH 09/19] Initialize KLU symbolic pointer --- src/Algorithm/LinearSolvers/IpKLUSolverInterface.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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; From be6d6961ed73f2404f3e5651248339ddd01d5543 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 20 Jul 2026 17:21:12 -0700 Subject: [PATCH 10/19] Improve ReSolve interface error handling --- .../IpReSolveSolverInterface.cpp | 52 ++++++++++++++++--- .../IpReSolveSolverInterface.hpp | 6 ++- 2 files changed, 50 insertions(+), 8 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index daef8b4a..26f8c849 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -208,6 +208,31 @@ bool ReSolveSolverInterface::InitializeImpl(const OptionsList& options, const st 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; } @@ -302,6 +327,8 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no delete[] val_; val_ = new Number[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) { Jnlst().Printf( @@ -380,7 +407,9 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index (void)ia; (void)ja; - if( new_matrix ) + // Ipopt updates val_ directly, so mark the host matrix current before + // synchronizing updated values to the GPU. + if (new_matrix) { if( A_->setUpdated(ReSolve::memory::HOST) != 0 ) { @@ -617,10 +646,14 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index return SYMSOLVER_FATAL_ERROR; } status = resolve_KLU_->solve(vec_rhs_, vec_x_); - //printf("Solving using KLU!\n"); if (status != 0) { - std::cout << "KLU solve status: " << status << std::endl; + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve KLU solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; } } @@ -793,7 +826,8 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } } - // Copy vec_x cuda to vec_x in cpu + // 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( @@ -854,7 +888,8 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } } - // Copy vec_x cuda to vec_x in cpu + // 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( @@ -905,7 +940,12 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index int status = resolve_KLU_->solve(vec_rhs_, vec_x_); if (status != 0) { - std::cout << "KLU solve status: " << status << std::endl; + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve KLU solve failed with status %d.\n", + status); + return SYMSOLVER_FATAL_ERROR; } } } diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp index c1ef0e58..a131e3b5 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp @@ -14,8 +14,10 @@ #include "IpoptConfig.h" #include -#include -#include +#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 From cb583a8c52ee905bc1abe2138f1ed888007e1656 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 20 Jul 2026 19:32:14 -0700 Subject: [PATCH 11/19] Support multiple RHS in ReSolve MultiSolve --- .../IpReSolveSolverInterface.cpp | 388 +++++++++++------- 1 file changed, 242 insertions(+), 146 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index 26f8c849..d553cf33 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -7,9 +7,11 @@ #include "IpReSolveSolverInterface.hpp" #include "IpoptConfig.h" +#include #include #include #include +#include namespace Ipopt { @@ -399,6 +401,14 @@ 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; @@ -636,24 +646,40 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } } - // Copy rhs_vals to vec_rhs - if (vec_rhs_->copyFromExternal(rhs_vals, 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) + for (Index irhs = 0; irhs < nrhs; ++irhs) { - Jnlst().Printf( - J_ERROR, - J_LINEAR_ALGEBRA, - "ReSolve KLU solve failed with status %d.\n", - status); - return SYMSOLVER_FATAL_ERROR; + 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; + } + } } } @@ -758,145 +784,186 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index # ifdef RESOLVE_USE_CUDA if (method_ == resolve_glu || method_ == resolve_rf || method_ == resolve_rf_fgmres) { - // Copy rhs_vals to vec_rhs cuda - if (vec_rhs_->copyFromExternal(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::DEVICE) != 0) + for (Index irhs = 0; irhs < nrhs; ++irhs) { - Jnlst().Printf( - J_ERROR, - J_LINEAR_ALGEBRA, - "Failed to copy rhs data to ReSolve device storage.\n"); - return SYMSOLVER_FATAL_ERROR; - } + 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; + } - if (method_ == resolve_glu) - { - status = resolve_GLU_->solve(vec_rhs_, vec_x_); - if (status != 0) + if (method_ == resolve_glu) { - std::cout << "GLU 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) - { - status = resolve_Rf_->solve(vec_rhs_, vec_x_); - if (status != 0) + else if (method_ == resolve_rf) { - Jnlst().Printf( - J_ERROR, - J_LINEAR_ALGEBRA, - "ReSolve CUDA RF solve failed with status %d.\n", - status); - return SYMSOLVER_FATAL_ERROR; + 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) + else if (method_ == resolve_rf_fgmres) { - Jnlst().Printf( - J_ERROR, - J_LINEAR_ALGEBRA, - "ReSolve CUDA RF initial solve failed with status %d.\n", - status); - return SYMSOLVER_FATAL_ERROR; + 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; + } + + 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; + } } - status = resolve_FGMRES_->resetMatrix(A_); - if (status != 0) + // 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 reset the ReSolve CUDA FGMRES matrix " - "with status %d.\n", - status); + J_ERROR, + J_LINEAR_ALGEBRA, + "Failed to synchronize the Resolve solution to the host.\n"); return SYMSOLVER_FATAL_ERROR; } - - status = resolve_FGMRES_->solve(vec_rhs_, vec_x_); - if (status != 0) + if (nrhs > 1) { - Jnlst().Printf( - J_ERROR, - J_LINEAR_ALGEBRA, - "ReSolve CUDA FGMRES solve failed with status %d.\n", - status); - return SYMSOLVER_FATAL_ERROR; + 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; + } } } - - // 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; - } matrix_handler_->setValuesChanged(true, ReSolve::memory::DEVICE); } # elif defined(RESOLVE_USE_HIP) if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - // Copy rhs_vals to vec_rhs cuda - if (vec_rhs_->copyFromExternal(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::DEVICE) != 0) + for (Index irhs = 0; irhs < nrhs; ++irhs) { - Jnlst().Printf( - J_ERROR, - J_LINEAR_ALGEBRA, - "Failed to copy rhs data to ReSolve device storage.\n"); - return SYMSOLVER_FATAL_ERROR; - } - - status = resolve_Rf_->solve(vec_rhs_, vec_x_); - if (status != 0) - { - Jnlst().Printf( - J_ERROR, - J_LINEAR_ALGEBRA, - "ReSolve HIP RF solve failed with status %d.\n", - status); - return SYMSOLVER_FATAL_ERROR; - } + 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; + } - if (method_ == resolve_rf_fgmres) - { - status = resolve_FGMRES_->resetMatrix(A_); + status = resolve_Rf_->solve(vec_rhs_, vec_x_); if (status != 0) { Jnlst().Printf( - J_ERROR, - J_LINEAR_ALGEBRA, - "Failed to reset the ReSolve HIP FGMRES matrix " - "with status %d.\n", - status); + J_ERROR, + J_LINEAR_ALGEBRA, + "ReSolve HIP RF solve failed with status %d.\n", + status); return SYMSOLVER_FATAL_ERROR; } - status = resolve_FGMRES_->solve(vec_rhs_, vec_x_); - if (status != 0) + 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; + } + } + + // 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, - "ReSolve HIP FGMRES solve failed with status %d.\n", - status); + 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_; - // 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 (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 @@ -928,36 +995,65 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } } - // Copy rhs_vals to vec_rhs cuda - if (vec_rhs_->copyFromExternal(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::HOST) != 0) + for (Index irhs = 0; irhs < nrhs; ++irhs) { - Jnlst().Printf( - J_ERROR, - J_LINEAR_ALGEBRA, - "Failed to copy rhs data to ReSolve host storage.\n"); - return SYMSOLVER_FATAL_ERROR; - } - int 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; + 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 - if (vec_x_->copyToExternal(rhs_vals, ReSolve::memory::HOST, ReSolve::memory::HOST) != 0) + if (nrhs == 1) { - Jnlst().Printf( - J_ERROR, - J_LINEAR_ALGEBRA, - "Failed to copy the Resolve solution to Ipopt.\n"); - return SYMSOLVER_FATAL_ERROR; + 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()) From cb5f6ee9069d277b2f7ac8512ccdd8a75dfe9b85 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Mon, 20 Jul 2026 20:06:19 -0700 Subject: [PATCH 12/19] Handle ReSolve refactorization failures --- .../IpReSolveSolverInterface.cpp | 28 +++++++++++++++++-- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index d553cf33..2cc7d729 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -170,7 +170,7 @@ void ReSolveSolverInterface::RegisterOptions(SmartPtr roption roptions->AddIntegerOption("resolve_n_skip_refactoring", // "How many iterations to skip refactoring", // 1, // - "Integer, Start Refactoring after k-th iteration", // + "Number of initial KLU iterations before switching to the selected refactorization method. Must be at least 1.", // false); roptions->AddBoolOption("resolve_use_rcond", // @@ -201,6 +201,18 @@ bool ReSolveSolverInterface::InitializeImpl(const OptionsList& options, const st 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. + if (n_skip_refactoring < 1) + { + Jnlst().Printf( + J_ERROR, + J_LINEAR_ALGEBRA, + "resolve_n_skip_refactoring must be at least 1.\n"); + return false; + } + k_ = n_skip_refactoring; options.GetBoolValue("resolve_halt_if_singular", halt_if_singular_, prefix); @@ -563,7 +575,12 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index 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) @@ -601,7 +618,12 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index 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; From a3e657f723ceefd19ab5d7479f89997844c51014 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Fri, 24 Jul 2026 02:39:22 -0700 Subject: [PATCH 13/19] Handle ReSolve GLU setup failures --- .../IpReSolveSolverInterface.cpp | 23 +++++++++++++++---- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index 2cc7d729..79b194a0 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -518,13 +518,26 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index 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 } From af602695dfe69dee50fdf6c7534d01a8c1d13dc7 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Fri, 24 Jul 2026 02:43:08 -0700 Subject: [PATCH 14/19] Remove unsupported ReSolve KLU options --- .../IpReSolveSolverInterface.cpp | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index 79b194a0..974dffd9 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -149,19 +149,6 @@ 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 @@ -193,12 +180,6 @@ bool ReSolveSolverInterface::InitializeImpl(const OptionsList& options, const st options.GetIntegerValue("resolve_ordering", ordering_, prefix); - Index btf; - options.GetIntegerValue("resolve_btf", btf, prefix); - - Index scale; - options.GetIntegerValue("resolve_scale", scale, prefix); - Index n_skip_refactoring; options.GetIntegerValue("resolve_n_skip_refactoring", n_skip_refactoring, prefix); From ff61d8b9f093f69b3fc4ab174078d3c8e06eda50 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Fri, 24 Jul 2026 03:27:11 -0700 Subject: [PATCH 15/19] Correct ReSolve quality and option handling --- .../IpReSolveSolverInterface.cpp | 35 ++++--------------- .../IpReSolveSolverInterface.hpp | 1 - 2 files changed, 6 insertions(+), 30 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index 974dffd9..85ec7ead 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -34,7 +34,6 @@ ReSolveSolverInterface::ReSolveSolverInterface() ndim_(0), val_(NULL), numneg_(0), - pivtol_changed_(false), re_factorize_(false), factorize_(false), method_(resolve_klu), @@ -154,10 +153,12 @@ void ReSolveSolverInterface::RegisterOptions(SmartPtr roption 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, // - "Number of initial KLU iterations before switching to the selected refactorization method. Must be at least 1.", // + 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", // @@ -185,15 +186,6 @@ bool ReSolveSolverInterface::InitializeImpl(const OptionsList& options, const st // ReSolve's GPU refactorization methods require an initial KLU solve // to construct the factors and permutations used during setup. - if (n_skip_refactoring < 1) - { - Jnlst().Printf( - J_ERROR, - J_LINEAR_ALGEBRA, - "resolve_n_skip_refactoring must be at least 1.\n"); - return false; - } - k_ = n_skip_refactoring; options.GetBoolValue("resolve_halt_if_singular", halt_if_singular_, prefix); @@ -385,7 +377,6 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no n_iteration_ = 0; initialized_ = true; - pivtol_changed_ = false; return retval; } @@ -535,21 +526,6 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // 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 // check if a re-factorization has to be done DBG_PRINT((1, "new_matrix = %d\n", new_matrix)); @@ -1097,7 +1073,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 a131e3b5..39ce0f11 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp @@ -139,7 +139,6 @@ 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_; From d52bce1390e0f827725f1c2e5ebc5cdfbfb10bf4 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Fri, 24 Jul 2026 03:39:34 -0700 Subject: [PATCH 16/19] Clean up ReSolve interface logging --- .../IpReSolveSolverInterface.cpp | 90 +++++++------------ .../IpReSolveSolverInterface.hpp | 3 - 2 files changed, 30 insertions(+), 63 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index 85ec7ead..55d43082 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -9,7 +9,6 @@ #include #include -#include #include #include @@ -44,7 +43,6 @@ ReSolveSolverInterface::ReSolveSolverInterface() halt_if_singular_(false), rcond_val_(1e-128), use_rcond_(false), - factor_by_t_(2), resolve_KLU_(NULL), workspace_CPU_(NULL), A_(NULL), @@ -64,16 +62,6 @@ ReSolveSolverInterface::ReSolveSolverInterface() vector_handler_(NULL) { DBG_START_METH("ReSolveSolverInterface::ReSolveSolverInterface()", dbg_verbosity); -#ifdef RESOLVE_USE_GPU - printf("Resolve with GPU.\n"); -# ifdef RESOLVE_USE_CUDA - printf("Resolve with CUDA.\n"); -# else - printf("Resolve with HIP.\n"); -# endif -#else - printf("Resolve with CPU. CUDA or HIP Unavailable.\n"); -#endif } ReSolveSolverInterface::~ReSolveSolverInterface() @@ -174,9 +162,6 @@ void ReSolveSolverInterface::RegisterOptions(SmartPtr roption bool ReSolveSolverInterface::InitializeImpl(const OptionsList& options, const std::string& prefix) { - - // printf("ReSolveSolverInterface::InitializeImpl is Called\n"); - options.GetNumericValue("resolve_tol", pivot_tol_, prefix); options.GetIntegerValue("resolve_ordering", ordering_, prefix); @@ -229,7 +214,6 @@ 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_) { @@ -240,7 +224,6 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no if (method_ == resolve_klu) { - printf("Resolve::KLU Setup\n"); workspace_CPU_ = new ReSolve::LinAlgWorkspaceCpu(); matrix_handler_ = new ReSolve::MatrixHandler(workspace_CPU_); vector_handler_ = new ReSolve::VectorHandler(workspace_CPU_); @@ -249,8 +232,6 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no #ifdef RESOLVE_USE_GPU else { - printf("Resolve::GPU Setup\n"); - workspace_GPU_ = new workspace_type(); workspace_GPU_->initializeHandles(); @@ -259,13 +240,11 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - printf("Resolve::RF Setup\n"); resolve_Rf_ = new rf_solver(workspace_GPU_); } # ifdef RESOLVE_USE_CUDA else if (method_ == resolve_glu) { - printf("Resolve::GLU Setup\n"); resolve_GLU_ = new ReSolve::LinSolverDirectCuSolverGLU(workspace_GPU_); } @@ -273,8 +252,6 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no if (method_ == resolve_rf_fgmres) { - printf("Resolve::FGMRES Setup\n"); - GS_ = new ReSolve::GramSchmidt( vector_handler_, ReSolve::GramSchmidt::CGS2); @@ -306,7 +283,6 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no // Store size for later use ndim_ = dim; nonzeros_ = nonzeros; - printf("Using Refactorization after %d iterations\n\n", k_); if (!initialized_) { @@ -428,17 +404,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index #endif } - // FACTORIZE - - // Every factor_by_t_ iteration do a Factorization!!! - if (n_iteration_ % factor_by_t_ == 0) - { - // factorize_ = true; - // re_factorize_ = true; - } - if( n_iteration_ == 0){ - printf("First Iteration: %d: Performing KLU Factorization\n", n_iteration_); // Symbolic Factorization if (HaveIpData()) { @@ -447,8 +413,11 @@ 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()) @@ -459,8 +428,6 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index if (factorize_ && (new_matrix || re_factorize_)) { - // printf("Iteration: %d: Performing KLU Factorization\n", n_iteration_); - // perform the factorization if (HaveIpData()) { @@ -468,17 +435,18 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } // First Factorization is always done by KLU - // std::cout << "%" << n_iteration_ << "%" << "KLU FULL FACTORIZATION" << 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) @@ -486,8 +454,6 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index #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(); ReSolve::index_type* P = resolve_KLU_->getPOrdering(); @@ -523,7 +489,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_); } // REFACTORIZE @@ -541,7 +506,6 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index // Actual Refactorize if (method_ == resolve_glu) { - //std::cout << "%" << n_iteration_ << "%" << "GLU->refactorize()" << std::endl; status = resolve_GLU_->refactorize(); if (status != 0) { @@ -584,7 +548,6 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index if (method_ == resolve_klu) { - //std::cout << "%" << n_iteration_ << "%" << "KLU->refactorize()" << std::endl; status = resolve_KLU_->refactorize(); if (status != 0) { @@ -617,12 +580,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 @@ -631,8 +596,11 @@ 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; } } @@ -681,8 +649,6 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index # ifdef RESOLVE_USE_CUDA if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - printf("CUDA: Iteration: %d: Setting up %s\n", n_iteration_, method_.c_str()); - ReSolve::matrix::Sparse* L = resolve_KLU_->getLFactor(); ReSolve::matrix::Sparse* U = resolve_KLU_->getUFactor(); ReSolve::index_type* P = resolve_KLU_->getPOrdering(); @@ -720,7 +686,6 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index # elif defined(RESOLVE_USE_HIP) if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - printf("HIP: Iteration: %d: Setting up %s\n", n_iteration_, method_.c_str()); ReSolve::matrix::Sparse* L = resolve_KLU_->getLFactor(); ReSolve::matrix::Sparse* U = resolve_KLU_->getUFactor(); ReSolve::index_type* P = resolve_KLU_->getPOrdering(); @@ -859,7 +824,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index Jnlst().Printf( J_ERROR, J_LINEAR_ALGEBRA, - "Failed to synchronize the Resolve solution to the host.\n"); + "Failed to synchronize the ReSolve solution to the host.\n"); return SYMSOLVER_FATAL_ERROR; } if (nrhs > 1) @@ -937,7 +902,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index Jnlst().Printf( J_ERROR, J_LINEAR_ALGEBRA, - "Failed to synchronize the Resolve solution to the host.\n"); + "Failed to synchronize the ReSolve solution to the host.\n"); return SYMSOLVER_FATAL_ERROR; } if (nrhs > 1) @@ -966,12 +931,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 @@ -980,8 +947,11 @@ 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; } } diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp index 39ce0f11..efbc91ef 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.hpp @@ -9,7 +9,6 @@ #include #include -#include #include "IpoptConfig.h" #include @@ -152,8 +151,6 @@ class ReSolveSolverInterface : public SparseSymLinearSolverInterface Number rcond_val_; bool use_rcond_; - int factor_by_t_; - ReSolve::LinSolverDirectKLU* resolve_KLU_; ReSolve::LinAlgWorkspaceCpu* workspace_CPU_; From e95574080df1a60ca1e56b2fce34429ab8cf7a23 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Fri, 24 Jul 2026 05:40:52 -0700 Subject: [PATCH 17/19] Add ReSolve solver interface regression coverage --- src/Makefile.am | 2 +- src/Makefile.in | 38 ++--- test/Makefile.am | 15 ++ test/Makefile.in | 48 ++++++- test/resolve_multirhs.cpp | 285 ++++++++++++++++++++++++++++++++++++++ 5 files changed, 362 insertions(+), 26 deletions(-) create mode 100644 test/resolve_multirhs.cpp 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; +} From d7acdd86592f2042913ecf5085f246e1ea1cfab3 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Fri, 24 Jul 2026 06:02:37 -0700 Subject: [PATCH 18/19] Initialize ReSolve RF factorization after setup --- .../IpReSolveSolverInterface.cpp | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index 55d43082..ed01a415 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -673,6 +673,17 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index 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 && resolve_FGMRES_->setup(A_) != 0) @@ -722,6 +733,17 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index 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 && resolve_FGMRES_->setup(A_) != 0) From 4778fdb007872e809e9ee074721fd9bbe5151689 Mon Sep 17 00:00:00 2001 From: Tamar DeWilde Date: Fri, 24 Jul 2026 09:29:26 -0700 Subject: [PATCH 19/19] Clean up ReSolve interface warnings --- src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp index ed01a415..1443d033 100644 --- a/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp +++ b/src/Algorithm/LinearSolvers/IpReSolveSolverInterface.cpp @@ -288,7 +288,7 @@ ESymSolverStatus ReSolveSolverInterface::InitializeStructure(Index dim, Index no { A_ = new ReSolve::matrix::Csr(dim, dim, nonzeros); delete[] val_; - val_ = new Number[nonzeros]; + 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_. @@ -370,12 +370,13 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } int status; - int status_refactor = 0; bool full_factor_done = false; (void)ia; (void)ja; + (void)check_NegEVals; + (void)numberOfNegEVals; // Ipopt updates val_ directly, so mark the host matrix current before // synchronizing updated values to the GPU. @@ -519,7 +520,7 @@ ESymSolverStatus ReSolveSolverInterface::MultiSolve(bool new_matrix, const Index } else if (method_ == resolve_rf || method_ == resolve_rf_fgmres) { - status_refactor = resolve_Rf_->refactorize(); + int status_refactor = resolve_Rf_->refactorize(); if (status_refactor != 0) { Jnlst().Printf(