[SciPy-User] Matrix indexing and updating

[Luca Manini]

    Hi everybody,

    I'm new to this list and here comes my first question about matrix
    indexing and updating (in scipy, of course).

    I'm starting writing some code (for my wife) to solve finite elements
    problems.  One typical action is to access (both for reading and
    writing) som "sparse" submatrices of a given BIG matrix.  The
    submatrices are defined by two lists of indices, one for the rows and
    one for the cols.

    I've two problems with that.

    The first one is how to ACCESS (read) the matrix entries from two
    lists of indices, in "pseudo code":

        mat = scipy.matrix( range(0,100) ) 
        mat( [1,3,4], [4,5,9] )

    in order to get all the elements mat[i,j] for i in [1,3,4] and j in
    [4,5,9].

    The second is how to CHANGE those elements.  I manage to solve the
    first problem (but the sintax seems to me "uglier than needed") but
    not the second.

    Here goes what I can do....

    Defining and filling the matrix is easy enough.

        import scipy
        mat = scipy.matrix( range(0,100) ) 
        mat.shape = (10,10)
        print mat

            [[ 0  1  2  3  4  5  6  7  8  9]
             [10 11 12 13 14 15 16 17 18 19]
             [20 21 22 23 24 25 26 27 28 29]
             [30 31 32 33 34 35 36 37 38 39]
             [40 41 42 43 44 45 46 47 48 49]
             [50 51 52 53 54 55 56 57 58 59]
             [60 61 62 63 64 65 66 67 68 69]
             [70 71 72 73 74 75 76 77 78 79]
             [80 81 82 83 84 85 86 87 88 89]
             [90 91 92 93 94 95 96 97 98 99]]

    Getting a submatrix with CONTIGUOUS indices is also quite easy:

        print mat[0:3,4:6]

            [[ 4  5]
             [14 15]
             [24 25]]

    Getting a submatrix with NOT contiguous indices is more difficult,
    partly because the indices come from a "vector" that is really scipy
    (1,N) matrix (and scipy "insist" in keeping two indices).

        rr = scipy.matrix('[2 4 5]')
        cc = scipy.matrix('[1,2]')
        print rr
        print cc
        print mat[[1,2,5],:][:,[4,8]]

        --- rr:
        [[2 4 5]]
        --- cc:
        [[1 2]]
        --- mat:
        [[14 18]
         [24 28]
         [54 58]]

    Now I want to add a given matrix to the submatrix I've "selected".

        off = scipy.matrix( [1000] * 6 )
        off.shape = (3,2)
        print '--- off\n', off

        print '--- mat + off'
        mat[[1,2,5],:][:,[4,8]] + off

        print --- 'try mat += off ... mat'
        mat[[1,2,5],:][:,[4,8]] += off
        print mat[[1,2,5],:][:,[4,8]] 

        --- off
        [[1000 1000]
         [1000 1000]
         [1000 1000]]
        --- mat + off
        [[1014 1018]
         [1024 1028]
         [1054 1058]]
        --- try mat += off ... mat
        [[14 18]
         [24 28]
         [54 58]]

    So it seems that the submatrix is a kind of "read-on view" that i can
    not update, but updating is just what I need!!!

        Any hint? Thanks in advance, Luca

    PS: the matrix (mat) that I'm "constructing" will be (in the future)
        BIG (1 million x 1 million, up to 100K x 100K is already feasible
        in MATLAB) SPARSE (some 20 or 30 not null elements per row) and
        COMPLEX, and the associated linear problem will need to be solved
        with some iterative method (GMRES or similar) .... but that's
        another story :)