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 ""