Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
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
# (dx, dy, dz) := -(1-eta) * (rx, ry, rz)
xscal(0.0, dx); xaxpy(rx, dx, alpha = -1.0 + eta)
yscal(0.0, dy); yaxpy(ry, dy, alpha = -1.0 + eta)
blas.scal(0.0, dz)
blas.axpy(rznl, dz, alpha = -1.0 + eta)
blas.axpy(rzl, dz, alpha = -1.0 + eta, offsety = mnl)
try: f4(dx, dy, dz, ds)
except ArithmeticError:
if iters == 0:
raise ValueError("Rank(A) < p or "\
"Rank([H(x); A; Df(x); G]) < n")
else:
sl, zl = s[mnl:], z[mnl:]
ind = dims['l'] + sum(dims['q'])
for m in dims['s']:
misc.symm(sl, m, ind)
misc.symm(zl, m, ind)
ind += m**2
ts = misc.max_step(s, dims, mnl)
tz = misc.max_step(z, dims, mnl)
# 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
# (dx, dy, dz) := -(1 - eta) * (rx, ry, rz)
dx *= 0.0
dx += (-1.0 + eta) * rx
dy *= 0.0
dy += (-1.0 + eta) * ry
blas.scal(0.0, dz)
blas.axpy(rz, dz, alpha=-1.0 + eta)
try:
f4(dx, dy, dz, ds)
except ArithmeticError:
# rx = c + A'*y + Df'*z[:mnl] + G'*z[mnl:]
xcopy(c, rx)
fA(y, rx, beta = 1.0, trans = 'T')
fDf(z[:mnl], rx, beta = 1.0, trans = 'T')
fG(z[mnl:], rx, beta = 1.0, trans = 'T')
resx = math.sqrt(xdot(rx, rx))
# ry = A*x - b
ycopy(b, ry)
fA(x, ry, alpha = 1.0, beta = -1.0)
resy = math.sqrt(ydot(ry, ry))
# rznl = s[:mnl] + f
blas.copy(s[:mnl], rznl)
blas.axpy(f, rznl)
resznl = blas.nrm2(rznl)
# rzl = s[mnl:] + G*x - h
blas.copy(s[mnl:], rzl)
blas.axpy(h, rzl, alpha = -1.0)
fG(x, rzl, beta = 1.0)
reszl = misc.snrm2(rzl, dims)
# Statistics for stopping criteria.
# pcost = c'*x
# dcost = c'*x + y'*(A*x-b) + znl'*f(x) + zl'*(G*x-h)
# = c'*x + y'*(A*x-b) + znl'*(f(x)+snl) + zl'*(G*x-h+sl)
# - z'*s
# = c'*x + y'*ry + znl'*rznl + zl'*rzl - gap
pcost = xdot(c,x)
# [ bx + GG' * W^{-1} * W^{-T} * bz ]
# = [ ]
# [ by ]
#
# and return x, y, W*z = W^{-T} * (GG*x - bz).
blas.copy(z, g)
scale(g, W, trans = 'T', inverse = 'I')
scale(g, W, inverse = 'I')
if mnl:
base.gemv(Df, g, u, trans = 'T')
beta = 1.0
else:
beta = 0.0
sgemv(G, g, u, dims, trans = 'T', offsetx = mnl, beta = beta)
blas.axpy(x, u)
blas.copy(y, u, offsety = n)
if p: lapack.sytrs(K, ipiv, u)
else: lapack.potrs(K, u)
blas.copy(u, x, n = n)
blas.copy(u, y, offsetx = n, n = p)
if mnl: base.gemv(Df, x, z, alpha = 1.0, beta = -1.0)
sgemv(G, x, z, dims, alpha = 1.0, beta = -1.0, offsety = mnl)
scale(z, W, trans = 'T', inverse = 'I')
blas.scal(0.0, u, offset = n-p)
# x[p:] := R3^{-1} * u[:n-p]
blas.copy(u, x, offsety = p, n = n-p)
lapack.trtrs(Gs, x, uplo='U', n = n-p, offsetA = Gs.size[0]*p,
offsetB = p)
# x is now [ R1^{-T}*by; R3^{-1}*u[:n-p] ]
# x := [Q1 Q2]*x
lapack.ormqr(QA, tauA, x)
# u := [Q3, Q4] * u - w
# = Q3 * u[:n-p] - w
lapack.ormqr(Gs, tauG, u, k = n-p, m = cdim_pckd, offsetA =
Gs.size[0]*p)
blas.axpy(w, u, alpha = -1.0)
# y := R1^{-1} * ( v[:p] - Gs1'*u )
# = R1^{-1} * ( Q1'*bx - Gs1'*u )
blas.copy(vv, y, n = p)
blas.gemv(Gs, u, y, m = cdim_pckd, n = p, trans = 'T', alpha =
-1.0, beta = 1.0)
lapack.trtrs(QA, y, uplo = 'U', n=p)
unpack(u, z, dims)
if refinement:
wx = matrix(np.copy(x))
wy = matrix(np.copy(y))
wz = matrix(np.copy(z))
ws = matrix(np.copy(s))
f4_no_ir(x, y, z, s)
for i in range(refinement):
wx2 = matrix(np.copy(wx))
wy2 = matrix(np.copy(wy))
wz2 = matrix(np.copy(wz))
ws2 = matrix(np.copy(ws))
res(x, y, z, s, wx2, wy2, wz2, ws2, W, lmbda)
f4_no_ir(wx2, wy2, wz2, ws2)
y += wx2
y += wy2
blas.axpy(wz2, z)
blas.axpy(ws2, s)
# s := lmbda o\ s
# = lmbda o\ bs
misc.sinv(s, lmbda, dims, mnl)
# z := z - W'*s
# = bz - W' * (lambda o\ bs)
blas.copy(s, ws3)
misc.scale(ws3, W, trans = 'T')
blas.axpy(ws3, z, alpha = -1.0)
# Solve for ux, uy, uz
f3(x, y, z)
# s := s - z
# = lambda o\ bs - z.
blas.axpy(z, s, alpha = -1.0)
# gradfi := Fi' * yi
# = Fi' * exp(Fi*x+gi) / sum(exp(Fi*x+gi))
base.gemv(F, y, Df, trans='T', m=stop-start, incy=mnl+1,
offsetA=start, offsetx=start, offsety=i)
if z is not None:
# Hi = Fi' * (diag(yi) - yi*yi') * Fi
# = Fisc' * Fisc
# where
# Fisc = diag(yi)^1/2 * (I - 1*yi') * Fi
# = diag(yi)^1/2 * (Fi - 1*gradfi')
Fsc[:K[i], :] = F[start:stop, :]
for k in range(start,stop):
blas.axpy(Df, Fsc, n=n, alpha=-1.0, incx=mnl+1,
incy=Fsc.size[0], offsetx=i, offsety=k-start)
blas.scal(math.sqrt(y[k]), Fsc, inc=Fsc.size[0],
offset=k-start)
# H += z[i]*Hi = z[i] * Fisc' * Fisc
blas.syrk(Fsc, H, trans='T', k=stop-start, alpha=z[i],
beta=1.0)
if z is None: return f, Df
else: return f, Df, H
wx = matrix(np.copy(x))
wy = matrix(np.copy(y))
wz = matrix(np.copy(z))
ws = matrix(np.copy(s))
f4_no_ir(x, y, z, s)
for i in range(refinement):
wx2 = matrix(np.copy(wx))
wy2 = matrix(np.copy(wy))
wz2 = matrix(np.copy(wz))
ws2 = matrix(np.copy(ws))
res(x, y, z, s, wx2, wy2, wz2, ws2, W, lmbda)
f4_no_ir(wx2, wy2, wz2, ws2)
y += wx2
y += wy2
blas.axpy(wz2, z)
blas.axpy(ws2, s)