Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80391815
sparse_solver_eigen.cc
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
Sat, Aug 31, 12:02
Size
2 KB
Mime Type
text/x-c
Expires
Mon, Sep 2, 12:02 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20359470
Attached To
rAKA akantu
sparse_solver_eigen.cc
View Options
#include "sparse_solver_eigen.hh"
#include "dof_manager_default.hh"
#include "solver_vector_default.hh"
#include "sparse_matrix_aij.hh"
namespace
akantu
{
/* -------------------------------------------------------------------------- */
SparseSolverEigen
::
SparseSolverEigen
(
DOFManagerDefault
&
dof_manager
,
const
ID
&
matrix_id
,
const
ID
&
id
)
:
SparseSolver
(
dof_manager
,
matrix_id
,
id
),
dof_manager
(
dof_manager
)
{
AKANTU_DEBUG_ASSERT
(
communicator
.
getNbProc
()
==
1
,
"This solver does not work in parallel"
);
}
/* -------------------------------------------------------------------------- */
void
SparseSolverEigen
::
setA
()
{
const
auto
&
Aij
=
aka
::
as_type
<
SparseMatrixAIJ
>
(
dof_manager
.
getMatrix
(
matrix_id
));
if
(
this
->
last_value_release
==
Aij
.
getRelease
()
and
this
->
last_value_release
!=
-
1
)
{
return
;
}
using
Triplet
=
Eigen
::
Triplet
<
Real
>
;
std
::
vector
<
Triplet
>
triplets
;
triplets
.
reserve
(
Aij
.
getNbNonZero
());
for
(
auto
&&
[
i
,
j
,
a
]
:
zip
(
Aij
.
getIRN
(),
Aij
.
getJCN
(),
Aij
.
getA
()))
{
triplets
.
emplace_back
(
i
-
1
,
j
-
1
,
a
);
if
(
Aij
.
getMatrixType
()
==
_symmetric
and
i
!=
j
)
{
triplets
.
emplace_back
(
j
-
1
,
i
-
1
,
a
);
}
}
A
.
resize
(
Aij
.
size
(),
Aij
.
size
());
A
.
setZero
();
A
.
setFromTriplets
(
triplets
.
begin
(),
triplets
.
end
());
}
/* -------------------------------------------------------------------------- */
void
SparseSolverEigen
::
solve
()
{
const
auto
&
Aij
=
aka
::
as_type
<
SparseMatrixAIJ
>
(
dof_manager
.
getMatrix
(
matrix_id
));
auto
&
b_
=
aka
::
as_type
<
SolverVectorDefault
>
(
this
->
dof_manager
.
getResidual
())
.
getGlobalVector
();
Array
<
Real
>
x_
(
b_
.
size
()
*
b_
.
getNbComponent
());
VectorProxy
<
Real
>
x
(
x_
.
data
(),
x_
.
size
()
*
x_
.
getNbComponent
());
VectorProxy
<
Real
>
b
(
b_
.
data
(),
b_
.
size
()
*
b_
.
getNbComponent
());
setA
();
if
(
this
->
last_profile_release
!=
Aij
.
getProfileRelease
())
{
solver
.
analyzePattern
(
A
);
this
->
last_profile_release
=
Aij
.
getProfileRelease
();
}
if
(
this
->
last_value_release
!=
Aij
.
getValueRelease
())
{
solver
.
factorize
(
A
);
this
->
last_value_release
=
Aij
.
getValueRelease
();
}
x
=
solver
.
solve
(
b
);
aka
::
as_type
<
SolverVectorDefault
>
(
this
->
dof_manager
.
getSolution
())
.
setGlobalVector
(
x_
);
this
->
dof_manager
.
splitSolutionPerDOFs
();
}
}
// namespace akantu
Event Timeline
Log In to Comment