From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 24906 invoked by alias); 5 Oct 2009 23:45:33 -0000 Received: (qmail 24884 invoked by uid 22791); 5 Oct 2009 23:45:32 -0000 X-SWARE-Spam-Status: No, hits=-2.0 required=5.0 tests=BAYES_00,SARE_MSGID_LONG40,SPF_PASS X-Spam-Check-By: sourceware.org Received: from mail-yw0-f177.google.com (HELO mail-yw0-f177.google.com) (209.85.211.177) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Mon, 05 Oct 2009 23:45:24 +0000 Received: by ywh7 with SMTP id 7so1595690ywh.24 for ; Mon, 05 Oct 2009 16:45:17 -0700 (PDT) MIME-Version: 1.0 Received: by 10.101.193.25 with SMTP id v25mr756330anp.132.1254786317572; Mon, 05 Oct 2009 16:45:17 -0700 (PDT) In-Reply-To: <1254783367.28192.98.camel@manticore.lanl.gov> References: <1254708349.18519.4.camel@ForbiddenPlanet> <7f1eaee30910050750l738876b1p41e6bd8ae5aa6d16@mail.gmail.com> <1254783367.28192.98.camel@manticore.lanl.gov> Date: Mon, 05 Oct 2009 23:45:00 -0000 Message-ID: <7f1eaee30910051645v5654353fj53416ed88b98844@mail.gmail.com> Subject: Re: containers tentative design summary From: James Bergstra To: Gerard Jungman Cc: gsl-discuss@sourceware.org Content-Type: text/plain; charset=ISO-8859-1 Content-Transfer-Encoding: quoted-printable Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2009-q4/txt/msg00005.txt.bz2 On Mon, Oct 5, 2009 at 6:56 PM, Gerard Jungman wrote: > On Mon, 2009-10-05 at 10:50 -0400, James Bergstra wrote: >> Two comments: >> >> I'm a bit rusty with my C structs... but you would need two distinct >> static classes to have const and non-const data pointers for your view >> right? > > Yes, sorry. I left that out. > > >> Also, it sounds like writing code that will work for a tensor of any >> rank (e.g. add two tensors together) might be either tedious or >> impossible. =A0I recognize that part of the problem is the lack of >> templating and polymorphism, but it would at least be comforting to >> see just how bad the situation is via a use case or two in the design >> documentation. =A0 I (naively?) fear that to get good performance will >> require a whole library of functions for even the most basic of >> operations. >> >> gsl_marray_add_1_0( gsl_marray_1, double ); >> gsl_marray_add_1_1( gsl_marray_1, gsl_marray_1); >> gsl_marray_add_1_2( gsl_marray_1, gsl_marray_2); >> gsl_marray_add_2_2(... ) >> ... >> >> gsl_marray_sub_1_0( ... ) >> >> Maybe a system of macros could be designed to help here, but it sounds >> like it will never be as easy as writing a couple of for-loops. > > > First, I want to be sure that we distinguish the multi-array > case from the vector/matrix/tensor case. Algebra for vectors and > matrices is well-defined and already on place; we only need to > re-write some wrappers, etc. It is handled by BLAS calls. As I > mentioned, this places a constraint on matrices, that the > fast traversals be of unit stride. It seems like we just > have to live with that constraint, for the matrix type. > > Addition, subtraction, scaling of multi-arrays is not hard > because it is only defined within the same rank. So there > is only a linear complexity, and no combinatorial explosion > in the interfaces, for these operations. Of course, there > are issues of shape conformance at run-time, but that is > also true for vectors and matrices. > > That leaves multi-array algebra as an open question. By algebra, > I mean technically the "cross-rank" operations, which form > some kind of general multi-linear algebra. Sounds pretty > hairy, as you suspect. > Another idea here would be to implement general-cases in terms of structures with a dynamic rank right? The rank does not need to be part of the struct, but can be a variable. I'm thinking of something like struct marray { int rank; int * dims; int * strides; dtype * base; }; marray_add_inplace( const marray a, const marray b, marray c); //N.B. the passing on the stack Inside that function, we can check for special cases of interest and optimize them as necessary. One disadvantage of such a structure I think you'll agree, is the internal pointers, which interfere with convenient stack allocation. A somewhat ugly alternative (which I am nonetheless suggesting) would be the following... { int rank; dtype * base; int dim0, dim1, dim2; int str0, str1, str2; int * extra_dims; int * extra_strides; }; For low-rank multiarrays, this structure can be manipulated on the stack, so it is possible to write expressions like marray_add_inplace(marray_reshape2(a, 5, 50), marray_slice(b, 17), marray_transpose(c)); // add a reshaped view of a to the 17th matrix-slice of b, storing result in a transposed view of c. At the same time, arbitrarily high ranks are supported if the user is willing to put up with the syntax of heap-based manipulation. (Of course we can haggle over how many extra ranks can be built into the static part of the structure.) James --=20 http://www-etud.iro.umontreal.ca/~bergstrj