sparse.src

没有说明

/*
**> sparseHConcat
**
**  Purpose:  concatenates a sparse matrix on the right of
**            another sparse matrix.
**
**  Format:  z = sparseHConcat(y,x);
**
**  Input:   y    MxN sparse matrix.
**
**           x    MxL sparse matrix.
**
**  Output:  z    Mx(L+N) sparse matrix.
*/

proc sparseHConcat(y,x);
    local n,c_y,c_x,c,z;

    if not isSparse(y);
        if not trapchk(4);
            errorlog "sparseHConcat:  first argument not a sparse matrix";
            end;
        else;
            retp(error(0));
        endif;
    endif;
    if not isSparse(x);
        if not trapchk(4);
            errorlog "sparseHConcat:  second argument not a sparse matrix";
            end;
        else;
            retp(error(0));
        endif;
    endif;

    c_y = sparseCols(y);
    c_x = sparseCols(x);

    if c_y /= c_x;
        if not trapchk(4);
            errorlog "sparseHConcat:  sparse matrices not conformable";
            end;
        else;
            retp(error(0));
        endif;
    endif;

    n = sparseNZE(y) + sparseNZE(x);
    z = zeros(floor((sparseRows(y)+n+5)/2)+n,1);

    dllcall sparseHConcatdll(z,y,x);
    retp(z);
endp;




/*
**> sparseRows
**
**  Purpose:  get number of rows in sparse matrix.
**
**  Format:  r = sparseRows(x);
**
**  Input:   x    MxN sparse matrix.
**
**  Output:  r    scalar, number of rows.
*/

proc sparseRows(x);
    local r;
    if not isSparse(x);
        if not trapchk(4);
            errorlog "sparseRows:  not a sparse matrix";
            end;
        else;
            retp(error(0));
        endif;
    endif;
    r = 0;
    dllcall sparseRowsdll(x,r);
    retp(r);
endp;


/*
**> sparseCols
**
**  Purpose:  get number of cols in sparse matrix.
**
**  Format:  c = sparseCols(x);
**
**  Input:   x    MxN sparse matrix.
**
**  Output:  c    scalar, number of cols.
*/

proc sparseCols(x);
    local c;
    if not isSparse(x);
        if not trapchk(4);
            errorlog "sparseCols:  not a sparse matrix";
            end;
        else;
            retp(error(0));
        endif;
    endif;
    c = 0;
    dllcall sparseColsdll(x,c);
    retp(c);
endp;


/*
**> sparseNZE
**
**  Purpose:  get number of nonzero elements.
**
**  Format:  r = sparseNZE(x);
**
**  Input:   x    MxN sparse matrix.
**
**  Output:  r    scalar, number of nonzero elements.
*/

proc sparseNZE(x);
    local r;
    if not isSparse(x);
        if not trapchk(4);
            errorlog "sparseRows:  not a sparse matrix";
            end;
        else;
            retp(error(0));
        endif;
    endif;
    r = 0;
    dllcall sparseNZEdll(x,r);
    retp(r);
endp;


/*
**> denseSubmat
**
**  Purpose:  get dense submatrix of sparse matrix.
**
**  Format:  e = denseSubmat(x,r,c);
**
**  Input:   x    MxN sparse matrix.
**
**           r    Kx1 vector, row indices.
**
**           c    Lx1 vector, column indices.
**
**  Output:  e    KxL dense matrix.
**
**  Remarks:  if r or c are scalar zeros, all rows or columns will
**            be returned.
*/


proc denseSubmat(x,r,c);
    local e, rows_r, rows_c;

    if not isSparse(x);
        if not trapchk(4);
            errorlog "denseSubmat:  first argument not a sparse matrix";
            end;
        else;
            retp(error(0));
        endif;
    endif;

    r=vec(r);
    c=vec(c);
    if r == 0;
        r = seqa(1,1,sparseRows(x));
    endif;
    if c == 0;
        c = seqa(1,1,sparseCols(x));
    endif;

    rows_r = rows(r);
    rows_c = rows(c);
    e = ones(rows(r),rows(c));

    dllcall denseSubmatdll(x,e,r,rows_r,c,rows_c);

    retp(e);
endp;


/*
**> sparseSubmat
**
**  Purpose:  get sparse submatrix of sparse matrix.
**
**  Format:  e = sparseSubmat(x,r,c);
**
**  Input:   x    MxN sparse matrix.
**
**           r    Kx1 vector, row indices.
**
**           c    Lx1 vector, column indices.
**
**  Output:  e    KxL sparse matrix.
**
**  Remarks:  if r or c are scalar zeros, all rows or columns will
**            be returned.
*/

proc sparseSubmat(x,r,c);
    local e, rows_r, rows_c, len;

    if not isSparse(x);
        if not trapchk(4);
            errorlog "sparseSubmat:  first argument not a sparse matrix";
            end;
        else;
            retp(error(0));
        endif;
    endif;

    r=vec(r);
    c=vec(c);
    len = 0;

    if r == 0;
        r = seqa(1,1,sparseRows(x));
    endif;
    if minc(r) < 0 or maxc(r) > sparseRows(x);
        if not trapchk(4);
            errorlog "sparseSubmat:  rows of submatrix exceed rows of matrix";
            end;
        else;
            retp(error(0));
        endif;
    endif;
    if c /= 0;
        if minc(c) < 0 or maxc(c) > sparseCols(x);
            if not trapchk(4);
                errorlog "sparseSubmat:  columns of submatrix exceed columns"\
                         " of matrix";
                end;
            else;
                retp(error(0));
            endif;
        endif;
    endif;

    rows_r = rows(r);

    if c == 0;

        e = zeros(floor((rows_r+sparseRows(x)+5)/2)+sparseRows(x),1);
        dllcall sparseSubmatRowsdll(x,e,r,rows_r,len);

    else;

        rows_c = rows(c);
        e = zeros(floor((rows_r+rows_r*rows_c+5)/2)+rows_r*rows_c,1);
        dllcall sparseSubmatdll(x,e,r,rows_r,c,rows_c,len);

    endif;
    retp(e[1:len]);
endp;




/*
**> sparseTD
**
**  Purpose:  multiply sparse matrix times dense matrix.
**
**  Format:  z = sparseTD(x,y);
**
**  Input:   x    MxN sparse matrix.
**
**           y    NxL matrix, dense matrix.
**
**  Output:  z    MxL dense matrix, the result of x times y.
**
*/


proc sparseTD(x,y);
    local r_x, c_x, r_y, c_y, z;

    if not isSparse(x);
        if not trapchk(4);
            errorlog "sparseTD:  first argument not a sparse matrix";
            end;
        else;
            retp(error(0));
        endif;
    endif;

    r_y = rows(y);
    c_x = sparseCols(x);

    if c_x /= r_y;
        if not trapchk(4);
            errorlog "sparseTimesdense:  matrices not conformable";
            end;
        else;
            retp(error(0));
        endif;
    endif;

    r_x = sparseRows(x);
    c_y = cols(y);

    z = zeros(r_x,c_y);

    dllcall sparseTDdll(x,y,z,c_y);

    retp(z);
endp;




/*
**> sparseTrTD
**
**  Purpose:  multiply sparse matrix transposed times dense matrix.
**
**  Format:  z = sparseTrTD(x,y);
**
**  Input:   x    NxM GAUSS sparse matrix.
**
**           y    NxL matrix, dense matrix.
**
**  Output:  z    MxL dense matrix, the result of x transpose times y.
**
*/


proc sparseTrTD(x,y);
    local r_x, c_x, r_y, c_y, z;

    if not isSparse(x);
        if not trapchk(4);
            errorlog "sparseTD:  first argument not a sparse matrix";
            end;
        else;
            retp(error(0));
        endif;
    endif;

    r_y = rows(y);
    r_x = sparseRows(x);

    if r_x /= r_y;
        if not trapchk(4);
            errorlog "sparseTransposeTimesdense:  matrices not conformable";
            end;
        else;
            retp(error(0));
        endif;
    endif;

    c_y = cols(y);
    c_x = sparseCols(x);

    z = zeros(c_x,c_y);

    dllcall sparseTrTDdll(x,y,z,c_y);

    retp(z);
endp;


/*
**  sparseSet  -  resets global matrices to default values
**
*/



proc(0) = sparseSet;

  _sparse_Damping = 0;
  _sparse_Atol = 0;
  _sparse_Btol = 0;
  _sparse_CondLimit = 0;
  _sparse_NumIters = 0;
  _sparse_RetCode = 0;
  _sparse_Anorm = 0;
  _sparse_Acond = 0;
  _sparse_Rnorm = 0;
  _sparse_ARnorm = 0;
  _sparse_Xnorm = 0;
endp;





#endif