Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
blas.copy(W0['d'], W['d'])
blas.copy(W0['di'], W['di'])
for k in range(len(dims['q'])):
blas.copy(W0['v'][k], W['v'][k])
W['beta'][k] = W0['beta'][k]
for k in range(len(dims['s'])):
blas.copy(W0['r'][k], W['r'][k])
blas.copy(W0['rti'][k], W['rti'][k])
xcopy(x0, x);
ycopy(y0, y);
blas.copy(s0, s); blas.copy(z0, z)
blas.copy(lmbda0, lmbda)
blas.copy(lmbdasq, lmbdasq0)
xcopy(rx0, rx); ycopy(ry0, ry)
resx = math.sqrt(xdot(rx, rx))
blas.copy(rznl0, rznl); blas.copy(rzl0, rzl);
resznl = blas.nrm2(rznl)
relaxed_iters = -1
try: f3 = kktsolver(x, z[:mnl], W)
except ArithmeticError:
singular_kkt_matrix = True
else:
singular_kkt_matrix = True
if singular_kkt_matrix:
sl, zl = s[mnl:], z[mnl:]
ind = dims['l'] + sum(dims['q'])
for m in dims['s']:
misc.symm(sl, m, ind)
def f4(x, y, z, s):
if refinement or DEBUG:
xcopy(x, wx)
ycopy(y, wy)
blas.copy(z, wz)
blas.copy(s, ws)
f4_no_ir(x, y, z, s)
for i in range(refinement):
xcopy(wx, wx2)
ycopy(wy, wy2)
blas.copy(wz, wz2)
blas.copy(ws, ws2)
res(x, y, z, s, wx2, wy2, wz2, ws2)
f4_no_ir(wx2, wy2, wz2, ws2)
xaxpy(wx2, x)
yaxpy(wy2, y)
blas.axpy(wz2, z)
blas.axpy(ws2, s)
if DEBUG:
res(x, y, z, s, wx, wy, wz, ws)
print("KKT residuals:")
print(" 'x': %e" %math.sqrt(xdot(wx, wx)))
print(" 'y': %e" %math.sqrt(ydot(wy, wy)))
print(" 'z': %e" %misc.snrm2(wz, dims, mnl))
print(" 's': %e" %misc.snrm2(ws, dims, mnl))
# Update lambda and scaling.
misc.update_scaling(W, lmbda, ds, dz)
# Unscale s, z (unscaled variables are used only to compute
# feasibility residuals).
blas.copy(lmbda, s, n = mnl + dims['l'] + sum(dims['q']))
ind = mnl + dims['l'] + sum(dims['q'])
ind2 = ind
for m in dims['s']:
blas.scal(0.0, s, offset = ind2)
blas.copy(lmbda, s, offsetx = ind, offsety = ind2, n = m,
incy = m+1)
ind += m
ind2 += m*m
misc.scale(s, W, trans = 'T')
blas.copy(lmbda, z, n = mnl + dims['l'] + sum(dims['q']))
ind = mnl + dims['l'] + sum(dims['q'])
ind2 = ind
for m in dims['s']:
blas.scal(0.0, z, offset = ind2)
blas.copy(lmbda, z, offsetx = ind, offsety = ind2, n = m,
incy = m+1)
ind += m
ind2 += m*m
misc.scale(z, W, inverse = 'I')
blas.scal(1.0/2.0/cc, v)
# v[k] = 1/sqrt(2*(vk0 + 1)) * ( vk + e ), e = [1; 0]
v[0] += 1.0
blas.scal(1.0/math.sqrt(2.0 * v[0]), v)
# To get the scaled variable lambda_k
#
# d = sk0/a + zk0/b + 2*c
# lambda_k = [ c;
# (c + zk0/b)/d * sk1/a + (c + sk0/a)/d * zk1/b ]
# lambda_k *= sqrt(a * b)
lmbda[ind] = cc
dd = 2*cc + s[ind]/aa + z[ind]/bb
blas.copy(s, lmbda, offsetx = ind+1, offsety = ind+1, n = m-1)
blas.scal((cc + z[ind]/bb)/dd/aa, lmbda, n = m-1, offset = ind+1)
blas.axpy(z, lmbda, (cc + s[ind]/aa)/dd/bb, n = m-1, offsetx =
ind+1, offsety = ind+1)
blas.scal(math.sqrt(aa*bb), lmbda, offset = ind, n = m)
ind += m
# For the 's' blocks: compute two lists 'r' and 'rti'.
#
# r[k]' * sk^{-1} * r[k] = diag(lambda_k)^{-1}
# r[k]' * zk * r[k] = diag(lambda_k)
#
# where sk and zk are the entries inds[k] : inds[k+1] of
# s and z, reshaped into symmetric matrices.
#
W['rti'] = [ matrix(0.0, (m,m)) for m in dims['s'] ]
work = matrix(0.0, (max( [0] + dims['s'] )**2, 1))
Ls = matrix(0.0, (max( [0] + dims['s'] )**2, 1))
Lz = matrix(0.0, (max( [0] + dims['s'] )**2, 1))
ind2 = ind
for k in range(len(dims['s'])):
m = dims['s'][k]
r, rti = W['r'][k], W['rti'][k]
# Factor sk = Ls*Ls'; store Ls in ds[inds[k]:inds[k+1]].
blas.copy(s, Ls, offsetx = ind2, n = m**2)
lapack.potrf(Ls, n = m, ldA = m)
# Factor zs[k] = Lz*Lz'; store Lz in dz[inds[k]:inds[k+1]].
blas.copy(z, Lz, offsetx = ind2, n = m**2)
lapack.potrf(Lz, n = m, ldA = m)
# SVD Lz'*Ls = U*diag(lambda_k)*V'. Keep U in work.
for i in range(m):
blas.scal(0.0, Ls, offset = i*m, n = i)
blas.copy(Ls, work, n = m**2)
blas.trmm(Lz, work, transA = 'T', ldA = m, ldB = m, n = m, m = m)
lapack.gesvd(work, lmbda, jobu = 'O', ldA = m, m = m, n = m,
offsetS = ind)
# r = Lz^{-T} * U
blas.copy(work, r, n = m*m)
blas.trsm(Lz, r, transA = 'T', m = m, n = m, ldA = m)
# rti = Lz * U
blas.copy(work, rti, n = m*m)
for k in range(len(dims['s'])):
blas.copy(W['r'][k], W0['r'][k])
blas.copy(W['rti'][k], W0['rti'][k])
xcopy(x, x0); xcopy(dx, dx0)
ycopy(y, y0); ycopy(dy, dy0)
blas.copy(s, s0); blas.copy(z, z0)
blas.copy(ds, ds0)
blas.copy(dz, dz0)
blas.copy(ds2, ds20)
blas.copy(dz2, dz20)
blas.copy(lmbda, lmbda0)
blas.copy(lmbdasq, lmbdasq0)
dsdz0 = dsdz
sigma0, eta0 = sigma, eta
xcopy(rx, rx0); ycopy(ry, ry0)
blas.copy(rznl, rznl0); blas.copy(rzl, rzl0)
relaxed_iters = 1
backtrack = False
elif 0 <= relaxed_iters < MAX_RELAXED_ITERS > 0:
if newphi <= phi0 + ALPHA * step0 * dphi0:
# Relaxed l.s. gives sufficient decrease.
relaxed_iters = 0
else:
# Relaxed line search
relaxed_iters += 1
backtrack = False
elif relaxed_iters == MAX_RELAXED_ITERS > 0:
blas.copy(W['dnl'], W0['dnl'])
blas.copy(W['dnli'], W0['dnli'])
blas.copy(W['d'], W0['d'])
blas.copy(W['di'], W0['di'])
for k in range(len(dims['q'])):
blas.copy(W['v'][k], W0['v'][k])
W0['beta'][k] = W['beta'][k]
for k in range(len(dims['s'])):
blas.copy(W['r'][k], W0['r'][k])
blas.copy(W['rti'][k], W0['rti'][k])
xcopy(x, x0); xcopy(dx, dx0)
ycopy(y, y0); ycopy(dy, dy0)
blas.copy(s, s0); blas.copy(z, z0)
blas.copy(ds, ds0)
blas.copy(dz, dz0)
blas.copy(ds2, ds20)
blas.copy(dz2, dz20)
blas.copy(lmbda, lmbda0)
blas.copy(lmbdasq, lmbdasq0)
dsdz0 = dsdz
sigma0, eta0 = sigma, eta
xcopy(rx, rx0); ycopy(ry, ry0)
blas.copy(rznl, rznl0); blas.copy(rzl, rzl0)
relaxed_iters = 1
backtrack = False
elif 0 <= relaxed_iters < MAX_RELAXED_ITERS > 0:
if newphi <= phi0 + ALPHA * step0 * dphi0:
# Relaxed l.s. gives sufficient decrease.
relaxed_iters = 0
def pack(x, y, dims, mnl = 0, offsetx = 0, offsety = 0):
"""
Copy x to y using packed storage.
The vector x is an element of S, with the 's' components stored in
unpacked storage. On return, x is copied to y with the 's' components
stored in packed storage and the off-diagonal entries scaled by
sqrt(2).
"""
nlq = mnl + dims['l'] + sum(dims['q'])
blas.copy(x, y, n = nlq, offsetx = offsetx, offsety = offsety)
iu, ip = offsetx + nlq, offsety + nlq
for n in dims['s']:
for k in range(n):
blas.copy(x, y, n = n-k, offsetx = iu + k*(n+1), offsety = ip)
y[ip] /= math.sqrt(2)
ip += n-k
iu += n**2
np = sum([ int(n*(n+1)/2) for n in dims['s'] ])
blas.scal(math.sqrt(2.0), y, n = np, offset = offsety+nlq)
# [ ux[0] ] [ bx[0] + bx[1] * gradf0 ]
# [ uy ] = K^-1 * [ by ].
# [ uz[1:] ] [ bz[1:] ]
#
# ux[1] = gradf0' * ux[0] - W['dnl'][0]**2 * uz[0] - bz[0]
# = gradf0' * ux[0] + W['dnl'][0]**2 * bx[1] - bz[0].
#
# Instead of uz we return the scaled solution W*uz.
a = z[0]
xcopy(x[0], ux)
xaxpy(gradf0, ux, alpha = x[1])
blas.copy(z, uz, offsetx = 1)
g(ux, y, uz)
z[0] = -x[1] * W['dnl'][0]
blas.copy(uz, z, offsety = 1)
xcopy(ux, x[0])
x[1] = xdot(gradf0, x[0]) + W['dnl'][0]**2 * x[1] - a