Re: [Pysparse-users] PysparseSpDiagsMatrix usage from spdiags in matlab

[Dominique O. <dom...@gm...>]

On Thu, Nov 5, 2009 at 4:00 AM, zipeppe <zi...@gm...> wrote:

> Dear all,
>
> I don't know it this is the right place to ask, but unfortunately the
> official Python-forum doesn't help.
>

This is the right place for Pysparse questions.


> I am quite new to Python and I am unable to correctly translate the
> following MATLAB script in Python code, specially because of the *spdiags*proprietary function:
>
> function bd=BLC(b,A,B,s)
>    L=length(b);
>    Bs=-B/s;  As=-A/s;
>    bd=ones(L,1)*median(b);bd0=b;nm=norm(bd-bd0);nm0=RealMax;
>    M0=-ones(L,1)/As;
>    e=ones(L,1);
>    D0=spdiags([1*e -4*e 6*e -4*e 1*e],-2:2,L,L);
>    D0(1,1)=2; D0(L,L)=2;
>    D0(2,1)=-4; D0(1,2)=-4; D0(L,L-1)=-4; D0(L-1,L)=-4;
>    D0(2,2)=10;D0(L-1,L-1)=10;
>    I=0;
>    while nm>10 & I<30
>        %& nm<nm0;
>        I=I+1;
>        M=M0;D=D0;bd0=bd;nm0=nm;
>        for i=1:L
>            if bd(i)>b(i)
>                M(i)=M(i)+2*Bs*b(i)/As;
>                D(i,i)=D(i,i)+2*Bs/As;
>            end
>        end
>        bd=D\M;
>        nm=norm(bd0-bd);
>    end
>
> The algorithm receives *b* as (Nx1 data vector), *A* (1x1), *B* (1x1) and
> *s* (1x1) as local parameters, and creates *bd* (Nx1data vector) as
> output.
>

The Matlab command that you highlighted creates an L x L matrix with the
vectors e, -4*e, 6*e, -4*e and e on and around the main diagonal. The
position argument is -2:2 which means that the first vector (e) goes on the
sub-sub-diagonal, the second vector (-4e) goes on the sub-diagonal, 6e goes
on the main diagonal, -4e goes one above the main diagonal and e goes two
above the main diagonal like so.

Here is a code fragment that performs the same operation with Pysparse:

In [1]: from pysparse.pysparseMatrix import PysparseSpDiagsMatrix as spdiags
In [2]: import numpy as np
In [3]: L = 5  # The length of b in your Matlab code
In [4]: e = np.ones(L)
In [5]: D0 = spdiags(L, (e, -4*e, 6*e, -4*e, e), (-2,-1,0,1,2))
In [6]: print D0
 6.000000  -4.000000   1.000000      ---        ---
-4.000000   6.000000  -4.000000   1.000000      ---
 1.000000  -4.000000   6.000000  -4.000000   1.000000
    ---     1.000000  -4.000000   6.000000  -4.000000
    ---        ---     1.000000  -4.000000   6.000000

The first argument (L) specifies the matrix size, the second is the set of
vectors that should be laid out on the diagonals and the third gives the
indices of the diagonals that are concerned.

For a more Matlab-like notation, the last argument can also be written
np.r_[-2:2] (see the Numpy function r_ ... note the underscore).

I hope this helps.

-- 
Dominique