--- Alexei Sheplyakov <varg@theor.jinr.ru> escribió:

Hello,

On Sun, Feb 18, 2007 at 12:02:25AM -0300, Charlls Quarra wrote:

...

Why am i trying to do this as an indexed object and not with a matrix?

well, because later, i would also want to make computations with an object like

D( f_{i} ) / dx_{j} dx_{k}

What about using vector (std::vector) of matrices (or GiNaC::lst of thereof)?

right now im trying to do something like this with this recursive function whose purpose is to produce an object that can be accessed like ex_to<matrix>( ex_to<matrix>(mat( i, j ))( k , l ) ) ( m , n ).. which at its time is wrapped with () operator of the correct size. However there is a problem when constructing a matrix of 3 indices and range 2. Concretely the assignment of the submatrices obtained from the tail recursion fail: matrix multi_index_matrix( int n , int range ) { cout << " constructing matrix of " << n << " indices of range " << range << endl; if (n == 1) { return matrix( 1 , range ); //return; } if (n == 2) { cout << " matrix of " << range << "x" << range << endl; return matrix( range , range ); //return; } if (n % 2 == 1) { // *put = matrix(1 , range); matrix put( 1 , range ); cout << " top odd matrix " << endl; for (int i=0; i < range; i++) { put( 1 , i ) = multi_index_matrix( n-1 , range ); cout << " done on index " << i << endl; } return put; } matrix put( range , range ); for (int i=0; i< range ; i++) { for (int j=0; j < range; j++) { put( i , j ) = multi_index_matrix( n - 2 , range ); } } return put; }; int main() { matrix foo = multi_index_matrix( 3 , 2 ); cout << "done! " << endl; } __________________________________________________ Preguntá. Respondé. Descubrí. Todo lo que querías saber, y lo que ni imaginabas, está en Yahoo! Respuestas (Beta). ¡Probalo ya! http://www.yahoo.com.ar/respuestas

I just noted that the following write access attempt will fail: matrix foo(1 , 2); for (int i=0 ; i< 2 ; i++) { foo(0,i) = matrix(2,2); } ex_to<matrix>( foo(0,1) )(1,1) = 34; cout << "it was wrote -> " << ex_to<matrix>( foo( 0 , 1 ) )(1,1) << endl; i get a compilation error: passing 'const GiNaC::ex' as 'this' argument of 'GiNaC::ex& GiNaC::ex::operator=(const GiNaC::ex&)' discards qualifiers any idea how to write to the entries of a submatrix in this way? __________________________________________________ Preguntá. Respondé. Descubrí. Todo lo que querías saber, y lo que ni imaginabas, está en Yahoo! Respuestas (Beta). ¡Probalo ya! http://www.yahoo.com.ar/respuestas

On Mon, Feb 19, 2007 at 11:38:20PM -0300, Charlls Quarra wrote:

ok i found that this awfully-looking hack does the job:

( *(matrix*) & ex_to<matrix>(foo(0,1)) )(1,1) = 34;

Ugh, this is plain EVIL.

im eager to find out what its the ginac-sanctioned way to do it though

As I explained in the previous mail, there is no way to assign anything to references returned by ex_to<T>. If you really need multidimensional arrays (and vector<matrix> or vector<vector<vector <ex> > > are not good enough) then GiNaC is not proper tool for your task. Best regards, Alexei -- All science is either physics or stamp collecting.

Hello! On Tue, Feb 20, 2007 at 09:44:54AM -0300, Charlls Quarra wrote:

This morning i found a non-hacky (less-hacky?) way to do it.

Could you please explain why vector<matrix> is not good enough? I admit there might be some reason for that, but I can't see any.

matrix foo(1 , 2); for (int i=0 ; i< 2 ; i++) { foo(0,i) = matrix(2,2); } matrix sub = ex_to<matrix>(foo( 0, 1));

[unnecessary] copy operation

sub( 1 , 1) = 34; foo( 0 , 1 ) = sub;

yet another [unncessary] copy operation

cout << "yes --> " << ex_to<matrix>( foo( 0 , 1 ) )(1,1) << endl;

if this is breaks refcounting too,

It does not. But it involves two copy operations...

then i guess one should ask for a revision why is that matrix entries behave like unidirectional membranes (you can put but cannot change)

Actually this is documented in the manual: see the section "Expressions are reference counted" in the the "Internal structures" chapter. Since I dislike RTFM as an answer (although I think reading manuals *does* help), I'll briefly summarise it here. GiNaC is designed to handle large expressions, so at minimum it needs to store them in efficient way. Reference counting and copy on write semantics are simple and powerful techniques to achive this. The price to pay is "undirectional membrane" behaviour (no direct write access to class data).

you make it sound quite final though, so i'll assume you meant that there was no way to do it strictly thru the ex_to<> interface, and no that there was no way to do it "at all"

It is impossible to do it in sane, efficient way. This is _feature_.

of course that this leaves unaswered how to perform simple assignments like t3_idx( qi , qj , qk ) = foo from the wrapper interface object... I guess ill have to think about that some more

1) Use some C++ container (with ex as an element type) 2) Write your own class which provides such an operator (using GiNaC::matrix as a starting point) 3) Use some other [template] library which provide multidimensional arrays (with ex as an element type). Everything else should be a bug (in GiNaC, not in your code) or evil hack (a la reinterpret_cast). Best regards, Alexei -- All science is either physics or stamp collecting.

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.

--- Chris Dams <Chris.Dams@mi.infn.it> escribió:

Why don't you simply invent a nice name for your multi-index object (say "M") and store the values of them in an exmap. As in

exmap m; m[indexed(M, idx(1,3), idx(2,3), idx(3,3))] = value_of_entry_1_2_3; After you have filled your exmap with all the indices you like, you can use the exmap in a call to .subs to replace tensor objects with numerical indices bij the appropriate value.

Thats a great way to do it when the object is quite sparse, i'll keep that way in mind for special cases, thanks

Does this help?

it did. Thanks __________________________________________________ Preguntá. Respondé. Descubrí. Todo lo que querías saber, y lo que ni imaginabas, está en Yahoo! Respuestas (Beta). ¡Probalo ya! http://www.yahoo.com.ar/respuestas