Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79127014
LinearSolver.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
Sat, Aug 24, 21:37
Size
15 KB
Mime Type
text/x-c
Expires
Mon, Aug 26, 21:37 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
20154113
Attached To
rLAMMPS lammps
LinearSolver.cpp
View Options
// Header file for this class
#include "LinearSolver.h"
#include <sstream>
using
std
::
stringstream
;
namespace
ATC
{
const
double
kPenalty
=
1.0e4
;
const
double
kTol
=
1.0e-8
;
const
int
kMaxDirect
=
1000
;
// ====================================================================
// LinearSolver
// ====================================================================
LinearSolver
::
LinearSolver
(
const
SPAR_MAT
&
A
,
const
BC_SET
&
bcs
,
const
int
solverType
,
const
int
constraintHandlerType
,
bool
parallel
)
:
solverType_
(
solverType
),
constraintHandlerType_
(
constraintHandlerType
),
nVariables_
(
0
),
initialized_
(
false
),
initializedMatrix_
(
false
),
initializedInverse_
(
false
),
matrixModified_
(
false
),
allowReinitialization_
(
false
),
homogeneousBCs_
(
false
),
bcs_
(
&
bcs
),
rhs_
(
NULL
),
rhsDense_
(),
b_
(
NULL
),
matrix_
(
A
),
matrixDense_
(),
matrixFreeFree_
(),
matrixFreeFixed_
(),
matrixInverse_
(),
penalty_
(
0
),
maxIterations_
(
0
),
maxRestarts_
(
0
),
tol_
(
0
),
parallel_
(
parallel
)
{
// deep copy
matrixCopy_
=
A
;
matrixSparse_
=
&
matrixCopy_
;
setup
();
}
LinearSolver
::
LinearSolver
(
const
SPAR_MAT
&
A
,
const
int
solverType
,
bool
parallel
)
:
solverType_
(
solverType
),
constraintHandlerType_
(
NO_CONSTRAINTS
),
nVariables_
(
0
),
initialized_
(
false
),
initializedMatrix_
(
true
),
initializedInverse_
(
false
),
matrixModified_
(
false
),
allowReinitialization_
(
false
),
homogeneousBCs_
(
false
),
bcs_
(
NULL
),
// null implies no contraints will be added later
rhs_
(
NULL
),
rhsDense_
(),
b_
(
NULL
),
matrix_
(
A
),
matrixDense_
(),
matrixFreeFree_
(),
matrixFreeFixed_
(),
matrixInverse_
(),
penalty_
(
0
),
maxIterations_
(
0
),
maxRestarts_
(
0
),
tol_
(
0
),
parallel_
(
parallel
)
{
// shallow copy
matrixSparse_
=
&
A
;
setup
();
}
// --------------------------------------------------------------------
// Setup
// --------------------------------------------------------------------
void
LinearSolver
::
setup
(
void
)
{
penalty_
=
kPenalty
;
// relative to matrix diagonal
tol_
=
kTol
;
nVariables_
=
matrix_
.
nRows
();
maxIterations_
=
2
*
nVariables_
;
maxRestarts_
=
nVariables_
;
// switch method based on size
if
(
solverType_
<
0
)
{
if
(
nVariables_
>
kMaxDirect
)
{
solverType_
=
ITERATIVE_SOLVE_SYMMETRIC
;
constraintHandlerType_
=
PENALIZE_CONSTRAINTS
;
}
else
{
solverType_
=
DIRECT_SOLVE
;
}
}
if
(
constraintHandlerType_
<
0
)
{
constraintHandlerType_
=
PENALIZE_CONSTRAINTS
;
if
(
solverType_
==
DIRECT_SOLVE
)
constraintHandlerType_
=
CONDENSE_CONSTRAINTS
;
}
if
(
solverType_
==
DIRECT_SOLVE
&&
constraintHandlerType_
==
CONDENSE_CONSTRAINTS
)
allowReinitialization_
=
true
;
if
(
solverType_
==
ITERATIVE_SOLVE_SYMMETRIC
&&
constraintHandlerType_
==
CONDENSE_CONSTRAINTS
)
throw
ATC_Error
(
"LinearSolver::unimplemented method"
);
}
// --------------------------------------------------------------------
// Initialize
// --------------------------------------------------------------------
void
LinearSolver
::
allow_reinitialization
(
void
)
{
if
(
constraintHandlerType_
==
PENALIZE_CONSTRAINTS
)
{
if
(
matrixModified_
)
throw
ATC_Error
(
"LinearSolver: can't allow reinitialization after matrix has been modified"
);
matrixOriginal_
=
*
matrixSparse_
;
}
allowReinitialization_
=
true
;
}
void
LinearSolver
::
initialize
(
const
BC_SET
*
bcs
)
{
if
(
bcs
)
{
if
(
!
allowReinitialization_
)
throw
ATC_Error
(
"LinearSolver: reinitialization not allowed"
);
//if (! bcs_ ) throw ATC_Error("LinearSolver: adding constraints after constructing without constraints is not allowed");
// shallow --> deep copy
if
(
!
bcs_
)
{
// constraintHandlerType_ == NO_CONSTRAINTS
if
(
matrixModified_
)
{
throw
ATC_Error
(
"LinearSolver: adding constraints after constructing without constraints is not allowed if matrix has been modified"
);
}
else
{
matrixCopy_
=
*
matrixSparse_
;
matrixSparse_
=
&
matrixCopy_
;
constraintHandlerType_
=
-
1
;
setup
();
}
}
bcs_
=
bcs
;
initializedMatrix_
=
false
;
initializedInverse_
=
false
;
if
(
matrixModified_
)
{
matrixCopy_
=
matrixOriginal_
;
matrixSparse_
=
&
matrixCopy_
;
}
}
initialize_matrix
();
initialize_inverse
();
initialize_rhs
();
initialized_
=
true
;
}
// --------------------------------------------------------------------
// initialize_matrix
// --------------------------------------------------------------------
void
LinearSolver
::
initialize_matrix
(
void
)
{
if
(
initializedMatrix_
)
return
;
if
(
constraintHandlerType_
==
PENALIZE_CONSTRAINTS
)
{
add_matrix_penalty
();
}
else
if
(
constraintHandlerType_
==
CONDENSE_CONSTRAINTS
)
{
partition_matrix
();
}
initializedMatrix_
=
true
;
}
// --------------------------------------------------------------------
// initialize_inverse
// --------------------------------------------------------------------
void
LinearSolver
::
initialize_inverse
(
void
)
{
if
(
initializedInverse_
)
return
;
if
(
solverType_
==
ITERATIVE_SOLVE_SYMMETRIC
||
solverType_
==
ITERATIVE_SOLVE
)
{
matrixDiagonal_
=
matrixSparse_
->
diag
();
// preconditioner
}
else
{
// DIRECT_SOLVE
if
(
constraintHandlerType_
==
CONDENSE_CONSTRAINTS
)
{
if
(
num_unknowns
()
>
0
)
{
matrixInverse_
=
inv
(
matrixFreeFree_
);
}
}
else
{
// NO_CONSTRAINTS || PENALIZE_CONSTRAINTS
matrixDense_
=
matrixSparse_
->
dense_copy
();
// need dense for lapack
matrixInverse_
=
inv
(
matrixDense_
);
}
}
initializedInverse_
=
true
;
}
// --------------------------------------------------------------------
// initialize_rhs
// --------------------------------------------------------------------
void
LinearSolver
::
initialize_rhs
(
void
)
{
if
(
!
rhs_
)
return
;
if
(
!
bcs_
)
{
b_
=
rhs_
;
return
;
}
if
(
constraintHandlerType_
==
PENALIZE_CONSTRAINTS
)
{
add_rhs_penalty
();
}
else
if
(
constraintHandlerType_
==
CONDENSE_CONSTRAINTS
)
{
add_rhs_influence
();
}
}
// --------------------------------------------------------------------
// add matrix penalty
// - change matrix for Dirichlet conditions: add penalty
// --------------------------------------------------------------------
void
LinearSolver
::
add_matrix_penalty
(
void
)
{
SPAR_MAT
&
A
=
matrixCopy_
;
penalty_
*=
(
A
.
diag
()).
maxabs
();
BC_SET
::
const_iterator
itr
;
for
(
itr
=
bcs_
->
begin
();
itr
!=
bcs_
->
end
();
itr
++
)
{
int
i
=
itr
->
first
;
A
.
add
(
i
,
i
,
penalty_
);
// modifies matrix
}
A
.
compress
();
matrixModified_
=
true
;
}
// --------------------------------------------------------------------
// partition matrix
// - partition matrix based on Dirichlet constraints
// --------------------------------------------------------------------
void
LinearSolver
::
partition_matrix
(
void
)
{
fixedSet_
.
clear
();
BC_SET
::
const_iterator
itr
;
for
(
itr
=
bcs_
->
begin
();
itr
!=
bcs_
->
end
();
itr
++
)
{
int
i
=
itr
->
first
;
fixedSet_
.
insert
(
i
);
}
freeSet_
.
clear
();
freeGlobalToCondensedMap_
.
clear
();
int
j
=
0
;
// local index
for
(
int
i
=
0
;
i
<
nVariables_
;
i
++
)
{
if
(
fixedSet_
.
find
(
i
)
==
fixedSet_
.
end
()
)
{
freeSet_
.
insert
(
i
);
freeGlobalToCondensedMap_
[
i
]
=
j
++
;
}
}
if
(
matrixDense_
.
nRows
()
==
0
)
matrixDense_
=
matrixSparse_
->
dense_copy
();
DENS_MAT
&
K
=
matrixDense_
;
K
.
row_partition
(
freeSet_
,
matrixFreeFree_
,
matrixFreeFixed_
);
}
// --------------------------------------------------------------------
// add_rhs_penalty
// --------------------------------------------------------------------
void
LinearSolver
::
add_rhs_penalty
()
{
// deep copy
VECTOR
&
b
=
rhsDense_
;
const
VECTOR
&
r
=
*
rhs_
;
int
size
=
r
.
nRows
();
b
.
reset
(
size
);
for
(
int
i
=
0
;
i
<
size
;
i
++
)
{
b
(
i
)
=
r
(
i
);
}
if
(
!
homogeneousBCs_
){
BC_SET
::
const_iterator
itr
;
for
(
itr
=
bcs_
->
begin
();
itr
!=
bcs_
->
end
();
itr
++
)
{
int
i
=
itr
->
first
;
double
v
=
itr
->
second
;
b
(
i
)
+=
penalty_
*
v
;
}
}
b_
=
&
rhsDense_
;
}
// --------------------------------------------------------------------
// add_rhs_influence
// --------------------------------------------------------------------
void
LinearSolver
::
add_rhs_influence
()
{
if
(
!
initializedMatrix_
)
partition_matrix
();
// rhs = rhs + K_free,fixed * x_fixed
int
nbcs
=
bcs_
->
size
();
if
(
nbcs
==
0
)
{
// no bcs to handle
b_
=
rhs_
;
}
else
{
DENS_VEC
&
b
=
rhsDense_
;
if
(
!
homogeneousBCs_
){
DENS_VEC
xFixed
(
nbcs
);
BC_SET
::
const_iterator
itr
;
int
i
=
0
;
for
(
itr
=
bcs_
->
begin
();
itr
!=
bcs_
->
end
();
itr
++
,
i
++
)
{
double
v
=
itr
->
second
;
xFixed
(
i
,
0
)
=
-
v
;
}
b
=
matrixFreeFixed_
*
xFixed
;
// matrix and bcs have same ordering
}
else
{
b
.
reset
(
matrixFreeFixed_
.
nRows
());
}
const
VECTOR
&
r
=
*
rhs_
;
set
<
int
>::
const_iterator
iter
;
int
i
=
0
;
for
(
iter
=
freeSet_
.
begin
();
iter
!=
freeSet_
.
end
();
iter
++
,
i
++
)
{
b
(
i
)
+=
r
(
*
iter
);
}
b_
=
&
rhsDense_
;
}
}
// --------------------------------------------------------------------
// set fixed values
// - {x_i = y_i}
// --------------------------------------------------------------------
void
LinearSolver
::
set_fixed_values
(
VECTOR
&
X
)
{
BC_SET
::
const_iterator
itr
;
for
(
itr
=
bcs_
->
begin
();
itr
!=
bcs_
->
end
();
itr
++
)
{
int
i
=
itr
->
first
;
double
v
=
0
;
if
(
!
homogeneousBCs_
)
v
=
itr
->
second
;
X
(
i
)
=
v
;
}
}
// --------------------------------------------------------------------
// Eigensystem
// --------------------------------------------------------------------
// calls lapack
void
LinearSolver
::
eigen_system
(
DENS_MAT
&
eigenvalues
,
DENS_MAT
&
eigenvectors
,
const
DENS_MAT
*
M
)
/* const */
{
initialize_matrix
();
// no inverse needed
const
DENS_MAT
*
Kp
=
NULL
;
const
DENS_MAT
*
Mp
=
M
;
DENS_MAT
MM
;
DENS_MAT
KM
;
if
(
constraintHandlerType_
==
CONDENSE_CONSTRAINTS
)
{
Kp
=
&
matrixFreeFree_
;
if
(
M
)
{
DENS_MAT
MfreeFixed
;
// not used
M
->
row_partition
(
freeSet_
,
MM
,
MfreeFixed
);
Mp
=
&
MM
;
}
}
else
{
if
(
matrixDense_
.
nRows
()
==
0
)
matrixDense_
=
matrixSparse_
->
dense_copy
();
Kp
=
&
matrixDense_
;
}
if
(
!
M
)
{
MM
.
identity
(
Kp
->
nRows
());
Mp
=
&
MM
;
}
DENS_MAT
eVecs
,
eVals
;
eVecs
=
eigensystem
(
*
Kp
,
*
Mp
,
eVals
);
eigenvalues
.
reset
(
nVariables_
,
1
);
eigenvectors
.
reset
(
nVariables_
,
nVariables_
);
set
<
int
>::
const_iterator
itr
;
for
(
int
i
=
0
;
i
<
Kp
->
nRows
();
i
++
)
{
// ordering is by energy not node
eigenvalues
(
i
,
0
)
=
eVals
(
i
,
0
);
int
j
=
0
;
for
(
itr
=
freeSet_
.
begin
();
itr
!=
freeSet_
.
end
();
itr
++
,
j
++
)
{
int
jj
=
*
itr
;
eigenvectors
(
jj
,
i
)
=
eVecs
(
j
,
i
);
// transpose
}
}
}
// --------------------------------------------------------------------
// solve
// - solves A x = b
// - if a "b" is provided it is used as the new rhs
// --------------------------------------------------------------------
bool
LinearSolver
::
solve
(
VECTOR
&
x
,
const
VECTOR
&
b
)
{
SPAR_MAT
*
A
=
NULL
;
rhs_
=
&
b
;
initialized_
=
false
;
initialize
();
if
(
num_unknowns
()
==
0
)
{
set_fixed_values
(
x
);
return
true
;
}
const
VECTOR
&
r
=
*
b_
;
if
(
solverType_
==
ITERATIVE_SOLVE_SYMMETRIC
)
{
if
(
parallel_
)
{
A
=
new
PAR_SPAR_MAT
(
LammpsInterface
::
instance
()
->
world
(),
*
matrixSparse_
);
}
else
{
A
=
new
SPAR_MAT
(
*
matrixSparse_
);
}
DIAG_MAT
&
PC
=
matrixDiagonal_
;
int
niter
=
maxIterations_
;
double
tol
=
tol_
;
int
convergence
=
CG
(
*
A
,
x
,
r
,
PC
,
niter
,
tol
);
// CG changes niter, tol
if
(
convergence
>
0
)
{
stringstream
ss
;
ss
<<
"CG solve did not converge,"
;
ss
<<
" iterations: "
<<
niter
;
ss
<<
" residual: "
<<
tol
;
throw
ATC_Error
(
ss
.
str
());
}
}
else
if
(
solverType_
==
ITERATIVE_SOLVE
)
{
if
(
parallel_
)
{
A
=
new
PAR_SPAR_MAT
(
LammpsInterface
::
instance
()
->
world
(),
*
matrixSparse_
);
}
else
{
A
=
new
SPAR_MAT
(
*
matrixSparse_
);
}
const
DIAG_MAT
&
PC
=
matrixDiagonal_
;
int
iterations
=
maxIterations_
;
int
restarts
=
maxRestarts_
;
double
tol
=
tol_
;
DENS_MAT
H
(
maxRestarts_
+
1
,
maxRestarts_
);
DENS_VEC
xx
(
nVariables_
);
DENS_VEC
bb
;
bb
=
b
;
int
convergence
=
GMRES
(
*
A
,
xx
,
bb
,
PC
,
H
,
restarts
,
iterations
,
tol
);
if
(
convergence
>
0
)
{
stringstream
ss
;
ss
<<
"GMRES greens_function solve did not converge,"
;
ss
<<
" iterations: "
<<
iterations
;
ss
<<
" residual: "
<<
tol
;
throw
ATC_Error
(
ss
.
str
());
}
x
.
copy
(
xx
.
ptr
(),
xx
.
nRows
());
}
else
{
// DIRECT_SOLVE
const
DENS_MAT
&
invA
=
matrixInverse_
;
if
(
constraintHandlerType_
==
CONDENSE_CONSTRAINTS
)
{
DENS_MAT
xx
=
invA
*
r
;
int
i
=
0
;
set
<
int
>::
const_iterator
itr
;
for
(
itr
=
freeSet_
.
begin
();
itr
!=
freeSet_
.
end
();
itr
++
,
i
++
)
{
int
ii
=
*
itr
;
x
(
ii
)
=
xx
(
i
,
0
);
}
set_fixed_values
(
x
);
}
else
{
DENS_VEC
xx
=
invA
*
r
;
for
(
int
i
=
0
;
i
<
xx
.
nRows
();
i
++
)
{
x
(
i
)
=
xx
(
i
);
}
}
}
delete
A
;
return
true
;
}
// --------------------------------------------------------------------
// greens function
// - returns the solution to a Kronecker delta rhs b = {0 0 .. 1 .. 0 0}
// and with homogeneous constraints {x_i = 0}
// --------------------------------------------------------------------
void
LinearSolver
::
greens_function
(
int
I
,
VECTOR
&
G_I
)
{
SPAR_MAT
*
A
=
NULL
;
initialize_matrix
();
initialize_inverse
();
G_I
.
reset
(
nVariables_
);
VECTOR
&
x
=
G_I
;
if
(
solverType_
==
ITERATIVE_SOLVE_SYMMETRIC
)
{
DENS_VEC
b
(
nVariables_
);
b
=
0.0
;
b
(
I
)
=
1.0
;
if
(
parallel_
)
{
A
=
new
PAR_SPAR_MAT
(
LammpsInterface
::
instance
()
->
world
(),
*
matrixSparse_
);
}
else
{
A
=
new
SPAR_MAT
(
*
matrixSparse_
);
}
const
DIAG_MAT
&
PC
=
matrixDiagonal_
;
int
niter
=
maxIterations_
;
double
tol
=
tol_
;
int
convergence
=
CG
(
*
A
,
x
,
b
,
PC
,
niter
,
tol
);
if
(
convergence
>
0
)
{
stringstream
ss
;
ss
<<
"CG greens_function solve did not converge,"
;
ss
<<
" iterations: "
<<
niter
;
ss
<<
" residual: "
<<
tol
;
throw
ATC_Error
(
ss
.
str
());
}
}
else
if
(
solverType_
==
ITERATIVE_SOLVE
)
{
DENS_VEC
b
(
nVariables_
);
b
=
0.0
;
b
(
I
)
=
1.0
;
// VECTOR & bb = b;
if
(
parallel_
)
{
A
=
new
PAR_SPAR_MAT
(
LammpsInterface
::
instance
()
->
world
(),
*
matrixSparse_
);
}
else
{
A
=
new
SPAR_MAT
(
*
matrixSparse_
);
}
// const DENS_MAT A = matrixSparse_->dense_copy();
const
DIAG_MAT
&
PC
=
matrixDiagonal_
;
int
iterations
=
maxIterations_
;
int
restarts
=
maxRestarts_
;
double
tol
=
tol_
;
DENS_MAT
H
(
maxRestarts_
+
1
,
maxRestarts_
);
DENS_VEC
xx
(
nVariables_
);
int
convergence
=
GMRES
(
*
A
,
xx
,
b
,
PC
,
H
,
restarts
,
iterations
,
tol
);
if
(
convergence
>
0
)
{
stringstream
ss
;
ss
<<
"GMRES greens_function solve did not converge,"
;
ss
<<
" iterations: "
<<
iterations
;
ss
<<
" residual: "
<<
tol
;
throw
ATC_Error
(
ss
.
str
());
}
x
.
copy
(
xx
.
ptr
(),
xx
.
nRows
());
}
else
{
const
DENS_MAT
&
invA
=
matrixInverse_
;
if
(
constraintHandlerType_
==
CONDENSE_CONSTRAINTS
)
{
set
<
int
>::
const_iterator
itr
;
for
(
itr
=
fixedSet_
.
begin
();
itr
!=
fixedSet_
.
end
();
itr
++
)
{
int
ii
=
*
itr
;
x
(
ii
)
=
0
;
}
itr
=
freeSet_
.
find
(
I
);
if
(
itr
!=
freeSet_
.
end
()
)
{
int
j
=
freeGlobalToCondensedMap_
[
I
];
int
i
=
0
;
for
(
itr
=
freeSet_
.
begin
();
itr
!=
freeSet_
.
end
();
itr
++
,
i
++
)
{
int
ii
=
*
itr
;
x
(
ii
)
=
invA
(
j
,
i
);
}
}
}
else
{
for
(
int
i
=
0
;
i
<
nVariables_
;
++
i
)
x
(
i
)
=
invA
(
I
,
i
);
}
}
delete
A
;
}
}
// namespace ATC
Event Timeline
Log In to Comment