For the MEX experts out there: I would like to factor a huge
matrix (passed in as a one-dimensional array) in-place in a
mex file.  Unfortunately, the documentation says there may
be side effects (and there are: the Matlab executive is
apparently getting all mixed up).

Is there no way I could modify an argument in-place? If that
is not possible, the penalty would be huge.  We are talking
about a matrix of several gigabytes size...

Thanks for your help.

Petr

In article <fp0fqv$3gd$1@fred.mathworks.com>,
Petr Krysl <pkryslNOSP@Mucsd.edu> wrote:
>For the MEX experts out there: I would like to factor a huge
>matrix (passed in as a one-dimensional array) in-place in a
>mex file.  Unfortunately, the documentation says there may
>be side effects (and there are: the Matlab executive is
>apparently getting all mixed up).

>Is there no way I could modify an argument in-place?

I know little about MEX, but generally speaking, the memory
that you are handed for any Matlab variable might be shared
with other Matlab variables due to the "lazy copy" behaviour.
The columns you see might be linked to columns of other variables.
When you change the memory in-place, you would have changed those
other variables too.
-- 
  "Let me live in my house by the side of the road --
   It's here the race of men go by.
   They are good, they are bad, they are weak, they are strong
   Wise, foolish -- so am I;"                 -- Sam Walter Foss

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
message <fp0gke$gqq$1@canopus.cc.umanitoba.ca>...
> I know little about MEX, but generally speaking, the memory
> that you are handed for any Matlab variable might be shared
> with other Matlab variables due to the "lazy copy" behaviour.
> The columns you see might be linked to columns of other
variables.
> When you change the memory in-place, you would have
changed those
> other variables too.

I am not sure I follow: I can certainly access the
individual elements of the array that I want to modify by
read-access (and I see the right data). Why can't I write
into the same memory locations? I understand that I might be
modifying a copy of the input argument, but that does not
seem to be the case. Actually, when I look at the matrix
after I've modified the input argument and returned from the
mex file, the matrix contains correctly modified data. 
Unfortunately something else seems to have been affected,
and the Matlab interpreter stops working correctly.

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
message <fp0gke$gqq$1@canopus.cc.umanitoba.ca>...
> 
> I know little about MEX, but generally speaking, the memory
> that you are handed for any Matlab variable might be shared
> with other Matlab variables due to the "lazy copy" behaviour.
> When you change the memory in-place, you would have
changed those
> other variables too.

That is the danger, of course.

> The columns you see might be linked to columns of other
variables.

I suppose that is possible, but I seriously doubt it. When
you create a variable, MATLAB grabs memory from the heap to
do it. It has to grab the data memory as one contiguous
chunk in order for all of the other things you do with
variables to work properly, like reshape & calls to the BLAS
etc. Real data & imag data are two separate memory chunks,
but each one is contiguous. When freeing this data memory
back to the heap, it would have to be done all at once, not
piecemeal. That would be impossible if other variables were
allowed to share only portions of the data block like
columns. You wouldn't be able to free it until all of the
variables that had pieces were cleared. This would be a
tremendous waste of resources and a nightmare to manage. For
example, do the following:

format debug
a = rand(3,3)
b = a(:,1)
c = a(:,:)
d = a

You can see that the real data pointer pr for b and c are
different from a, indicating that a copy of the data was
made, whereas d and a share the same data pointer. In this
case a copy of the data was made to create c even though it
was the same size and could have used the same data pointer
as a. Now do this:

e = b(:)

Suddenly MATLAB realizes that all of the data is being
copied so an actual copy isn't made ... e simply uses the
same data pointer as b.  I don't know what the rules for
data sharing are, but you can see that it can be tricky to
try and guess what MATLAB will do in every case.

But getting back to OP's original question. Can you modify
input variable data in place and get away with it? The
answer is yes, but only if you are very careful. And even
then you will have to accept some risk. If you create a huge
2GB variable in MATLAB and don't do anything that would
cause data sharing (like the d = a statement above), then
you can get away with modifying the data in place inside a
mex routineand not screw up other variables. The problem, of
course, is that there is no official list (that I am aware
of) of things you  should avoid in order to prevent data
sharing. You should avoid assignments like the d = a above,
of course, but what else?

My advice to OP is to go ahead and do it (at risk) if the
memory/speed penalty is too great otherwise, but never do
MATLAB assignments and function calls that might result in
data sharing with this variable before your mex call.

None of this would be recommended by the Mathworks, of course.

James Tursa

"Petr Krysl" <pkryslNOSP@Mucsd.edu> schrieb im Newsbeitrag 
news:fp0fqv$3gd$1@fred.mathworks.com...
> For the MEX experts out there: I would like to factor a huge
> matrix (passed in as a one-dimensional array) in-place in a
> mex file.  Unfortunately, the documentation says there may
> be side effects (and there are: the Matlab executive is
> apparently getting all mixed up).
>
> Is there no way I could modify an argument in-place? If that
> is not possible, the penalty would be huge.  We are talking
> about a matrix of several gigabytes size...
>
> Thanks for your help.
>
> Petr

Hi,

besides the warnings given by others, I have another question:
do you really need to do this in a mex file? And why?
MATLAB provides some sort of in-place operation for
variables, e.g.,

function X = somefunction(X)
% modify X

and you call this function by

foo = somefunction(foo);

then MATLAB will actually not copy your matrix foo...

Titus 

"Titus" <titus.edelhofer@mathworks.de> wrote in message
<fp0utr$4gr$1@fred.mathworks.com>...
> 
> "Petr Krysl" <pkryslNOSP@Mucsd.edu> schrieb im Newsbeitrag 
> news:fp0fqv$3gd$1@fred.mathworks.com...
> > For the MEX experts out there: I would like to factor a huge
> > matrix (passed in as a one-dimensional array) in-place in a
> > mex file.  Unfortunately, the documentation says there may
> > be side effects (and there are: the Matlab executive is
> > apparently getting all mixed up).
> >
> > Is there no way I could modify an argument in-place? If that
> > is not possible, the penalty would be huge.  We are talking
> > about a matrix of several gigabytes size...
> >
> > Thanks for your help.
> >
> > Petr
> 
> Hi,
> 
> besides the warnings given by others, I have another question:
> do you really need to do this in a mex file? And why?
> MATLAB provides some sort of in-place operation for
> variables, e.g.,
> 
> function X = somefunction(X)
> % modify X
> 
> and you call this function by
> 
> foo = somefunction(foo);
> 
> then MATLAB will actually not copy your matrix foo...
> 
> Titus 
> 
> 

James' advice is perfect.

Titus:  there are lots of things you can do in mexFunctions
that are too slow in M, or take too much memory.  So
sometimes you just can avoid going to C.

The new feature where MATLAB can avoid a memory copy in
functions like x=f(x) is really handy, but unfortunately
it's not made available to us mex authors.  At least not
yet.  It would be great if there could be query like:

if (mxIsThisVariableAShallowCopy (prhs [0]))
{
    oops, make a copy of prhs[0], and modify the copy
}
else
{
    modify prhs[0] in place
}
plhs [0] = prhs [0] ;

There would be lots of cool things a mexFunction could do if
it could modify its inputs (like a sparse cholupdate with an
O(1) best-case run time, for xample), but it's very dangerous.

The safest thing to do for a mexFunction that modifies its
inputs is to wrap it in an M-file that guarantees the
parameters to the hazardous mexFunction aren't a shallow
copy of another variable.

"Tim Davis" <davis@cise.ufl.edu> wrote in message 
<fp1g61$elr$1@fred.mathworks.com>...
 
> The new feature where MATLAB can avoid a memory copy in
> functions like x=f(x) is really handy, but unfortunately
> it's not made available to us mex authors.

Not to mention that, as you recalled, it is a relatively 
new functionality that is not available in older ML 
versions and in compiled code.
MEX in place modifications are very usefull for memory 
savings. I have done it a lot without regrets.

J. Luis 

"Tim Davis" <davis@cise.ufl.edu> writes:

> The new feature where MATLAB can avoid a memory copy in
> functions like x=f(x) is really handy, but unfortunately
> it's not made available to us mex authors.  At least not
> yet.  It would be great if there could be query like:
>
> if (mxIsThisVariableAShallowCopy (prhs [0]))
> {
>     oops, make a copy of prhs[0], and modify the copy

This stuff is getting to be a bigger and bigger problem.  More and more
features are being added to MATLAB without access methods for users, or
without proper documentation.  In many cases, it is almost as if MATLAB
is aiming to be an analysis package rather than a programming
environment.  Whenever a MATLAB employee can produce a chunk of code
that users cannot, the value of the software as a programming
environment has gone way down.

There's a bunch of new graphics stuff (like lineseries objects) that are
poorly (if at all) documented.  There are opaque objects that users
cannot create on their own.  There are features like you mention above
which are just "magic".

It's either intentional, or evidence of poor design, or ...

-Peter

"Petr Krysl" <pkryslNOSP@Mucsd.edu> writes:

> roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
> message <fp0gke$gqq$1@canopus.cc.umanitoba.ca>...
>> I know little about MEX, but generally speaking, the memory
>> that you are handed for any Matlab variable might be shared
>> with other Matlab variables due to the "lazy copy" behaviour.
>> The columns you see might be linked to columns of other
> variables.
>> When you change the memory in-place, you would have
> changed those
>> other variables too.
>
> I am not sure I follow: I can certainly access the
> individual elements of the array that I want to modify by
> read-access (and I see the right data). Why can't I write
> into the same memory locations? I understand that I might be
> modifying a copy of the input argument, but that does not
> seem to be the case. Actually, when I look at the matrix
> after I've modified the input argument and returned from the
> mex file, the matrix contains correctly modified data. 
> Unfortunately something else seems to have been affected,
> and the Matlab interpreter stops working correctly.

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. 

This can happen in strange places if you're not careful.

The answer is to use the undocumented mxUnshareArray() before modifying
the input.  We think the prototype is:

int mxUnshareArray(mxArray *, int);

We don't know what the second argument does.  Zero seems to do the trick.


-Peter

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

"Tim Davis" <davis@cise.ufl.edu> schrieb im Newsbeitrag 
news:fp1g61$elr$1@fred.mathworks.com...
> "Titus" <titus.edelhofer@mathworks.de> wrote in message
> <fp0utr$4gr$1@fred.mathworks.com>...
>>
>> "Petr Krysl" <pkryslNOSP@Mucsd.edu> schrieb im Newsbeitrag
>> news:fp0fqv$3gd$1@fred.mathworks.com...
>> > For the MEX experts out there: I would like to factor a huge
>> > matrix (passed in as a one-dimensional array) in-place in a
>> > mex file.  Unfortunately, the documentation says there may
>> > be side effects (and there are: the Matlab executive is
>> > apparently getting all mixed up).
>> >
>> > Is there no way I could modify an argument in-place? If that
>> > is not possible, the penalty would be huge.  We are talking
>> > about a matrix of several gigabytes size...
>> >
>> > Thanks for your help.
>> >
>> > Petr
>>
>> Hi,
>>
>> besides the warnings given by others, I have another question:
>> do you really need to do this in a mex file? And why?
>> MATLAB provides some sort of in-place operation for
>> variables, e.g.,
>>
>> function X = somefunction(X)
>> % modify X
>>
>> and you call this function by
>>
>> foo = somefunction(foo);
>>
>> then MATLAB will actually not copy your matrix foo...
>>
>> Titus
>>
>>
>
> James' advice is perfect.
>
> Titus:  there are lots of things you can do in mexFunctions
> that are too slow in M, or take too much memory.  So
> sometimes you just can avoid going to C.
>
> The new feature where MATLAB can avoid a memory copy in
> functions like x=f(x) is really handy, but unfortunately
> it's not made available to us mex authors.  At least not
> yet.  It would be great if there could be query like:
>
> if (mxIsThisVariableAShallowCopy (prhs [0]))
> {
>    oops, make a copy of prhs[0], and modify the copy
> }
> else
> {
>    modify prhs[0] in place
> }
> plhs [0] = prhs [0] ;
>
> There would be lots of cool things a mexFunction could do if
> it could modify its inputs (like a sparse cholupdate with an
> O(1) best-case run time, for xample), but it's very dangerous.
>
> The safest thing to do for a mexFunction that modifies its
> inputs is to wrap it in an M-file that guarantees the
> parameters to the hazardous mexFunction aren't a shallow
> copy of another variable.

Hi Tim,
no doubt, MEX files are often quite handy, it's not that I want
to put this in question. The only thing I try to keep pointing at,
is, that I've seen quite a few MEX files of people who could
have avoided the effort when using (o.k., in newer versions)
the capabilities e.g. JIT has added to the MATLAB programming
language.

Titus 

"Titus" <titus.edelhofer@mathworks.de> wrote in message
<fp1kf6$pgq$1@fred.mathworks.com>...
> 
> "Tim Davis" <davis@cise.ufl.edu> schrieb im Newsbeitrag 
> news:fp1g61$elr$1@fred.mathworks.com...
> > "Titus" <titus.edelhofer@mathworks.de> wrote in message
> no doubt, MEX files are often quite handy, it's not that I
want
> to put this in question. The only thing I try to keep
pointing at,
> is, that I've seen quite a few MEX files of people who could
> have avoided the effort when using (o.k., in newer versions)
> the capabilities e.g. JIT has added to the MATLAB programming
> language.
> 
> Titus 

Titus,

Actually, my solver already modifies the matrix in place
when I use pure Matlab code.  The factorization is about an
order of magnitude slower than the Matlab Choleski sparse
solver (which is a compiled mex code). Therefore, I see a
potential to speed up while at the same time preserving the
in-place aspect...

Petr

"Petr Krysl" <pkryslNOSP@Mucsd.edu> writes:

> 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.

I'm trying to explain to you the problem with in-place modification.
The deal is that b=a does NOT copy the data, but instead shares it,
which leads to strange behavior. Then I give the solution to it.

If you can avoid it altogether, because you are sure it is never
assigned (this includes passing into functions, for example), then
ignore the whole thing!  Modify your data at will!

-Peter

Peter Boettcher <boettcher@ll.mit.edu> wrote in message 
<muyir0r1v1p.fsf@G99-Boettcher.llan.ll.mit.edu>...
> 
> The answer is to use the undocumented mxUnshareArray() 
before modifying
> the input.  We think the prototype is:
> 
> int mxUnshareArray(mxArray *, int);
> 
> We don't know what the second argument does.  Zero seems 
to do the trick.
> 
> 
> -Peter

If mxUnshareArray does what I think it does, which is make 
a copy of the data (if needed) and reset the data pointers 
to the new copy, then I disagree that this is *the* answer 
to OP's problem. OP is trying to avoid multiple copies of 
his data, not make new ones just so his mex routine can 
modify the data in place. For example, suppose A was a 2GB 
array and then you did B = A, and then you called the mex 
routine mymexroutine(B). What is the difference, memory 
usage wise, between calling mxUnshareArray(prhs[0]) and 
then modifying prhs[0] in place vs just calling plhs[0]
=mxDuplicateArray(prhs[0]) and then modifying plhs[0]? 
Nothing. You still get two copies of the 2GB data either 
way. I still think the best thing to do is be very careful 
on the MATLAB side to never do assignments or function 
calls that might result in data sharing *before* you call 
the mex routine. Then you can call mxUnshareArray inside 
the mex routine just to be safe. In those cases where you 
screwed up on the MATLAB side and the data was in fact 
shared, you would wind up paying the memory/speed penalty 
or maybe even exhaust the heap. But at least you would 
return with an error instead of screwing up other 
variables had you not called mxUnshareArray.

Do you know if mxUnshareArray behaves like other mx 
functions that allocate memory? i.e., if it fails in a mex 
routine does it free all memory allocated inside the mex 
routine and return control to MATLAB? And in an engine 
application if it fails does it return 1?

James Tursa

Just did a quick check. mxUnshareArray does indeed return 
control to MATLAB (with an out-of-memory error message) if 
there is an out of memory error inside the call. I don't 
know about engine applications yet.

James Tursa

"Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message
<fp1lbe$8v7$1@fred.mathworks.com>...
> "Titus" <titus.edelhofer@mathworks.de> wrote in message
> <fp1kf6$pgq$1@fred.mathworks.com>...
> > 
> > "Tim Davis" <davis@cise.ufl.edu> schrieb im Newsbeitrag 
> > news:fp1g61$elr$1@fred.mathworks.com...
> > > "Titus" <titus.edelhofer@mathworks.de> wrote in message
> > no doubt, MEX files are often quite handy, it's not that I
> want
> > to put this in question. The only thing I try to keep
> pointing at,
> > is, that I've seen quite a few MEX files of people who could
> > have avoided the effort when using (o.k., in newer versions)
> > the capabilities e.g. JIT has added to the MATLAB
programming
> > language.
> > 
> > Titus 
> 
> Titus,
> 
> Actually, my solver already modifies the matrix in place
> when I use pure Matlab code.  The factorization is about an
> order of magnitude slower than the Matlab Choleski sparse
> solver (which is a compiled mex code). Therefore, I see a
> potential to speed up while at the same time preserving the
> in-place aspect...
> 
> Petr
> 

<Advance notice of particular bias>
I'm the author of sparse Cholesky in MATLAB
</end notice>

I wouldn't recommend using a skyline code.  Even written in
C, theres no way any skyline code can beat a true sparse
Cholesky code.   For more details, see:

http://doi.acm.org/10.1145/1236463.1236465

The question about how to modify a mxArray in place is still
a good one, but you shouldn't try to do it for this case.

"Tim Davis" <davis@cise.ufl.edu> wrote in message
<fp2r4p$s2u$1@fred.mathworks.com>...
> 
> <Advance notice of particular bias>
> I'm the author of sparse Cholesky in MATLAB
> </end notice>
> 
> I wouldn't recommend using a skyline code.  Even written in
> C, theres no way any skyline code can beat a true sparse
> Cholesky code.   For more details, see:
> 
> http://doi.acm.org/10.1145/1236463.1236465
> 
> The question about how to modify a mxArray in place is still
> a good one, but you shouldn't try to do it for this case.
> 

Tim,

I hear you, and I know that the Cholesky solver in Matlab is
lightning fast.  However, I'm up against its memory usage
here: the largest matrix I want to deal with has on the
order of  20 GB of nonzeros, and that is just the upper
triangle. 

The skyline is slower, but its memory usage is way down
because I can do the factorization in place.

Petr

> 
> Tim,
> 
> I hear you, and I know that the Cholesky solver in Matlab is
> lightning fast.  However, I'm up against its memory usage
> here: the largest matrix I want to deal with has on the
> order of  20 GB of nonzeros, and that is just the upper
> triangle. 
> 
> The skyline is slower, but its memory usage is way down
> because I can do the factorization in place.
> 
> Petr

What about an iterative solver then, like something from the
preconditioned conjugate gradient family. You only need to
store the matrix plus some additional vectors, so it should
be much lighter in terms of memory use than trying to form L*L'.

You could probably get away with only storing the upper
triangle as well, as you only need to implement a
matrix-vector multiply.

Skyline solvers are generally really, really slow, and they
should also probably use more memory to store the matrix in
skyline format as opposed to a true sparse format.


Darren

"Titus" <titus.edelhofer@mathworks.de> wrote in message ...
> Hi Tim,
> no doubt, MEX files are often quite handy, it's not that I
want
> to put this in question. The only thing I try to keep
pointing at,
> is, that I've seen quite a few MEX files of people who could
> have avoided the effort when using (o.k., in newer versions)
> the capabilities e.g. JIT has added to the MATLAB programming
> language.
> 
> Titus 
> 
> 

Yes, I agree with you.  Sometimes using mex is not needed,
when the latest JIT can do just as well.  mex should be the
last resort.

"Darren Engwirda" <d_engwirda.nospam@hotmail.com> wrote in
message <fp3p66$4vu$1@fred.mathworks.com>...
> 
> > 
> > Tim,
> > 
> > I hear you, and I know that the Cholesky solver in Matlab is
> > lightning fast.  However, I'm up against its memory usage
> > here: the largest matrix I want to deal with has on the
> > order of  20 GB of nonzeros, and that is just the upper
> > triangle. 
> > 
> > The skyline is slower, but its memory usage is way down
> > because I can do the factorization in place.
> > 
> > Petr
> 
> What about an iterative solver then, like something from the
> preconditioned conjugate gradient family. You only need to
> store the matrix plus some additional vectors, so it should
> be much lighter in terms of memory use than trying to form
L*L'.
> 
> You could probably get away with only storing the upper
> triangle as well, as you only need to implement a
> matrix-vector multiply.
> 
> Skyline solvers are generally really, really slow, and they
> should also probably use more memory to store the matrix in
> skyline format as opposed to a true sparse format.
> 
> 
> Darren
> 

If chol is running out of memory, then iterative methods are
definitely worth trying.

Regarding the use of chol:  internally, chol uses a
different data structure for L than the standard MATLAB
sparse matrix.  So x=A\b when A is symmetric positive
definite takes less memory than what is required to hold a
MATLAB L.  However, L=chol(A,'lower') is worse; it must
create a MATLAB sparse matrix L from its (more compact)
supernodal L.   This is not quite as bad as you may imagine,
because some of the memory converting the supernodal L into
the MATLAB L is done in-place.

What is still worse, R=chol(A) must make another copy,
because the matrix computed by CHOLMOD is lower triangular,
not upper.

So R=chol(A) first creates a supernodal L, and then modifies
it partly in place into the MATLAB sparse L.  Then R=L' is
performed.  that last step is a pure copy, and thus the
total memory usage of R=chol(A) is about twice that of
L=chol(A,'lower').

So ... the moral is ... if you're running out of memory in
R=chol(A), use L=chol(A,'lower') instead.

Of course, in both cases, you  should actually use
[R,q,q] = chol(A) or [L,p,q] = chol(A,'lower'), since
otherwise you don't get the fill-reducing ordering, q.

You can read the code for "chol" itself in
SuiteSparse/CHOLMOD/MATLAB/chol2.c (for the R=chol(A) usage)
or SuiteSparse/CHOLMOD/MATLAB/lchol.c for the
L=chol(A,'lower') usage.  The "chol" code itself is not
there; but it's pretty much an amalgam of chol2.c and
lchol.c.  Note the extra transpose in chol2.c.  SuiteSparse
is on the File Exchange.

"James Tursa" <aclassyguywithaknotac@hotmail.com> writes:

> Peter Boettcher <boettcher@ll.mit.edu> wrote in message 
> <muyir0r1v1p.fsf@G99-Boettcher.llan.ll.mit.edu>...
>> 
>> The answer is to use the undocumented mxUnshareArray() before
>> modifying the input.  We think the prototype is:
>> 
>> int mxUnshareArray(mxArray *, int);
>> 
>> We don't know what the second argument does.  Zero seems to do the
>> trick.
>
> If mxUnshareArray does what I think it does, which is make 
> a copy of the data (if needed) and reset the data pointers 
> to the new copy, then I disagree that this is *the* answer 
> to OP's problem. OP is trying to avoid multiple copies of 
> his data, not make new ones just so his mex routine can 
> modify the data in place. For example, suppose A was a 2GB 
> array and then you did B = A, and then you called the mex 
> routine mymexroutine(B). What is the difference, memory 
> usage wise, between calling mxUnshareArray(prhs[0]) and 
> then modifying prhs[0] in place vs just calling plhs[0]
> =mxDuplicateArray(prhs[0]) and then modifying plhs[0]? 
> Nothing. You still get two copies of the 2GB data either 
> way.

Well... yes, in that case.  The difference is that mxUnshareArray()
copies the data only when necessary, while mxDuplicateArray will do it
every time, all the time.

> I still think the best thing to do is be very careful 
> on the MATLAB side to never do assignments or function 
> calls that might result in data sharing *before* you call 
> the mex routine. Then you can call mxUnshareArray inside 
> the mex routine just to be safe. In those cases where you 
> screwed up on the MATLAB side and the data was in fact 
> shared, you would wind up paying the memory/speed penalty 
> or maybe even exhaust the heap. But at least you would 
> return with an error instead of screwing up other 
> variables had you not called mxUnshareArray.

Right.

> Do you know if mxUnshareArray behaves like other mx 
> functions that allocate memory? i.e., if it fails in a mex 
> routine does it free all memory allocated inside the mex 
> routine and return control to MATLAB? And in an engine 
> application if it fails does it return 1?

I assume this to be the case...

-Peter

"Tim Davis" <davis@cise.ufl.edu> wrote in message
<fp3uvq$75i$1@fred.mathworks.com>...
> "Darren Engwirda" <d_engwirda.nospam@hotmail.com> wrote in
> message <fp3p66$4vu$1@fred.mathworks.com>...
> > 
> > > 

Thank you very much, Tim.   The information below is worth
its weight in gold (figuratively ;))...

Petr

> Regarding the use of chol:  internally, chol uses a
> different data structure for L than the standard MATLAB
> sparse matrix.  So x=A\b when A is symmetric positive
> definite takes less memory than what is required to hold a
> MATLAB L.  However, L=chol(A,'lower') is worse; it must
> create a MATLAB sparse matrix L from its (more compact)
> supernodal L.   This is not quite as bad as you may imagine,
> because some of the memory converting the supernodal L into
> the MATLAB L is done in-place.
> 
> What is still worse, R=chol(A) must make another copy,
> because the matrix computed by CHOLMOD is lower triangular,
> not upper.
> 
> So R=chol(A) first creates a supernodal L, and then modifies
> it partly in place into the MATLAB sparse L.  Then R=L' is
> performed.  that last step is a pure copy, and thus the
> total memory usage of R=chol(A) is about twice that of
> L=chol(A,'lower').
> 
> So ... the moral is ... if you're running out of memory in
> R=chol(A), use L=chol(A,'lower') instead.
> 
> Of course, in both cases, you  should actually use
> [R,q,q] = chol(A) or [L,p,q] = chol(A,'lower'), since
> otherwise you don't get the fill-reducing ordering, q.
> 
> You can read the code for "chol" itself in
> SuiteSparse/CHOLMOD/MATLAB/chol2.c (for the R=chol(A) usage)
> or SuiteSparse/CHOLMOD/MATLAB/lchol.c for the
> L=chol(A,'lower') usage.  The "chol" code itself is not
> there; but it's pretty much an amalgam of chol2.c and
> lchol.c.  Note the extra transpose in chol2.c.  SuiteSparse
> is on the File Exchange.

"Darren Engwirda" <d_engwirda.nospam@hotmail.com> wrote in
message <fp3p66$4vu$1@fred.mathworks.com>...
> What about an iterative solver then, like something from the
> preconditioned conjugate gradient family. You only need to
> store the matrix plus some additional vectors, so it should
> be much lighter in terms of memory use than trying to form
L*L'.
> 
> You could probably get away with only storing the upper
> triangle as well, as you only need to implement a
> matrix-vector multiply.


I agree with you there Darren, for a one-shot solution that
will be a good alternative.  However, I need to perform on
the order of a few dozen or a hundred solves (the linear
equation solver is called from within an eigenvalue solver
repeatedly). Then a single factorization may be amortized
quite quickly.


> 
> Skyline solvers are generally really, really slow, and they
> should also probably use more memory to store the matrix in
> skyline format as opposed to a true sparse format.

For finite element applications, with fairly limited number
of non-zeros and reasonable fill in, the skyline solver is
not too bad.  Furthermore, I'm hoping to amortize the cost
of the factorization as indicated above.

Petr

"Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message 

> Thank you very much, Tim.   The information below is worth
> its weight in gold (figuratively ;))...
> 
> Petr
> 

Considering that bits don't weigh very much at all, maybe
it's worth its weight in gold after all :-)

"Joaquim Luis" <jluis@--ualg--.pt> wrote in message
<fp1h86$f8o$1@fred.mathworks.com>...
> "Tim Davis" <davis@cise.ufl.edu> wrote in message 
> <fp1g61$elr$1@fred.mathworks.com>...
>  
> > The new feature where MATLAB can avoid a memory copy in
> > functions like x=f(x) is really handy, but unfortunately
> > it's not made available to us mex authors.
> 
> Not to mention that, as you recalled, it is a relatively 
> new functionality that is not available in older ML 
> versions and in compiled code.

From which version of MATLAB this kind of memory saving is
implemented?

Do you know why it is not implemented in compiled code?

Thanks,

Bruno

"Bruno Luong" <b.luong@fogale.fr> wrote in message 
<fp69vp$4vh$1@fred.mathworks.com>...
> "Joaquim Luis" <jluis@--ualg--.pt> wrote in message
> <fp1h86$f8o$1@fred.mathworks.com>...
> 
> > Not to mention that, as you recalled, it is a 
relatively 
> > new functionality that is not available in older ML 
> > versions and in compiled code.
> 
> From which version of MATLAB this kind of memory saving 
is
> implemented?

I don't know. Maybe at 2006b. I learnt about it in this 
Loren's blog

http://blogs.mathworks.com/loren/2007/03/22/in-place-
operations-on-data/
 
> Do you know why it is not implemented in compiled code?

I was refering to *real* compiler which ended its existance 
at R13.

L. Luis

"Tim Davis" <davis@cise.ufl.edu> wrote in message
<fp4mlb$fb4$1@fred.mathworks.com>...
> "Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message 
> 
> > Thank you very much, Tim.   The information below is worth
> > its weight in gold (figuratively ;))...

Tim,

I would like to avoid the transpose of L =chol(A,'lower')
when solving by forward/backward substitution (this is where
I run out of memory).  Is there support in the cholmod
package for this operation?

Thanks very much,

Petr

"Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message 

> I would like to avoid the transpose of L =chol(A,'lower')
> when solving by forward/backward substitution (this is where
> I run out of memory).  Is there support in the cholmod
> package for this operation?

Actually, I have just had some sort of inspiration ;). I
would do

c=A\V; x = (c'/A)';

instead of

x = A'\(A\V);

That seems to do the trick.

Petr

"Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message
<fp8402$n2u$1@fred.mathworks.com>...
> "Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message 
> 
> > I would like to avoid the transpose of L =chol(A,'lower')
> > when solving by forward/backward substitution (this is where
> > I run out of memory).  Is there support in the cholmod
> > package for this operation?
> 
> Actually, I have just had some sort of inspiration ;). I
> would do
> 
> c=A\V; x = (c'/A)';
> 
> instead of
> 
> x = A'\(A\V);
> 
> That seems to do the trick.

It seemed a good idea, but the implementation got in the
way.  Apparently, the / is implemented with A'\c, which
means that the matrix is first transposed.  Exactly what I
wanted to avoid.  So I'm still looking for a solution. Ideas
and suggestions are welcome.

Petr

PS: Here is the test that shows that the memory is increased
during the backward substitution phase (at least with 2007b).


clear all
n = 6500;
disp(' Starting in Five seconds')
pause(5);
% A = spalloc(n,n,nzmax);
% for i=1:n
%     v = randn(i-prof(i)+1,1);
%     A(prof(i):i,i) = v;
%     A(i,prof(i):i) = v';
%     A(i,i) = abs(A(i,i));
%     i
% end
disp(' Starting assembly')
A= sprandsym(n,0.1);
A = A + diag(sum(abs(A)));
V =rand(n,1);
disp(' Finished assembling the matrix')
% p=symamd(A);
% A = A(p,p);
% spy(A)

pause(5);
disp('  Starting factorization')
A = chol(A,'lower');
disp(' Finished  factorization of the matrix')
pause(5);
disp('  Starting Forward substitution')
c=A\V;
disp(' Finished forward Substitution')
pause(5);
disp('  Starting backward substitution')
x=(c'/A)';
disp(' Finished backward Substitution')

"Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message
<fp812a$qe2$1@fred.mathworks.com>...
> "Tim Davis" <davis@cise.ufl.edu> wrote in message
> <fp4mlb$fb4$1@fred.mathworks.com>...
> > "Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message 
> > 
> > > Thank you very much, Tim.   The information below is worth
> > > its weight in gold (figuratively ;))...
> 
> Tim,
> 
> I would like to avoid the transpose of L =chol(A,'lower')
> when solving by forward/backward substitution (this is where
> I run out of memory).  Is there support in the cholmod
> package for this operation?
> 
> Thanks very much,
> 
> Petr

You want to do the forward/backsolve after getting L, right?

Try cs_lsolve and cs_ltsolve, in CSparse/CXSparse.  I can't
recall off the top of my head if I have that same code in
CHOLMOD ... I'd have to look it up.  But I know it's in C*Sparse

"Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message
<fpd2c7$eqt$1@fred.mathworks.com>...
> "Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message
> <fp8402$n2u$1@fred.mathworks.com>...
> > "Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message 
> > 
> > > I would like to avoid the transpose of L =chol(A,'lower')
> > > when solving by forward/backward substitution (this is
where
> > > I run out of memory).  Is there support in the cholmod
> > > package for this operation?
> > 
> > Actually, I have just had some sort of inspiration ;). I
> > would do
> > 
> > c=A\V; x = (c'/A)';
> > 
> > instead of
> > 
> > x = A'\(A\V);
> > 
> > That seems to do the trick.
> 
> It seemed a good idea, but the implementation got in the
> way.  Apparently, the / is implemented with A'\c, which
> means that the matrix is first transposed.  Exactly what I
> wanted to avoid.  So I'm still looking for a solution. Ideas
> and suggestions are welcome.
> 
> Petr
> 
> PS: Here is the test that shows that the memory is increased
> during the backward substitution phase (at least with 2007b).
> 
> 
> clear all
> n = 6500;
> disp(' Starting in Five seconds')
> pause(5);
> % A = spalloc(n,n,nzmax);
> % for i=1:n
> %     v = randn(i-prof(i)+1,1);
> %     A(prof(i):i,i) = v;
> %     A(i,prof(i):i) = v';
> %     A(i,i) = abs(A(i,i));
> %     i
> % end
> disp(' Starting assembly')
> A= sprandsym(n,0.1);
> A = A + diag(sum(abs(A)));
> V =rand(n,1);
> disp(' Finished assembling the matrix')
> % p=symamd(A);
> % A = A(p,p);
> % spy(A)
> 
> pause(5);
> disp('  Starting factorization')
> A = chol(A,'lower');
> disp(' Finished  factorization of the matrix')
> pause(5);
> disp('  Starting Forward substitution')
> c=A\V;
> disp(' Finished forward Substitution')
> pause(5);
> disp('  Starting backward substitution')
> x=(c'/A)';
> disp(' Finished backward Substitution')


I see ... I was confused by your use of "A" for "L".  A is
the original matrix, so c=A\v looked like a full chol to me,
not a forward solve.

Don't use symamd ... let chol do the work (using amd
internally):

[L,ignore,p] = chol(A, 'lower') ;
x = cs_ltsolve (L, cs_lsolve (L, b(p))) ;
x (p) = x ;

I might be off on permutations ... I always have to rethink
them from scratch - I'm a bit dyslexic and get rows and
columns mixed up sometimes :-)

Also - you've picked an exceedingly slow and memory
intensive method for creating A.  See Loren's March 1st,
2007, blog for more details on how to do it better.

"Tim Davis" <davis@cise.ufl.edu> wrote in message 
> [L,ignore,p] = chol(A, 'lower') ;
> x = cs_ltsolve (L, cs_lsolve (L, b(p))) ;
> x (p) = x ;
> 

Thank you Tim, this works like a charm.  I do have one more
wish though: is this possible to do this with cholmod?
Otherwise I'd have to distribute csparse with my toolbox and
also I'd have to make sure that csparse is installed... Not
that it is a big deal, the install script works very well
indeed.  (By the way, I just ordered the companion SIAM
book, looks like a good read.)

Thanks again for your help,

Petr

cs_lsolve and cs_ltsolve are included in Csparse,
but not in CHOLMOD or LDL, as long as I can find. 
cs_lsolve and cs_ltsolve in Csparse are pretty slow,
compared to the triangular solvers in 
CHOLMOD/LDL/mldivide. (For the triangular solve,
mldivide was 3 times faster than cs_lsolve and cs_ltsolve 
in Csparse, for my problem.)

"J Mathgeek" <jmathgeek@gmail.com> wrote in message
<fpeliu$2h$1@fred.mathworks.com>...
> cs_lsolve and cs_ltsolve are included in Csparse,
> but not in CHOLMOD or LDL, as long as I can find. 
> cs_lsolve and cs_ltsolve in Csparse are pretty slow,
> compared to the triangular solvers in 
> CHOLMOD/LDL/mldivide. (For the triangular solve,
> mldivide was 3 times faster than cs_lsolve and cs_ltsolve 
> in Csparse, for my problem.)
> 

If b is sparse, as I think you're trying to do, CSparse will
be slow.  Its triangular solvers are designed for dense
right-hand-sides.  CSparse has a solver, cs_spsolve, for
sparse right-hand-side, but it doesn't have a MATLAB interface.

What if there are multiple sparse right-hand sides, will 
cs_spsolve still be faster than mldivide?

(Are you planning on giving cs_spsolve a MATLAB interface?)

"J Mathgeek" <jmathgeek@gmail.com> wrote in message
<fpmdhu$kng$1@fred.mathworks.com>...
> What if there are multiple sparse right-hand sides, will 
> cs_spsolve still be faster than mldivide?
> 
> (Are you planning on giving cs_spsolve a MATLAB interface?)
> 

In my book, I think that's an exercise left to the reader
.... :-)

It's trickier than it looks at first glance.  In part, it
requires an O(k*log(k)) time sort, where k = nnz(x) and
x=L\b is the result.  That sort is "useless",
algorithmically speaking.  It's there only because MATLAB
requires sparse matrices with sorted columns.  However, that
sort can dominate the total run time.  The cs_spsolve time
can be as low as Omega(k).

The mexFunction must not do any O(n) work, either, so it
can't do any "calloc" of workspaces ("malloc" instead).

"Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message 
<fp0fqv$3gd$1@fred.mathworks.com>...
> For the MEX experts out there: I would like to factor a 
huge
> matrix (passed in as a one-dimensional array) in-place in 
a
> mex file.  Unfortunately, the documentation says there may
> be side effects (and there are: the Matlab executive is
> apparently getting all mixed up).
> 
> Is there no way I could modify an argument in-place? If 
that
> is not possible, the penalty would be huge.  We are 
talking
> about a matrix of several gigabytes size...
> 
> Thanks for your help.
> 
> Petr



Back on this original point, I've haven't had any problems 
abusing the const nature of the data of a mex input.  You 
can unintentionally change other matlab variables, but this 
hasn't caused Matlab internal errors.  For example the mex 
file MyFun.c:

#include "mex.h"
#include "matrix.h"

void mexFunction(int nlhs, mxArray *plhs[], 
     int nrhs, const mxArray *prhs[])
{
    double *pData;
    int i;
    pData = mxGetPr(prhs[0]);  //Grab the data pointer

    for (i = 0; i < mxGetNumberOfElements(prhs[0]); i++)
        pData[i]++;
}

And then:

a = [1,2]; b = a;
MyFun(a);

a =

     2     3
b =

     2     3

b = b + 1;

a =

     2     3

b =

     3     4


No Matlab internal errors... If you try to change the data 
pointer of the input it spits the dummy though.

Cheers,
Andrew

"Andrew " <awbsmith@itee.uq.edu.au> wrote in message 
....
> Back on this original point, I've haven't had any problems 
> abusing the const nature of the data of a mex input.  You 
> can unintentionally change other matlab variables, but this 
> hasn't caused Matlab internal errors.  For example the mex 
> file MyFun.c:
> 

Danger Will Robinson, Danger!

:-)

Modifying the inputs in a mexFunction is a Really Bad Idea -
unless and until it's supported by The MathWorks.  You may
not think you're breaking something when you really are, if
you were to test the use of your mexFunction carefully.

Give me any mexFunction that can modify its inputs, and I
can merrily destroy MATLAB for you.  Or at least I can
destroy things you don't expect to be destroyed.

If I changed my example function to check the input is a 
non-sparse double precision matrix what could possibli go 
wrong? I'm not seeing how modifying a data value could be 
anything but opaque to Matlab.

Andrew  



"Tim Davis" <davis@cise.ufl.edu> wrote in message 
<fqcmvt$ss7$1@fred.mathworks.com>...
> "Andrew " <awbsmith@itee.uq.edu.au> wrote in message 
> ...
> > Back on this original point, I've haven't had any 
problems 
> > abusing the const nature of the data of a mex input.  
You 
> > can unintentionally change other matlab variables, but 
this 
> > hasn't caused Matlab internal errors.  For example the 
mex 
> > file MyFun.c:
> > 
> 
> Danger Will Robinson, Danger!
> 
> :-)
> 
> Modifying the inputs in a mexFunction is a Really Bad 
Idea -
> unless and until it's supported by The MathWorks.  You may
> not think you're breaking something when you really are, 
if
> you were to test the use of your mexFunction carefully.
> 
> Give me any mexFunction that can modify its inputs, and I
> can merrily destroy MATLAB for you.  Or at least I can
> destroy things you don't expect to be destroyed.

"Andrew " <awbsmith@itee.uq.edu.au> wrote in message
<fqfk55$t1a$1@fred.mathworks.com>...
> If I changed my example function to check the input is a 
> non-sparse double precision matrix what could possibli go 
> wrong? I'm not seeing how modifying a data value could be 
> anything but opaque to Matlab.
> 
> Andrew  
> 

It can fail.  It has nothing to do with sparsity.  What if
you do something like

x = fun(a+b)

where "fun" is your mangling mexFunction?  What if you do

a = b ;
x = fun (a) ;

then look at b?  You mexFunction has modified *both* a and b.

Modifying the inputs to a mexFunction is far from opaque to
MATLAB.

"Tim Davis" <davis@cise.ufl.edu> wrote in message 
> > 
> 
> It can fail.  It has nothing to do with sparsity.  What if
> you do something like

There are actually even worse failings. Here is an example
that completely destroys Matlab's representation of the
matrix that is being modified in place.  And still, all I'm
doing is over writing the elements of the array in which
Matlab holds the input matrix. So not only the variable
referencing is in danger here.

For details refer to
http://www.mathworks.com/matlabcentral/newsreader/view_thread/163748#415068

Petr

We're trying to make an inplace function, so we don't need 
to assign outputs.  Think of it as a C function with the 
form:

void myFun(double * pIn, int* pSize);

Making the inplace function return its input is really 
bad.  So put in an "if (nlhs > 0) then error" line if you 
like.  

fun(a+b) works fine, but since a+b makes a temporary 
mxArray you can't access the result.

As I mentioned before:
a = b ;
fun(a);

will modify a and b, so the programmers discretion is 
required.  But I think this thread started with someone 
wanting an inplace funtion because they couldn't hold two 
copies of the input in memory anyway, i.e.

a = not_inplace_fun(b) generates an out of memory error 
anyway...





"Tim Davis" <davis@cise.ufl.edu> wrote in message 
<fqrbve$jol$1@fred.mathworks.com>...
> "Andrew " <awbsmith@itee.uq.edu.au> wrote in message
> <fqfk55$t1a$1@fred.mathworks.com>...
> > If I changed my example function to check the input is 
a 
> > non-sparse double precision matrix what could possibli 
go 
> > wrong? I'm not seeing how modifying a data value could 
be 
> > anything but opaque to Matlab.
> > 
> > Andrew  
> > 
> 
> It can fail.  It has nothing to do with sparsity.  What if
> you do something like
> 
> x = fun(a+b)
> 
> where "fun" is your mangling mexFunction?  What if you do
> 
> a = b ;
> x = fun (a) ;
> 
> then look at b?  You mexFunction has modified *both* a 
and b.
> 
> Modifying the inputs to a mexFunction is far from opaque 
to
> MATLAB.

I took a quick look at this Petr and couldn't see where its 
doing anything in place.  You can typecast the pointers to 
the inputs to const and it still compiles...

const COLSOL_Index *MAXA = (const COLSOL_Index *)mxGetData
(pargin[1]);
const double *Ain= (const double*)mxGetPr(pargin[0]) ;


"Petr Krysl" <pkryslNOSP@Mucsd.edu> wrote in message 
<fqsv4a$npj$1@fred.mathworks.com>...
> "Tim Davis" <davis@cise.ufl.edu> wrote in message 
> > > 
> > 
> > It can fail.  It has nothing to do with sparsity.  What 
if
> > you do something like
> 
> There are actually even worse failings. Here is an example
> that completely destroys Matlab's representation of the
> matrix that is being modified in place.  And still, all 
I'm
> doing is over writing the elements of the array in which
> Matlab holds the input matrix. So not only the variable
> referencing is in danger here.
> 
> For details refer to
> 
http://www.mathworks.com/matlabcentral/newsreader/view_threa
d/163748#415068
> 
> Petr

"Andrew " <awbsmith@itee.uq.edu.au> wrote in message
<fr26sg$sr9$1@fred.mathworks.com>...
> I took a quick look at this Petr and couldn't see where its 
> doing anything in place.  You can typecast the pointers to 
> the inputs to const and it still compiles...
> 

Apologies.  I posted code which was tweaked in my attempts
to understand what was going on.  The correct code is

....
if (mxIsClass(pargin [0], "double")) {
        double *A= mxGetPr (pargin[0]) ;
        int n;
//         for (n = 0; n < MAXA [neq]; ++n) { printf ("A
[%d] =%g\n",n,A[n]); }
        if (!  (decomp_double(A, MAXA, neq)==0)) {
            mexErrMsgTxt ("colsoldecomp: A factorization
failed") ;
        }
//         for (n = 0; n < MAXA [neq]; ++n) { printf ("A
[%d] =%g\n",n,A[n]); }
   } else if (mxIsClass(pargin [0], "single")) {
     ...
    } 
....
This is a call to an in-place factorization routine. The
doubles in A are modified in place.  The correct mex routine
and the test data are posted in the archive


http://hogwarts.ucsd.edu/~pkrysl/public/colsol.zip


When compiled using colsol_make.m, the test reports 
decomposed matrix of the form

A =

  0.0e+000 *

    1.#INF
    1.#INF
   -1.#INF
 ...
   -1.#INF
         0
    1.#INF

I have even difficulty understanding what this printout means!

Petr

"Andrew " <awbsmith@itee.uq.edu.au> wrote in message
<fr2642$ltq$1@fred.mathworks.com>...
> We're trying to make an inplace function, so we don't need 
> to assign outputs.  Think of it as a C function with the 
> form:
> 
> void myFun(double * pIn, int* pSize);
> 
> Making the inplace function return its input is really 
> bad.  So put in an "if (nlhs > 0) then error" line if you 
> like.  
> 
> fun(a+b) works fine, but since a+b makes a temporary 
> mxArray you can't access the result.
> 
> As I mentioned before:
> a = b ;
> fun(a);
> 
> will modify a and b, so the programmers discretion is 
> required.  But I think this thread started with someone 
> wanting an inplace funtion because they couldn't hold two 
> copies of the input in memory anyway, i.e.
> 
> a = not_inplace_fun(b) generates an out of memory error 
> anyway...
> 
> 
> 

That's not my point.  A mexFunction that modifies its input
is free to return (another) output, or not.  The problem has
nothing to do with the outputs.  so

a = b ;
x = fun (b) ;

is just as bad (or good) as

a = b ;
fun (b) ;

Your test for "if (nargin > 0) mexErrMsgTxt ( )." won't
really help (nor will it make anything worse), unless you
are also trying to do

b = fun(b)

for example (or if inside the mexFunction you try to return
your input, as in plhs[0] = prhs[0]).

My point is that modifying your inputs is dangerous.  It's
not supported by The MathWorks.  Future versions of MATLAB
would be free to crash in a most merrily manner if told to
modify the inputs of a mexFunctions; your code would break
and you can't complain to The MathWorks for breaking your code.

Having said all that ... if you get it to work - that's
fine.  Just be aware that you are living outside the letter
of the law, so to speak, and have to be willing to see codes
die on you, if not now then maybe later on...

"Andrew " <awbsmith@itee.uq.edu.au> wrote in message
<fr2642$ltq$1@fred.mathworks.com>...
> We're trying to make an inplace function, so we don't need 
> to assign outputs.  Think of it as a C function with the 
> form:
> 
> void myFun(double * pIn, int* pSize);
> 
> Making the inplace function return its input is really 
> bad.  So put in an "if (nlhs > 0) then error" line if you 
> like.  
> 
> fun(a+b) works fine, but since a+b makes a temporary 
> mxArray you can't access the result.
> 
> As I mentioned before:
> a = b ;
> fun(a);
> 
> will modify a and b, so the programmers discretion is 
> required.  But I think this thread started with someone 
> wanting an inplace funtion because they couldn't hold two 
> copies of the input in memory anyway, i.e.
> 
> a = not_inplace_fun(b) generates an out of memory error 
> anyway...
> 
> 
> 

That's not my point.  A mexFunction that modifies its input
is free to return (another) output, or not.  The problem has
nothing to do with the outputs.  so

a = b ;
x = fun (b) ;

is just as bad (or good) as

a = b ;
fun (b) ;

Your test for "if (nargin > 0) mexErrMsgTxt ( )." won't
really help (nor will it make anything worse), unless you
are also trying to do

b = fun(b)

for example (or if inside the mexFunction you try to return
your input, as in plhs[0] = prhs[0]).

My point is that modifying your inputs is dangerous.  It's
not supported by The MathWorks.  Future versions of MATLAB
would be free to crash in a most merrily manner if told to
modify the inputs of a mexFunctions; your code would break
and you can't complain to The MathWorks for breaking your code.

Having said all that ... if you get it to work - that's
fine.  Just be aware that you are living outside the letter
of the law, so to speak, and have to be willing to see codes
die on you, if not now then maybe later on...

"Andrew " <awbsmith@itee.uq.edu.au> wrote in message
<fr2642$ltq$1@fred.mathworks.com>...
> We're trying to make an inplace function, so we don't need 
> to assign outputs.  Think of it as a C function with the 
> form:
> 
> void myFun(double * pIn, int* pSize);
> 
> Making the inplace function return its input is really 
> bad.  So put in an "if (nlhs > 0) then error" line if you 
> like.  
> 
> fun(a+b) works fine, but since a+b makes a temporary 
> mxArray you can't access the result.
> 
> As I mentioned before:
> a = b ;
> fun(a);
> 
> will modify a and b, so the programmers discretion is 
> required.  But I think this thread started with someone 
> wanting an inplace funtion because they couldn't hold two 
> copies of the input in memory anyway, i.e.
> 
> a = not_inplace_fun(b) generates an out of memory error 
> anyway...
> 
> 
> 

That's not my point.  A mexFunction that modifies its input
is free to return (another) output, or not.  The problem has
nothing to do with the outputs.  so

a = b ;
x = fun (b) ;

is just as bad (or good) as

a = b ;
fun (b) ;

Your test for "if (nargin > 0) mexErrMsgTxt ( )." won't
really help (nor will it make anything worse), unless you
are also trying to do

b = fun(b)

for example (or if inside the mexFunction you try to return
your input, as in plhs[0] = prhs[0]).

My point is that modifying your inputs is dangerous.  It's
not supported by The MathWorks.  Future versions of MATLAB
would be free to crash in a most merrily manner if told to
modify the inputs of a mexFunctions; your code would break
and you can't complain to The MathWorks for breaking your code.

Having said all that ... if you get it to work - that's
fine.  Just be aware that you are living outside the letter
of the law, so to speak, and have to be willing to see codes
die on you, if not now then maybe later on...

"Andrew " <awbsmith@itee.uq.edu.au> wrote in message
<fr2642$ltq$1@fred.mathworks.com>...
> We're trying to make an inplace function, so we don't need 
> to assign outputs.  Think of it as a C function with the 
> form:
> 
> void myFun(double * pIn, int* pSize);
> 
> Making the inplace function return its input is really 
> bad.  So put in an "if (nlhs > 0) then error" line if you 
> like.  
> 
> fun(a+b) works fine, but since a+b makes a temporary 
> mxArray you can't access the result.
> 
> As I mentioned before:
> a = b ;
> fun(a);
> 
> will modify a and b, so the programmers discretion is 
> required.  But I think this thread started with someone 
> wanting an inplace funtion because they couldn't hold two 
> copies of the input in memory anyway, i.e.
> 
> a = not_inplace_fun(b) generates an out of memory error 
> anyway...
> 
> 
> 

That's not my point.  A mexFunction that modifies its input
is free to return (another) output, or not.  The problem has
nothing to do with the outputs.  so

a = b ;
x = fun (b) ;

is just as bad (or good) as

a = b ;
fun (b) ;

Your test for "if (nargin > 0) mexErrMsgTxt ( )." won't
really help (nor will it make anything worse), unless you
are also trying to do

b = fun(b)

for example (or if inside the mexFunction you try to return
your input, as in plhs[0] = prhs[0]).

My point is that modifying your inputs is dangerous.  It's
not supported by The MathWorks.  Future versions of MATLAB
would be free to crash in a most merrily manner if told to
modify the inputs of a mexFunctions; your code would break
and you can't complain to The MathWorks for breaking your code.

Having said all that ... if you get it to work - that's
fine.  Just be aware that you are living outside the letter
of the law, so to speak, and have to be willing to see codes
die on you, if not now then maybe later on...

I learned about mex functions recently. It seems that mex 
function input arguments can indeed be modified in place.


This is s simple mex function:

// begin of source code
#include "Mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const 
mxArray *prhs[])
{
 double *pElement = mxGetPr(prhs[0]);

 for(unsigned int n = 0; n < mxGetNumberOfElements(prhs
[0]); n++)
 {
  pElement[n] = pElement[n] + 2;
 };

 plhs[0] = (mxArray*) prhs[0];

 return;
}
// end of source code


In Matlab you can do things like these:

>> a = [1 2]

>> a = 
       [1 2]

>> a = MexInline(a)

>> a = 
       [3 4]

>> a = magic(3)

a =

     8     1     6
     3     5     7
     4     9     2

>> a = 'Everything is fine.'

a =

Everything is fine.

Indeed everything has been fine until now. However, things 
can get messy if you type something like this:

>> a = MexInplace(a + 1)

a =

    4.0000    5.0000

Of course you are not going to use an inplace mex function 
this way. For the statement a + 1 at least one temporary 
variable must be created and you loose the advantages of 
your in place function. Moreover, the memory for a + 1 is 
released when the function returns and so is the memory of 
a. Matlab can crash if you access a in any way.

But what if you have forgotten how MexInplace works or 
someone else who isn't aware of the implementation of 
MexInplace is trying to use this function? You can't catch 
an expression like a + 1 inside the mex function. So you 
should document your mex function very well and give it a 
name that reminds the users of its effects.

I learned about mex functions recently. It seems that mex 
function input arguments can indeed be modified in place.


This is s simple mex function:

// begin of source code
#include "Mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const 
mxArray *prhs[])
{
 double *pElement = mxGetPr(prhs[0]);

 for(unsigned int n = 0; n < mxGetNumberOfElements(prhs
[0]); n++)
 {
  pElement[n] = pElement[n] + 2;
 };

 plhs[0] = (mxArray*) prhs[0];

 return;
}
// end of source code


In Matlab you can do things like these:

>> a = [1 2]

>> a = 
       [1 2]

>> a = MexInline(a)

>> a = 
       [3 4]

>> a = magic(3)

a =

     8     1     6
     3     5     7
     4     9     2

>> a = 'Everything is fine.'

a =

Everything is fine.

Indeed everything has been fine until now. However, things 
can get messy if you type something like this:

>> a = MexInplace(a + 1)

a =

    4.0000    5.0000

Of course you are not going to use an inplace mex function 
this way. For the statement a + 1 at least one temporary 
variable must be created and you loose the advantages of 
your in place function. Moreover, the memory for a + 1 is 
released when the function returns and so is the memory of 
a. Matlab can crash if you access a in any way.

But what if you have forgotten how MexInplace works or 
someone else who isn't aware of the implementation of 
MexInplace is trying to use this function? You can't catch 
an expression like a + 1 inside the mex function. So you 
should document your mex function very well and give it a 
name that reminds the users of its effects.

> Is there no way I could modify an argument in-place? If 
that
> is not possible, the penalty would be huge.  We are 
talking
> about a matrix of several gigabytes size...


What if you make your variable global, and ask for it's 
pointer from your mex file using mexGetArrayPtr ? (pass 
only the name of the var in the mex)
That looks more safe, although still not supported.
Gergely Katona

"Gergely Katona" <asdasd15@freemail.hu> wrote in message
<ft2jps$8oq$1@fred.mathworks.com>...
> 
> > Is there no way I could modify an argument in-place? If 
> that
> > is not possible, the penalty would be huge.  We are 
> talking
> > about a matrix of several gigabytes size...
> 
> 
> What if you make your variable global, and ask for it's 
> pointer from your mex file using mexGetArrayPtr ? (pass 
> only the name of the var in the mex)
> That looks more safe, although still not supported.
> Gergely Katona

That doesn't make any difference.

This thread is dragging on far past its useful life.

The final word is, yes, you "can" modify your inputs, but
doing so can have indeterminate side effects.  It is
unsupported.  Do so only at your own peril.  Modifying the
inputs to your mexFunction is  very dangerous, but sometimes
useful.  It is never safe.  It should never be done in a
mexFunction that you intend for others to use.

> What if you make your variable global, and ask for it's 
> pointer from your mex file using mexGetArrayPtr ? (pass 
> only the name of the var in the mex)
> That looks more safe, although still not supported.
> Gergely Katona

That doesn't work. Global Variables can't be accessed by 
mex files, because they are unknown there.