Skip to content

Updated cvx work vectors #3179

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
victorfors wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Updated cvx work vectors #3179

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4082e95
Added test cases for convexify with larger strongly connecter matrice…
victorfors May 5, 2023
c838922
created work vectors for householder values
victorfors May 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 7 additions & 1 deletion casadi/core/convexify.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -271,7 +271,13 @@ namespace casadi {

d.sz_w = 0;
if (d.config.strategy==CVX_EIGEN_REFLECT || d.config.strategy==CVX_EIGEN_CLIP) {
d.sz_w = std::max(block_size, 2*(block_size-1)*d.config.max_iter_eig);
// Work vector size needed for Givens rotations (c and s)
casadi_int sz_w_givens = 2*(block_size-1)*d.config.max_iter_eig;
// Work vector size needed for house holder transformations (beta and p)
casadi_int sz_w_house = block_size-2 + block_size;
// Work vector size needed for casadi_cvx
d.sz_w = sz_w_givens + sz_w_house;
// Additional work vectors
if (d.config.Hsp_project) d.sz_w = std::max(d.sz_w, Hsp.size1());
if (d.config.scc_transform) d.sz_w += block_size*block_size;
if (inplace) d.sz_w = std::max(d.sz_w, Hsp.size1()+d.Hrsp.nnz());
Expand Down
33 changes: 20 additions & 13 deletions casadi/core/runtime/casadi_cvx.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,12 +118,12 @@ void casadi_cvx_house_apply_symm(casadi_int n, casadi_int k, T1* A, T1* p, T1* v
// Upper triangular part contains compact Housholder factorisations
//
// A: n-by-n dense
// beta: work vector; length n-2
// p: work vector; length n
//
// Reference: Golub & Van Loan, Alg. 8.3.1
template<typename T1>
void casadi_cvx_tri(T1* A, casadi_int n, T1* beta, T1* p) {
T1 pp[1000];
casadi_int k, N;
T1 *A_base, *v;
for (k=0;k<n-2;++k) {
Expand All @@ -138,7 +138,7 @@ void casadi_cvx_tri(T1* A, casadi_int n, T1* beta, T1* p) {
// Assign 2-norm
*A_base = casadi_cvx_house(v, &beta[k], N);

casadi_cvx_house_apply_symm(n, k, A, pp, v, beta[k]);
casadi_cvx_house_apply_symm(n, k, A, p, v, beta[k]);

}
}
Expand Down Expand Up @@ -222,8 +222,8 @@ void casadi_cvx_implicit_qr(casadi_int n, T1* t_diag, T1* t_off, T1* cs) {
//
// tolerance greater than machine precision
//
// trace_meta: length 1+3*n_iter
// trace: length 2*(n-1)*n_iter
// trace_meta: length 1+3*max_iter
// trace: length 2*(n-1)*max_iter
//
/// Golub & Van Loan Alg. 8.3.3
template<typename T1>
Expand Down Expand Up @@ -375,8 +375,8 @@ T1 casadi_cvx_scalar(T1 epsilon, casadi_int reflect, T1 eig) {
// SYMBOL "cvx"
// Convexify a dense symmetric Hessian
//
// w real work vector: length max(n,2*(n-1)*n_iter)
// iw integer work vector: 1+3*n_iter
// w real work vector: length 2*(n-1)*max_iter + n + n-2
// iw integer work vector: 1+3*max_iter
//
// tol: tolerance for symmetric schur
// epsilon: minimum magnitude of eigenvalues
Expand All @@ -387,8 +387,15 @@ int casadi_cvx(casadi_int n, T1 *A, T1 epsilon, T1 tol, casadi_int reflect, casa
casadi_int i, j, k, n_iter, nn, p, trace_offset;
casadi_int *t_meta;
T1 c, s, t_off0;
T1 *cs, *t_diag, *t_off;
T1 beta[100];
T1 *cs, *w_trace, *w_beta, *w_p, *t_diag, *t_off;

// Work vector for Givens rotations: length 2*(n-1)*max_iter
w_trace = w;
// Work vectors for householder transform
// Length n
w_p = w_trace + 2*(n-1)*max_iter;
// Length n-2
w_beta = w_p + n;

// Short-circuit for empty matrices
if (n==0) return 0;
Expand All @@ -399,7 +406,7 @@ int casadi_cvx(casadi_int n, T1 *A, T1 epsilon, T1 tol, casadi_int reflect, casa
return 0;
}

casadi_cvx_tri(A, n, beta, w);
casadi_cvx_tri(A, n, w_beta, w_p);

for (i=0;i<n;++i) {
for (j=0;j<n;++j) {
Expand All @@ -422,7 +429,7 @@ int casadi_cvx(casadi_int n, T1 *A, T1 epsilon, T1 tol, casadi_int reflect, casa
}

// Diagonalize matrix by Symmetric QR
if (casadi_cvx_symm_schur(n, t_diag, t_off, tol, max_iter, iw, w)) return 1;
if (casadi_cvx_symm_schur(n, t_diag, t_off, tol, max_iter, iw, w_trace)) return 1;

// Retain diagonals (eigenvalues)
for (i=0;i<n;++i) {
Expand All @@ -442,7 +449,7 @@ int casadi_cvx(casadi_int n, T1 *A, T1 epsilon, T1 tol, casadi_int reflect, casa
nn = *t_meta++;
p = *t_meta++;
trace_offset = *t_meta++;
cs = w+trace_offset;
cs = w_trace+trace_offset;
t_meta-= 6;
for (j=0;j<nn-1;j++) {
s = *--cs;
Expand All @@ -455,8 +462,8 @@ int casadi_cvx(casadi_int n, T1 *A, T1 epsilon, T1 tol, casadi_int reflect, casa
for (k = n-3; k>=0; --k) {
casadi_int N = n-k-1;
T1 *v = A+N*n;
casadi_cvx_house_apply_symm(n, k, A, w, v, beta[k]);
casadi_cvx_house_apply(k+1, N, n, A+k+1, w, v, beta[k]);
casadi_cvx_house_apply_symm(n, k, A, w_p, v, w_beta[k]);
casadi_cvx_house_apply(k+1, N, n, A+k+1, w_p, v, w_beta[k]);
}

return 0;
Expand Down
44 changes: 44 additions & 0 deletions test/python/mx.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3032,6 +3032,50 @@ def test_convexify_bugs(self):
Ac = evalf(convexify(A,{"strategy":"eigen-reflect"}))
self.checkarray(A,Ac,digits=8)

def test_convexify_buffer_overflow(self):
# Size of matrix to convexify
ns = [15, 105]
for n in ns:
# Use hessian function to create a symmetric connected nxn matrix
z = MX.sym("z", n)
f = z[2:]*z[:-2]*z[1:-1]
j = sumsqr(z)
# Keep almost diagonal to reduce number of iterations required for convegence
H, _ = hessian(j+1e-10*sum1(f),z)
H_convex = convexify(H,{"strategy":"eigen-reflect"})

# Setup casadi function and input values
func = Function("convexify_hessian", [z], [H_convex])
z_vals = np.sin(np.linspace(0,n,n))
# Evaluate to check for runtime errors
func(z_vals)

# Code generate casadi function with -fstack-protector-all and -D_FORTIFY_SOURCE=2 to increase chance to catch some overflow and memory bugs
options = ["-fstack-protector-all", "-D_FORTIFY_SOURCE=2"]
self.check_codegen(func,inputs=[z_vals], extra_options=options)

@slow()
@memory_heavy()
@known_bug()
def test_code_generage_convexify_with_main(self):
# Size of matrix to convexify
n = 1005
# Use hessian function to create a symmetric connected nxn matrix
z = MX.sym("z", n)
f = z[2:]*z[:-2]*z[1:-1]
j = sumsqr(z)
# Keep almost diagonal to reduce number of iterations required for convegence
H, _ = hessian(j+1e-10*sum1(f),z)
H_convex = convexify(H,{"strategy":"eigen-reflect","max_iter_eig":1000})

# Setup casadi function and input values
func = Function("convexify_hessian", [z], [H_convex])
z_vals = np.sin(np.linspace(0,n,n))

# Code generate casadi function with -fstack-protector-all and -D_FORTIFY_SOURCE=2 to increase chance to catch some overflow and memory bugs
options = ["-fstack-protector-all", "-D_FORTIFY_SOURCE=2"]
self.check_codegen(func,inputs=[z_vals], extra_options=options) # Works but is slow
self.check_codegen(func,inputs=[z_vals], main=True, extra_options=options) # main=True results in segmentation fault

def test_logsumexp(self):
x = MX.sym("x",3)
Expand Down