4ce66b6cbb49700e093f1b9c25a433eb4df1e634
[stefan/supreme.git] / supreme / resolve / tests / test_lsqr.py
1 import numpy as np
2
3 from supreme.resolve.lsqr import lsqr
4 from time import time
5
6 # Set up a test problem
7 n = 35
8 G = np.eye(n)
9 normal = np.random.normal
10 norm = np.linalg.norm
11
12 for jj in range(5):
13     gg = normal(size=n)
14     hh = gg * gg.T
15     G += (hh + hh.T) * 0.5
16     G += normal(size=n) * normal(size=n)
17
18 b = normal(size=n)
19
20 tol = 1e-10
21 show = False
22 maxit = None
23
24 def test_basic():
25     svx = np.linalg.solve(G, b)
26     X = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit)
27     xo = X[0]
28     assert norm(svx - xo) < 1e-5
29
30 if __name__ == "__main__":
31     svx = np.linalg.solve(G, b)
32
33     tic = time()
34     X = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit)
35     xo = X[0]
36     phio = X[3]
37     psio = X[7]
38     k = X[2]
39     chio = X[8]
40     mg = np.amax(G - G.T)
41     if mg > 1e-14:
42         sym='No'
43     else:
44         sym='Yes'
45
46     print 'LSQR'
47     print "Is linear operator symmetric? " + sym
48     print "n: %3g  iterations:   %3g" % (n, k)
49     print "Norms computed in %.2fs by LSQR" % (time() - tic)
50     print " ||x||  %9.4e  ||r|| %9.4e  ||Ar||  %9.4e " %( chio, phio, psio)
51     print "Residual norms computed directly:"
52     print " ||x||  %9.4e  ||r|| %9.4e  ||Ar||  %9.4e" %  (norm(xo),
53                                                           norm(G*xo - b),
54                                                           norm(G.T*(G*xo-b)))
55     print "Direct solution norms:"
56     print " ||x||  %9.4e  ||r|| %9.4e " %  (norm(svx), norm(G*svx -b))
57     print ""
58     print " || x_{direct} - x_{LSQR}|| %9.4e " % norm(svx-xo)
59     print ""