Commit 7971041f authored by Niklas Eiling's avatar Niklas Eiling
Browse files

add error handling and logging to MnaSolverGpuSparse


Signed-off-by: Niklas Eiling's avatarNiklas Eiling <niklas.eiling@eonerc.rwth-aachen.de>
parent 2ffe6924
......@@ -99,10 +99,10 @@ build:linux-profiling:
- cmake -DWITH_PROFILING=ON -DWITH_ASAN=ON -DWITH_CUDA=OFF -DWITH_SPDLOG_SUBMODULE=ON ..
- make -j 32
image: ${DOCKER_IMAGE_DEV}-centos:${DOCKER_TAG}
cache:
paths:
- build
key: build-linux-profiling
#cache:
# paths:
# - build
# key: build-linux-profiling
artifacts:
paths:
- build
......
......@@ -26,6 +26,8 @@ namespace DPsim {
/// Intermediate Vector
cuda::Vector<double> mGpuIntermediateVec;
using Solver::mSLog;
/// Initialize cuSparse-library
void initialize();
/// ILU factorization
......@@ -41,6 +43,8 @@ namespace DPsim {
csrsv2Info_t info_L = nullptr;
csrsv2Info_t info_U = nullptr;
void checkCusparseStatus(cusparseStatus_t status, std::string additionalInfo="cuSparse Error:");
public:
MnaSolverGpuSparse(String name,
CPS::Domain domain = CPS::Domain::DP,
......
#include "dpsim/Definitions.h"
#include <dpsim/MNASolverGpuSparse.h>
#include <dpsim/SequentialScheduler.h>
#include <Eigen/Eigen>
......@@ -16,21 +17,45 @@ MnaSolverGpuSparse<VarType>::MnaSolverGpuSparse(String name,
mGpuRhsVec(0), mGpuLhsVec(0), mGpuIntermediateVec(0),
pBuffer(0) {
//TODO Error-Handling
cusparseCreate(&mCusparsehandle);
//TODO: Error-Handling
cusolverSpCreate(&mCusolverhandle);
}
template <typename VarType>
MnaSolverGpuSparse<VarType>::~MnaSolverGpuSparse() {
//TODO Deconstructor
if (mCusparsehandle != nullptr) {
cusparseDestroyMatDescr(descr_L);
cusparseDestroyMatDescr(descr_U);
cusparseDestroyCsrsv2Info(info_L);
cusparseDestroyCsrsv2Info(info_U);
cusparseDestroy(mCusparsehandle);
}
if (mCusolverhandle != nullptr) {
cusolverSpDestroy(mCusolverhandle);
}
}
template <typename VarType>
inline void MnaSolverGpuSparse<VarType>::checkCusparseStatus(cusparseStatus_t status, std::string additionalInfo)
{
if (status != CUSPARSE_STATUS_SUCCESS) {
mSLog->error("{} {}", additionalInfo, cusparseGetErrorString(status));
throw SolverException();
}
}
template <typename VarType>
void MnaSolverGpuSparse<VarType>::initialize() {
MnaSolver<VarType>::initialize();
cusparseStatus_t csp_status;
cusolverStatus_t cso_status;
csp_status = cusparseCreate(&mCusparsehandle);
checkCusparseStatus(csp_status, "cuSparse initialization failed:");
if ((cso_status = cusolverSpCreate(&mCusolverhandle)) != CUSOLVER_STATUS_SUCCESS) {
mSLog->error("cuSolver initialization failed: Error code {}", cso_status);
throw SolverException();
}
size_t N = this->mRightSideVector.rows();
auto hMat = this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)];
......@@ -53,50 +78,47 @@ void MnaSolverGpuSparse<VarType>::initialize() {
// - matrix U is base-1
// - matrix U is upper triangular
// - matrix U has non-unit diagonal
cusparseCreateMatDescr(&descr_M);
cusparseSetMatIndexBase(descr_M, CUSPARSE_INDEX_BASE_ZERO);
cusparseSetMatType(descr_M, CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseCreateMatDescr(&descr_L);
cusparseSetMatIndexBase(descr_L, CUSPARSE_INDEX_BASE_ZERO);
cusparseSetMatType(descr_L, CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatFillMode(descr_L, CUSPARSE_FILL_MODE_LOWER);
cusparseSetMatDiagType(descr_L, CUSPARSE_DIAG_TYPE_UNIT);
cusparseCreateMatDescr(&descr_U);
cusparseSetMatIndexBase(descr_U, CUSPARSE_INDEX_BASE_ZERO);
cusparseSetMatType(descr_U, CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatFillMode(descr_U, CUSPARSE_FILL_MODE_UPPER);
cusparseSetMatDiagType(descr_U, CUSPARSE_DIAG_TYPE_NON_UNIT);
checkCusparseStatus(cusparseCreateMatDescr(&descr_M));
checkCusparseStatus(cusparseSetMatIndexBase(descr_M, CUSPARSE_INDEX_BASE_ZERO));
checkCusparseStatus(cusparseSetMatType(descr_M, CUSPARSE_MATRIX_TYPE_GENERAL));
checkCusparseStatus(cusparseCreateMatDescr(&descr_L));
checkCusparseStatus(cusparseSetMatIndexBase(descr_L, CUSPARSE_INDEX_BASE_ZERO));
checkCusparseStatus(cusparseSetMatType(descr_L, CUSPARSE_MATRIX_TYPE_GENERAL));
checkCusparseStatus(cusparseSetMatFillMode(descr_L, CUSPARSE_FILL_MODE_LOWER));
checkCusparseStatus(cusparseSetMatDiagType(descr_L, CUSPARSE_DIAG_TYPE_UNIT));
checkCusparseStatus(cusparseCreateMatDescr(&descr_U));
checkCusparseStatus(cusparseSetMatIndexBase(descr_U, CUSPARSE_INDEX_BASE_ZERO));
checkCusparseStatus(cusparseSetMatType(descr_U, CUSPARSE_MATRIX_TYPE_GENERAL));
checkCusparseStatus(cusparseSetMatFillMode(descr_U, CUSPARSE_FILL_MODE_UPPER));
checkCusparseStatus(cusparseSetMatDiagType(descr_U, CUSPARSE_DIAG_TYPE_NON_UNIT));
// step 2: create a empty info structure
// we need one info for csrilu02 and two info's for csrsv2
cusparseCreateCsrilu02Info(&info_M);
cusparseCreateCsrsv2Info(&info_L);
cusparseCreateCsrsv2Info(&info_U);
checkCusparseStatus(cusparseCreateCsrilu02Info(&info_M));
checkCusparseStatus(cusparseCreateCsrsv2Info(&info_L));
checkCusparseStatus(cusparseCreateCsrsv2Info(&info_U));
// step 2a: permutate Matrix M' = P*M
int p_nnz = 0;
int p[N];
cusolverStatus_t ret = cusolverSpDcsrzfdHost(
cso_status = cusolverSpDcsrzfdHost(
mCusolverhandle, N,nnz, descr_M,
hMat.valuePtr(),
hMat.outerIndexPtr(),
hMat.innerIndexPtr(),
p, &p_nnz);
if (ret != CUSOLVER_STATUS_SUCCESS) {
this->mSLog->error("cusolverSpDcsrzfdHost returend an error");
if (cso_status != CUSOLVER_STATUS_SUCCESS) {
mSLog->error("cusolverSpDcsrzfdHost returend an error");
}
// create Eigen::PermutationMatrix from the p
mTransp = std::unique_ptr<Eigen::PermutationMatrix<Eigen::Dynamic> >(
new Eigen::PermutationMatrix<Eigen::Dynamic>(
Eigen::Map< Eigen::Matrix<int, Eigen::Dynamic, 1> >(p, N, 1)));
//std::cout << "Transposition:" << std::endl;
//std::cout << mTransp->indices() << std::endl;
// apply permutation
hMat = *mTransp * hMat;
......@@ -111,9 +133,20 @@ void MnaSolverGpuSparse<VarType>::initialize() {
// step 3: query how much memory used in csrilu02 and csrsv2, and allocate the buffer
int pBufferSize_M, pBufferSize_L, pBufferSize_U;
int pBufferSize;
cusparseDcsrilu02_bufferSize(mCusparsehandle, N, nnz, descr_M, d_csrVal, d_csrRowPtr, d_csrColInd, info_M, &pBufferSize_M);
cusparseDcsrsv2_bufferSize(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, nnz, descr_L, d_csrVal, d_csrRowPtr, d_csrColInd, info_L, &pBufferSize_L);
cusparseDcsrsv2_bufferSize(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, nnz, descr_U, d_csrVal, d_csrRowPtr, d_csrColInd, info_U, &pBufferSize_U);
csp_status = cusparseDcsrilu02_bufferSize(mCusparsehandle, N, nnz, descr_M,
d_csrVal, d_csrRowPtr, d_csrColInd,
info_M, &pBufferSize_M);
checkCusparseStatus(csp_status, "failed to get cusparse bufferSize:");
csp_status = cusparseDcsrsv2_bufferSize(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE,
N, nnz, descr_L, d_csrVal, d_csrRowPtr,
d_csrColInd, info_L, &pBufferSize_L);
checkCusparseStatus(csp_status, "failed to get cusparse bufferSize:");
csp_status = cusparseDcsrsv2_bufferSize(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE,
N, nnz, descr_U, d_csrVal, d_csrRowPtr,
d_csrColInd, info_U, &pBufferSize_U);
checkCusparseStatus(csp_status, "failed to get cusparse bufferSize:");
// Buffer
pBufferSize = std::max({pBufferSize_M, pBufferSize_L, pBufferSize_U});
......@@ -124,23 +157,51 @@ void MnaSolverGpuSparse<VarType>::initialize() {
// perform analysis of triangular solve on U
// The lower(upper) triangular part of M has the same sparsity pattern as L(U),
// we can do analysis of csrilu0 and csrsv2 simultaneously.
cusparseDcsrilu02_analysis(mCusparsehandle, N, nnz, descr_M, d_csrVal, d_csrRowPtr, d_csrColInd, info_M, CUSPARSE_SOLVE_POLICY_NO_LEVEL, pBuffer.data());
auto status = cusparseXcsrilu02_zeroPivot(mCusparsehandle, info_M, &structural_zero);
if (CUSPARSE_STATUS_ZERO_PIVOT == status){
std::cout << "A(" << structural_zero << ',' << structural_zero << ") is missing\n" << std::endl;
csp_status = cusparseDcsrilu02_analysis(mCusparsehandle, N, nnz, descr_M,
d_csrVal, d_csrRowPtr, d_csrColInd,
info_M, CUSPARSE_SOLVE_POLICY_NO_LEVEL,
pBuffer.data());
checkCusparseStatus(csp_status, "failed to analyse cusparse problem:");
csp_status = cusparseXcsrilu02_zeroPivot(mCusparsehandle, info_M, &structural_zero);
if (csp_status == CUSPARSE_STATUS_ZERO_PIVOT){
mSLog->error("A({},{}) is missing", structural_zero, structural_zero);
checkCusparseStatus(csp_status);
throw SolverException();
}
cusparseDcsrsv2_analysis(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, nnz, descr_L, d_csrVal, d_csrRowPtr, d_csrColInd, info_L, CUSPARSE_SOLVE_POLICY_NO_LEVEL, pBuffer.data());
cusparseDcsrsv2_analysis(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, nnz, descr_U, d_csrVal, d_csrRowPtr, d_csrColInd, info_U, CUSPARSE_SOLVE_POLICY_USE_LEVEL, pBuffer.data());
csp_status = cusparseDcsrsv2_analysis(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE,
N, nnz, descr_L, d_csrVal, d_csrRowPtr,
d_csrColInd, info_L, CUSPARSE_SOLVE_POLICY_NO_LEVEL,
pBuffer.data());
checkCusparseStatus(csp_status, "failed to analyse cusparse problem:");
csp_status = cusparseDcsrsv2_analysis(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE,
N, nnz, descr_U, d_csrVal, d_csrRowPtr,
d_csrColInd, info_U, CUSPARSE_SOLVE_POLICY_USE_LEVEL,
pBuffer.data());
checkCusparseStatus(csp_status, "failed to analyse cusparse problem:");
// step 5: M = L * U
cusparseDcsrilu02(mCusparsehandle, N, nnz, descr_M, d_csrVal, d_csrRowPtr, d_csrColInd, info_M, CUSPARSE_SOLVE_POLICY_NO_LEVEL, pBuffer.data());
status = cusparseXcsrilu02_zeroPivot(mCusparsehandle, info_M, &numerical_zero);
if (CUSPARSE_STATUS_ZERO_PIVOT == status){
printf("U(%d,%d) is zero\n", numerical_zero, numerical_zero);
csp_status = cusparseDcsrilu02(mCusparsehandle, N, nnz, descr_M, d_csrVal,
d_csrRowPtr, d_csrColInd, info_M, CUSPARSE_SOLVE_POLICY_NO_LEVEL,
pBuffer.data());
checkCusparseStatus(csp_status, "failed to perform cusparse ILU:");
csp_status = cusparseXcsrilu02_zeroPivot(mCusparsehandle, info_M, &numerical_zero);
if (csp_status == CUSPARSE_STATUS_ZERO_PIVOT){
mSLog->error("U({},{}) is zero\n", numerical_zero, numerical_zero);
checkCusparseStatus(csp_status);
throw SolverException();
}
//std::cout << *mSysMat.get() << std::endl;
csp_status = cusparseDestroyCsrilu02Info(info_M);
checkCusparseStatus(csp_status, "failed to destroy info:");
csp_status = cusparseDestroyMatDescr(descr_M);
checkCusparseStatus(csp_status, "failed to destroy MatDescr:");
}
template <typename VarType>
......@@ -169,6 +230,8 @@ Task::List MnaSolverGpuSparse<VarType>::getTasks() {
template <typename VarType>
void MnaSolverGpuSparse<VarType>::solve(Real time, Int timeStepCount) {
cudaError_t status;
cusparseStatus_t csp_status;
int size = this->mRightSideVector.rows();
// Reset source vector
this->mRightSideVector.setZero();
......@@ -181,22 +244,38 @@ void MnaSolverGpuSparse<VarType>::solve(Real time, Int timeStepCount) {
//Copy right vector to device
//Permutate right side: R' = P * R
this->mRightSideVector = *mTransp * this->mRightSideVector;
cudaMemcpy(mGpuRhsVec.data(), &this->mRightSideVector(0), size * sizeof(Real), cudaMemcpyHostToDevice);
status = cudaMemcpy(mGpuRhsVec.data(), &this->mRightSideVector(0), size * sizeof(Real), cudaMemcpyHostToDevice);
if (status != cudaSuccess) {
mSLog->error("Cuda Error: {}", cudaGetErrorString(status));
throw SolverException();
}
const double alpha = 1.;
// Solve
// step 6: solve L*z = x
cusparseDcsrsv2_solve(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE, size, mSysMat->non_zero, &alpha, descr_L,
mSysMat->val.data(), mSysMat->row.data(), mSysMat->col.data(), info_L,
mGpuRhsVec.data(), mGpuIntermediateVec.data(), CUSPARSE_SOLVE_POLICY_NO_LEVEL, pBuffer.data());
csp_status = cusparseDcsrsv2_solve(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE,
size, mSysMat->non_zero, &alpha, descr_L,
mSysMat->val.data(), mSysMat->row.data(),
mSysMat->col.data(), info_L, mGpuRhsVec.data(),
mGpuIntermediateVec.data(), CUSPARSE_SOLVE_POLICY_NO_LEVEL,
pBuffer.data());
checkCusparseStatus(csp_status, "failed to solve L*z=x:");
// step 7: solve U*y = z
cusparseDcsrsv2_solve(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE, size, mSysMat->non_zero, &alpha, descr_U,
mSysMat->val.data(), mSysMat->row.data(), mSysMat->col.data(), info_U,
mGpuIntermediateVec.data(), mGpuLhsVec.data(), CUSPARSE_SOLVE_POLICY_USE_LEVEL, pBuffer.data());
csp_status = cusparseDcsrsv2_solve(mCusparsehandle, CUSPARSE_OPERATION_NON_TRANSPOSE,
size, mSysMat->non_zero, &alpha, descr_U,
mSysMat->val.data(), mSysMat->row.data(),
mSysMat->col.data(), info_U, mGpuIntermediateVec.data(),
mGpuLhsVec.data(), CUSPARSE_SOLVE_POLICY_USE_LEVEL,
pBuffer.data());
checkCusparseStatus(csp_status, "failed to solve U*y=z:");
//Copy Solution back
cudaMemcpy(&this->mLeftSideVector(0), mGpuLhsVec.data(), size * sizeof(Real), cudaMemcpyDeviceToHost);
status = cudaMemcpy(&this->mLeftSideVector(0), mGpuLhsVec.data(), size * sizeof(Real), cudaMemcpyDeviceToHost);
if (status != cudaSuccess) {
mSLog->error("Cuda Error: {}", cudaGetErrorString(status));
throw SolverException();
}
//Apply inverse Permutation L = P' * L'
this->mLeftSideVector = mTransp->inverse() * this->mLeftSideVector;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment