from cvxopt import solvers, matrix, mul, normal, setseed from cvxopt.lapack import * from cvxopt.blas import * from symmlq import * import time as t import cvxopt.misc as misc import numpy as np setseed(2) n=35 G=matrix(np.eye(n), tc='d') for jj in range(5): gg=normal(n,1) hh=gg*gg.T G+=(hh+hh.T)*0.5 G+=normal(n,1)*normal(1,n) G=(G+G.T)/2 b=normal(n,1) def Gfun(x,y,trans='N'): gemv(G,x,y,trans) svx=+b gesv(G,svx) tol=1e-10 show=False maxit=None t1=t.time() X = symmlq( b, Gfun, show=show, rtol=tol, maxit=maxit) xo =X[0] phio=X[5] k =X[2] chio=X[6] mg=max(G-G.T) if mg>1e-14:sym='No' else: sym='Yes' alg='SYMMLQ' print alg print "Is linear operator symmetric? " + sym print "n: %3g iterations: %3g" % (n, k) print " norms computed in ", alg print " ||x|| %9.4e ||r|| %9.4e " %( chio, phio ) print "Residual norms computed directly:" print " ||x|| %9.4e ||r|| %9.4e" % (nrm2(xo), nrm2(G*xo -b)) print "Direct solution norms:" print " ||x|| %9.4e ||r|| %9.4e " % (nrm2(svx), nrm2(G*svx -b)) print "" print " || x_{direct} - x_{"+alg+"}|| %9.4e " % nrm2(svx-xo) print ""