Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
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)
blas.trmm(Lz, rti, m = m, n = m, ldA = m)
# r := r * diag(sqrt(lambda_k))
# rti := rti * diag(1 ./ sqrt(lambda_k))
for i in range(m):
a = math.sqrt( lmbda[ind+i] )
blas.scal(a, r, offset = m*i, n = m)
blas.scal(1.0/a, rti, offset = m*i, n = m)
ind += m
ind2 += m*m
return W
else:
fA = A
if not customy:
if b is None: b = matrix(0.0, (0,1))
if type(b) is not matrix or b.typecode != 'd' or b.size[1] != 1:
raise TypeError("'b' must be a 'd' matrix with one column")
if not operatorA and b.size[0] != A.size[0]:
raise TypeError("'b' must have length %d" %A.size[0])
if b is None and customy:
raise ValueEror("use of non vector type for y requires b")
if xnewcopy is None: xnewcopy = matrix
if xdot is None: xdot = blas.dot
if xaxpy is None: xaxpy = blas.axpy
if xscal is None: xscal = blas.scal
def xcopy(x, y):
xscal(0.0, y)
xaxpy(x, y)
if ynewcopy is None: ynewcopy = matrix
if ydot is None: ydot = blas.dot
if yaxpy is None: yaxpy = blas.axpy
if yscal is None: yscal = blas.scal
def ycopy(x, y):
yscal(0.0, y)
yaxpy(x, y)
# The problem is solved by applying cpl() to the epigraph form
#
# minimize t
# subject to f0(x) - t <= 0
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)
blas.scal(step, sigs)
blas.scal(step, sigz)
sigs += 1.0
sigz += 1.0
blas.tbsv(lmbda, sigs, n = sum(dims['s']), k = 0, ldA = 1,
offsetA = mnl + dims['l'] + sum(dims['q']) )
blas.tbsv(lmbda, sigz, n = sum(dims['s']), k = 0, ldA = 1,
offsetA = mnl + dims['l'] + sum(dims['q']) )
# dsk := Ls = dsk * sqrt(sigs).
# dzk := Lz = dzk * sqrt(sigz).
ind2, ind3 = mnl + dims['l'] + sum(dims['q']), 0
for k in range(len(dims['s'])):
m = dims['s'][k]
for i in range(m):
blas.scal(math.sqrt(sigs[ind3+i]), ds, offset = ind2 + m*i,
n = m)
blas.scal(math.sqrt(sigz[ind3+i]), dz, offset = ind2 + m*i,
n = m)
ind2 += m*m
ind3 += m
# 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']))
# where beta = W['beta'][k], v = W['v'][k], J = [1, 0; 0, -I].
#
# Inverse scaling is
#
# xk := 1/beta * (2*J*v*v'*J - J) * xk
# = 1/beta * (-J) * (2*v*((-J*xk)'*v)' + xk).
w = matrix(0.0, (x.size[1], 1))
for k in range(len(W['v'])):
v = W['v'][k]
m = v.size[0]
if inverse == 'I':
blas.scal(-1.0, x, offset = ind, inc = x.size[0])
blas.gemv(x, v, w, trans = 'T', m = m, n = x.size[1], offsetA =
ind, ldA = x.size[0])
blas.scal(-1.0, x, offset = ind, inc = x.size[0])
blas.ger(v, w, x, alpha = 2.0, m = m, n = x.size[1], ldA =
x.size[0], offsetA = ind)
if inverse == 'I':
blas.scal(-1.0, x, offset = ind, inc = x.size[0])
a = 1.0 / W['beta'][k]
else:
a = W['beta'][k]
for i in range(x.size[1]):
blas.scal(a, x, n = m, offset = ind + i*x.size[0])
ind += m
# Scaling for 's' component xk is
#
# xk := vec( r' * mat(xk) * r ) if trans = 'N'
# xk := vec( r * mat(xk) * r' ) if trans = 'T'.
# xk := beta * (2*v*v' - J) * xk
# = beta * (2*v*(xk'*v)' - J*xk)
#
# where beta = W['beta'][k], v = W['v'][k], J = [1, 0; 0, -I].
#
# Inverse scaling is
#
# xk := 1/beta * (2*J*v*v'*J - J) * xk
# = 1/beta * (-J) * (2*v*((-J*xk)'*v)' + xk).
w = matrix(0.0, (x.size[1], 1))
for k in range(len(W['v'])):
v = W['v'][k]
m = v.size[0]
if inverse == 'I':
blas.scal(-1.0, x, offset=ind, inc=x.size[0])
blas.gemv(x, v, w, trans='T', m=m, n=x.size[1], offsetA=ind, ldA=x.size[0])
blas.scal(-1.0, x, offset=ind, inc=x.size[0])
blas.ger(v, w, x, alpha=2.0, m=m, n=x.size[1], ldA=x.size[0], offsetA=ind)
if inverse == 'I':
blas.scal(-1.0, x, offset=ind, inc=x.size[0])
a = 1.0 / W['beta'][k]
else:
a = W['beta'][k]
for i in range(x.size[1]):
blas.scal(a, x, n=m, offset=ind + i * x.size[0])
ind += m
# Scaling for 's' component xk is
#
# xk := vec( r' * mat(xk) * r ) if trans = 'N'
# xk := vec( r * mat(xk) * r' ) if trans = 'T'.
blas.scal(step, sigs)
blas.scal(step, sigz)
sigs += 1.0
sigz += 1.0
blas.tbsv(lmbda, sigs, n=sum(dims['s']), k=0, ldA=1, offsetA=dims['l'] + sum(dims['q']))
blas.tbsv(lmbda, sigz, n=sum(dims['s']), k=0, ldA=1, offsetA=dims['l'] + sum(dims['q']))
# dsk := Ls = dsk * sqrt(sigs).
# dzk := Lz = dzk * sqrt(sigz).
ind2, ind3 = dims['l'] + sum(dims['q']), 0
for k in range(len(dims['s'])):
m = dims['s'][k]
for i in range(m):
blas.scal(math.sqrt(sigs[ind3 + i]), ds, offset=ind2 + m * i,
n=m)
blas.scal(math.sqrt(sigz[ind3 + i]), dz, offset=ind2 + m * i,
n=m)
ind2 += m * m
ind3 += m
# 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=dims['l'] + sum(dims['q']))
ind = 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,
in unpacked storage and off-diagonal entries scaled by sqrt(2).
On return, x is copied to y with the 's' components stored in
unpacked storage.
"""
nlq = mnl + dims['l'] + sum(dims['q'])
blas.copy(x, y, n = nlq, offsetx = offsetx, offsety = offsety)
iu, ip = offsety+nlq, offsetx+nlq
for n in dims['s']:
for k in range(n):
blas.copy(x, y, n = n-k, offsetx = ip, offsety = iu+k*(n+1))
y[iu+k*(n+1)] *= math.sqrt(2)
ip += n-k
iu += n**2
nu = sum([ n**2 for n in dims['s'] ])
blas.scal(1.0/math.sqrt(2.0), y, n = nu, offset = offsety+nlq)
def G(u, v, alpha = 1.0, beta = 0.0, trans = 'N'):
"""
v := alpha*[I, -I; -I, -I] * u + beta * v (trans = 'N' or 'T')
"""
blas.scal(beta, v)
blas.axpy(u, v, n = n, alpha = alpha)
blas.axpy(u, v, n = n, alpha = -alpha, offsetx = n)
blas.axpy(u, v, n = n, alpha = -alpha, offsety = n)
blas.axpy(u, v, n = n, alpha = -alpha, offsetx = n, offsety = n)
# Solve
#
# [ 0 ] [ P A' G' ] [ dx ]
# [ 0 ] + [ A 0 0 ] * [ dy ] = -(1 - eta) * r
# [ W'*ds ] [ G 0 0 ] [ W^{-1}*dz ]
#
# lmbda o (dz + ds) = -lmbda o lmbda + sigma*mu*e (i=0)
# lmbda o (dz + ds) = -lmbda o lmbda - dsa o dza
# + sigma*mu*e (i=1) where dsa, dza
# are the solution for i=0.
# ds = -lmbdasq + sigma * mu * e (if i is 0)
# = -lmbdasq - dsa o dza + sigma * mu * e (if i is 1),
# where ds, dz are solution for i is 0.
blas.scal(0.0, ds)
if i == 1:
blas.axpy(ws3, ds, alpha=-1.0)
blas.axpy(lmbdasq, ds, n=dims['l'] + sum(dims['q']),
alpha=-1.0)
ds[:dims['l']] += sigma * mu
ind = dims['l']
for m in dims['q']:
ds[ind] += sigma * mu
ind += m
ind2 = ind
for m in dims['s']:
blas.axpy(lmbdasq, ds, n=m, offsetx=ind2, offsety=ind, incy=m + 1, alpha=-1.0)
ds[ind: ind + m * m: m + 1] += sigma * mu
ind += m * m
ind2 += m