from cvxopt import solvers, matrix, mul, normal, setseed from cvxopt.lapack import * from cvxopt.blas import * from lsqr 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) 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 = lsqr(Gfun, n, b, show=show, atol=tol, btol=tol,itnlim=maxit) xo =X[0] phio=X[3] psio=X[7] k =X[2] chio=X[8] mg=max(G-G.T) if mg>1e-14:sym='No' else: sym='Yes' alg='LSQR' 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 ||Ar|| %9.4e " %( chio, phio, psio) print "Residual norms computed directly:" print " ||x|| %9.4e ||r|| %9.4e ||Ar|| %9.4e" % (nrm2(xo), nrm2(G*xo -b), nrm2(G.T*(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 ""