Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Menu
Open sidebar
Guilhem Saurel
soth
Commits
3e1ae6c8
Commit
3e1ae6c8
authored
Dec 30, 2015
by
Olivier Stasse
Browse files
Add eigen 3.2.x support + Fix some warnings.
Keep backward compatibility.
parent
2ac34d4f
Changes
23
Hide whitespace changes
Inline
Side-by-side
exe/conic_simplification.cpp
View file @
3e1ae6c8
...
...
@@ -37,7 +37,7 @@ bool checkColumn( const Eigen::MatrixBase<D>& A,
const
Eigen
::
MatrixBase
<
D
>&
B
,
const
int
index
,
int
sign
=
+
1
)
{
const
int
NC
=
A
.
rows
(),
NA
=
A
.
cols
(),
NB
=
B
.
cols
(),
NX
=
NA
+
NB
;
const
int
NC
=
(
int
)
A
.
rows
(),
NA
=
(
int
)
A
.
cols
(),
NB
=
(
int
)
B
.
cols
(),
NX
=
(
int
)(
NA
+
NB
)
;
std
::
vector
<
Eigen
::
MatrixXd
>
J
(
3
);
std
::
vector
<
soth
::
VectorBound
>
b
(
3
);
...
...
@@ -98,7 +98,7 @@ int checkAndModify( SubMatrix<MatrixXd,ColPermutation> &A,
}
if
(
checkColumn
(
A
,
B
,
index
,
-
1
)
)
{
int
ref
=
A
.
removeCol
(
index
);
int
ref
=
(
int
)
A
.
removeCol
(
index
);
B
.
pushColFront
(
ref
);
cout
<<
"transfert ... "
<<
endl
;
return
0
;
...
...
@@ -116,7 +116,7 @@ bool checkOut( const MatrixXd &X,
/* shoot ab = rand, with a>0 and target=AB*ab
* solve X*x=target, with x>0. */
const
int
NC
=
A
.
rows
(),
NA
=
A
.
cols
(),
NB
=
B
.
cols
(),
NX
=
X
.
cols
();
const
int
NC
=
(
int
)
A
.
rows
(),
NA
=
(
int
)
A
.
cols
(),
NB
=
(
int
)
B
.
cols
(),
NX
=
(
int
)
X
.
cols
();
VectorXd
ab
=
VectorXd
::
Random
(
NA
+
NB
);
ab
.
head
(
NA
)
+=
VectorXd
::
Ones
(
NA
);
ab
.
head
(
NA
)
/=
2
;
VectorXd
target
=
A
*
ab
.
head
(
NA
)
+
B
*
ab
.
tail
(
NB
);
cout
<<
"ab = "
<<
(
MATLAB
)
ab
<<
endl
;
...
...
@@ -162,7 +162,7 @@ bool checkIn( const MatrixXd &X,
/* shoot x = rand>0, and target=X*x
* solve Aa+Bb = x, with a>0. */
const
int
NC
=
A
.
rows
(),
NA
=
A
.
cols
(),
NB
=
B
.
cols
(),
NX
=
X
.
cols
();
const
int
NC
=
(
int
)
A
.
rows
(),
NA
=
(
int
)
A
.
cols
(),
NB
=
(
int
)
B
.
cols
(),
NX
=
(
int
)
X
.
cols
();
VectorXd
x
=
VectorXd
::
Random
(
NX
);
x
+=
VectorXd
::
Ones
(
NX
);
x
/=
2
;
VectorXd
target
=
X
*
x
;
cout
<<
"x = "
<<
(
MATLAB
)
x
<<
endl
;
...
...
@@ -201,7 +201,7 @@ bool checkIn( const MatrixXd &X,
int
main
(
int
argc
,
char
**
argv
)
int
main
(
int
,
char
**
)
{
# ifndef NDEBUG
sotDebugTrace
::
openFile
();
...
...
exe/simple.cpp
View file @
3e1ae6c8
...
...
@@ -55,7 +55,7 @@ int main ()
std
::
vector
<
soth
::
VectorBound
>
b
;
soth
::
Random
::
setSeed
(
704819
);
const
int
size
=
100
;
//
const int size = 100;
//generateFixedSizeRandomProfile(size,
// 1,0.99,0.99,NB_STAGE,RANKFREE,RANKLINKED,NR,NC);
...
...
@@ -79,7 +79,7 @@ int main ()
{
RANKFREE
+=
6
,
3
,
7
,
2
,
5
,
5
,
44
,
6
,
2
,
4
,
5
,
3
,
7
,
1
,
50
;
RANKLINKED
.
resize
(
RANKFREE
.
size
(),
0
);
NB_STAGE
=
RANKFREE
.
size
();
NB_STAGE
=
(
unsigned
int
)
RANKFREE
.
size
();
NR
=
RANKFREE
;
NC
=
150
;
}
...
...
@@ -87,7 +87,7 @@ int main ()
std
::
cout
<<
"nVar
\t
= "
<<
NC
<<
std
::
endl
;
std
::
cout
<<
"nLevels
\t
= "
<<
NB_STAGE
<<
std
::
endl
;
std
::
cout
<<
"LevelDim
\t
= [ "
;
for
(
int
i
=
0
;
i
<
NB_STAGE
;
++
i
)
std
::
cout
<<
NR
[
i
]
<<
" "
;
for
(
int
i
=
0
;
i
<
(
int
)
NB_STAGE
;
++
i
)
std
::
cout
<<
NR
[
i
]
<<
" "
;
std
::
cout
<<
"]"
<<
std
::
endl
;
generateDeficientDataSet
(
J
,
b
,
NB_STAGE
,
RANKFREE
,
RANKLINKED
,
NR
,
NC
);
...
...
@@ -109,7 +109,7 @@ int main ()
gettimeofday
(
&
t0
,
NULL
);
for
(
int
i
=
0
;
i
<
1000
;
++
i
)
ehqp
(
hsolver
);
gettimeofday
(
&
t1
,
NULL
);
time
=
(
t1
.
tv_sec
-
t0
.
tv_sec
)
+
(
t1
.
tv_usec
-
t0
.
tv_usec
)
/
1.0e6
;
time
=
(
double
)
(
t1
.
tv_sec
-
t0
.
tv_sec
)
+
(
double
)(
t1
.
tv_usec
-
t0
.
tv_usec
)
/
1.0e6
;
cout
<<
"ehqp = "
<<
time
<<
" ms "
<<
std
::
endl
;
...
...
@@ -121,7 +121,7 @@ int main ()
Eigen
::
DestructiveColPivQR
<
MatrixXd
,
MatrixXd
>
mQR
(
A
,
Y
,
1e-6
);
//Eigen::DestructiveColPivQR<SubMatrixXd,MatrixXd > mQR(Ar,Y, 1e-6);
gettimeofday
(
&
t1
,
NULL
);
time
=
(
t1
.
tv_sec
-
t0
.
tv_sec
)
+
(
t1
.
tv_usec
-
t0
.
tv_usec
)
/
1.0e6
;
time
=
(
double
)
(
t1
.
tv_sec
-
t0
.
tv_sec
)
+
(
double
)(
t1
.
tv_usec
-
t0
.
tv_usec
)
/
1.0e6
;
cout
<<
"qr = "
<<
time
<<
" ms "
<<
std
::
endl
;
...
...
src/ActiveSet.cpp
View file @
3e1ae6c8
...
...
@@ -65,7 +65,7 @@ namespace soth
activeRow
(
unsigned
int
ref
,
Bound
::
bound_t
type
)
{
const
unsigned
int
row
=
getAFreeRow
();
assert
(
(
row
>=
0
)
&&
(
row
<
size
()
)
);
assert
(
row
<
size
()
);
active
(
ref
,
type
,
row
);
return
row
;
}
...
...
src/ActiveSet.hpp
View file @
3e1ae6c8
...
...
@@ -44,7 +44,7 @@ namespace soth
public:
/* --- Accessors --- */
/* Return the number of active constraint. */
inline
unsigned
int
nbActive
(
void
)
const
{
return
nba
;
}
inline
unsigned
int
size
(
void
)
const
{
return
cstMap
.
size
();
}
inline
unsigned
int
size
(
void
)
const
{
return
(
unsigned
int
)
cstMap
.
size
();
}
bool
isFreezed
(
unsigned
int
ref
)
const
;
bool
isActive
(
unsigned
int
ref
)
const
;
bool
wasActive
(
unsigned
int
ref
,
const
Bound
::
bound_t
type
)
const
;
...
...
src/Algebra.hpp
View file @
3e1ae6c8
...
...
@@ -66,8 +66,8 @@ namespace soth
template
<
typename
Derived
>
void
MATLAB
::
genericInit
(
const
MatrixBase
<
Derived
>
&
m1
)
{
if
(
m1
.
rows
()
==
0
)
initMatrixRowNull
(
m1
.
cols
()
);
else
if
(
m1
.
cols
()
==
0
)
initMatrixColNull
(
m1
.
rows
()
);
if
(
m1
.
rows
()
==
0
)
initMatrixRowNull
(
(
unsigned
int
)
m1
.
cols
()
);
else
if
(
m1
.
cols
()
==
0
)
initMatrixColNull
(
(
unsigned
int
)
m1
.
rows
()
);
else
if
(
m1
.
IsVectorAtCompileTime
)
{
if
(
m1
.
cols
()
==
1
)
initVector
(
m1
.
col
(
0
)
);
...
...
@@ -91,7 +91,7 @@ namespace soth
if
(
m1
.
size
()
!=
i
+
1
)
{
ostmp
<<
","
;
const
int
size
=
ostmp
.
str
().
length
();
const
int
size
=
(
int
)
ostmp
.
str
().
length
();
for
(
unsigned
int
i
=
size
;
i
<
8
;
++
i
)
ostmp
<<
" "
;
ostmp
<<
"
\t
"
;
}
...
...
@@ -138,7 +138,7 @@ namespace soth
if
(
m1
.
cols
()
!=
j
+
1
)
{
ostmp
<<
","
;
const
int
size
=
ostmp
.
str
().
length
();
const
int
size
=
(
int
)
ostmp
.
str
().
length
();
for
(
unsigned
int
i
=
size
;
i
<
8
;
++
i
)
ostmp
<<
" "
;
ostmp
<<
"
\t
"
;
}
...
...
src/BasicStage.cpp
View file @
3e1ae6c8
...
...
@@ -54,7 +54,7 @@ namespace soth
,
boundsMap
(
inbounds
.
data
(),
inbounds
.
size
(),
1
)
,
J
(
Jmap
),
bounds
(
boundsMap
)
,
nr
(
inJ
.
rows
()
),
nc
(
inJ
.
cols
()
)
,
nr
(
(
unsigned
int
)
inJ
.
rows
()
),
nc
(
(
unsigned
int
)
inJ
.
cols
()
)
,
Y
(
inY
)
{
...
...
src/DestructiveColPivQR.hpp
View file @
3e1ae6c8
...
...
@@ -5,6 +5,18 @@
namespace
Eigen
{
#if EIGEN_VERSION_AT_LEAST(3,2,0)
#define LOCAL_ABS(x) std::abs(x)
#define LOCAL_ABS2(x) numext::abs2(x)
#else
#define LOCAL_ABS(x) internal::abs(x)
#define LOCAL_ABS2(x) internal::abs2(x)
#endif
/** This class is a modification of Eigen's ColPivHouseholdeQR to perform a rank-revealing QR with column pivoting
* MP = QR with R*P' directly stored in the input matrix M and the householder vectors essential parts stored in the
* column of a different matrix given by the user.
...
...
@@ -41,14 +53,14 @@ class DestructiveColPivQR
m_q
(
householderEssentialStorage
),
//m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
m_hCoeffs
(
householderEssentialStorage
.
diagonal
()),
m_colsPermutation
(
matrix
.
cols
()),
m_colsTranspositions
(
matrix
.
cols
()),
m_colsIntTranspositions
(
matrix
.
cols
()),
m_temp
(
matrix
.
cols
()),
m_colSqNorms
(
matrix
.
cols
()),
m_colsPermutation
(
(
int
)
matrix
.
cols
()),
m_colsTranspositions
(
(
int
)
matrix
.
cols
()),
m_colsIntTranspositions
(
(
int
)
matrix
.
cols
()),
m_temp
(
(
int
)
matrix
.
cols
()),
m_colSqNorms
(
(
int
)
matrix
.
cols
()),
m_isInitialized
(
false
),
m_usePrescribedEpsilon
((
epsilon
==
0.
)
?
false
:
true
),
m_prescribedEpsilon
(
epsilon
)
m_prescribedEpsilon
(
epsilon
)
{
eigen_assert
(
epsilon
>=
0.
);
compute
();
...
...
@@ -262,7 +274,7 @@ typename MatrixType::RealScalar DestructiveColPivQR<MatrixType, HouseholderStror
{
eigen_assert
(
m_isInitialized
&&
"DestructiveColPivQR is not initialized."
);
eigen_assert
(
m_r
.
rows
()
==
m_r
.
cols
()
&&
"You can't take the determinant of a non-square matrix!"
);
return
internal
::
abs
(
m_r
.
diagonal
().
prod
());
return
LOCAL_ABS
(
m_r
.
diagonal
().
prod
());
}
template
<
typename
MatrixType
,
typename
HouseholderStrorageType
>
...
...
@@ -305,8 +317,7 @@ DestructiveColPivQR<MatrixType, HouseholderStrorageType>& DestructiveColPivQR<Ma
// The threshold should be decided wrt. to the norm of ML, while this class only consider
// the norm of L -> TODO: add an initialization of threshold by EPS*norm(L) ... is it really
// necessary, 'cos it is time consuming.
RealScalar
threshold_helper
=
internal
::
abs2
(
epsilon
());
RealScalar
threshold_helper
=
LOCAL_ABS2
(
epsilon
());
m_nonzero_pivots
=
size
;
// the generic case is that in which all pivots are nonzero (invertible case)
m_maxpivot
=
RealScalar
(
0
);
...
...
@@ -362,7 +373,7 @@ DestructiveColPivQR<MatrixType, HouseholderStrorageType>& DestructiveColPivQR<Ma
m_r
.
col
(
m_colsIntTranspositions
[
k
]).
tail
(
rows
-
k
-
1
).
setZero
();
// remember the maximum absolute value of diagonal coefficients
if
(
internal
::
abs
(
beta
)
>
m_maxpivot
)
m_maxpivot
=
internal
::
abs
(
beta
);
if
(
LOCAL_ABS
(
beta
)
>
m_maxpivot
)
m_maxpivot
=
LOCAL_ABS
(
beta
);
// apply the householder transformation
for
(
Index
l
=
k
+
1
;
l
<
cols
;
++
l
)
...
...
@@ -383,9 +394,9 @@ DestructiveColPivQR<MatrixType, HouseholderStrorageType>& DestructiveColPivQR<Ma
//std::cout << "norms=" << std::endl << m_colSqNorms << std::endl;
}
m_colsPermutation
.
setIdentity
(
cols
);
m_colsPermutation
.
setIdentity
(
(
int
)
cols
);
for
(
Index
k
=
0
;
k
<
m_nonzero_pivots
;
++
k
)
m_colsPermutation
.
applyTranspositionOnTheRight
(
k
,
m_colsTranspositions
.
coeff
(
k
));
m_colsPermutation
.
applyTranspositionOnTheRight
(
(
int
)
k
,
(
int
)
m_colsTranspositions
.
coeff
(
k
));
m_det_pq
=
(
number_of_transpositions
%
2
)
?
-
1
:
1
;
m_isInitialized
=
true
;
...
...
src/HCOD.cpp
View file @
3e1ae6c8
...
...
@@ -130,7 +130,7 @@ namespace soth
{
int
r
=
0
;
for
(
size_t
i
=
0
;
i
<
stages
.
size
();
++
i
)
r
+=
stages
[
i
]
->
rank
();
r
+=
(
int
)
stages
[
i
]
->
rank
();
return
r
;
}
...
...
src/HCOD.hpp
View file @
3e1ae6c8
...
...
@@ -47,7 +47,7 @@ namespace soth
//sizes
int
sizeA
()
const
;
int
rank
()
const
;
unsigned
int
nbStages
()
const
{
return
stages
.
size
();
}
unsigned
int
nbStages
()
const
{
return
(
unsigned
int
)
stages
.
size
();
}
/* --- Decomposition --- */
public:
...
...
src/Stage.cpp
View file @
3e1ae6c8
...
...
@@ -133,7 +133,7 @@ namespace soth
}
cstref_vector_t
Stage
::
getOptimalActiveSet
(
bool
withTwin
)
getOptimalActiveSet
(
bool
)
{
cstref_vector_t
res
(
sizeA
());
int
loop
=
0
;
...
...
@@ -235,7 +235,7 @@ namespace soth
/* Compute ML=J(initIr,:)*Y. */
computeInitialJY
();
assert
(
Yinit
.
getRank
()
>=
0
&&
Yinit
.
getRank
()
<=
(
int
)
nc
);
const
unsigned
int
previousRank
=
Y
.
getRank
();
const
unsigned
int
previousRank
=
(
unsigned
int
)
Y
.
getRank
();
isWIdenty
=
true
;
...
...
@@ -292,7 +292,7 @@ namespace soth
* sizeL = mQR.rank();
*/
assert
(
mQR
.
rank
()
>=
0
&&
mQR
.
rank
()
<=
int
(
sizeL
)
);
const
unsigned
int
rank
=
mQR
.
rank
();
const
unsigned
int
rank
=
(
unsigned
int
)
mQR
.
rank
();
conditionalWinit
(
sizeL
==
rank
);
while
(
sizeL
>
rank
)
...
...
@@ -344,13 +344,13 @@ namespace soth
const
Index
r
=
(
in_r
<
0
)
?
row
-
1
:
in_r
;
for
(
Index
i
=
r
-
1
;
i
>=
0
;
--
i
)
// PSEUDOZEROS
{
Givens
G1
(
L
.
col
(
i
),
i
,
row
);
Givens
G1
(
L
.
col
(
i
),
(
int
)
i
,(
int
)
row
);
G1
.
applyTransposeOnTheRight
(
Mr
);
G1
.
applyTransposeOnTheRight
(
L
,
r
);
G1
.
applyTransposeOnTheRight
(
L
,
(
int
)
r
);
G1
.
applyThisOnTheLeft
(
Wr
);
}
removeARowFromL
(
row
);
removeARowFromL
(
(
unsigned
int
)
row
);
sotDEBUG
(
15
)
<<
"W = "
<<
MATLAB
(
W
,
isWIdenty
)
<<
std
::
endl
;
sotDEBUG
(
15
)
<<
"L = "
<<
(
MATLAB
)
L
<<
std
::
endl
;
...
...
@@ -387,7 +387,8 @@ namespace soth
{
if
(
std
::
abs
(
lambda
[
i
])
>
EPSILON
)
{
const
unsigned
int
cstref
=
activeSet
.
mapInv
(
i
);
const
unsigned
int
cstref
=
(
unsigned
int
)
activeSet
.
mapInv
((
unsigned
int
)
i
);
if
(
!
activeSet
.
isFreezed
(
cstref
)
)
{
activeSet
.
freeze
(
cstref
);
...
...
@@ -529,7 +530,7 @@ namespace soth
sotDEBUG
(
5
)
<<
"W0 = "
<<
MATLAB
(
W
,
isWIdenty
)
<<
std
::
endl
;
sotDEBUG
(
5
)
<<
"M0 = "
<<
(
MATLAB
)
M
<<
std
::
endl
;
sotDEBUG
(
5
)
<<
"L0 = "
<<
(
MATLAB
)
L
<<
std
::
endl
;
L
.
pushColFront
(
M
.
popColBack
()
);
L
.
pushColFront
(
(
int
)
M
.
popColBack
()
);
sizeM
--
;
/* Check if one of the M's grown. */
...
...
@@ -544,7 +545,7 @@ namespace soth
Block
<
MatrixXd
>
ML
(
ML_
,
0
,
0
,
nr
,
sizeM
+
1
);
for
(
Index
j
=
i
+
1
;
j
<
sizeN
();
++
j
)
// PSEUDOZERO
{
Givens
G1
(
Ln
,
i
,
j
);
Givens
G1
(
Ln
,
(
int
)
i
,(
int
)
j
);
G1
.
transpose
()
>>
M
;
G1
.
transpose
()
>>
Ln
;
W
<<
G1
;
...
...
@@ -631,7 +632,7 @@ namespace soth
const
int
rowRankInL
=
col
-
sizeN
();
activeSet
.
unactiveRow
(
row
);
const
unsigned
int
wcoldown
=
W
.
removeCol
(
col
);
const
unsigned
int
wcoldown
=
(
unsigned
int
)
W
.
removeCol
(
col
);
if
(
rowRankInL
>=
0
)
{
L
.
removeRow
(
rowRankInL
);
sizeL
--
;
}
freeML
.
put
(
wcoldown
);
...
...
@@ -711,7 +712,7 @@ namespace soth
/* Choose the first available row in ML. */
const
unsigned
int
wcolup
=
freeML
.
get
();
assert
(
(
wcolup
>=
0
)
&&
(
wcolup
<
nr
)
);
assert
(
wcolup
<
nr
);
sotDEBUG
(
5
)
<<
" wc="
<<
wcolup
<<
endl
;
/* Add a line to ML. */
...
...
@@ -728,7 +729,7 @@ namespace soth
if
(
norm2
>
EPSILON
*
EPSILON
)
{
sotDEBUG
(
45
)
<<
"Oversize value is x"
<<
i
<<
" = "
<<
JupY
(
i
)
<<
std
::
endl
;
rankJ
=
i
+
1
;
rankJ
=
(
unsigned
int
)
i
+
1
;
break
;
}
}
...
...
@@ -738,7 +739,7 @@ namespace soth
{
/* Rank increase: remove the tail of JuY. */
for
(
Index
i
=
rankJ
-
1
;
i
>
int
(
sizeM
+
sizeL
);
--
i
)
{
Givens
G1
(
JupY
,
i
-
1
,
i
,
true
);
Givens
G1
(
JupY
,
(
int
)
i
-
1
,(
int
)
i
,
true
);
Yup
.
push
(
G1
);
}
addARow
(
wcolup
);
...
...
@@ -1394,7 +1395,8 @@ namespace soth
bool
res
=
false
;
EI_FOREACH
(
i
,
lambda
)
{
const
unsigned
int
cstref
=
activeSet
.
mapInv
(
i
);
const
unsigned
int
cstref
=
(
unsigned
int
)
activeSet
.
mapInv
((
unsigned
int
)
i
);
Bound
::
bound_t
btype
=
activeSet
.
whichBound
(
cstref
);
if
(
activeSet
.
isFreezed
(
cstref
)
)
continue
;
...
...
@@ -1403,7 +1405,8 @@ namespace soth
case
Bound
::
BOUND_TWIN
:
break
;
// Nothing to do.
case
Bound
::
BOUND_SUP
:
sotDEBUG
(
5
)
<<
name
<<
": row"
<<
i
<<
", cst"
<<
which
(
i
)
<<
": l="
<<
lambda
(
i
,
0
)
<<
endl
;
sotDEBUG
(
5
)
<<
name
<<
": row"
<<
i
<<
", cst"
<<
which
((
unsigned
int
)
i
)
<<
": l="
<<
lambda
(
i
,
0
)
<<
endl
;
if
(
-
lambda
[
i
]
>
lmax
)
{
// double Ju = J.row(cstref)*u;
...
...
@@ -1411,12 +1414,13 @@ namespace soth
{
res
=
true
;
lmax
=-
lambda
[
i
];
row
=
i
;
row
=
(
unsigned
int
)
i
;
}
}
break
;
case
Bound
::
BOUND_INF
:
sotDEBUG
(
5
)
<<
name
<<
": row"
<<
i
<<
", cst"
<<
which
(
i
)
<<
": l="
<<
lambda
(
i
,
0
)
<<
endl
;
sotDEBUG
(
5
)
<<
name
<<
": row"
<<
i
<<
", cst"
<<
which
((
unsigned
int
)
i
)
<<
": l="
<<
lambda
(
i
,
0
)
<<
endl
;
if
(
-
lambda
[
i
]
>
lmax
)
// TODO: change the sign of the bound-inf cst.
{
// double Ju = J.row(cstref)*u;
...
...
@@ -1424,7 +1428,7 @@ namespace soth
{
res
=
true
;
lmax
=-
lambda
[
i
];
row
=
i
;
row
=
(
unsigned
int
)
i
;
}
}
break
;
...
...
src/Stage.hpp
View file @
3e1ae6c8
...
...
@@ -243,7 +243,7 @@ namespace soth
inline
unsigned
int
nbConstraints
(
void
)
const
{
return
nr
;
}
inline
unsigned
int
sizeA
(
void
)
const
{
return
activeSet
.
nbActive
();
}
// sizeN = card(In) = sizeA-sizeL.
inline
int
sizeN
(
void
)
const
{
assert
((
int
)
sizeA
()
-
sizeL
>=
0
);
return
sizeA
()
-
sizeL
;
}
inline
int
sizeN
(
void
)
const
{
assert
((
int
)
(
sizeA
()
-
sizeL
)
>=
0
);
return
sizeA
()
-
sizeL
;
}
inline
Index
rank
()
const
{
return
sizeL
;}
inline
int
getSizeM
()
const
{
return
sizeM
;
}
...
...
src/SubMatrix.hpp
View file @
3e1ae6c8
...
...
@@ -56,11 +56,11 @@ namespace Eigen
{
Index
size
=
high
-
low
;
if
(
size
>
1
)
return
VectorType
::
LinSpaced
(
size
,
low
,
high
-
1
);
return
VectorType
::
LinSpaced
(
(
int
)
size
,
(
int
)
low
,
(
int
)
high
-
1
);
else
if
(
size
==
0
)
return
VectorType
();
else
return
VectorType
::
Constant
(
1
,
low
);
return
VectorType
::
Constant
(
1
,
(
int
)
low
);
}
};
...
...
@@ -218,7 +218,7 @@ namespace Eigen
eigen_assert
(
j
>=
0
&&
j
<
rows
());
Index
tmp
=
rowIndices
[
i
];
rowIndices
[
i
]
=
rowIndices
[
j
];
rowIndices
[
j
]
=
tmp
;
rowIndices
[
j
]
=
(
int
)
tmp
;
}
void
pushRowFront
(
Index
index
)
{
...
...
@@ -227,14 +227,14 @@ namespace Eigen
rowIndices
.
conservativeResize
(
s
+
1
);
/* Cannot use block for this operation. */
for
(
Index
i
=
s
;
i
>
0
;
--
i
)
{
rowIndices
[
i
]
=
rowIndices
[
i
-
1
];}
rowIndices
[
0
]
=
index
;
rowIndices
[
0
]
=
(
int
)
index
;
}
void
pushRowBack
(
Index
index
)
{
eigen_assert
(
inMRange
(
index
)
);
const
Index
s
=
rows
();
rowIndices
.
conservativeResize
(
s
+
1
);
rowIndices
[
s
]
=
index
;
rowIndices
[
s
]
=
(
int
)
index
;
}
Index
removeRow
(
Index
index
)
{
...
...
@@ -363,7 +363,7 @@ namespace Eigen
eigen_assert
(
j
>=
0
&&
j
<
cols
());
Index
tmp
=
colIndices
[
i
];
colIndices
[
i
]
=
colIndices
[
j
];
colIndices
[
j
]
=
tmp
;
colIndices
[
j
]
=
(
int
)
tmp
;
}
void
pushColFront
(
Index
index
)
{
...
...
@@ -371,14 +371,14 @@ namespace Eigen
const
Index
s
=
cols
();
colIndices
.
conservativeResize
(
s
+
1
);
for
(
Index
i
=
s
;
i
>
0
;
--
i
)
{
colIndices
[
i
]
=
colIndices
[
i
-
1
];}
colIndices
[
0
]
=
index
;
colIndices
[
0
]
=
(
int
)
index
;
}
void
pushColBack
(
Index
index
)
{
eigen_assert
(
inMRange
(
index
)
);
const
Index
s
=
cols
();
colIndices
.
conservativeResize
(
s
+
1
);
colIndices
[
s
]
=
index
;
colIndices
[
s
]
=
(
int
)
index
;
}
Index
removeCol
(
Index
index
)
{
...
...
src/solvers.hpp
View file @
3e1ae6c8
...
...
@@ -22,7 +22,7 @@ namespace soth
assert
(
lhs
.
rows
()
==
lhs
.
cols
());
assert
(
lhs
.
rows
()
>
0
);
assert
(
rhs
.
size
()
==
lhs
.
rows
());
const
int
n
=
lhs
.
rows
();
const
int
n
=
(
int
)
lhs
.
rows
();
for
(
int
i
=
0
;
i
<
n
-
1
;
++
i
)
{
rhs
[
i
]
/=
lhs
(
i
,
i
);
...
...
@@ -47,7 +47,7 @@ namespace soth
assert
(
lhs
.
rows
()
==
lhs
.
cols
());
assert
(
lhs
.
rows
()
>
0
);
assert
(
rhs
.
size
()
==
lhs
.
rows
());
const
int
n
=
lhs
.
rows
();
const
int
n
=
(
int
)
lhs
.
rows
();
for
(
int
i
=
n
-
1
;
i
>
0
;
--
i
)
{
rhs
[
i
]
/=
lhs
(
i
,
i
);
...
...
unitTesting/COD.hpp
View file @
3e1ae6c8
...
...
@@ -78,7 +78,7 @@ namespace soth
void
compute
(
const
MatrixXd
&
A
,
const
int
rank_
=-
1
,
const
double
EPSILON
=
1e-8
)
{
NC
=
A
.
cols
();
NR
=
A
.
rows
();
NC
=
(
int
)
A
.
cols
();
NR
=
(
int
)
A
.
rows
();
#ifdef DEBUG
Ainit
=
A
;
#endif
...
...
unitTesting/tbasic.cpp
View file @
3e1ae6c8
...
...
@@ -20,7 +20,7 @@ void testBasicStage()
{
MatrixXd
m1
(
5
,
4
);
Map
<
MatrixXd
>
map1
(
m1
.
data
(),
m1
.
size
(),
1
);
map1
=
VectorXd
::
LinSpaced
(
m1
.
size
(),
0
,
m1
.
size
()
-
1
);
map1
=
VectorXd
::
LinSpaced
(
(
long
int
)
m1
.
size
(),
0
.0
,
(
double
)
m1
.
size
()
-
1
);
std
::
cout
<<
"m1 = "
<<
m1
<<
endl
;
VectorBound
b1
(
5
);
...
...
unitTesting/tchrono.cpp
View file @
3e1ae6c8
...
...
@@ -37,7 +37,8 @@ namespace soth
};
double
operator
-
(
const
Now
&
t1
,
const
Now
&
t0
)
{
return
(
t1
.
sec
-
t0
.
sec
)
*
1000.0
+
(
t1
.
usec
-
t0
.
usec
+
0.0
)
/
1000.0
;
return
((
double
)(
t1
.
sec
-
t0
.
sec
))
*
1000.0
+
((
double
)(
t1
.
usec
-
t0
.
usec
)
+
0.0
)
/
1000.0
;
}
std
::
ostream
&
operator
<<
(
std
::
ostream
&
os
,
const
Now
&
now
)
{
return
os
<<
now
.
sec
<<
"' "
<<
now
.
usec
;
}
...
...
@@ -95,7 +96,7 @@ int main (int argc, char** argv)
/* Initialize the seed. */
struct
timeval
tv
;
gettimeofday
(
&
tv
,
NULL
);
int
seed
=
tv
.
tv_usec
%
7919
;
//= 7594;
int
seed
=
(
int
)(
tv
.
tv_usec
%
7919
)
;
//= 7594;
if
(
argc
==
2
)
{
seed
=
atoi
(
argv
[
1
]);
}
std
::
cout
<<
"seed = "
<<
seed
<<
std
::
endl
;
...
...
unitTesting/tcod.cpp
View file @
3e1ae6c8
...
...
@@ -60,14 +60,15 @@ int main ( void )
gettimeofday
(
&
t0
,
NULL
);
COD
Acod
;
Acod
.
compute
(
A
,
RANK
);
gettimeofday
(
&
t1
,
NULL
);
double
time
=
(
t1
.
tv_sec
-
t0
.
tv_sec
)
+
(
t1
.
tv_usec
-
t0
.
tv_usec
)
/
1.0e6
;
double
time
=
(
double
)
(
t1
.
tv_sec
-
t0
.
tv_sec
)
+
(
double
)(
t1
.
tv_usec
-
t0
.
tv_usec
)
/
1.0e6
;
totalTime
+=
time
;
gettimeofday
(
&
t0
,
NULL
);
//JacobiSVD<MatrixXd> Asvd(A, ComputeThinU | ComputeThinV);
gettimeofday
(
&
t1
,
NULL
);
double
timeSVD
=
(
t1
.
tv_sec
-
t0
.
tv_sec
)
+
(
t1
.
tv_usec
-
t0
.
tv_usec
)
/
1.0e6
;
double
timeSVD
=
(
double
)(
t1
.
tv_sec
-
t0
.
tv_sec
)
+
(
double
)(
t1
.
tv_usec
-
t0
.
tv_usec
)
/
1.0e6
;
totalTimeSVD
+=
timeSVD
;
if
(
!
(
shoot
%
100
))
std
::
cout
<<
time
*
1000
<<
"
\t
"
<<
timeSVD
*
1000
<<
std
::
endl
;
}
...
...
unitTesting/tdump_problem.cpp
View file @
3e1ae6c8
...
...
@@ -48,7 +48,7 @@ int main (int argc, char** argv)
struct
timeval
tv
;
gettimeofday
(
&
tv
,
NULL
);
int
seed
=
tv
.
tv_usec
%
7919
;
//= 7594;
int
seed
=
(
int
)(
tv
.
tv_usec
%
7919
)
;
//= 7594;
if
(
argc
==
2
)
{
seed
=
atoi
(
argv
[
1
]);
}
std
::
cout
<<
"seed = "
<<
seed
<<
std
::
endl
;
...
...
unitTesting/tgivens.cpp
View file @
3e1ae6c8
...
...
@@ -43,7 +43,7 @@ void testSimpleGivens()
const
double
n1
=
M
.
row
(
0
).
tail
(
M
.
cols
()
-
1
).
norm
();
UNUSED
(
n1
);
sotDEBUG
(
1
)
<<
endl
<<
"Nullify row 0 ..."
<<
endl
;
for
(
int
i
=
M
.
cols
()
-
1
;
i
>
1
;
--
i
)
for
(
int
i
=
(
int
)
M
.
cols
()
-
1
;
i
>
1
;
--
i
)
{
//Givens G(M.row(0), i-1, i);
M
<<
Givens
(
M
.
row
(
0
),
i
-
1
,
i
);
...
...
@@ -90,7 +90,7 @@ void testSubMatrixGivens()
const
double
n1
=
M
.
row
(
0
).
tail
(
M
.
cols
()
-
1
).
norm
();
UNUSED
(
n1
);
sotDEBUG
(
1
)
<<
endl
<<
"Nullify row 0 ..."
<<
endl
;
for
(
int
i
=
M
.
cols
()
-
1
;
i
>
1
;
--
i
)
for
(
int
i
=
(
int
)
M
.
cols
()
-
1
;
i
>
1
;
--
i
)
{
M
<<
Givens
(
M
.
row
(
0
),
i
-
1
,
i
);
//Givens (M.row(0), i-1, i).applyThisOnTheLeft(M);
...
...
@@ -134,7 +134,7 @@ void testSequenceSub()
U
.
push
(
G3
);
Givens
Gr
;
for
(
int
i
=
M
.
cols
()
-
1
;
i
>
1
;
--
i
)
for
(
int
i
=
(
int
)
M
.
cols
()
-
1
;
i
>
1
;
--
i
)
{
Gr
.
makeGivens
(
M
.
row
(
0
),
i
-
1
,
i
);
Gr
.
applyThisOnTheLeft
(
M
);
V
.
push
(
Gr
);
...
...
unitTesting/thcod_full.cpp
View file @
3e1ae6c8
...
...
@@ -104,7 +104,7 @@ clearIteralively( soth::HCOD& hcod )
int
main
(
int
argc
,
char
**
argv
)
int
main
(
int
,
char
**
)
{
bool
exitOk
=
true
;
const
int
executeAll
=
1
;
...
...
@@ -118,7 +118,7 @@ int main (int argc, char** argv)
{
struct
timeval
tv
;
gettimeofday
(
&
tv
,
NULL
);
int
seed
=
tv
.
tv_usec
%
7919
;
int
seed
=
(
int
)(
tv
.
tv_usec
%
7919
)
;
std
::
cout
<<
"seed = "
<<
seed
<<
std
::
endl
;