[GiNaC-list] write access to entries of submatrices [Was: multi
index storage of expressions]
Sheplyakov Alexei
varg at theor.jinr.ru
Wed Feb 21 12:12:32 CET 2007
Hello!
On Tue, Feb 20, 2007 at 04:16:08PM -0300, Charlls Quarra wrote:
> I know i'm new to the list and probably you feel you have to be so
> overly cynical and aggresive to make some sort of social standpoint.
I'm sorry, no offence was intended... My point was that GiNaCs design
makes some operations inefficient (and/or ugly), so GiNaC is not proper
tool to perform them. But apparently such operations are not necessary
for your task. So I wonder what is exactly aggressive and/or cynical in
my posts...
> Oh well, first you said that it wasn't possible to do
> it in ginac and i should look somewhere else, now you
> speak about it being inefficient.
Sure it is possible _in principle_, as in "the determined Real Programmer
can write FORTRAN programs in any language". As for me "being unable to do
something in a natural and efficient way" and "impossible" is essentially
the same thing.
> I happen to have read that chapter that you RTFM me
> before asking about write access, and precisely on
> that page it says that copy is done lazily, so in fact
> the copy operation:
>
> matrix sub = ex_to<matrix>( mat(1,2) );
>
> should be "light-weight", at least if the matrix is
> really just a handle like the ex.
Matrices are unshareable by default (because they provide write access
to data):
/** Very common ctor. Initializes to r x c-dimensional zero-matrix.
*
* @param r number of rows
* @param c number of cols */
matrix::matrix(unsigned r, unsigned c)
: inherited(TINFO_matrix), row(r), col(c), m(r*c, _ex0)
{
setflag(status_flags::not_shareable);
}
so the manual is somewhat inaccurate...
> now, when i do
> sub(3,4) = 15;
>
> this is arguably doing a copy (again, according to the
> RTFMed chapter) because both sub and mat(1,2) point to
> the same object,
They don't. For illustration consider the following code:
// matmat.cpp
#include <iostream>
#include <string>
#include <ginac/ginac.h>
#include <stdexcept>
using namespace std;
using namespace GiNaC;
void test_ex()
{
symbol x("x"), y("y");
ex foo = pow(x,2)*y + pow(y,2)*x + pow(y,3) + pow(x,3);
cout << foo << endl;
ex bar = foo;
bar += pow(x,2)*pow(y,2);
cout << bar << endl;
}
void test_matr()
{
matrix A = ex_to<matrix>(symbolic_matrix(2, 2, string("a")));
cout << A << endl;
matrix B = A;
B(0, 0) = ex(1);
cout << B << endl;
}
int main(int argc, char** argv)
{
test_ex();
test_matr();
return 0;
}
Let's run it in the debugger:
$ cat matmat.gdb
break matmat.cpp:12
break matmat.cpp:14
break matmat.cpp:15
break matmat.cpp:21
break matmat.cpp:23
break matmat.cpp:24
commands 1
print foo.bp.p.refcount
end
commands 2
print foo.bp.p.refcount
print bar.bp.p.refcount
end
commands 3
print foo.bp.p.refcount
print bar.bp.p.refcount
end
commands 4
print A.refcount
end
commands 5
print A.refcount
print B.refcount
print A.m
print B.m
end
commands 6
print A.refcount
print B.refcount
print A.m
print B.m
end
$ echo $CXXFLAGS
-Wall -O2 -g -pipe -funroll-loops -ffast-math -fomit-frame-pointer -finline-limit=1200 -m32 -march=pentium4 -mmmx -msse2 -mfpmath=sse
$ g++-4.1 $CXXFLAGS -o matmat matmat.cpp -lginac
$ gdb --silent ./matmat
Using host libthread_db library "/lib/tls/i686/cmov/libthread_db.so.1".
(gdb) source matmat.gdb
Breakpoint 1 at 0x80491fb: file matmat.cpp, line 12.
Breakpoint 2 at 0x8049225: file matmat.cpp, line 14.
Breakpoint 3 at 0x804932a: file matmat.cpp, line 15.
Breakpoint 4 at 0x8049a25: file matmat.cpp, line 21.
Breakpoint 5 at 0x8049c47: file matmat.cpp, line 23.
Breakpoint 6 at 0x8049ca7: file matmat.cpp, line 24.
(gdb) run
Starting program: /afs/diastp.jinr.ru/user/varg/tmp/matmat
Breakpoint 1, test_ex () at matmat.cpp:12
12 cout << foo << endl;
$1 = 1
(gdb) cont
Continuing.
x^3+y*x^2+y^3+y^2*x
Breakpoint 2, test_ex () at matmat.cpp:14
14 bar += pow(x,2)*pow(y,2);
$2 = 2
$3 = 2
// there is only one copy with reference count 2
(gdb)
Continuing.
Breakpoint 3, test_ex () at matmat.cpp:15
15 cout << bar << endl;
$4 = 1
$5 = 1
// now there are two different expressions
(gdb)
Continuing.
x^3+y*x^2+y^3+y^2*x^2+y^2*x
Breakpoint 4, test_matr () at matmat.cpp:21
21 cout << A << endl;
$6 = 0
// matrices are not sharable, so reference counting is pointless
(gdb)
Continuing.
[[a00,a01],[a10,a11]]
Breakpoint 5, test_matr () at matmat.cpp:23
23 B(0, 0) = ex(1);
$7 = 0
$8 = 0
$9 = {<std::_Vector_base<GiNaC::ex,std::allocator<GiNaC::ex> >> = {
_M_impl = {<std::allocator<GiNaC::ex>> = {<__gnu_cxx::new_allocator<GiNaC::ex>> = {<No data fields>}, <No data fields>},
_M_start = 0x805d710, _M_finish = 0x805d720, _M_end_of_storage = 0x805d720}}, <No data fields>}
$10 = {<std::_Vector_base<GiNaC::ex,std::allocator<GiNaC::ex> >> = {
_M_impl = {<std::allocator<GiNaC::ex>> = {<__gnu_cxx::new_allocator<GiNaC::ex>> = {<No data fields>}, <No data fields>},
_M_start = 0x805a9f0, _M_finish = 0x805aa00, _M_end_of_storage = 0x805aa00}}, <No data fields>}
// Adresses of actual storage (_M_start, _M_finish) are different for
// A and B/ Both A and B are _different_ objects
(gdb)
Continuing.
Breakpoint 6, test_matr () at matmat.cpp:24
24 cout << B << endl;
$11 = 0
$12 = 0
$13 = {<std::_Vector_base<GiNaC::ex,std::allocator<GiNaC::ex> >> = {
_M_impl = {<std::allocator<GiNaC::ex>> = {<__gnu_cxx::new_allocator<GiNaC::ex>> = {<No data fields>}, <No data fields>},
_M_start = 0x805d710, _M_finish = 0x805d720, _M_end_of_storage = 0x805d720}}, <No data fields>}
$14 = {<std::_Vector_base<GiNaC::ex,std::allocator<GiNaC::ex> >> = {
_M_impl = {<std::allocator<GiNaC::ex>> = {<__gnu_cxx::new_allocator<GiNaC::ex>> = {<No data fields>}, <No data fields>},
_M_start = 0x805a9f0, _M_finish = 0x805aa00, _M_end_of_storage = 0x805aa00}}, <No data fields>}
// The address of storage did not change when we assigned something to B
// (_M_start and _M_finish are the same is in the previous breakpoint)
(gdb)
Continuing.
[[1,a01],[a10,a11]]
Program exited normally.
(gdb)
> but the write operation on a handle triggers the copy.
> However, now that we are talking about _efficiency_,
> this is not unavoidable, at least from what
> refcounting refers. If one would do:
>
> matrix sub = ex_to<matrix>( mat( 1 , 2 ) );
> mat( 1 , 2 ) = 0; // <-- this kills the first
> reference
> sub(3, 4) = 15; // no copy should be involved anymore,
> there is a single reference to the underlying matrix
> mat( 1 , 2 ) = sub;
>
> so in this case, Apparently no deep copy should be
> involved.
Unfortunately, you are wrong: for matrices there is a copy operation
in every assignment. Usually this is not a big deal: the actual elements
of matrices are ex (which are usually shareable), so copying is not that
expansive (but still there are nxm handles to copy). Storing non-shareable
objects (like matrices) in matrix is _very_ inefficient.
> I actually wasnt trying to be fixed on doing this with matrix or ginac
> objects only other than esthetical/completeness reasons: I even could
> come around to consider doing it with std::vector after all (or even
> better, blitz), but that was off the point by several miles.
> I was really just trying to learn ginac by doing something that i
> would need commonly.
And I was trying to explain that one don't have to treat every object
as an expression (like in CASes). If it is natural to store (and manipulate)
your data as an array of expressions (or whatever else) -- just use that.
There is no need to pack <your data type> it into yet another expression.
> I was also intrigued if others don't seem to require multi indexing
> (i use it all the time on maple using lists of lists, but it doesnt
> seem to be efficient with deep stacking).
I think because STL containers provide this functionality.
Best regards,
Alexei
--
All science is either physics or stamp collecting.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 827 bytes
Desc: Digital signature
Url : http://www.cebix.net/pipermail/ginac-list/attachments/20070221/8a03926a/attachment.pgp
More information about the GiNaC-list
mailing list