Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Modifying MEX arguments in place?
Date: Thu, 14 Feb 2008 14:48:02 +0000 (UTC)
Organization: University of California San Diego
Lines: 201
Message-ID: <fp1kb2$nrv$1@fred.mathworks.com>
References: <fp0fqv$3gd$1@fred.mathworks.com> <muyir0r1v1p.fsf@G99-Boettcher.llan.ll.mit.edu>
Reply-To: <HIDDEN>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1203000482 24447 172.30.248.35 (14 Feb 2008 14:48:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 14 Feb 2008 14:48:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 163259
Xref: news.mathworks.com comp.soft-sys.matlab:451397



Peter Boettcher <boettcher@ll.mit.edu> wrote in message
<muyir0r1v1p.fsf@G99-Boettcher.llan.ll.mit.edu>...
> "Petr Krysl" <pkryslNOSP@Mucsd.edu> writes:
> 
> First generate a dataset.  Then do b=a;  (Copy the
dataset).  Then run
> your in-place code on "a".  Then look at b.  That's the
problem. 
> 

Actually, I'm trying carefully to avoid having two copies of
the data (because of the tremendous size). The variable I
want to modify is an attribute of a class (variable A): see the 
solve method.

% Constructor of a sparse system matrix. This class provides
methods for
% constructing a skyline-storage system matrix, In-place
factorization into
% L*D*L^T, and solution via forward and backward substitution.
%
% The factorization and forward/backward substitution code
is adapted from
% a Fortran 77 kernel by Bathe/1976 (see Bittnar, Rericha 1983).
%
% function retobj = sparse_colsol_sysmat ()
%
% See discussion of constructors in <a href="matlab:helpwin
'sofea/classes/Contents'">classes/Contents</a>.
%
function self = sparse_colsol_sysmat (varargin)
    class_name='sparse_colsol_sysmat';
    parent=sysmat;

    Factorized =false;
    Single = false;
    invalid_eqnum = 0;
    A = [];
    MAXA= [];
    Column_height = [];
    NEQ= 0;
    MEPSIL = 1.0e-307;
    Factorization_progress = [];

    function assemble_method(ems)
        ne =length(ems);
        %         Compute the height of the columns
        Column_height =zeros(1,NEQ);
        for m = 1:ne
            eqnums=get(ems(m), 'eqnums');
            minen=min(eqnums(eqnums~= 0));
            for I=1:length (eqnums);
                II=eqnums(I);
                if (II~=invalid_eqnum)
                    ME=II-minen;
                    if (ME >Column_height(II))
                        Column_height(II) = ME;
                    end
                end
            end
        end
        %         Compute the addresses of the Diagonal elements
        MAXA=zeros(1,NEQ+1);
        MAXA(1:2) =(1:2);
        for I = 2:NEQ
            MAXA(I+1)=MAXA(I)+Column_height(I)+1;
        end
        %         Allocate the Matrix
        if (Single)
            disp ([' Allocating ' num2str(MAXA(end)) '
single-precision matrix elements'])
            A =zeros(1,MAXA(end),'single');
        else
            disp ([' Allocating ' num2str(MAXA(end)) '
double-precision matrix elements'])
            A =zeros(1,MAXA(end),'double');
        end
        %         Assemble the element Matrices
        for m = 1:ne
            eqnums = get(ems(m), 'eqnums');
            mat    = get(ems(m), 'mat');
            ND =length(eqnums);
            for I = 1:ND %DO 200 I = 1,ND
                II    = eqnums(I);
                if (II ~=invalid_eqnum)
                    MI = MAXA(II);
                    for J = 1:ND %DO 220 J = 1,ND
                        JJ = eqnums(J);
                        if (JJ ~=invalid_eqnum)
                            if (II >= JJ)
                                KK      = MI+II-JJ;
                                A( KK ) = A( KK ) + mat( I,J );
                            end
                        end
                    end%220       CONTINUE
                end
            end%200 CONTINUE
        end
    end
    self.assemble_method=@assemble_method;

    function x=solve_method (V)
%         This is a MEX code
        if (~Factorized)
            MAXA =uint32(MAXA);
            colsoldecomp(A,  MAXA);
            Factorized = true;
        end
        if (Single)
            V =single(V);
        end
        x=REDBAK_precision_sensitive(V);
        if (Single)
            x =double(x);
        end
%         This is pure-Matlab code
%         if (~Factorized)
%             DECOMP_precision_sensitive;
%             Factorized = true;
%         end
%         if (Single)
%             V =single(V);
%         end
%         x=REDBAK_precision_sensitive(V);
%         if (Single)
%             x =double(x);
%         end
    end
    self.solve_method=@solve_method;

    function retobj=start_method (dim)
        NEQ =dim;
        Factorized = ~true;
        retobj =self;
    end
    self.start_method=@start_method;

    function val = get_method(varargin)
        if (nargin == 0)
            val = {{{'dim'}, {'dimension, array '}};
                {{'mat'}, {'sparse data matrix, array
[dim,dim]'}};
                {{'MAXA'}, {'diagonal element address array'}}};
            return;
        else
            prop_name=varargin{1};
            switch prop_name
                case 'mat'
                    val = sparse(NEQ,NEQ);
                    for m = 1:NEQ
                        k = MAXA(m);
                        r=m;
                        while (k <MAXA(m+1))
                            val(m,r)=A(k);
                            val(r,m)=A(k);
                            r=r+1; k=k+1;
                        end
                    end
                case 'dim'
                    val = NEQ;
                case 'MAXA'
                    val=MAXA;
                case 'Column_height'
                    val=Column_height;
                otherwise
                    error(['Unknown property name '''
prop_name '''!']);
            end
        end
        return;
    end
    self.get_method=@get_method;

    function retval=set_method(varargin)
        if (nargin == 1)
            retval = {{{'Factorization_progress'},
{'progress report function Handle'}};
                {{'Single'}, {'use storage in single
precision, true or false'}};};
            return;
        end
        prop_name=varargin{1};
        val=varargin{2};
        switch prop_name
            case 'Factorization_progress'
                Factorization_progress = val;
            case 'Single'
                Single = (val==true);
            otherwise
                error(['Unknown property ''' prop_name '''!'])
        end
        retval=self; % return self
    end
    self.set_method=@set_method;

    self = class(self, class_name);
end