Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Praetorius, Simon
mtl4
Commits
0efe17d3
Commit
0efe17d3
authored
Oct 29, 2016
by
Praetorius, Simon
Browse files
replaced boost methods by std methods
parent
d196221a
Changes
75
Hide whitespace changes
Inline
Side-by-side
boost/numeric/itl/krylov/cg.hpp
View file @
0efe17d3
// Software License for MTL
//
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_CG_INCLUDE
...
...
@@ -16,7 +16,6 @@
#include <cmath>
#include <cassert>
#include <iostream>
#include <boost/mpl/bool.hpp>
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/itl/itl_fwd.hpp>
...
...
@@ -35,9 +34,9 @@
namespace
itl
{
/// Conjugate Gradients without preconditioning
template
<
typename
LinearOperator
,
typename
HilbertSpaceX
,
typename
HilbertSpaceB
,
template
<
typename
LinearOperator
,
typename
HilbertSpaceX
,
typename
HilbertSpaceB
,
typename
Iteration
>
int
cg
(
const
LinearOperator
&
A
,
HilbertSpaceX
&
x
,
const
HilbertSpaceB
&
b
,
int
cg
(
const
LinearOperator
&
A
,
HilbertSpaceX
&
x
,
const
HilbertSpaceB
&
b
,
Iteration
&
iter
)
{
mtl
::
vampir_trace
<
7001
>
tracer
;
...
...
@@ -48,20 +47,20 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
Scalar
rho
(
0
),
rho_1
(
0
),
alpha
(
0
),
alpha_1
(
0
);
Vector
p
(
resource
(
x
)),
q
(
resource
(
x
)),
r
(
resource
(
x
)),
z
(
resource
(
x
));
r
=
b
-
A
*
x
;
rho
=
dot
(
r
,
r
);
while
(
!
iter
.
finished
(
Real
(
sqrt
(
abs
(
rho
)))))
{
++
iter
;
if
(
iter
.
first
())
p
=
r
;
else
p
=
r
+
(
rho
/
rho_1
)
*
p
;
else
p
=
r
+
(
rho
/
rho_1
)
*
p
;
// q = A * p; alpha = rho / dot(p, q);
(
lazy
(
q
)
=
A
*
p
)
||
(
lazy
(
alpha_1
)
=
lazy_dot
(
p
,
q
));
alpha
=
rho
/
alpha_1
;
x
+=
alpha
*
p
;
rho_1
=
rho
;
(
lazy
(
r
)
-=
alpha
*
q
)
||
(
lazy
(
rho
)
=
lazy_unary_dot
(
r
));
...
...
@@ -71,9 +70,9 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
}
/// Conjugate Gradients
template
<
typename
LinearOperator
,
typename
HilbertSpaceX
,
typename
HilbertSpaceB
,
template
<
typename
LinearOperator
,
typename
HilbertSpaceX
,
typename
HilbertSpaceB
,
typename
Preconditioner
,
typename
Iteration
>
int
cg
(
const
LinearOperator
&
A
,
HilbertSpaceX
&
x
,
const
HilbertSpaceB
&
b
,
int
cg
(
const
LinearOperator
&
A
,
HilbertSpaceX
&
x
,
const
HilbertSpaceB
&
b
,
const
Preconditioner
&
L
,
Iteration
&
iter
)
{
using
pc
::
is_identity
;
...
...
@@ -88,7 +87,7 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
Scalar
rho
(
0
),
rho_1
(
0
),
rr
,
alpha
(
0
),
alpha_1
;
Vector
p
(
resource
(
x
)),
q
(
resource
(
x
)),
r
(
resource
(
x
)),
z
(
resource
(
x
));
r
=
b
-
A
*
x
;
rr
=
dot
(
r
,
r
);
while
(
!
iter
.
finished
(
Real
(
sqrt
(
abs
(
rr
)))))
{
...
...
@@ -97,12 +96,12 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
if
(
iter
.
first
())
p
=
z
;
else
else
p
=
z
+
(
rho
/
rho_1
)
*
p
;
(
lazy
(
q
)
=
A
*
p
)
||
(
lazy
(
alpha_1
)
=
lazy_dot
(
p
,
q
));
alpha
=
rho
/
alpha_1
;
x
+=
alpha
*
p
;
rho_1
=
rho
;
(
lazy
(
r
)
-=
alpha
*
q
)
||
(
lazy
(
rr
)
=
lazy_unary_dot
(
r
));
...
...
@@ -111,9 +110,9 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
}
/// Conjugate Gradients with ignored right preconditioner to unify interface
template
<
typename
LinearOperator
,
typename
HilbertSpaceX
,
typename
HilbertSpaceB
,
template
<
typename
LinearOperator
,
typename
HilbertSpaceX
,
typename
HilbertSpaceB
,
typename
Preconditioner
,
typename
RightPreconditioner
,
typename
Iteration
>
int
cg
(
const
LinearOperator
&
A
,
HilbertSpaceX
&
x
,
const
HilbertSpaceB
&
b
,
int
cg
(
const
LinearOperator
&
A
,
HilbertSpaceX
&
x
,
const
HilbertSpaceB
&
b
,
const
Preconditioner
&
L
,
const
RightPreconditioner
&
,
Iteration
&
iter
)
{
return
cg
(
A
,
x
,
b
,
L
,
iter
);
...
...
@@ -121,7 +120,7 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
/// Solver class for CG method; right preconditioner ignored (prints warning if not identity)
/** Methods inherited from \ref base_solver. **/
template
<
typename
LinearOperator
,
typename
Preconditioner
,
template
<
typename
LinearOperator
,
typename
Preconditioner
,
typename
RightPreconditioner
>
class
cg_solver
:
public
base_solver
<
cg_solver
<
LinearOperator
,
Preconditioner
,
RightPreconditioner
>
,
LinearOperator
>
...
...
@@ -129,7 +128,7 @@ class cg_solver
typedef
base_solver
<
cg_solver
<
LinearOperator
,
Preconditioner
,
RightPreconditioner
>
,
LinearOperator
>
base
;
public:
/// Construct solver from a linear operator; generate (left) preconditioner from it
explicit
cg_solver
(
const
LinearOperator
&
A
)
:
base
(
A
),
L
(
A
)
explicit
cg_solver
(
const
LinearOperator
&
A
)
:
base
(
A
),
L
(
A
)
{
// MTL_STATIC_ASSERT((!pc::static_is_identity<RightPreconditioner>::value),
// "Right preconditioner must be identity!");
...
...
@@ -138,7 +137,7 @@ class cg_solver
}
/// Construct solver from a linear operator and (left) preconditioner
cg_solver
(
const
LinearOperator
&
A
,
const
Preconditioner
&
L
)
:
base
(
A
),
L
(
L
)
cg_solver
(
const
LinearOperator
&
A
,
const
Preconditioner
&
L
)
:
base
(
A
),
L
(
L
)
{
if
(
!
pc
::
static_is_identity
<
RightPreconditioner
>::
value
)
std
::
cerr
<<
"Right Preconditioner ignored!"
<<
std
::
endl
;
...
...
boost/numeric/itl/krylov/repeating_solver.hpp
View file @
0efe17d3
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG, www.simunova.com.
// 2010 SimuNova UG, www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also tools/license/license.mtl.txt in the distribution.
#ifndef ITL_REPEATING_SOLVER_INCLUDE
#define ITL_REPEATING_SOLVER_INCLUDE
#include <
boost/mpl/if.hpp
>
#include <boost/static_assert.hpp>
#include <
type_traits
>
#include <boost/numeric/itl/itl_fwd.hpp>
namespace
itl
{
...
...
@@ -25,14 +25,14 @@ namespace itl {
template
<
typename
Solver
,
unsigned
N
,
bool
Stored
>
class
repeating_solver
{
typedef
typename
boo
st
::
mpl
::
if_c
<
Stored
,
Solver
,
const
Solver
&>::
type
solver_type
;
typedef
typename
st
d
::
conditional
<
Stored
,
Solver
,
const
Solver
&>::
type
solver_type
;
public:
explicit
repeating_solver
(
const
Solver
&
s
)
:
s
(
s
)
{}
template
<
typename
Matrix
>
explicit
repeating_solver
(
const
Matrix
&
A
)
:
s
(
A
)
{
BOOST_STATIC_ASSERT
(
(
Stored
)
);
// if matrix is passed class must own solver
static_assert
(
Stored
,
""
);
// if matrix is passed class must own solver
}
/// Perform N iterations of the referred solver
...
...
boost/numeric/itl/pc/concat.hpp
View file @
0efe17d3
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG, www.simunova.com.
// 2010 SimuNova UG, www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also tools/license/license.mtl.txt in the distribution.
#ifndef ITL_PC_CONCAT_INCLUDE
#define ITL_PC_CONCAT_INCLUDE
#include <boost/mpl/if.hpp>
#include <type_traits>
#include <boost/numeric/mtl/operation/resource.hpp>
#include <boost/numeric/mtl/interface/vpt.hpp>
#include <boost/numeric/itl/pc/solver.hpp>
...
...
@@ -24,14 +25,14 @@ namespace itl { namespace pc {
template
<
typename
PC1
,
typename
PC2
,
typename
Matrix
,
bool
Store1
=
true
,
bool
Store2
=
true
>
class
concat
{
typedef
typename
boo
st
::
mpl
::
if_c
<
Store1
,
PC1
,
const
PC1
&>::
type
pc1_type
;
typedef
typename
boo
st
::
mpl
::
if_c
<
Store2
,
PC2
,
const
PC2
&>::
type
pc2_type
;
typedef
typename
st
d
::
conditional
<
Store1
,
PC1
,
const
PC1
&>::
type
pc1_type
;
typedef
typename
st
d
::
conditional
<
Store2
,
PC2
,
const
PC2
&>::
type
pc2_type
;
public:
/// Construct both preconditioners from matrix \p A
explicit
concat
(
const
Matrix
&
A
)
:
A
(
A
),
pc1
(
A
),
pc2
(
A
)
{
BOOST_STATIC_ASSERT
(
(
Store1
&&
Store2
)
);
static_assert
(
Store1
&&
Store2
,
""
);
}
/// Both preconditioners are already constructed and passed as arguments
...
...
boost/numeric/itl/pc/ic_0.hpp
View file @
0efe17d3
// Software License for MTL
//
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_PC_IC_0_INCLUDE
#define ITL_PC_IC_0_INCLUDE
#include <
boost/mpl/bool.hpp
>
#include <
type_traits
>
#include <boost/numeric/linear_algebra/identity.hpp>
#include <boost/numeric/linear_algebra/inverse.hpp>
...
...
@@ -99,7 +99,7 @@ class ic_0
void
adjoint_solve
(
const
VectorIn
&
x
,
VectorOut
&
y
)
const
{
mtl
::
vampir_trace
<
5044
>
tracer
;
solve
(
x
,
y
);
solve
(
x
,
y
);
}
...
...
@@ -109,46 +109,46 @@ class ic_0
protected:
template
<
typename
VectorOut
,
typename
Solver
>
friend
struct
ic_0_evaluator
;
// Dummy type to perform factorization in initializer list not in
// Dummy type to perform factorization in initializer list not in
struct
factorizer
{
factorizer
(
const
Matrix
&
A
,
U_type
&
U
)
{
factorize
(
A
,
U
,
mtl
::
traits
::
is_sparse
<
Matrix
>
(),
boo
st
::
is_same
<
Value
,
typename
mtl
::
Collection
<
Matrix
>::
value_type
>
());
}
{
factorize
(
A
,
U
,
mtl
::
traits
::
is_sparse
<
Matrix
>
(),
st
d
::
is_same
<
Value
,
typename
mtl
::
Collection
<
Matrix
>::
value_type
>
());
}
template
<
typename
T
>
void
factorize
(
const
Matrix
&
,
U_type
&
,
boost
::
mpl
::
false_
,
T
)
void
factorize
(
const
Matrix
&
,
U_type
&
,
std
::
false_
type
,
T
)
{
MTL_THROW_IF
(
true
,
mtl
::
logic_error
(
"IC(0) is not suited for dense matrices"
));
}
// When we change the value_type then the factorization is still performed with that of A
template
<
typename
UF
>
void
factorize
(
const
Matrix
&
A
,
UF
&
U
,
boost
::
mpl
::
true_
,
boost
::
mpl
::
false_
)
void
factorize
(
const
Matrix
&
A
,
UF
&
U
,
std
::
true_type
,
std
::
false_
type
)
{
typedef
mtl
::
mat
::
compressed2D
<
typename
mtl
::
Collection
<
Matrix
>::
value_type
,
para
>
tmp_type
;
tmp_type
U_tmp
;
factorize
(
A
,
U_tmp
,
boost
::
mpl
::
true_
(),
boost
::
mpl
::
true_
());
factorize
(
A
,
U_tmp
,
std
::
true_type
(),
std
::
true_
type
());
U
=
U_tmp
;
}
// Factorization adapted from Saad
// Undefined (runtime) behavior if matrix is not symmetric
// Undefined (runtime) behavior if matrix is not symmetric
// UF is type for the factorization
template
<
typename
UF
>
void
factorize
(
const
Matrix
&
A
,
UF
&
U
,
boost
::
mpl
::
true_
,
boost
::
mpl
::
true_
)
template
<
typename
UF
>
void
factorize
(
const
Matrix
&
A
,
UF
&
U
,
std
::
true_type
,
std
::
true_
type
)
{
using
namespace
mtl
;
using
namespace
mtl
::
tag
;
using
mtl
::
traits
::
range_generator
;
using
namespace
mtl
;
using
namespace
mtl
::
tag
;
using
mtl
::
traits
::
range_generator
;
using
math
::
reciprocal
;
using
mtl
::
mat
::
upper
;
mtl
::
vampir_trace
<
5035
>
tracer
;
// For the factorization we take still the value_type of A and later we copy it maybe to another value_type
typedef
typename
mtl
::
Collection
<
Matrix
>::
value_type
value_type
;
typedef
typename
range_generator
<
row
,
UF
>::
type
cur_type
;
typedef
typename
range_generator
<
nz
,
cur_type
>::
type
icur_type
;
typedef
typename
range_generator
<
row
,
UF
>::
type
cur_type
;
typedef
typename
range_generator
<
nz
,
cur_type
>::
type
icur_type
;
MTL_THROW_IF
(
num_rows
(
A
)
!=
num_cols
(
A
),
mtl
::
matrix_not_square
());
U
=
upper
(
A
);
typename
mtl
::
traits
::
col
<
UF
>::
type
col
(
U
);
typename
mtl
::
traits
::
value
<
UF
>::
type
value
(
U
);
typename
mtl
::
traits
::
value
<
UF
>::
type
value
(
U
);
cur_type
kc
=
begin
<
row
>
(
U
),
kend
=
end
<
row
>
(
U
);
for
(
size_type
k
=
0
;
kc
!=
kend
;
++
kc
,
++
k
)
{
...
...
@@ -159,7 +159,7 @@ class ic_0
// U[k][k]= 1.0 / sqrt(U[k][k]);
value_type
inv_dia
=
reciprocal
(
sqrt
(
value
(
*
ic
)));
value
(
*
ic
,
inv_dia
);
// icur_type jbegin=
// icur_type jbegin=
++
ic
;
for
(;
ic
!=
iend
;
++
ic
)
{
// U[k][i] *= U[k][k]
...
...
@@ -174,8 +174,8 @@ class ic_0
icur_type
jc
=
begin
<
nz
>
(
irow
),
jend
=
end
<
nz
>
(
irow
);
while
(
col
(
*
jc
)
<=
k
)
++
jc
;
while
(
col
(
*--
jend
)
>
i
)
;
++
jend
;
++
jend
;
for
(;
jc
!=
jend
;
++
jc
)
{
size_type
j
=
col
(
*
jc
);
U
.
lvalue
(
j
,
i
)
-=
d
*
U
[
k
][
j
];
...
...
@@ -191,7 +191,7 @@ class ic_0
L_type
L
;
lower_solver_t
lower_solver
;
upper_solver_t
upper_solver
;
};
};
#if 0
template <typename Matrix, typename Value, typename Vector>
...
...
@@ -204,9 +204,9 @@ struct ic_0_solver
template <typename VectorOut>
void assign_to(VectorOut& y) const
{ P.solve(x, y); }
{ P.solve(x, y); }
const ic_0<Matrix, Value>& P;
const ic_0<Matrix, Value>& P;
const Vector& x;
};
#endif
...
...
@@ -219,7 +219,7 @@ struct ic_0_evaluator
typedef
typename
mtl
::
Collection
<
VectorOut
>::
value_type
out_value_type
;
ic_0_evaluator
(
VectorOut
&
y
,
const
Solver
&
s
)
ic_0_evaluator
(
VectorOut
&
y
,
const
Solver
&
s
)
:
y
(
y
),
s
(
s
),
U
(
s
.
P
.
U
),
y0
(
s
.
P
.
solve_lower
(
s
.
x
,
y
))
{
MTL_DEBUG_ARG
(
lr
=
99999999
);
}
...
...
boost/numeric/itl/pc/ilu_0.hpp
View file @
0efe17d3
// Software License for MTL
//
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_PC_ILU_0_INCLUDE
#define ITL_PC_ILU_0_INCLUDE
#include <boost/mpl/bool.hpp>
#include <type_traits>
#include <boost/numeric/linear_algebra/identity.hpp>
#include <boost/numeric/mtl/concept/collection.hpp>
...
...
@@ -30,7 +31,7 @@
namespace
itl
{
namespace
pc
{
// Dummy type to perform factorization in initializer list not in
// Dummy type to perform factorization in initializer list not in
struct
ilu_0_factorizer
{
template
<
typename
Matrix
,
typename
L_type
,
typename
U_type
>
...
...
@@ -38,14 +39,14 @@ struct ilu_0_factorizer
{
factorize
(
A
,
L
,
U
,
mtl
::
traits
::
is_sparse
<
Matrix
>
());
}
template
<
typename
Matrix
,
typename
L_type
,
typename
U_type
>
void
factorize
(
const
Matrix
&
,
L_type
&
,
U_type
&
,
boost
::
mpl
::
false_
)
void
factorize
(
const
Matrix
&
,
L_type
&
,
U_type
&
,
std
::
false_
type
)
{
MTL_THROW_IF
(
true
,
mtl
::
logic_error
(
"ILU is not intended for dense matrices"
));
}
template
<
typename
Matrix
,
typename
L_type
,
typename
U_type
>
void
factorize
(
const
Matrix
&
A
,
L_type
&
L
,
U_type
&
U
,
boost
::
mpl
::
true_
)
void
factorize
(
const
Matrix
&
A
,
L_type
&
L
,
U_type
&
U
,
std
::
true_
type
)
{
using
namespace
mtl
;
using
namespace
mtl
::
tag
;
using
mtl
::
traits
::
range_generator
;
using
math
::
reciprocal
;
using
namespace
mtl
;
using
namespace
mtl
::
tag
;
using
mtl
::
traits
::
range_generator
;
using
math
::
reciprocal
;
MTL_THROW_IF
(
num_rows
(
A
)
!=
num_cols
(
A
),
mtl
::
matrix_not_square
());
mtl
::
vampir_trace
<
5038
>
tracer
;
...
...
@@ -55,10 +56,10 @@ struct ilu_0_factorizer
typedef
mtl
::
mat
::
compressed2D
<
value_type
,
para
>
LU_type
;
LU_type
LU
(
A
);
typedef
typename
range_generator
<
row
,
LU_type
>::
type
cur_type
;
typedef
typename
range_generator
<
nz
,
cur_type
>::
type
icur_type
;
typedef
typename
range_generator
<
row
,
LU_type
>::
type
cur_type
;
typedef
typename
range_generator
<
nz
,
cur_type
>::
type
icur_type
;
typename
mtl
::
traits
::
col
<
LU_type
>::
type
col
(
LU
);
typename
mtl
::
traits
::
value
<
LU_type
>::
type
value
(
LU
);
typename
mtl
::
traits
::
value
<
LU_type
>::
type
value
(
LU
);
mtl
::
vec
::
dense_vector
<
value_type
,
mtl
::
vec
::
parameters
<>
>
inv_dia
(
num_rows
(
A
));
cur_type
ic
=
begin
<
row
>
(
LU
),
iend
=
end
<
row
>
(
LU
);
for
(
size_type
i
=
0
;
ic
!=
iend
;
++
ic
,
++
i
)
{
...
...
@@ -72,14 +73,14 @@ struct ilu_0_factorizer
for
(
icur_type
jc
=
kc
+
1
;
jc
!=
kend
;
++
jc
)
value
(
*
jc
,
value
(
*
jc
)
-
aik
*
LU
[
k
][
col
(
*
jc
)]);
// std::cout << "LU after eliminating A[" << i << "][" << k << "] =\n" << LU;
// std::cout << "LU after eliminating A[" << i << "][" << k << "] =\n" << LU;
}
inv_dia
[
i
]
=
reciprocal
(
LU
[
i
][
i
]);
}
invert_diagonal
(
LU
);
invert_diagonal
(
LU
);
L
=
strict_lower
(
LU
);
U
=
upper
(
LU
);
}
}
};
template
<
typename
Matrix
,
typename
Value
=
typename
mtl
::
Collection
<
Matrix
>
::
value_type
>
...
...
boost/numeric/itl/pc/ilut.hpp
View file @
0efe17d3
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_PC_ILUT_INCLUDE
#define ITL_PC_ILUT_INCLUDE
#include <type_traits>
#include <boost/numeric/mtl/vector/sparse_vector.hpp>
#include <boost/numeric/mtl/operation/invert_diagonal.hpp>
#include <boost/numeric/mtl/interface/vpt.hpp>
...
...
@@ -29,26 +31,26 @@ struct ilut_factorizer
// column-major matrices are copied first
template
<
typename
Matrix
,
typename
Para
,
typename
L_type
,
typename
U_type
,
bool
B
>
void
factorize
(
const
Matrix
&
A
,
const
Para
&
p
,
L_type
&
L
,
U_type
&
U
,
boo
st
::
mpl
::
bool
_
<
B
>
)
void
factorize
(
const
Matrix
&
A
,
const
Para
&
p
,
L_type
&
L
,
U_type
&
U
,
st
d
::
integral_constant
<
bool
,
B
>
)
{
typedef
typename
mtl
::
Collection
<
Matrix
>::
value_type
value_type
;
typedef
typename
mtl
::
Collection
<
Matrix
>::
size_type
size_type
;
typedef
mtl
::
mat
::
parameters
<
mtl
::
row_major
,
mtl
::
index
::
c_index
,
mtl
::
non_fixed
::
dimensions
,
false
,
size_type
>
para
;
typedef
mtl
::
mat
::
compressed2D
<
value_type
,
para
>
LU_type
;
LU_type
LU
(
A
);
factorize
(
LU
,
p
,
L
,
U
,
boost
::
mpl
::
true_
());
factorize
(
LU
,
p
,
L
,
U
,
std
::
true_
type
());
}
// According Yousef Saad: ILUT, NLAA, Vol 1(4), 387-402 (1994)
#if 0
template <typename Value, typename MPara, typename Para, typename L_type, typename U_type>
factorize(const mtl::mat::compressed2D<Value, MPara>& A, const Para& p, L_type& L, U_type& U,
boost::mpl
::true_)
factorize(const mtl::mat::compressed2D<Value, MPara>& A, const Para& p, L_type& L, U_type& U,
std
::true_
type
)
#endif
template
<
typename
Matrix
,
typename
Para
,
typename
L_type
,
typename
U_type
>
void
factorize
(
const
Matrix
&
A
,
const
Para
&
p
,
L_type
&
L
,
U_type
&
U
,
boost
::
mpl
::
true_
)
void
factorize
(
const
Matrix
&
A
,
const
Para
&
p
,
L_type
&
L
,
U_type
&
U
,
std
::
true_
type
)
{
{
mtl
::
vampir_trace
<
5049
>
tracer
;
using
std
::
abs
;
using
mtl
::
traits
::
range_generator
;
using
mtl
::
begin
;
using
mtl
::
end
;
using
namespace
mtl
::
tag
;
...
...
@@ -56,22 +58,22 @@ struct ilut_factorizer
typedef
typename
mtl
::
Collection
<
Matrix
>::
value_type
value_type
;
typedef
typename
mtl
::
Collection
<
Matrix
>::
size_type
size_type
;
typedef
typename
range_generator
<
row
,
Matrix
>::
type
cur_type
;
typedef
typename
range_generator
<
nz
,
cur_type
>::
type
icur_type
;
typedef
typename
range_generator
<
row
,
Matrix
>::
type
cur_type
;
typedef
typename
range_generator
<
nz
,
cur_type
>::
type
icur_type
;
typename
mtl
::
traits
::
col
<
Matrix
>::
type
col
(
A
);
typename
mtl
::
traits
::
const_value
<
Matrix
>::
type
value
(
A
);
typename
mtl
::
traits
::
const_value
<
Matrix
>::
type
value
(
A
);
size_type
n
=
num_rows
(
A
);
L
.
change_dim
(
n
,
n
);
L
.
change_dim
(
n
,
n
);
U
.
change_dim
(
n
,
n
);
{
mtl
::
mat
::
inserter
<
L_type
>
L_ins
(
L
,
p
.
first
);
mtl
::
mat
::
inserter
<
U_type
>
U_ins
(
U
,
p
.
first
+
1
);
// plus one for diagonal
mtl
::
sparse_vector
<
value_type
>
vec
(
n
);
// corr. row in paper
cur_type
ic
=
begin
<
row
>
(
A
);
// , iend= end<row>(A);
for
(
size_type
i
=
0
;
i
<
n
;
++
i
,
++
ic
)
{
for
(
icur_type
kc
=
begin
<
nz
>
(
ic
),
kend
=
end
<
nz
>
(
ic
);
kc
!=
kend
;
++
kc
)
// row= A[i][*]
vec
.
insert
(
col
(
*
kc
),
value
(
*
kc
));
// std::cerr << "vec_" << i << " = " << vec << std::endl;
...
...
@@ -96,7 +98,7 @@ struct ilut_factorizer
// std::cerr << "vec_" << i << " = " << vec << std::endl;
vec
.
sort_on_data
();
// std::cerr << "vec_" << i << " sorted on data = " << vec << std::endl;
// std::cout << "vec at " << i << " is " << vec << '\n';
// mtl::vampir_trace<9904> tracer2;
bool
diag_found
=
false
;
...
...
@@ -111,7 +113,7 @@ struct ilut_factorizer
U_ins
[
i
][
k
]
<<
v
;
}
else
// i > k
if
(
cntl
++
<
p
.
first
)
L_ins
[
i
][
k
]
<<
v
;
L_ins
[
i
][
k
]
<<
v
;
}
if
(
!
diag_found
)
std
::
cerr
<<
"Deleted diagonal!!!!
\n
"
;
vec
.
make_empty
();
...
...
@@ -128,7 +130,7 @@ class ilut
{
typedef
ilu
<
Matrix
,
ilut_factorizer
,
Value
>
base
;
public:
ilut
(
const
Matrix
&
A
,
std
::
size_t
p
,
typename
mtl
::
Collection
<
Matrix
>::
value_type
tau
)
ilut
(
const
Matrix
&
A
,
std
::
size_t
p
,
typename
mtl
::
Collection
<
Matrix
>::
value_type
tau
)
:
base
(
A
,
std
::
make_pair
(
p
,
tau
))
{}
};
...
...
boost/numeric/itl/pc/is_identity.hpp
View file @
0efe17d3
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_PC_IS_IDENTITY_INCLUDE
#define ITL_PC_IS_IDENTITY_INCLUDE
#include <boost/mpl/bool.hpp>
#include <type_traits>
#include <boost/numeric/itl/itl_fwd.hpp>
namespace
itl
{
namespace
pc
{
template
<
typename
PC
>
bool
is_identity
(
const
PC
&
)
bool
is_identity
(
const
PC
&
)
{
return
false
;
}
template
<
typename
Matrix
,
typename
Value
>
...
...
@@ -28,12 +29,12 @@ bool is_identity(const itl::pc::identity<Matrix, Value>&)
template
<
typename
PC
>
struct
static_is_identity
:
boost
::
mpl
::
false_
:
std
::
false_
type
{};
template
<
typename
Matrix
,
typename
Value
>
struct
static_is_identity
<
itl
::
pc
::
identity
<
Matrix
,
Value
>
>
:
boost
::
mpl
::
true_
:
std
::
true_
type