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
CoMISo
CoMISo
Commits
dea8257e
Commit
dea8257e
authored
Dec 12, 2016
by
Martin Marinov
Browse files
Merge remote-tracking branch 'origin/warning_fixes'
parents
957aeb32
54dc54db
Changes
33
Hide whitespace changes
Inline
Side-by-side
Solver/CholmodSolver.cc
View file @
dea8257e
...
...
@@ -82,7 +82,7 @@ bool CholmodSolver::calc_system( const std::vector<int>& _colptr,
rowind_
=
_rowind
;
values_
=
_values
;
in
t
n
=
colptr_
.
size
()
-
1
;
size_
t
n
=
colptr_
.
size
()
-
1
;
cholmod_sparse
matA
;
...
...
@@ -167,7 +167,7 @@ bool CholmodSolver::calc_system_prepare_pattern( const std::vector<int>& _col
rowind_
=
_rowind
;
values_
=
_values
;
in
t
n
=
colptr_
.
size
()
-
1
;
size_
t
n
=
colptr_
.
size
()
-
1
;
// setup matrix matA
cholmod_sparse
matA
;
...
...
@@ -282,7 +282,7 @@ bool CholmodSolver::update_system( const std::vector<int>& _colptr,
colptr_
=
_colptr
;
rowind_
=
_rowind
;
values_
=
_values
;
in
t
n
=
colptr_
.
size
()
-
1
;
size_
t
n
=
colptr_
.
size
()
-
1
;
cholmod_sparse
matA
;
...
...
@@ -329,7 +329,7 @@ bool CholmodSolver::update_downdate_factor( const std::vector<int>& _colptr,
rowind_
=
_rowind
;
values_
=
_values
;
in
t
n
=
colptr_
.
size
()
-
1
;
size_
t
n
=
colptr_
.
size
()
-
1
;
cholmod_sparse
matA
;
...
...
@@ -387,7 +387,7 @@ bool CholmodSolver::update_downdate_factor( const std::vector<int>& _colptr,
bool
CholmodSolver
::
solve
(
double
*
_x
,
double
*
_b
)
{
const
unsigned
in
t
n
=
mp_L
->
n
;
const
size_
t
n
=
mp_L
->
n
;
cholmod_dense
*
x
,
b
;
...
...
Solver/ConstrainedSolver.hh
View file @
dea8257e
...
...
@@ -303,13 +303,13 @@ public:
private:
template
<
class
RowT
,
class
MatrixT
>
void
add_row
(
int
_row_i
,
double
_coeff
,
RowT
_row
,
MatrixT
&
_mat
);
void
add_row
(
gmm
::
size_type
_row_i
,
double
_coeff
,
RowT
_row
,
MatrixT
&
_mat
);
template
<
class
RowT
,
class
RMatrixT
,
class
CMatrixT
>
void
add_row_simultaneously
(
int
_row_i
,
void
add_row_simultaneously
(
gmm
::
size_type
_row_i
,
double
_coeff
,
RowT
_row
,
RMatrixT
&
_rmat
,
...
...
Solver/ConstrainedSolverT.cc
View file @
dea8257e
...
...
@@ -217,10 +217,9 @@ solve(
if
(
_show_miso_settings
)
miso_
.
show_options_dialog
();
int
nrows
=
gmm
::
mat_nrows
(
_A
);
int
ncols
=
gmm
::
mat_ncols
(
_A
);
int
ncons
=
gmm
::
mat_nrows
(
_constraints
);
gmm
::
size_type
nrows
=
gmm
::
mat_nrows
(
_A
);
gmm
::
size_type
ncols
=
gmm
::
mat_ncols
(
_A
);
gmm
::
size_type
ncons
=
gmm
::
mat_nrows
(
_constraints
);
DEB_out_if
(
_show_timings
,
1
,
"Initital dimension: "
<<
nrows
<<
" x "
<<
ncols
<<
", number of constraints: "
<<
ncons
...
...
@@ -290,8 +289,8 @@ resolve(
// COMISO_GMM::factored_to_quadratic(_B, A, rhs);
//TODO only compute rhs, not complete A for efficiency
unsigned
int
m
=
gmm
::
mat_nrows
(
_B
);
unsigned
int
n
=
gmm
::
mat_ncols
(
_B
);
gmm
::
size_type
m
=
gmm
::
mat_nrows
(
_B
);
gmm
::
size_type
n
=
gmm
::
mat_ncols
(
_B
);
typedef
typename
gmm
::
linalg_traits
<
RMatrixT
>::
const_sub_row_type
CRowT
;
typedef
typename
gmm
::
linalg_traits
<
RMatrixT
>::
sub_row_type
RowT
;
...
...
@@ -349,8 +348,8 @@ resolve(
gmm
::
mult
(
rhs_update_table_
.
D_
,
*
_constraint_rhs
,
rhs_update_table_
.
cur_constraint_rhs_
);
// update rhs of stored constraints
unsigned
int
nc
=
gmm
::
mat_ncols
(
rhs_update_table_
.
constraints_p_
);
for
(
unsigned
int
i
=
0
;
i
<
rhs_update_table_
.
cur_constraint_rhs_
.
size
();
++
i
)
gmm
::
size_type
nc
=
gmm
::
mat_ncols
(
rhs_update_table_
.
constraints_p_
);
for
(
gmm
::
size_type
i
=
0
;
i
<
rhs_update_table_
.
cur_constraint_rhs_
.
size
();
++
i
)
rhs_update_table_
.
constraints_p_
(
i
,
nc
-
1
)
=
-
rhs_update_table_
.
cur_constraint_rhs_
[
i
];
}
if
(
_rhs
)
...
...
@@ -397,14 +396,14 @@ make_constraints_independent(
{
DEB_enter_func
;
// setup linear transformation for rhs, start with identity
unsigned
int
nr
=
gmm
::
mat_nrows
(
_constraints
);
gmm
::
size_type
nr
=
gmm
::
mat_nrows
(
_constraints
);
gmm
::
resize
(
rhs_update_table_
.
D_
,
nr
,
nr
);
gmm
::
clear
(
rhs_update_table_
.
D_
);
for
(
unsigned
int
i
=
0
;
i
<
nr
;
++
i
)
rhs_update_table_
.
D_
(
i
,
i
)
=
1.0
;
for
(
gmm
::
size_type
i
=
0
;
i
<
nr
;
++
i
)
rhs_update_table_
.
D_
(
i
,
i
)
=
1.0
;
// Base::StopWatch sw;
// number of variables
int
n_vars
=
gmm
::
mat_ncols
(
_constraints
);
const
gmm
::
size_type
n_vars
=
gmm
::
mat_ncols
(
_constraints
);
// TODO Check: HZ added 14.08.09
_c_elim
.
clear
();
...
...
@@ -453,7 +452,7 @@ make_constraints_independent(
for
(;
row_it
!=
row_end
;
++
row_it
)
{
int
cur_j
=
row_it
.
index
();
int
cur_j
=
static_cast
<
int
>
(
row_it
.
index
()
)
;
// do not use the constant part
if
(
cur_j
!=
n_vars
-
1
)
{
...
...
@@ -462,7 +461,7 @@ make_constraints_independent(
{
if
(
fabs
(
*
row_it
)
>
max_elim_val
)
{
elim_j
=
cur_j
;
elim_j
=
(
int
)
cur_j
;
max_elim_val
=
fabs
(
*
row_it
);
}
//break;
...
...
@@ -473,19 +472,18 @@ make_constraints_independent(
// gcd
// if the coefficient of an integer variable is not an integer, then
// the variable most problably will not be (expect if all coeffs are the same, e.g. 0.5)
if
(
(
double
(
int
(
cur_row_val
))
-
cur_row_val
)
!=
0.0
)
{
// std::cerr << __FUNCTION__ << " Warning: coefficient of integer variable is NOT integer: "
// << cur_row_val << std::endl;
gcd_update_valid
=
false
;
}
if
((
double
(
int
(
cur_row_val
))
-
cur_row_val
)
!=
0.0
)
{
DEB_warning
(
2
,
"coefficient of integer variable is NOT integer : "
<<
cur_row_val
)
gcd_update_valid
=
false
;
}
v_gcd
[
n_ints
]
=
cur_row_val
;
v_gcd
[
n_ints
]
=
static_cast
<
int
>
(
cur_row_val
)
;
++
n_ints
;
// store integer closest to 1, must be greater than epsilon_
if
(
fabs
(
cur_row_val
-
1.0
)
<
elim_val
&&
cur_row_val
>
epsilon_
)
{
{
elim_int_j
=
cur_j
;
elim_val
=
fabs
(
cur_row_val
-
1.0
);
}
...
...
@@ -558,7 +556,7 @@ make_constraints_independent(
// sw.start();
double
val
=
-
(
*
c_it
)
/
elim_val_cur
;
add_row_simultaneously
(
c_it
.
index
(),
val
,
gmm
::
mat_row
(
_constraints
,
i
),
_constraints
,
constraints_c
);
add_row_simultaneously
(
(
int
)
c_it
.
index
(),
val
,
gmm
::
mat_row
(
_constraints
,
i
),
_constraints
,
constraints_c
);
// make sure the eliminated entry is 0 on all other rows and not 1e-17
_constraints
(
c_it
.
index
(),
elim_j
)
=
0
;
constraints_c
(
c_it
.
index
(),
elim_j
)
=
0
;
...
...
@@ -586,14 +584,17 @@ make_constraints_independent_reordering(
{
DEB_enter_func
;
// setup linear transformation for rhs, start with identity
unsigned
int
nr
=
gmm
::
mat_nrows
(
_constraints
);
gmm
::
size_type
nr
=
gmm
::
mat_nrows
(
_constraints
);
gmm
::
resize
(
rhs_update_table_
.
D_
,
nr
,
nr
);
gmm
::
clear
(
rhs_update_table_
.
D_
);
for
(
unsigned
int
i
=
0
;
i
<
nr
;
++
i
)
rhs_update_table_
.
D_
(
i
,
i
)
=
1.0
;
for
(
gmm
::
size_type
i
=
0
;
i
<
nr
;
++
i
)
rhs_update_table_
.
D_
(
i
,
i
)
=
1.0
;
// Base::StopWatch sw;
// number of variables
int
n_vars
=
gmm
::
mat_ncols
(
_constraints
);
// AF: Why was n_vars signed? Can it be zero? Later we subtract 1
const
gmm
::
size_type
n_vars
=
gmm
::
mat_ncols
(
_constraints
);
// TODO Check: HZ added 14.08.09
_c_elim
.
clear
();
...
...
@@ -601,7 +602,7 @@ make_constraints_independent_reordering(
// build round map
std
::
vector
<
bool
>
roundmap
(
n_vars
,
false
);
for
(
unsigned
in
t
i
=
0
;
i
<
_idx_to_round
.
size
();
++
i
)
for
(
size_
t
i
=
0
;
i
<
_idx_to_round
.
size
();
++
i
)
roundmap
[
_idx_to_round
[
i
]]
=
true
;
// copy constraints into column matrix (for faster update via iterators)
...
...
@@ -613,17 +614,19 @@ make_constraints_independent_reordering(
// init priority queue
MutablePriorityQueueT
<
unsigned
int
,
unsigned
int
>
queue
;
queue
.
clear
(
nr
);
queue
.
clear
(
static_cast
<
int
>
(
nr
)
);
for
(
unsigned
int
i
=
0
;
i
<
nr
;
++
i
)
{
int
cur_nnz
=
gmm
::
nnz
(
gmm
::
mat_row
(
_constraints
,
i
));
if
(
_constraints
(
i
,
n_vars
-
1
)
!=
0.0
)
--
cur_nnz
;
gmm
::
size_type
cur_nnz
=
gmm
::
nnz
(
gmm
::
mat_row
(
_constraints
,
i
));
if
(
_constraints
(
i
,
n_vars
-
1
)
!=
0.0
)
--
cur_nnz
;
queue
.
update
(
i
,
cur_nnz
);
queue
.
update
(
i
,
static_cast
<
int
>
(
cur_nnz
)
)
;
}
std
::
vector
<
bool
>
row_visited
(
nr
,
false
);
std
::
vector
<
unsigned
int
>
row_ordering
;
row_ordering
.
reserve
(
nr
);
std
::
vector
<
gmm
::
size_type
>
row_ordering
;
row_ordering
.
reserve
(
nr
);
// for all conditions
...
...
@@ -663,16 +666,16 @@ make_constraints_independent_reordering(
for
(;
row_it
!=
row_end
;
++
row_it
)
{
int
cur_j
=
row_it
.
index
();
int
cur_j
=
static_cast
<
int
>
(
row_it
.
index
()
)
;
// do not use the constant part
if
(
cur_j
!=
n_vars
-
1
)
if
(
cur_j
!=
n_vars
-
1
)
{
// found real valued var? -> finished (UPDATE: no not any more, find biggest real value to avoid x/1e-13)
if
(
!
roundmap
[
cur_j
])
if
(
!
roundmap
[
cur_j
])
{
if
(
fabs
(
*
row_it
)
>
max_elim_val
)
if
(
fabs
(
*
row_it
)
>
max_elim_val
)
{
elim_j
=
cur_j
;
elim_j
=
(
int
)
cur_j
;
max_elim_val
=
fabs
(
*
row_it
);
}
//break;
...
...
@@ -683,20 +686,19 @@ make_constraints_independent_reordering(
// gcd
// if the coefficient of an integer variable is not an integer, then
// the variable most problably will not be (expect if all coeffs are the same, e.g. 0.5)
if
(
(
double
(
int
(
cur_row_val
))
-
cur_row_val
)
!=
0.0
)
{
// std::cerr << __FUNCTION__ << " Warning: coefficient of integer variable is NOT integer: "
// << cur_row_val << std::endl;
gcd_update_valid
=
false
;
}
if
((
double
(
int
(
cur_row_val
))
-
cur_row_val
)
!=
0.0
)
{
DEB_warning
(
2
,
"coefficient of integer variable is NOT integer : "
<<
cur_row_val
);
gcd_update_valid
=
false
;
}
v_gcd
[
n_ints
]
=
cur_row_val
;
v_gcd
[
n_ints
]
=
static_cast
<
int
>
(
cur_row_val
)
;
++
n_ints
;
// store integer closest to 1, must be greater than epsilon_
if
(
fabs
(
cur_row_val
-
1.0
)
<
elim_val
&&
cur_row_val
>
epsilon_
)
{
elim_int_j
=
cur_j
;
elim_int_j
=
(
int
)
cur_j
;
elim_val
=
fabs
(
cur_row_val
-
1.0
);
}
}
...
...
@@ -734,7 +736,7 @@ make_constraints_independent_reordering(
if
(
do_gcd_
&&
gcd_update_valid
)
{
// perform gcd update
bool
gcd_ok
=
update_constraint_gcd
(
_constraints
,
i
,
elim_j
,
v_gcd
,
n_ints
);
bool
gcd_ok
=
update_constraint_gcd
(
_constraints
,
(
int
)
i
,
elim_j
,
v_gcd
,
n_ints
);
DEB_warning_if
(
!
gcd_ok
&&
(
noisy_
>
0
),
1
,
" GCD update failed! "
<<
DEB_os_str
(
gmm
::
mat_const_row
(
_constraints
,
i
))
)
}
...
...
@@ -755,7 +757,7 @@ make_constraints_independent_reordering(
// is this condition dependent?
if
(
elim_j
!=
-
1
)
if
(
elim_j
!=
-
1
)
{
// get elim variable value
double
elim_val_cur
=
_constraints
(
i
,
elim_j
);
...
...
@@ -764,30 +766,34 @@ make_constraints_independent_reordering(
CVector
col
=
constraints_c
.
col
(
elim_j
);
// iterate over column
typename
gmm
::
linalg_traits
<
CVector
>::
const_iterator
c_it
=
gmm
::
vect_const_begin
(
col
);
typename
gmm
::
linalg_traits
<
CVector
>::
const_iterator
c_end
=
gmm
::
vect_const_end
(
col
);
typename
gmm
::
linalg_traits
<
CVector
>::
const_iterator
c_it
=
gmm
::
vect_const_begin
(
col
);
typename
gmm
::
linalg_traits
<
CVector
>::
const_iterator
c_end
=
gmm
::
vect_const_end
(
col
);
for
(;
c_it
!=
c_end
;
++
c_it
)
// if( c_it.index() > i)
if
(
!
row_visited
[
c_it
.
index
()])
for
(;
c_it
!=
c_end
;
++
c_it
)
{
// if( c_it.index() > i)
if
(
!
row_visited
[
c_it
.
index
()])
{
// sw.start();
double
val
=
-
(
*
c_it
)
/
elim_val_cur
;
add_row_simultaneously
(
c_it
.
index
(),
val
,
gmm
::
mat_row
(
_constraints
,
i
),
_constraints
,
constraints_c
);
// sw.start();
double
val
=
-
(
*
c_it
)
/
elim_val_cur
;
add_row_simultaneously
(
(
int
)
c_it
.
index
(),
val
,
gmm
::
mat_row
(
_constraints
,
i
),
_constraints
,
constraints_c
);
// make sure the eliminated entry is 0 on all other rows and not 1e-17
_constraints
(
c_it
.
index
(),
elim_j
)
=
0
;
_constraints
(
c_it
.
index
(),
elim_j
)
=
0
;
constraints_c
(
c_it
.
index
(),
elim_j
)
=
0
;
int
cur_idx
=
c_it
.
index
();
int
cur_nnz
=
gmm
::
nnz
(
gmm
::
mat_row
(
_constraints
,
cur_idx
));
if
(
_constraints
(
cur_idx
,
n_vars
-
1
)
!=
0.0
)
--
cur_nnz
;
gmm
::
size_type
cur_idx
=
c_it
.
index
();
gmm
::
size_type
cur_nnz
=
gmm
::
nnz
(
gmm
::
mat_row
(
_constraints
,
cur_idx
));
if
(
_constraints
(
cur_idx
,
n_vars
-
1
)
!=
0.0
)
--
cur_nnz
;
queue
.
update
(
cur_idx
,
cur_nnz
);
queue
.
update
(
static_cast
<
int
>
(
cur_idx
),
static_cast
<
int
>
(
cur_nnz
));
// update linear transition of rhs
gmm
::
add
(
gmm
::
scaled
(
gmm
::
mat_row
(
rhs_update_table_
.
D_
,
i
),
val
),
gmm
::
mat_row
(
rhs_update_table_
.
D_
,
c_it
.
index
()));
gmm
::
mat_row
(
rhs_update_table_
.
D_
,
c_it
.
index
()));
}
}
}
}
// // check result
...
...
@@ -857,10 +863,10 @@ update_constraint_gcd( RMatrixT& _constraints,
for
(
;
row_it
!=
row_end
;
++
row_it
)
{
int
cur_j
=
row_it
.
index
();
gmm
::
size_type
cur_j
=
row_it
.
index
();
_constraints
(
_row_i
,
cur_j
)
=
(
*
row_it
)
/
i_gcd
;
}
int
elim_coeff
=
abs
(
_constraints
(
_row_i
,
_elim_j
));
gmm
::
size_type
elim_coeff
=
static_cast
<
gmm
::
size_type
>
(
std
::
abs
(
_constraints
(
_row_i
,
_elim_j
))
)
;
DEB_error_if
(
elim_coeff
!=
1
,
"elimination coefficient "
<<
elim_coeff
<<
" will (most probably) NOT lead to an integer solution!"
)
return
true
;
...
...
@@ -1026,7 +1032,7 @@ eliminate_constraints(
//_rhs[con_it.index()] -= cur_rhs * (( *con_it )/cur_val);
// rhs_update_table_.append(con_it.index(), -1.0 * (( *con_it )/cur_val), cur_j);
rhs_update_table_
.
append
(
con_it
.
index
(),
-
1.0
*
((
*
con_it
)
/
cur_val
),
cur_j
,
false
);
rhs_update_table_
.
append
(
static_cast
<
int
>
(
con_it
.
index
()
)
,
-
1.0
*
((
*
con_it
)
/
cur_val
),
static_cast
<
int
>
(
cur_j
)
,
false
);
//std::cerr << con_it.index() << " += " << -1.0*(( *con_it )/cur_val) << " * " << cur_rhs << " (["<<cur_j<<"] = "<<_rhs[cur_j]<<") " << std::endl;
}
...
...
@@ -1047,7 +1053,7 @@ eliminate_constraints(
_A
(
col_it
.
index
(),
con_it
.
index
()
)
-=
(
*
col_it
)
*
((
*
con_it
)
/
cur_val
);
//_rhs[col_it.index()] += constraint_rhs*( *col_it )/cur_val;
// rhs_update_table_.append(col_it.index(), constraint_rhs*( *col_it )/cur_val);
rhs_update_table_
.
append
(
col_it
.
index
(),
-
(
*
col_it
)
/
cur_val
,
i
,
true
);
rhs_update_table_
.
append
(
static_cast
<
int
>
(
col_it
.
index
()
)
,
-
(
*
col_it
)
/
cur_val
,
i
,
true
);
//std::cerr << col_it.index() << " += " << constraint_rhs*( *col_it )/cur_val << std::endl;
}
...
...
@@ -1108,7 +1114,7 @@ eliminate_constraints(
template
<
class
RowT
,
class
MatrixT
>
void
ConstrainedSolver
::
add_row
(
int
_row_i
,
add_row
(
gmm
::
size_type
_row_i
,
double
_coeff
,
RowT
_row
,
MatrixT
&
_mat
)
...
...
@@ -1128,11 +1134,11 @@ add_row( int _row_i,
template
<
class
RowT
,
class
RMatrixT
,
class
CMatrixT
>
void
ConstrainedSolver
::
add_row_simultaneously
(
int
_row_i
,
double
_coeff
,
RowT
_row
,
RMatrixT
&
_rmat
,
CMatrixT
&
_cmat
)
add_row_simultaneously
(
gmm
::
size_type
_row_i
,
double
_coeff
,
RowT
_row
,
RMatrixT
&
_rmat
,
CMatrixT
&
_cmat
)
{
typedef
typename
gmm
::
linalg_traits
<
RowT
>::
const_iterator
RIter
;
RIter
r_it
=
gmm
::
vect_const_begin
(
_row
);
...
...
@@ -1248,7 +1254,7 @@ restore_eliminated_vars( RMatrixT& _constraints,
_x
.
back
()
=
1.0
;
// reverse iterate from prelast element
for
(
int
i
=
_new_idx
.
size
()
-
2
;
i
>=
0
;
--
i
)
for
(
int
i
=
static_cast
<
int
>
(
_new_idx
.
size
()
)
-
2
;
i
>=
0
;
--
i
)
// AF: Can this be negative?
{
if
(
_new_idx
[
i
]
!=
-
1
)
{
...
...
@@ -1259,7 +1265,7 @@ restore_eliminated_vars( RMatrixT& _constraints,
}
// reverse iterate
for
(
int
i
=
_c_elim
.
size
()
-
1
;
i
>=
0
;
--
i
)
for
(
int
i
=
static_cast
<
int
>
(
_c_elim
.
size
()
)
-
1
;
i
>=
0
;
--
i
)
// AF: Can this be negative?
{
int
cur_var
=
_c_elim
[
i
];
...
...
Solver/EigenLDLTSolver.hh
View file @
dea8257e
...
...
@@ -45,9 +45,12 @@
#include <vector>
#include <Base/Code/Quality.hh>
LOW_CODE_QUALITY_SECTION_BEGIN
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
LOW_CODE_QUALITY_SECTION_END
//== NAMESPACES ===============================================================
...
...
@@ -55,44 +58,46 @@
namespace
COMISO
{
//== CLASS DEFINITION =========================================================
class
COMISODLLEXPORT
EigenLDLTSolver
class
EigenLDLTSolver
{
public:
// _size is maximal size this instance can handle (smaller problems are possible!!!)
EigenLDLTSolver
();
~
EigenLDLTSolver
();
COMISODLLEXPORT
EigenLDLTSolver
();
COMISODLLEXPORT
~
EigenLDLTSolver
();
COMISODLLEXPORT
bool
calc_system
(
const
std
::
vector
<
int
>&
_colptr
,
const
std
::
vector
<
int
>&
_rowind
,
const
std
::
vector
<
double
>&
_values
);
template
<
class
GMM_MatrixT
>
bool
calc_system_gmm
(
const
GMM_MatrixT
&
_mat
);
template
<
class
Eigen_MatrixT
>
bool
calc_system_eigen
(
const
Eigen_MatrixT
&
_mat
);
COMISODLLEXPORT
bool
update_system
(
const
std
::
vector
<
int
>&
_colptr
,
const
std
::
vector
<
int
>&
_rowind
,
const
std
::
vector
<
double
>&
_values
);
template
<
class
GMM_MatrixT
>
bool
update_system_gmm
(
const
GMM_MatrixT
&
_mat
);
template
<
class
Eigen_MatrixT
>
bool
update_system_eigen
(
const
Eigen_MatrixT
&
_mat
);
COMISODLLEXPORT
bool
solve
(
double
*
_x0
,
double
*
_b
);
COMISODLLEXPORT
bool
solve
(
std
::
vector
<
double
>&
_x0
,
std
::
vector
<
double
>&
_b
);
COMISODLLEXPORT
bool
&
show_timings
();
COMISODLLEXPORT
int
dimension
();
private:
...
...
Solver/Eigen_Tools.hh
View file @
dea8257e
...
...
@@ -33,16 +33,18 @@
//== INCLUDES =================================================================
#include <Base/Code/Quality.hh>
#include <iostream>
#include <vector>
#include <algorithm>
#include <limits>
#include <cmath>
LOW_CODE_QUALITY_SECTION_BEGIN
#include <Eigen/Dense>
#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
#include <Eigen/Sparse>
LOW_CODE_QUALITY_SECTION_END
#if COMISO_SUITESPARSE_AVAILABLE
#include <cholmod.h>
...
...
Solver/Eigen_ToolsT.cc
View file @
dea8257e
...
...
@@ -340,7 +340,7 @@ void cholmod_to_eigen( const cholmod_sparse& _AC, MatrixT& _A)
typedef
Eigen
::
Triplet
<
Scalar
>
Triplet
;
size_t
nzmax
(
_AC
.
nzmax
);
std
::
cerr
<<
__FUNCTION__
<<
" row "
<<
_AC
.
nrow
<<
" col "
<<
_AC
.
ncol
<<
" stype "
<<
_AC
.
stype
<<
std
::
endl
;
_A
=
MatrixT
(
_AC
.
nrow
,
_AC
.
ncol
);
_A
=
MatrixT
(
(
int
)
_AC
.
nrow
,
(
int
)
_AC
.
ncol
);
std
::
vector
<
Triplet
>
triplets
;
triplets
.
reserve
(
nzmax
);
...
...
@@ -365,7 +365,7 @@ void cholmod_to_eigen( const cholmod_sparse& _AC, MatrixT& _A)
for
(
SuiteSparse_long
i
=
0
;
i
<
(
SuiteSparse_long
)
_AC
.
ncol
;
++
i
)
for
(
SuiteSparse_long
j
=
P
[
i
];
j
<
P
[
i
+
1
];
++
j
)
//_A( I[j], i) += X[j]; // += really needed?
triplets
.
push_back
(
Triplet
(
I
[
j
],
i
,
X
[
j
]));
triplets
.
push_back
(
Triplet
(
(
int
)
I
[
j
],
(
int
)
i
,
X
[
j
]));
}
else
{
...
...
@@ -374,7 +374,7 @@ void cholmod_to_eigen( const cholmod_sparse& _AC, MatrixT& _A)
for
(
int
i
=
0
;
i
<
(
int
)
_AC
.
ncol
;
++
i
)
for
(
int
j
=
P
[
i
];
j
<
P
[
i
+
1
];
++
j
)
triplets
.
push_back
(
Triplet
(
I
[
j
],
i
,
X
[
j
]));
triplets
.
push_back
(
Triplet
(
(
int
)
I
[
j
],
(
int
)
i
,
X
[
j
]));
//_A( I[j], i) += X[j];
}
...
...
@@ -391,11 +391,11 @@ void cholmod_to_eigen( const cholmod_sparse& _AC, MatrixT& _A)
for
(
SuiteSparse_long
j
=
P
[
i
];
j
<
P
[
i
+
1
];
++
j
)
{
//_A(I[j], i) += X[j];
triplets
.
push_back
(
Triplet
(
I
[
j
],
i
,
X
[
j
]));
triplets
.
push_back
(
Triplet
(
(
int
)
I
[
j
],
(
int
)
i
,
X
[
j
]));
// add up symmetric part
if
(
I
[
j
]
!=
i
)
triplets
.
push_back
(
Triplet
(
i
,
I
[
j
],
X
[
j
]));
triplets
.
push_back
(
Triplet
(
(
int
)
i
,
(
int
)
I
[
j
],
X
[
j
]));
//_A(i,I[j]) += X[j];
}
}
...
...
@@ -731,19 +731,21 @@ void f(const gmm::row_matrix<GMM_VectorT>& _G, EIGEN_MatrixT& _E)
std
::
vector
<
Triplet
>
triplets
;
triplets
.
reserve
(
gmm
::
nnz
(
_G
));
for
(
unsigned
int
i
=
0
;
i
<
gmm
::
mat_nrows
(
_G
);
++
i
)
for
(
gmm
::
size_type
i
=
0
;
i
<
gmm
::
mat_nrows
(
_G
);
++
i
)
{
RowT
row
=
mat_const_row
(
_G
,
i
);
CIter
it
=
gmm
::
vect_const_begin
(
row
);
CIter
ite
=
gmm
::
vect_const_end
(
row
);
for
(
;
it
!=
ite
;
++
it
)
triplets
.
push_back
(
Triplet
(
i
,
it
.
index
(),
*
it
));
triplets
.
push_back
(
Triplet
(
static_cast
<
int
>
(
i
),
static_cast
<
int
>
(
it
.
index
()),
*
it
));
}
// generate eigen matrix
_E
=
EIGEN_MatrixT
(
gmm
::
mat_nrows
(
_G
),
gmm
::
mat_ncols
(
_G
));
_E
=
EIGEN_MatrixT
(
static_cast
<
int
>
(
gmm
::
mat_nrows
(
_G
)),
static_cast
<
int
>
(
gmm
::
mat_ncols
(
_G
)));
_E
.
setFromTriplets
(
triplets
.
begin
(),
triplets
.
end
());
}
...
...
@@ -769,11 +771,12 @@ void f(const gmm::csc_matrix<GMM_RealT,0>& _G, EIGEN_MatrixT& _E)
CIter
it
=
gmm
::
vect_const_begin
(
col
);
CIter
ite
=
gmm
::
vect_const_end
(
col
);
for
(
;
it
!=
ite
;
++
it
)
triplets
.
push_back
(
Triplet
(
it
.
index
(),
i
,
*
it
));
triplets
.
push_back
(
Triplet
(
static_cast
<
int
>
(
it
.
index
()
)
,
i
,
*
it
));
}
// generate eigen matrix
_E
=
EIGEN_MatrixT
(
gmm
::
mat_nrows
(
_G
),
gmm
::
mat_ncols
(
_G
));
_E
=
EIGEN_MatrixT
(
static_cast
<
int
>
(
gmm
::
mat_nrows
(
_G
)),
static_cast
<
int
>
(
gmm
::
mat_ncols
(
_G
)));
_E
.
setFromTriplets
(
triplets
.
begin
(),
triplets
.
end
());
}
...
...
Solver/GMM_ToolsT.cc
View file @
dea8257e
...
...
@@ -65,13 +65,13 @@ void eliminate_csc_vars(
typedef
unsigned
int
uint
;
typedef
typename
gmm
::
csc_matrix
<
RealT
>
MatrixT
;
unsigned
int
nc
=
_A
.
nc
;
unsigned
int
nr
=
_A
.
nr
;
gmm
::
size_type
nc
=
_A
.
nc
;
gmm
::
size_type
nr
=
_A
.
nr
;
unsigned
int
n_new
=
nc
-
_evar
.
size
();
gmm
::
size_type
n_new
=
nc
-
_evar
.
size
();
// modify rhs
for
(
unsigned
int
k
=
0
;
k
<
_evar
.
size
();
++
k
)
for
(
gmm
::
size_type
k
=
0
;
k
<
_evar
.
size
();
++
k
)
{
IntegerT
i
=
_evar
[
k
];
...
...
@@ -95,7 +95,7 @@ void eliminate_csc_vars(
// build subindex set and update rhs
std
::
vector
<
size_t
>
si
(
n_new
);
int
cur_evar_idx
=
0
;
for
(
unsigned
int
i
=
0
;
i
<
nc
;
++
i
)
for
(
gmm
::
size_type
i
=
0
;
i
<
nc
;
++
i
)
{
unsigned
int
next_i
=
evar
[
cur_evar_idx
];
if
(
i
!=
next_i
)
...
...
@@ -119,7 +119,7 @@ void eliminate_csc_vars(