[frames] | no frames]

# Source Code for Module gmisclib.gpk_lapack

```  1
2  import Numeric
3
4 -class NoSolutionError(ValueError):
5 -        def __init__(self, s):
6                  ValueError.__init__(self, s)
7
8 -def solve(a, b0):
9          """ solves a*x = b."""
10          b = Numeric.asarray(b0, Numeric.Float)
11          if len(b.shape) == 1:
12                  b = Numeric.array((b,), Numeric.Float)
13          # print "b.shape=", b.shape
14          assert a.n == b.shape[1]
15          import lapack_dpb
16          # print "type(a.d)==", type(a.d)
17          # print "a.d=", a.d
18          # print "ravel(a.d)=", Numeric.ravel(a.d)
19          # print "a.ldab=", a.ldab
20          result = lapack_dpb.dpbsv('L', a.n, a.kd, b.shape[0], a.d,
21                                          a.ldab, b, max(1, a.n), 0)
22          # print "result=", result
23          if result['info'] != 0:
24                  raise NoSolutionError, 'Linear system has no solution. Lapack_dpb.dpbsv info code=%d' % result['info']
25          # print "b=", b
26          return b
27
28
29 -def test1():
30          x = sbd(100, 10)
31          x[0,0] = 1
32          assert x[0,0] == 1
33          x[99,99] = 2
34          assert x[99,99] == 2
35          x[10,0] = 10
36          assert x[10,0] == 10
37          assert x[0,10] == 10
38          try:
39                  x[12,1] = 2
40          except IndexError:
41                  pass
42          else:
43                  assert None
44
45
46 -def err(a, b):
47          return Numeric.sum(Numeric.square(a-b))
48
49 -def test2():
50          x = sbd(6,2)
51          b = Numeric.zeros((6,), Numeric.Float) + 10
52          for i in range(6):
53                  x[i,i] = 2
54          x[2,1] = 0.5
55          # print "x=", x
56          y = solve(x, b)
57          y_t = Numeric.array(([5, 4, 4, 5, 5, 5],))
58          # print "y=", y, "y_t=", y_t
59          # print "x=", x
60          assert err(y_t, y) < 1e-9
61
62
63 -def testa():
64          x = sbd(6, 2)
65          x[4,3] = 1
66          x[4,4] = 2
67          x[5,3] = 3
68          x[5,5] = -1
69
70          y = x[:]
71          # print "y=", y
72          assert y[4,3]==1
73          assert y[4,4]==2
74          assert y[5,3]==3
75          assert y[4,5]==0
76          assert y[3,4]==1
77          assert y[3,3]==0
78          assert y[4,5]==0
79          assert y[5,5] == -1
80          y = x[5:6]
81          # print "y=", y
82          assert y[0,5] == -1
83          assert y[0,3] == 3
84          assert y[0,4] == 0
85
86 -def testbdi():
87          x = sbd(6,2)
88          y = Numeric.array([[1, 2], [2, 3]])
89          x.bd_increment(1, y)
90          assert x[1,1]==1
91          assert x[2,2]==3
92          assert x[2,1]==2
93          assert x[1,2]==2
94          assert x[0,0]==0
95          assert x[3,3]==0
96          assert x[3,2]==0
97          z = Numeric.array([[-1, 0, 1], [0, 0, 0], [0, 0, 0]])
98          x.bd_increment(0, z)
99          assert x[0,0]==-1
100          assert x[0, 2]==1
101          assert x[1,1]==1
102          assert x[1, 2]==2
103
104  if __name__ == '__main__':
105          test1()
106          test2()
107          testa()
108          testbdi()
109
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Thu Sep 22 04:25:13 2011 http://epydoc.sourceforge.net