Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79499882
sparse_solvers.cpp
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
Mon, Aug 26, 18:25
Size
9 KB
Mime Type
text/x-c
Expires
Wed, Aug 28, 18:25 (2 d)
Engine
blob
Format
Raw Data
Handle
20194990
Attached To
rSPECMICP SpecMiCP / ReactMiCP
sparse_solvers.cpp
View Options
#include "catch.hpp"
#include "specmicp_common/sparse_solvers/sparse_solver.hpp"
#include <vector>
using
namespace
specmicp
;
constexpr
int
size
=
10
;
using
MatrixT
=
Eigen
::
SparseMatrix
<
scalar_t
>
;
// simple diffusion problem
// BC :
// left : fixed value
// right : no flow
void
fill_residuals
(
Vector
&
x
,
Vector
&
residuals
)
{
residuals
.
resize
(
size
);
for
(
int
i
=
0
;
i
<
size
-
1
;
++
i
)
{
residuals
(
i
)
=
(
-
2
*
x
(
i
+
1
)
+
x
(
i
)
+
x
(
i
+
2
));
}
residuals
(
9
)
=
x
(
9
)
-
x
(
10
);
}
void
update_solution
(
Vector
&
xp
,
const
Vector
&
x
,
const
Vector
solution
)
{
for
(
int
i
=
0
;
i
<
size
;
++
i
)
{
xp
(
i
+
1
)
=
x
(
i
+
1
)
+
solution
(
i
);
}
}
void
fill_jacobian
(
MatrixT
&
jacobian
)
{
using
Triplet
=
Eigen
::
Triplet
<
scalar_t
,
index_t
>
;
std
::
vector
<
Triplet
>
list_triplets
;
list_triplets
.
emplace_back
(
Triplet
(
0
,
0
,
-
2
));
list_triplets
.
emplace_back
(
Triplet
(
9
,
9
,
-
1
));
for
(
auto
i
=
1
;
i
<
size
-
1
;
++
i
)
{
list_triplets
.
emplace_back
(
Triplet
(
i
,
i
,
-
2
));
}
for
(
auto
i
=
0
;
i
<
size
-
1
;
++
i
)
{
list_triplets
.
emplace_back
(
Triplet
(
i
,
i
+
1
,
1
));
list_triplets
.
emplace_back
(
Triplet
(
i
+
1
,
i
,
1
));
}
jacobian
=
MatrixT
(
size
,
size
);
jacobian
.
setFromTriplets
(
list_triplets
.
begin
(),
list_triplets
.
end
());
}
TEST_CASE
(
"Sparse solvers"
,
"[Eigen],[Sparse],[Solver]"
)
{
Vector
x
(
size
+
1
);
x
<<
1
,
2
,
2
,
2
,
3
,
2
,
2
,
2
,
2
,
2
,
2
;
Vector
scaling
(
size
);
scaling
.
setConstant
(
1.0e3
);
SECTION
(
"LU solver"
,
"[LU]"
)
{
Vector
residuals
;
fill_residuals
(
x
,
residuals
);
MatrixT
jacobian
;
fill_jacobian
(
jacobian
);
auto
lu_solver
=
sparse_solvers
::
get_sparse_solver
<
MatrixT
,
Vector
,
Vector
>
(
sparse_solvers
::
SparseSolver
::
SparseLU
);
lu_solver
->
analyse_pattern
(
jacobian
);
auto
retcode
=
lu_solver
->
decompose
(
jacobian
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
solution
(
size
);
retcode
=
lu_solver
->
solve
(
residuals
,
solution
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
xp
(
size
+
1
);
xp
(
0
)
=
1
;
update_solution
(
xp
,
x
,
solution
);
fill_residuals
(
xp
,
residuals
);
REQUIRE
(
residuals
.
norm
()
<
1e-10
);
}
SECTION
(
"LU solver with scaling"
,
"[LU],[scaling]"
)
{
Vector
residuals
;
fill_residuals
(
x
,
residuals
);
MatrixT
jacobian
;
fill_jacobian
(
jacobian
);
jacobian
=
jacobian
*
scaling
.
asDiagonal
();
jacobian
.
makeCompressed
();
auto
lu_solver
=
sparse_solvers
::
get_sparse_solver
<
MatrixT
,
Vector
,
Vector
>
(
sparse_solvers
::
SparseSolver
::
SparseLU
);
lu_solver
->
analyse_pattern
(
jacobian
);
auto
retcode
=
lu_solver
->
decompose
(
jacobian
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
solution
(
size
);
retcode
=
lu_solver
->
solve_scaling
(
residuals
,
scaling
,
solution
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
xp
(
size
+
1
);
xp
(
0
)
=
1
;
update_solution
(
xp
,
x
,
solution
);
fill_residuals
(
xp
,
residuals
);
REQUIRE
(
residuals
.
norm
()
<
1e-10
);
}
SECTION
(
"QR solver"
,
"[QR]"
)
{
Vector
residuals
;
fill_residuals
(
x
,
residuals
);
MatrixT
jacobian
;
fill_jacobian
(
jacobian
);
auto
qr_solver
=
sparse_solvers
::
get_sparse_solver
<
MatrixT
,
Vector
,
Vector
>
(
sparse_solvers
::
SparseSolver
::
SparseQR
);
qr_solver
->
analyse_pattern
(
jacobian
);
auto
retcode
=
qr_solver
->
decompose
(
jacobian
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
solution
(
size
);
retcode
=
qr_solver
->
solve
(
residuals
,
solution
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
xp
(
size
+
1
);
xp
(
0
)
=
1
;
update_solution
(
xp
,
x
,
solution
);
fill_residuals
(
xp
,
residuals
);
REQUIRE
(
residuals
.
norm
()
<
1e-10
);
}
SECTION
(
"QR solver - pivots failed"
,
"[QR]"
)
{
Vector
residuals
;
fill_residuals
(
x
,
residuals
);
MatrixT
jacobian
;
fill_jacobian
(
jacobian
);
auto
qr_solver
=
sparse_solvers
::
get_sparse_solver
<
MatrixT
,
Vector
,
Vector
>
(
sparse_solvers
::
SparseSolver
::
SparseQR
,
1e8
);
qr_solver
->
analyse_pattern
(
jacobian
);
auto
retcode
=
qr_solver
->
decompose
(
jacobian
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
FailedDecomposition
);
}
SECTION
(
"QR solver - pivots work"
,
"[QR]"
)
{
Vector
residuals
;
fill_residuals
(
x
,
residuals
);
MatrixT
jacobian
;
fill_jacobian
(
jacobian
);
auto
qr_solver
=
sparse_solvers
::
get_sparse_solver
<
MatrixT
,
Vector
,
Vector
>
(
sparse_solvers
::
SparseSolver
::
SparseQR
,
1e-8
);
qr_solver
->
analyse_pattern
(
jacobian
);
auto
retcode
=
qr_solver
->
decompose
(
jacobian
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
}
SECTION
(
"QR solver with scaling"
,
"[QR],[scaling]"
)
{
Vector
residuals
;
fill_residuals
(
x
,
residuals
);
MatrixT
jacobian
;
fill_jacobian
(
jacobian
);
jacobian
=
jacobian
*
scaling
.
asDiagonal
();
jacobian
.
makeCompressed
();
auto
qr_solver
=
sparse_solvers
::
get_sparse_solver
<
MatrixT
,
Vector
,
Vector
>
(
sparse_solvers
::
SparseSolver
::
SparseQR
);
qr_solver
->
analyse_pattern
(
jacobian
);
auto
retcode
=
qr_solver
->
decompose
(
jacobian
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
solution
(
size
);
retcode
=
qr_solver
->
solve_scaling
(
residuals
,
scaling
,
solution
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
xp
(
size
+
1
);
xp
(
0
)
=
1
;
update_solution
(
xp
,
x
,
solution
);
fill_residuals
(
xp
,
residuals
);
REQUIRE
(
residuals
.
norm
()
<
1e-10
);
}
SECTION
(
"BiCGSTAB solver"
,
"[BiCGSTAB]"
)
{
Vector
residuals
;
fill_residuals
(
x
,
residuals
);
MatrixT
jacobian
;
fill_jacobian
(
jacobian
);
auto
bicgstab_solver
=
sparse_solvers
::
get_sparse_solver
<
MatrixT
,
Vector
,
Vector
>
(
sparse_solvers
::
SparseSolver
::
BiCGSTAB
);
bicgstab_solver
->
analyse_pattern
(
jacobian
);
auto
retcode
=
bicgstab_solver
->
decompose
(
jacobian
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
solution
(
size
);
retcode
=
bicgstab_solver
->
solve
(
residuals
,
solution
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
xp
(
size
+
1
);
xp
(
0
)
=
1
;
update_solution
(
xp
,
x
,
solution
);
fill_residuals
(
xp
,
residuals
);
REQUIRE
(
residuals
.
norm
()
<
1e-10
);
}
SECTION
(
"BiCGSTAB solver with scaling"
,
"[BiCGSTAB], [scaling]"
)
{
Vector
residuals
;
fill_residuals
(
x
,
residuals
);
MatrixT
jacobian
;
fill_jacobian
(
jacobian
);
jacobian
=
jacobian
*
scaling
.
asDiagonal
();
jacobian
.
makeCompressed
();
auto
bicgstab_solver
=
sparse_solvers
::
get_sparse_solver
<
MatrixT
,
Vector
,
Vector
>
(
sparse_solvers
::
SparseSolver
::
BiCGSTAB
);
bicgstab_solver
->
analyse_pattern
(
jacobian
);
auto
retcode
=
bicgstab_solver
->
decompose
(
jacobian
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
solution
(
size
);
retcode
=
bicgstab_solver
->
solve_scaling
(
residuals
,
scaling
,
solution
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
xp
(
size
+
1
);
xp
(
0
)
=
1
;
update_solution
(
xp
,
x
,
solution
);
fill_residuals
(
xp
,
residuals
);
REQUIRE
(
residuals
.
norm
()
<
1e-10
);
}
#ifdef EIGEN_UNSUPPORTED_FOUND
SECTION
(
"GMRES solver"
,
"[GMRES]"
)
{
Vector
residuals
;
fill_residuals
(
x
,
residuals
);
MatrixT
jacobian
;
fill_jacobian
(
jacobian
);
auto
gmres_solver
=
sparse_solvers
::
get_sparse_solver
<
MatrixT
,
Vector
,
Vector
>
(
sparse_solvers
::
SparseSolver
::
GMRES
);
gmres_solver
->
analyse_pattern
(
jacobian
);
auto
retcode
=
gmres_solver
->
decompose
(
jacobian
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
solution
(
size
);
retcode
=
gmres_solver
->
solve
(
residuals
,
solution
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
xp
(
size
+
1
);
xp
(
0
)
=
1
;
update_solution
(
xp
,
x
,
solution
);
fill_residuals
(
xp
,
residuals
);
REQUIRE
(
residuals
.
norm
()
<
1e-10
);
}
SECTION
(
"GMRES solver with scaling"
,
"[GMRES],[scaling]"
)
{
Vector
residuals
;
fill_residuals
(
x
,
residuals
);
MatrixT
jacobian
;
fill_jacobian
(
jacobian
);
jacobian
=
jacobian
*
scaling
.
asDiagonal
();
jacobian
.
makeCompressed
();
auto
gmres_solver
=
sparse_solvers
::
get_sparse_solver
<
MatrixT
,
Vector
,
Vector
>
(
sparse_solvers
::
SparseSolver
::
GMRES
);
gmres_solver
->
analyse_pattern
(
jacobian
);
auto
retcode
=
gmres_solver
->
decompose
(
jacobian
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
solution
(
size
);
retcode
=
gmres_solver
->
solve_scaling
(
residuals
,
scaling
,
solution
);
REQUIRE
(
retcode
==
sparse_solvers
::
SparseSolverReturnCode
::
Success
);
Vector
xp
(
size
+
1
);
xp
(
0
)
=
1
;
update_solution
(
xp
,
x
,
solution
);
fill_residuals
(
xp
,
residuals
);
REQUIRE
(
residuals
.
norm
()
<
1e-10
);
}
#endif
}
Event Timeline
Log In to Comment