From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 9150 invoked by alias); 6 Jan 2010 15:47:52 -0000 Received: (qmail 9130 invoked by uid 22791); 6 Jan 2010 15:47:50 -0000 X-SWARE-Spam-Status: No, hits=-2.5 required=5.0 tests=AWL,BAYES_05,RCVD_IN_DNSWL_LOW X-Spam-Check-By: sourceware.org Received: from mail.phy.duke.edu (HELO mail.phy.duke.edu) (152.3.182.2) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Wed, 06 Jan 2010 15:47:45 +0000 Received: from localhost (localhost [127.0.0.1]) by mail.phy.duke.edu (Postfix) with ESMTP id 32AB2780B6; Wed, 6 Jan 2010 10:47:37 -0500 (EST) Received: from mail.phy.duke.edu ([127.0.0.1]) by localhost (mail.phy.duke.edu [127.0.0.1]) (amavisd-new, port 10026) with LMTP id ITBin0-h4RI3; Wed, 6 Jan 2010 10:47:37 -0500 (EST) Received: from lilith.rgb.private.net (cpe-069-134-079-008.nc.res.rr.com [69.134.79.8]) by mail.phy.duke.edu (Postfix) with ESMTP id EA85A780B2; Wed, 6 Jan 2010 10:47:36 -0500 (EST) Date: Wed, 06 Jan 2010 15:47:00 -0000 From: "Robert G. Brown" To: Tuomo Keskitalo cc: Brian Gough , Gerard Jungman , GSL Discuss Mailing List Subject: Re: gsl container designs In-Reply-To: <4B4477B8.50305@iki.fi> Message-ID: References: <1259110486.3028.69.camel@manticore.lanl.gov> <4B4477B8.50305@iki.fi> User-Agent: Alpine 2.00 (LFD 1167 2008-08-23) MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed 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: 2010-q1/txt/msg00003.txt.bz2 On Wed, 6 Jan 2010, Tuomo Keskitalo wrote: > Here is a crazy idea: If we take Gerards design 1, would it be too absurd to > use a horrible big wrapper struct like > > typedef struct { > gsl_marray *m; > gsl_const_marray *cm; > gsl_marray_1 *m1; > gsl_const_marray_1 *cm1; > gsl_marray_2 *m2; > gsl_const_marray_2 *cm2; > gsl_marray_3 *m3; > gsl_const_marray_3 *cm3; > gsl_vector *vec; > gsl_const_vector *cvec; > gsl_matrix *mat; > gsl_const_matrix *cmat; > ... > } gsl_container; > > and make everything a gsl_container, and then always use them like > &a.mat or &b.vec? Maybe the const types could be in another > gsl_const_container struct.. > > One downside is that this would effectively mask the type of a or b in > program code (gsl_container a; // is it a vector, matrix, or marray?), which > might make reading code a pain.. IMO the primary reason for "fixing" vector/matrix objects in the GSL is to enable the single vector that represents the actual data in the object to be presented to the user in a matrix-intuitive way, even at some cost in code efficiency (and to try to minimize that cost or provide further benefits to counterbalance it by enabling efficient linear algebra routines to be matched to it, for example). Any of us who have coded for a decade or two know perfectly well how to address a vector via explicit displacement arithmetic, burying horrible constructs in one's code to compute the displacement of an i,j,k element in a vector containin aa "3d matrix". Or, one can use C pointer constructs that are equally unreadable for the same general purpose. We have also had to DEBUG this sort of thing, which is extremely unpleasant, or understand somebody else's displacement code, which is often unpleasant to the point of being nightmarish. What I want is extremely simple: * A way to construct/allocate a matrix of arbitrary dimension (out to some (un)reasonable max, e.g. 10, likely to be higher than anybody needs for coding purposes at least at this point, and deconstruct/free the matrix when I'm done with it, ideally with a single call in each case, ideally for any type of variable. I'm perfectly comfortable with requiring a pointer and cast to establish the type, but I want to be able to use it just like m[i][j][k]...[r][s] when all is said and done. * The specific dimension bounds need to be user specified. If I want to make a 6 dimensional triple triangular array to hold e.g. 3j symbols, I don't want to have to allocate a rectangular array most of which is empty wasted space. I'm not talking about sparse matrix here, just non-rectangular arrays. Angular momentum indices are perfect (and common) examples. * Access to the underlying actual data vector (non-opaque, in other words) so it can be passed to e.g. ODE solvers and the like that can only be sensibly written in a general way for a vector. I want to be able to pass the ODE solver the vector and its dimension (and perhaps an auxiliary similar vector for the derivatives with similar packing and indexing), but I want to be able to write the ODE deriv routine using the readable matrix form. None of this is particularly difficult to accomplish. I do it all the time in my own code. In the end, one can integrate sets of ODEs where the derivatives are indexed by angular momentum l,m, where the deriviative expression contains sums over l',m', l",m" with precomputed tables of e.g. Gaunt numbers, and where the derivative expression is HUMAN READABLE and directly comparable to an expression in a physics paper or textbook. My problem in this regard with the GSL is that so far, the matrix and vector constructs are pretty much useless -- they are far away from the standard way one usually does this sort of thing in C using pointers or just allocating a vector or matrix. One of the ADVANTAGES of C is that it permits one to work pointer magic to improve the readability (and sometimes execution efficiency) of code. Opaque types and set/get calls are dark evil in both regards. I just want to be able to do (double) I[l][m][l'][m'][l"][m"]aand/or (int) epsilon[i][j][k] transparently and reasonably efficiently, because both of these occur all the time in real numerical code tht I've written. I can do it myself and still use the GSL, and have in the past done so. For ODE solvers this isn't much of a problem. For linear algebra, however, it starts being a problem. Then issues like column-major become important and have to be solved in a standardized way. This would make even my own code more efficient, as well, as most of what I do with the I() matrix above is basically matrix multiplication or contraction -- summing out shared indices -- and linear algebra routines that operate on correctly allocated structures can be 2x-3x more efficient than ones that sum across big steps in the underlying vectors. rgb > > -- > Tuomo.Keskitalo@iki.fi > http://iki.fi/tuomo.keskitalo > Robert G. Brown http://www.phy.duke.edu/~rgb/ Duke University Dept. of Physics, Box 90305 Durham, N.C. 27708-0305 Phone: 1-919-660-2567 Fax: 919-660-2525 email:rgb@phy.duke.edu