Skip to content

Commit d6a3deb

Browse files
adhamsishakedregev
adhamsi
and
shakedregev
committed
Implement the HyKKT Cholesky module (#350)
* initial solver and cuda header * cuda implementation * minimal test on cuda passes * cpu implementation * tell cmake to find cholmod * add template hip implementation * dependencies * cpu: memory freeing and unused fields * Apply pre-commmit fixes * error handling and correct memory in cuda * randomized testing: different sparsity patterns * Apply pre-commmit fixes * fix for cuda: resetting workspace * Apply pre-commmit fixes * test reusing sparsity pattern * Apply pre-commmit fixes * Apply pre-commmit fixes * draft: rocsolver * Apply pre-commmit fixes * fix build issue * documentation * Apply pre-commmit fixes * error message when no gpu support Co-authored-by: Shaked Regev <35384901+shakedregev@users.noreply.github.com> * address minor review comments * Apply pre-commmit fixes * remove commented code * more detail on documentation * Apply pre-commmit fixes --------- Co-authored-by: adhamsi <adhamsi@users.noreply.github.com> Co-authored-by: Shaked Regev <35384901+shakedregev@users.noreply.github.com>
1 parent 6f9f69c commit d6a3deb

16 files changed

+1233
-3
lines changed

‎cmake/FindSuiteSparse.cmake‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,12 @@ Author(s):
99
- Cameron Rutherford <cameron.rutherford@pnnl.gov>
1010
1111
]]
12-
set(SUITESPARSE_MODULES amd colamd klu suitesparseconfig)
12+
set(SUITESPARSE_MODULES
13+
amd
14+
colamd
15+
klu
16+
cholmod
17+
suitesparseconfig)
1318

1419
find_library(SUITESPARSE_LIBRARY
1520
NAMES

‎resolve/CMakeLists.txt‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -149,7 +149,6 @@ endif()
149149
if(RESOLVE_USE_KLU)
150150
add_subdirectory(hykkt)
151151
endif()
152-
#list(APPEND ReSolve_Targets_List resolve_hykkt)
153152

154153
# Set installable targets
155154
install(TARGETS ${ReSolve_Targets_List} EXPORT ReSolveTargets)

‎resolve/hykkt/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,3 @@
11
add_subdirectory(permutation)
22
add_subdirectory(ruiz)
3+
add_subdirectory(cholesky)

‎resolve/hykkt/cholesky/CMakeLists.txt‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
#[[
2+
3+
@brief Build ReSolve matrix module
4+
5+
@author Slaven Peles <peless@ornl.gov>
6+
7+
]]
8+
9+
# C++ code
10+
set(HyKKT_CHOL_SRC
11+
CholeskySolver.cpp
12+
CholeskySolverCpu.cpp
13+
)
14+
15+
# C++ code that depends on CUDA SDK libraries
16+
set(HyKKT_CUDASDK_SRC
17+
CholeskySolverCuda.cpp
18+
)
19+
20+
# and on HIP
21+
set(HyKKT_ROCM_SRC
22+
CholeskySolverHip.cpp
23+
)
24+
25+
# Header files to be installed
26+
set(HyKKT_CHOL_HEADER_INSTALL
27+
CholeskySolver.hpp
28+
CholeskySolverImpl.hpp
29+
CholeskySolverCpu.hpp
30+
CholeskySolverCuda.hpp
31+
CholeskySolverHip.hpp
32+
)
33+
34+
add_library(resolve_hykkt_chol SHARED ${HyKKT_CHOL_SRC})
35+
36+
target_link_libraries(resolve_hykkt_chol PUBLIC ${suitesparse_cholmod})
37+
target_include_directories(resolve_hykkt_chol PUBLIC ${SUITESPARSE_INCLUDE_DIR})
38+
39+
# Link to CUDA ReSolve backend if CUDA is support enabled
40+
if (RESOLVE_USE_CUDA)
41+
target_sources(resolve_hykkt_chol PRIVATE ${HyKKT_CUDASDK_SRC})
42+
target_link_libraries(resolve_hykkt_chol PUBLIC resolve_backend_cuda)
43+
endif()
44+
45+
if (RESOLVE_USE_HIP)
46+
target_sources(resolve_hykkt_chol PRIVATE ${HyKKT_ROCM_SRC})
47+
target_link_libraries(resolve_hykkt_chol PUBLIC resolve_backend_hip)
48+
endif()
49+
50+
# Link to dummy device backend if GPU support is not enabled
51+
if (NOT RESOLVE_USE_GPU)
52+
target_link_libraries(resolve_hykkt_chol PUBLIC resolve_backend_cpu)
53+
endif()
54+
55+
target_link_libraries(resolve_hykkt_chol PUBLIC resolve_workspace
56+
resolve_vector
57+
resolve_matrix
58+
resolve_logger)
59+
60+
target_include_directories(resolve_hykkt_chol INTERFACE
61+
$<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}>
62+
$<INSTALL_INTERFACE:include>
63+
)
64+
65+
install(FILES ${HyKKT_CHOL_HEADER_INSTALL} DESTINATION include/resolve/hykkt)

‎resolve/hykkt/cholesky/CholeskySolver.cpp‎

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
/**
2+
* @file CholeskySolver.cpp
3+
* @author Adham Ibrahim (ibrahimas@ornl.gov)
4+
* @brief Cholesky decomposition solver implementation
5+
*/
6+
7+
#include "CholeskySolver.hpp"
8+
9+
#include "CholeskySolverCpu.hpp"
10+
#ifdef RESOLVE_USE_CUDA
11+
#include "CholeskySolverCuda.hpp"
12+
#elif defined(RESOLVE_USE_HIP)
13+
#include "CholeskySolverHip.hpp"
14+
#endif
15+
16+
namespace ReSolve
17+
{
18+
using real_type = ReSolve::real_type;
19+
using out = ReSolve::io::Logger;
20+
21+
namespace hykkt
22+
{
23+
/**
24+
* @brief Cholesky Solver constructor
25+
* @param[in] memspace - memory space to use for computations
26+
*/
27+
CholeskySolver::CholeskySolver(memory::MemorySpace memspace)
28+
: memspace_(memspace)
29+
{
30+
if (memspace_ == memory::HOST)
31+
{
32+
impl_ = new CholeskySolverCpu();
33+
}
34+
else
35+
{
36+
#ifdef RESOLVE_USE_CUDA
37+
impl_ = new CholeskySolverCuda();
38+
#elif defined(RESOLVE_USE_HIP)
39+
impl_ = new CholeskySolverHip();
40+
#else
41+
out::error() << "Memory space set to DEVICE, though no GPU support enabled. Must enable RESOLVE_USE_CUDA or RESOLVE_USE_HIP.\n";
42+
exit(1);
43+
#endif
44+
}
45+
}
46+
47+
/**
48+
* @brief Cholesky Solver destructor
49+
*/
50+
CholeskySolver::~CholeskySolver()
51+
{
52+
delete impl_;
53+
}
54+
55+
/**
56+
* @brief Loads or reloads matrix pointer to the solver
57+
* @param[in] A - pointer to the matrix in CSR format
58+
*/
59+
void CholeskySolver::addMatrixInfo(matrix::Csr* A)
60+
{
61+
A_ = A;
62+
impl_->addMatrixInfo(A);
63+
}
64+
65+
/**
66+
* @brief Performs symbolic analysis. Determines the sparsity pattern of
67+
* the factor L. Values will be computed by numerical analysis.
68+
* This need only be called once as long as the sparsity pattern does not change.
69+
*/
70+
void CholeskySolver::symbolicAnalysis()
71+
{
72+
impl_->symbolicAnalysis();
73+
}
74+
75+
/**
76+
* @brief Sets the pivot tolerance for the solver.
77+
*
78+
* This is only used in the CUDA implementation. For other backends,
79+
* it is ignored.
80+
*
81+
* @param[in] tol - pivot tolerance value
82+
*/
83+
void CholeskySolver::setPivotTolerance(real_type tol)
84+
{
85+
tol_ = tol;
86+
}
87+
88+
/**
89+
* @brief Performs numerical factorization. Fills in the values of the factor L such
90+
* that LL^T = A.
91+
*/
92+
void CholeskySolver::numericalFactorization()
93+
{
94+
impl_->numericalFactorization(tol_);
95+
}
96+
97+
/**
98+
* @brief Solves the linear system Ax = b and stores the result in x.
99+
*
100+
* @pre The vector x is allocated in the given memspace.
101+
*
102+
* @param[out] x - pointer to the solution vector
103+
* @param[in] b - pointer to the right-hand side vector
104+
*/
105+
void CholeskySolver::solve(vector::Vector* x, vector::Vector* b)
106+
{
107+
impl_->solve(x, b);
108+
x->setDataUpdated(memspace_);
109+
}
110+
} // namespace hykkt
111+
} // namespace ReSolve

‎resolve/hykkt/cholesky/CholeskySolver.hpp‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
/**
2+
* @file CholeskySolver.hpp
3+
* @author Adham Ibrahim (ibrahimas@ornl.gov)
4+
* @brief Cholesky decomposition solver header
5+
*/
6+
7+
#pragma once
8+
#include "CholeskySolverImpl.hpp"
9+
#include <resolve/MemoryUtils.hpp>
10+
#include <resolve/matrix/Csr.hpp>
11+
#include <resolve/vector/Vector.hpp>
12+
13+
namespace ReSolve
14+
{
15+
namespace hykkt
16+
{
17+
class CholeskySolver
18+
{
19+
public:
20+
CholeskySolver(memory::MemorySpace memspace);
21+
~CholeskySolver();
22+
23+
void addMatrixInfo(matrix::Csr* A);
24+
void symbolicAnalysis();
25+
void setPivotTolerance(real_type tol);
26+
void numericalFactorization();
27+
void solve(vector::Vector* x, vector::Vector* b);
28+
29+
private:
30+
memory::MemorySpace memspace_;
31+
32+
matrix::Csr* A_;
33+
real_type tol_ = 1e-12;
34+
CholeskySolverImpl* impl_;
35+
};
36+
} // namespace hykkt
37+
} // namespace ReSolve

‎resolve/hykkt/cholesky/CholeskySolverCpu.cpp‎

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
/**
2+
* @file CholeskySolverCpu.cpp
3+
* @author Adham Ibrahim (ibrahimas@ornl.gov)
4+
* @brief CPU implementation of Cholesky Solver
5+
*/
6+
7+
#include "CholeskySolverCpu.hpp"
8+
9+
namespace ReSolve
10+
{
11+
using real_type = ReSolve::real_type;
12+
using out = ReSolve::io::Logger;
13+
14+
namespace hykkt
15+
{
16+
CholeskySolverCpu::CholeskySolverCpu()
17+
{
18+
Common_.nmethods = 1;
19+
// Use natural ordering
20+
Common_.method[0].ordering = CHOLMOD_NATURAL;
21+
22+
A_chol_ = nullptr;
23+
factorization_ = nullptr;
24+
cholmod_start(&Common_);
25+
}
26+
27+
CholeskySolverCpu::~CholeskySolverCpu()
28+
{
29+
if (A_chol_)
30+
{
31+
cholmod_free_sparse(&A_chol_, &Common_);
32+
}
33+
if (factorization_)
34+
{
35+
cholmod_free_factor(&factorization_, &Common_);
36+
}
37+
cholmod_finish(&Common_);
38+
}
39+
40+
void CholeskySolverCpu::addMatrixInfo(matrix::Csr* A)
41+
{
42+
if (A_chol_)
43+
{
44+
cholmod_free_sparse(&A_chol_, &Common_);
45+
}
46+
A_chol_ = convertToCholmod(A);
47+
}
48+
49+
/**
50+
* @brief Performs symbolic analysis
51+
*
52+
* Uses the `cholmod_analyze` routine.
53+
*/
54+
void CholeskySolverCpu::symbolicAnalysis()
55+
{
56+
factorization_ = cholmod_analyze(A_chol_, &Common_);
57+
if (Common_.status < 0)
58+
{
59+
out::error() << "Cholesky symbolic analysis failed with status: " << Common_.status << "\n";
60+
}
61+
}
62+
63+
/**
64+
* @brief Performs symbolic analysis
65+
* @brief Performs numerical factorization
66+
*
67+
* Uses the `cholmod_factorize` routine.
68+
*
69+
* @param[in] tol - This is ignored in the CPU implementation
70+
*/
71+
void CholeskySolverCpu::numericalFactorization(real_type tol)
72+
{
73+
(void) tol; // Mark tol as unused
74+
75+
cholmod_factorize(A_chol_, factorization_, &Common_);
76+
if (Common_.status < 0)
77+
{
78+
out::error() << "Cholesky factorization failed with status: " << Common_.status << "\n";
79+
}
80+
}
81+
82+
/**
83+
* @brief Solves the system Ax = b and stores the result in x.
84+
*
85+
* Uses the `cholmod_solve` routine.
86+
*
87+
* @param[out] x - The solution vector.
88+
* @param[in] b - The right-hand side vector.
89+
*/
90+
void CholeskySolverCpu::solve(vector::Vector* x, vector::Vector* b)
91+
{
92+
cholmod_dense* b_chol = convertToCholmod(b);
93+
cholmod_dense* x_chol = cholmod_solve(CHOLMOD_A, factorization_, b_chol, &Common_);
94+
if (Common_.status < 0)
95+
{
96+
out::error() << "Cholesky solve failed with status: " << Common_.status << "\n";
97+
}
98+
x->copyDataFrom(static_cast<real_type*>(x_chol->x), memory::HOST, memory::HOST);
99+
}
100+
101+
/**
102+
* @brief Converts a ReSolve type CSR matrix to CHOLMOD format.
103+
*
104+
* @param[in] A - pointer to the CSR matrix
105+
* @return Pointer to the equivalent CHOLMOD sparse matrix
106+
*/
107+
cholmod_sparse* CholeskySolverCpu::convertToCholmod(matrix::Csr* A)
108+
{
109+
A_chol_ = cholmod_allocate_sparse((size_t) A->getNumRows(),
110+
(size_t) A->getNumColumns(),
111+
(size_t) A->getNnz(),
112+
1,
113+
1,
114+
1,
115+
CHOLMOD_REAL,
116+
&Common_);
117+
mem_.copyArrayHostToHost(
118+
static_cast<int*>(A_chol_->p), A->getRowData(memory::HOST), A->getNumRows() + 1);
119+
mem_.copyArrayHostToHost(
120+
static_cast<int*>(A_chol_->i), A->getColData(memory::HOST), A->getNnz());
121+
mem_.copyArrayHostToHost(
122+
static_cast<double*>(A_chol_->x), A->getValues(memory::HOST), A->getNnz());
123+
124+
return A_chol_;
125+
}
126+
127+
/**
128+
* @brief Converts a ReSolve type vector to CHOLMOD format.
129+
*
130+
* @param[in] v - pointer to the vector
131+
* @return Pointer to the equivalent CHOLMOD dense matrix
132+
*/
133+
cholmod_dense* CholeskySolverCpu::convertToCholmod(vector::Vector* v)
134+
{
135+
cholmod_dense* v_chol = cholmod_allocate_dense((size_t) v->getSize(),
136+
1,
137+
(size_t) v->getSize(),
138+
CHOLMOD_REAL,
139+
&Common_);
140+
mem_.copyArrayHostToHost(
141+
static_cast<double*>(v_chol->x), v->getData(memory::HOST), v->getSize());
142+
return v_chol;
143+
}
144+
} // namespace hykkt
145+
} // namespace ReSolve

0 commit comments

Comments
 (0)