[Scipy-tickets] [SciPy] #1766: Unify scipy.linalg.solve_banded and scipy.sparse.dia_matrix

SciPy Trac scipy-tickets@scipy....
Sun Nov 25 15:07:07 CST 2012


#1766: Unify scipy.linalg.solve_banded and scipy.sparse.dia_matrix
---------------------------------------+------------------------------------
 Reporter:  johntyree                  |       Owner:  somebody   
     Type:  enhancement                |      Status:  new        
 Priority:  normal                     |   Milestone:  Unscheduled
Component:  scipy.sparse               |     Version:  0.11.0     
 Keywords:  sparse, linalg,dia_matrix  |  
---------------------------------------+------------------------------------
Changes (by johntyree):

  * keywords:  sparse, linalg => sparse, linalg,dia_matrix
  * component:  Other => scipy.sparse


Comment:

 To make this even worse, it looks like if you contruct your dia_matrix
 from a dense matrix, it puts the data in the wrong order for solve_banded.
 It appears that it automatically sorts the offsets field, but solve banded
 is expected it to be sorted in reverse.



 {{{
 c
 Out[175]:
 array([[ 0.,  0.,  1.,  1.,  1.,  1.,  1.,  1.],
        [ 0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.],
        [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.],
        [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.],
        [ 1.,  1.,  1.,  1.,  1.,  1.,  0.,  0.]])

 cd = sps.dia_matrix((c, [2,1,0,-1,-2]), shape=(8,8))

 cd.todense()
 Out[177]:
 matrix([[ 1.,  1.,  1.,  0.,  0.,  0.,  0.,  0.],
         [ 1.,  1.,  1.,  1.,  0.,  0.,  0.,  0.],
         [ 1.,  1.,  1.,  1.,  1.,  0.,  0.,  0.],
         [ 0.,  1.,  1.,  1.,  1.,  1.,  0.,  0.],
         [ 0.,  0.,  1.,  1.,  1.,  1.,  1.,  0.],
         [ 0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.],
         [ 0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.],
         [ 0.,  0.,  0.,  0.,  0.,  1.,  1.,  1.]])

 sps.dia_matrix(cd.todense()).offsets
 Out[178]: array([-2, -1,  0,  1,  2], dtype=int32)

 sps.dia_matrix(cd.todense()).data
 Out[179]:
 array([[ 1.,  1.,  1.,  1.,  1.,  1.,  0.,  0.],
        [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.],
        [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.],
        [ 0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.],
        [ 0.,  0.,  1.,  1.,  1.,  1.,  1.,  1.]])
 }}}

 Then from scipy.linalg.solve_banded's docstring we have


 {{{
 The matrix a is stored in ab using the matrix diagonal ordered form::

     ab[u + i - j, j] == a[i,j]

 Example of ab (shape of a is (6,6), u=1, l=2)::

     *    a01  a12  a23  a34  a45
     a00  a11  a22  a33  a44  a55
     a10  a21  a32  a43  a54   *
     a20  a31  a42  a53   *    *
 }}}

 So the diagonals from dia_matrix are stored in reverse order.

-- 
Ticket URL: <http://projects.scipy.org/scipy/ticket/1766#comment:1>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.


More information about the Scipy-tickets mailing list