Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F106896191
sparse_solver.hpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Apr 2, 05:20
Size
6 KB
Mime Type
text/x-c++
Expires
Fri, Apr 4, 05:20 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
25300494
Attached To
rSPECMICP SpecMiCP / ReactMiCP
sparse_solver.hpp
View Options
#ifndef SPECMICP_UTILS_SPARSESOLVER_HPP
#define SPECMICP_UTILS_SPARSESOLVER_HPP
#include <Eigen/Sparse>
#ifdef EIGEN_UNSUPPORTED_FOUND
#include <iostream>
#include <Eigen/IterativeSolvers>
#endif
#include "sparse_solver_structs.hpp"
namespace
specmicp
{
// header
// ======
//! \brief Solve a sparse linear system
template
<
class
DerivedV1
,
class
DerivedV2
,
class
DerivedM
>
SparseSolverReturnCode
solve_sparse_linear_system
(
DerivedV1
&
residual
,
DerivedM
&
jacobian
,
DerivedV2
&
update
,
SparseSolver
solver
);
//! \brief Solve a sparse linear system using the QR decomposition
template
<
class
DerivedV1
,
class
DerivedV2
,
class
DerivedM
>
SparseSolverReturnCode
solve_linear_system_sparseqr
(
DerivedV1
&
residual
,
DerivedM
&
jacobian
,
DerivedV2
&
update
);
//! \brief Solve a sparse linear system using the LU decomposition
template
<
class
DerivedV1
,
class
DerivedV2
,
class
DerivedM
>
SparseSolverReturnCode
solve_linear_system_sparselu
(
DerivedV1
&
residual
,
DerivedM
&
jacobian
,
DerivedV2
&
update
);
//! \brief Solve a sparse linear system using the BiCGSTAB algorithm
template
<
class
DerivedV1
,
class
DerivedV2
,
class
DerivedM
>
SparseSolverReturnCode
solve_linear_system_bicgstab
(
DerivedV1
&
residual
,
DerivedM
&
jacobian
,
DerivedV2
&
update
);
//! \brief Solve a sparse linear system using the BiCGSTAB algorithm
template
<
class
DerivedV1
,
class
DerivedV2
,
class
DerivedM
>
SparseSolverReturnCode
solve_linear_system_gmres
(
DerivedV1
&
residual
,
DerivedM
&
jacobian
,
DerivedV2
&
update
);
//! \brief A sparse system solver
template
<
class
SolverType
>
class
SparseSystemSolver
{
public
:
//! \brief Type of the matrix;
using
MatrixType
=
typename
SolverType
::
MatrixType
;
//! \brief Analyse the sparsity matrix of the jacobian
void
analyse_sparsity
()
{
m_solver
.
analysePattern
(
m_jacobian
);
m_structure_is_known
=
true
;
}
//! \brief Compute the decomposition of the jacobian (if available)
SparseSolverReturnCode
compute
();
//! \brief Solve the problem, using 'rhs' and storing solution in 'solution'
template
<
class
RHSType
,
class
UpdateType
>
SparseSolverReturnCode
solve
(
const
RHSType
&
rhs
,
UpdateType
&
solution
);
//! \brief Return a reference to the jacobian, so it can be updated
MatrixType
&
jacobian
()
{
m_jacobian_has_changed
=
true
;
return
m_jacobian
;
}
//! \brief Const reference to the jacobian
const
MatrixType
&
jacobian
()
const
{
return
m_jacobian
;}
private
:
bool
m_structure_is_known
;
//!< True if the structure of the jacobian is known
bool
m_jacobian_has_changed
;
//!< True if the jacobian has changed
MatrixType
m_jacobian
;
//!< The jacobian
SolverType
m_solver
;
//!< The solver
};
// Implementation
// ==============
template
<
class
SolverType
>
SparseSolverReturnCode
SparseSystemSolver
<
SolverType
>::
compute
()
{
// first compress the matrix to avoid copy
if
(
not
m_jacobian
.
isCompressed
())
m_jacobian
.
makeCompressed
();
if
(
m_structure_is_known
)
m_solver
.
factorize
(
m_jacobian
);
else
m_solver
.
compute
(
m_jacobian
);
if
(
m_solver
.
info
()
!=
Eigen
::
Success
)
return
SparseSolverReturnCode
::
FailedDecomposition
;
return
SparseSolverReturnCode
::
Success
;
}
template
<
class
SolverType
>
template
<
class
RHSType
,
class
UpdateType
>
SparseSolverReturnCode
SparseSystemSolver
<
SolverType
>::
solve
(
const
RHSType
&
rhs
,
UpdateType
&
solution
)
{
if
(
m_jacobian_has_changed
)
compute
();
solution
=
m_solver
.
solve
(
rhs
);
if
(
m_solver
.
info
()
!=
Eigen
::
Success
)
return
SparseSolverReturnCode
::
FailedSystemSolving
;
return
SparseSolverReturnCode
::
Success
;
}
// Implementation
// ==============
template
<
class
DerivedV1
,
class
DerivedV2
,
class
DerivedM
>
SparseSolverReturnCode
solve_sparse_linear_system
(
DerivedV1
&
residual
,
DerivedM
&
jacobian
,
DerivedV2
&
update
,
SparseSolver
solver
)
{
if
(
not
jacobian
.
isCompressed
())
jacobian
.
makeCompressed
();
if
(
solver
==
SparseSolver
::
SparseQR
)
{
return
solve_linear_system_sparseqr
(
residual
,
jacobian
,
update
);
}
else
if
(
solver
==
SparseSolver
::
SparseLU
)
{
return
solve_linear_system_sparselu
(
residual
,
jacobian
,
update
);
#ifdef EIGEN_UNSUPPORTED_FOUND
}
else
if
(
solver
==
SparseSolver
::
GMRES
)
{
return
solve_linear_system_gmres
(
residual
,
jacobian
,
update
);
#endif
}
else
{
return
solve_linear_system_bicgstab
(
residual
,
jacobian
,
update
);
}
}
template
<
class
DerivedV1
,
class
DerivedV2
,
class
DerivedM
>
SparseSolverReturnCode
solve_linear_system_sparseqr
(
DerivedV1
&
residual
,
DerivedM
&
jacobian
,
DerivedV2
&
update
)
{
Eigen
::
SparseQR
<
DerivedM
,
Eigen
::
COLAMDOrdering
<
int
>>
solver
(
jacobian
);
if
(
solver
.
info
()
!=
Eigen
::
Success
)
{
return
SparseSolverReturnCode
::
FailedDecomposition
;
}
update
=
solver
.
solve
(
-
residual
);
if
(
solver
.
info
()
!=
Eigen
::
Success
)
{
return
SparseSolverReturnCode
::
FailedSystemSolving
;
}
return
SparseSolverReturnCode
::
Success
;
}
template
<
class
DerivedV1
,
class
DerivedV2
,
class
DerivedM
>
SparseSolverReturnCode
solve_linear_system_sparselu
(
DerivedV1
&
residual
,
DerivedM
&
jacobian
,
DerivedV2
&
update
)
{
Eigen
::
SparseLU
<
DerivedM
>
solver
(
jacobian
);
if
(
solver
.
info
()
!=
Eigen
::
Success
)
{
return
SparseSolverReturnCode
::
FailedDecomposition
;;
}
update
=
solver
.
solve
(
-
residual
);
if
(
solver
.
info
()
!=
Eigen
::
Success
)
{
return
SparseSolverReturnCode
::
FailedSystemSolving
;
}
return
SparseSolverReturnCode
::
Success
;
}
template
<
class
DerivedV1
,
class
DerivedV2
,
class
DerivedM
>
SparseSolverReturnCode
solve_linear_system_bicgstab
(
DerivedV1
&
residual
,
DerivedM
&
jacobian
,
DerivedV2
&
update
)
{
Eigen
::
BiCGSTAB
<
DerivedM
,
Eigen
::
IncompleteLUT
<
typename
DerivedM
::
Scalar
>>
solver
(
jacobian
);
update
=
solver
.
solve
(
-
residual
);
if
(
solver
.
info
()
!=
Eigen
::
Success
)
{
return
SparseSolverReturnCode
::
FailedSystemSolving
;
}
return
SparseSolverReturnCode
::
Success
;
}
template
<
class
DerivedV1
,
class
DerivedV2
,
class
DerivedM
>
SparseSolverReturnCode
solve_linear_system_gmres
(
DerivedV1
&
residual
,
DerivedM
&
jacobian
,
DerivedV2
&
update
)
{
Eigen
::
GMRES
<
DerivedM
,
Eigen
::
IncompleteLUT
<
typename
DerivedM
::
Scalar
>>
solver
(
jacobian
);
update
=
solver
.
solve
(
-
residual
);
if
(
solver
.
info
()
!=
Eigen
::
Success
)
{
return
SparseSolverReturnCode
::
FailedSystemSolving
;
}
return
SparseSolverReturnCode
::
Success
;
}
}
// end namespace specmicp
#endif
//SPECMICP_UTILS_SPARSESOLVER_HPP
Event Timeline
Log In to Comment