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