Skip to content

Commit

Permalink
Implement system solver class (#86)
Browse files Browse the repository at this point in the history
* First stab at system solver.

* Adjust direct and iterative solver interfaces.

* Clean-up KLU interface, avoid false positives in functionality tests.

* Add accessors to solver components; use more consistent names in SystemSolver class.

* SystemSolver lets its modules manage their own parameter settings.

* Set KLU defaults at only one place.

---------

Co-authored-by: kswirydo <[email protected]>
  • Loading branch information
pelesh and kswirydo authored Dec 5, 2023
1 parent 2d2b21a commit 5129c8c
Show file tree
Hide file tree
Showing 46 changed files with 2,042 additions and 312 deletions.
18 changes: 17 additions & 1 deletion examples/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,10 @@ target_link_libraries(klu_klu.exe PRIVATE ReSolve)
add_executable(klu_klu_standalone.exe r_KLU_KLU_standalone.cpp)
target_link_libraries(klu_klu_standalone.exe PRIVATE ReSolve)

# Build example with a configurable system solver
add_executable(system.exe r_SysSolver.cpp)
target_link_libraries(system.exe PRIVATE ReSolve)

# Create CUDA examples
if(RESOLVE_USE_CUDA)

Expand All @@ -40,6 +44,11 @@ if(RESOLVE_USE_CUDA)
#rand solver
add_executable(gmres_cusparse_rand.exe r_randGMRES_CUDA.cpp)
target_link_libraries(gmres_cusparse_rand.exe PRIVATE ReSolve)

# Build example with a configurable system solver
add_executable(system_cuda.exe r_SysSolverCuda.cpp)
target_link_libraries(system_cuda.exe PRIVATE ReSolve)

endif(RESOLVE_USE_CUDA)

# Create HIP examples
Expand All @@ -60,10 +69,17 @@ if(RESOLVE_USE_HIP)
# Rand GMRES test with rocsparse
add_executable(gmres_rocsparse_rand.exe r_randGMRES.cpp)
target_link_libraries(gmres_rocsparse_rand.exe PRIVATE ReSolve)
# Build example with a configurable system solver
add_executable(system_hip.exe r_SysSolverHip.cpp)
target_link_libraries(system_hip.exe PRIVATE ReSolve)

# Build example with a configurable system solver
add_executable(system_hip_fgmres.exe r_SysSolverHipRefine.cpp)
target_link_libraries(system_hip_fgmres.exe PRIVATE ReSolve)
endif(RESOLVE_USE_HIP)

# Install all examples in bin directory
set(installable_executables klu_klu.exe klu_klu_standalone.exe)
set(installable_executables klu_klu.exe klu_klu_standalone.exe system.exe)

if(RESOLVE_USE_CUDA)
set(installable_executables ${installable_executables} klu_glu.exe klu_rf.exe klu_rf_fgmres.exe klu_glu_values_update.exe gmres_cusparse_rand.exe)
Expand Down
13 changes: 4 additions & 9 deletions examples/r_KLU_GLU.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -115,9 +115,6 @@ int main(int argc, char *argv[])
}
std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<<std::endl;
//Now call direct solver
if (i == 0) {
KLU->setupParameters(1, 0.1, false);
}
int status;
if (i < 1) {
KLU->setup(A);
Expand All @@ -133,25 +130,23 @@ int main(int argc, char *argv[])
GLU->setup(A, L, U, P, Q);
status = GLU->solve(vec_rhs, vec_x);
std::cout<<"GLU solve status: "<<status<<std::endl;
// status = KLU->solve(vec_rhs, vec_x);
// std::cout<<"KLU solve status: "<<status<<std::endl;
} else {
//status = KLU->refactorize();
std::cout<<"Using CUSOLVER GLU"<<std::endl;
status = GLU->refactorize();
std::cout<<"CUSOLVER GLU refactorization status: "<<status<<std::endl;
status = GLU->solve(vec_rhs, vec_x);
std::cout<<"CUSOLVER GLU solve status: "<<status<<std::endl;
}
vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);


// Estimate solution error
vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
real_type bnorm = sqrt(vector_handler->dot(vec_r, vec_r, "cuda"));
matrix_handler->setValuesChanged(true, "cuda");
matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda");

std::cout << "\t 2-Norm of the residual: "
<< std::scientific << std::setprecision(16)
<< sqrt(vector_handler->dot(vec_r, vec_r, "cuda")) << "\n";
<< sqrt(vector_handler->dot(vec_r, vec_r, "cuda"))/bnorm << "\n";

} // for (int i = 0; i < numSystems; ++i)

Expand Down
3 changes: 0 additions & 3 deletions examples/r_KLU_GLU_matrix_values_update.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -125,9 +125,6 @@ int main(int argc, char *argv[])
}
std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<<std::endl;
//Now call direct solver
if (i == 0) {
KLU->setupParameters(1, 0.1, false);
}
int status;
if (i < 1){
KLU->setup(A);
Expand Down
3 changes: 0 additions & 3 deletions examples/r_KLU_KLU.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -116,9 +116,6 @@ int main(int argc, char *argv[])
}
std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<<std::endl;
//Now call direct solver
if (i == 0) {
KLU->setupParameters(1, 0.1, false);
}
int status;
if (i < 2){
KLU->setup(A);
Expand Down
1 change: 0 additions & 1 deletion examples/r_KLU_KLU_standalone.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,6 @@ int main(int argc, char *argv[])
vec_rhs->setDataUpdated(ReSolve::memory::HOST);
std::cout << "COO to CSR completed. Expanded NNZ: " << A->getNnzExpanded() << std::endl;
//Now call direct solver
KLU->setupParameters(1, 0.1, false);
int status;
KLU->setup(A);
status = KLU->analyze();
Expand Down
3 changes: 0 additions & 3 deletions examples/r_KLU_rf.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -115,9 +115,6 @@ int main(int argc, char *argv[] )
}
std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<<std::endl;
//Now call direct solver
if (i == 0) {
KLU->setupParameters(1, 0.1, false);
}
int status;
if (i < 2){
KLU->setup(A);
Expand Down
3 changes: 0 additions & 3 deletions examples/r_KLU_rf_FGMRES.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -119,9 +119,6 @@ int main(int argc, char *argv[])
}
std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<<std::endl;
//Now call direct solver
if (i == 0) {
KLU->setupParameters(1, 0.1, false);
}
int status;
real_type norm_b;
if (i < 2){
Expand Down
3 changes: 0 additions & 3 deletions examples/r_KLU_rf_FGMRES_reuse_factorization.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -121,9 +121,6 @@ int main(int argc, char *argv[])
}
std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<<std::endl;
//Now call direct solver
if (i == 0) {
KLU->setupParameters(1, 0.1, false);
}
int status;
real_type norm_b;
if (i < 2){
Expand Down
26 changes: 11 additions & 15 deletions examples/r_KLU_rocSolverRf_FGMRES.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -118,10 +118,6 @@ int main(int argc, char *argv[])
vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
}
std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<<std::endl;
//Now call direct solver
if (i == 0) {
KLU->setupParameters(1, 0.1, false);
}
int status;
real_type norm_b;
if (i < 2){
Expand Down Expand Up @@ -175,17 +171,17 @@ int main(int argc, char *argv[])
<< rnrm/norm_b << "\n";

vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
if(!std::isnan(rnrm) && !std::isinf(rnrm)) {
FGMRES->solve(vec_rhs, vec_x);

std::cout << "FGMRES: init nrm: "
<< std::scientific << std::setprecision(16)
<< FGMRES->getInitResidualNorm()/norm_b
<< " final nrm: "
<< FGMRES->getFinalResidualNorm()/norm_b
<< " iter: " << FGMRES->getNumIter() << "\n";
}
}
if(!std::isnan(rnrm) && !std::isinf(rnrm)) {
FGMRES->solve(vec_rhs, vec_x);

std::cout << "FGMRES: init nrm: "
<< std::scientific << std::setprecision(16)
<< FGMRES->getInitResidualNorm()/norm_b
<< " final nrm: "
<< FGMRES->getFinalResidualNorm()/norm_b
<< " iter: " << FGMRES->getNumIter() << "\n";
}
}

} // for (int i = 0; i < numSystems; ++i)

Expand Down
12 changes: 4 additions & 8 deletions examples/r_KLU_rocsolverrf.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -115,11 +115,8 @@ int main(int argc, char *argv[] )
}
std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<<std::endl;
//Now call direct solver
if (i == 0) {
KLU->setupParameters(1, 0.1, false);
}
int status;
if (i < 2){
if (i < 2) {
KLU->setup(A);
status = KLU->analyze();
std::cout<<"KLU analysis status: "<<status<<std::endl;
Expand All @@ -134,7 +131,6 @@ int main(int argc, char *argv[] )
index_type* Q = KLU->getQOrdering();
vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
Rf->setup(A, L, U, P, Q, vec_rhs);
Rf->refactorize();
}
} else {
std::cout<<"Using rocsolver rf"<<std::endl;
Expand All @@ -143,15 +139,15 @@ int main(int argc, char *argv[] )
status = Rf->solve(vec_rhs, vec_x);
std::cout<<"rocsolver rf solve status: "<<status<<std::endl;
}
vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);

vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
real_type bnorm = sqrt(vector_handler->dot(vec_r, vec_r, "hip"));
matrix_handler->setValuesChanged(true, "hip");

matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "hip");

std::cout << "\t 2-Norm of the residual: "
<< std::scientific << std::setprecision(16)
<< sqrt(vector_handler->dot(vec_r, vec_r, "hip")) << "\n";
<< sqrt(vector_handler->dot(vec_r, vec_r, "hip"))/bnorm << "\n";

} // for (int i = 0; i < numSystems; ++i)

Expand Down
3 changes: 0 additions & 3 deletions examples/r_KLU_rocsolverrf_redo_factorization.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -118,9 +118,6 @@ int main(int argc, char *argv[] )
}
std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<<std::endl;
//Now call direct solver
if (i == 0) {
KLU->setupParameters(1, 0.1, false);
}
int status;
if (i < 2){
KLU->setup(A);
Expand Down
160 changes: 160 additions & 0 deletions examples/r_SysSolver.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,160 @@
#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>

#include <resolve/matrix/Coo.hpp>
#include <resolve/matrix/Csr.hpp>
#include <resolve/matrix/Csc.hpp>
#include <resolve/vector/Vector.hpp>
#include <resolve/matrix/io.hpp>
#include <resolve/matrix/MatrixHandler.hpp>
#include <resolve/vector/VectorHandler.hpp>
#include <resolve/LinSolverDirectKLU.hpp>
#include <resolve/workspace/LinAlgWorkspace.hpp>
#include <resolve/SystemSolver.hpp>

using namespace ReSolve::constants;

int main(int argc, char *argv[])
{
// Use the same data types as those you specified in ReSolve build.
using index_type = ReSolve::index_type;
using real_type = ReSolve::real_type;
using vector_type = ReSolve::vector::Vector;

(void) argc; // TODO: Check if the number of input parameters is correct.
std::string matrixFileName = argv[1];
std::string rhsFileName = argv[2];

index_type numSystems = atoi(argv[3]);
std::cout<<"Family mtx file name: "<< matrixFileName << ", total number of matrices: "<<numSystems<<std::endl;
std::cout<<"Family rhs file name: "<< rhsFileName << ", total number of RHSes: " << numSystems<<std::endl;

std::string fileId;
std::string rhsId;
std::string matrixFileNameFull;
std::string rhsFileNameFull;

ReSolve::matrix::Coo* A_coo;
ReSolve::matrix::Csr* A;
ReSolve::LinAlgWorkspaceCpu* workspace = new ReSolve::LinAlgWorkspaceCpu();
ReSolve::MatrixHandler* matrix_handler = new ReSolve::MatrixHandler(workspace);
ReSolve::VectorHandler* vector_handler = new ReSolve::VectorHandler(workspace);
real_type* rhs = nullptr;
real_type* x = nullptr;

vector_type* vec_rhs;
vector_type* vec_x;
vector_type* vec_r;

ReSolve::SystemSolver* solver = new ReSolve::SystemSolver();

for (int i = 0; i < numSystems; ++i)
{
index_type j = 4 + i * 2;
fileId = argv[j];
rhsId = argv[j + 1];

matrixFileNameFull = "";
rhsFileNameFull = "";

// Read matrix first
matrixFileNameFull = matrixFileName + fileId + ".mtx";
rhsFileNameFull = rhsFileName + rhsId + ".mtx";
std::cout << std::endl << std::endl << std::endl;
std::cout << "========================================================================================================================"<<std::endl;
std::cout << "Reading: " << matrixFileNameFull << std::endl;
std::cout << "========================================================================================================================"<<std::endl;
std::cout << std::endl;
// Read first matrix
std::ifstream mat_file(matrixFileNameFull);
if(!mat_file.is_open())
{
std::cout << "Failed to open file " << matrixFileNameFull << "\n";
return -1;
}
std::ifstream rhs_file(rhsFileNameFull);
if(!rhs_file.is_open())
{
std::cout << "Failed to open file " << rhsFileNameFull << "\n";
return -1;
}
if (i == 0) {
A_coo = ReSolve::io::readMatrixFromFile(mat_file);
A = new ReSolve::matrix::Csr(A_coo->getNumRows(),
A_coo->getNumColumns(),
A_coo->getNnz(),
A_coo->symmetric(),
A_coo->expanded());

rhs = ReSolve::io::readRhsFromFile(rhs_file);
x = new real_type[A->getNumRows()];
vec_rhs = new vector_type(A->getNumRows());
vec_x = new vector_type(A->getNumRows());
vec_r = new vector_type(A->getNumRows());
} else {
ReSolve::io::readAndUpdateMatrix(mat_file, A_coo);
ReSolve::io::readAndUpdateRhs(rhs_file, &rhs);
}
std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows()
<< " x " << A->getNumColumns()
<< ", nnz: " << A->getNnz()
<< ", symmetric? " << A->symmetric()
<< ", Expanded? " << A->expanded() << std::endl;
mat_file.close();
rhs_file.close();

//Now convert to CSR.
if (i < 2) {
matrix_handler->coo2csr(A_coo, A, "cpu");
vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST);
vec_rhs->setDataUpdated(ReSolve::memory::HOST);
} else {
matrix_handler->coo2csr(A_coo, A, "cpu");
vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST);
}
std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<<std::endl;
//Now call direct solver
solver->setMatrix(A);
// if (i == 0) {
// solver->setupParameters(1, 0.1, false);
// }
int status;
if (i < 2){
// solver->setup(A);
status = solver->analyze();
std::cout<<"solver analysis status: "<<status<<std::endl;
status = solver->factorize();
std::cout<<"solver factorization status: "<<status<<std::endl;
status = solver->solve(vec_rhs, vec_x);
std::cout<<"solver solve status: "<<status<<std::endl;
} else {
status = solver->refactorize();
std::cout<<"solver re-factorization status: "<<status<<std::endl;
status = solver->solve(vec_rhs, vec_x);
std::cout<<"solver solve status: "<<status<<std::endl;
}
vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST);

matrix_handler->setValuesChanged(true, "cpu");

matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "cpu");

std::cout << "\t 2-Norm of the residual: "
<< std::scientific << std::setprecision(16)
<< sqrt(vector_handler->dot(vec_r, vec_r, "cpu")) << "\n";
}

//now DELETE
delete A;
delete solver;
delete [] x;
delete [] rhs;
delete vec_r;
delete vec_x;
delete matrix_handler;
delete vector_handler;

return 0;
}
Loading

0 comments on commit 5129c8c

Please sign in to comment.