From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 874 invoked by alias); 27 Sep 2009 08:03:19 -0000 Received: (qmail 671 invoked by uid 22791); 27 Sep 2009 08:03:17 -0000 X-SWARE-Spam-Status: No, hits=-2.1 required=5.0 tests=AWL,BAYES_00 X-Spam-Check-By: sourceware.org Received: from smtp4.welho.com (HELO smtp4.welho.com) (213.243.153.38) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Sun, 27 Sep 2009 08:03:12 +0000 Received: from [10.0.0.133] (cs181229105.pp.htv.fi [82.181.229.105]) by smtp4.welho.com (Postfix) with ESMTP id 239E15BC092; Sun, 27 Sep 2009 11:03:09 +0300 (EEST) Message-ID: <4ABF1C3C.6070801@iki.fi> Date: Sun, 27 Sep 2009 08:03:00 -0000 From: Tuomo Keskitalo User-Agent: Mozilla-Thunderbird 2.0.0.22 (X11/20090707) MIME-Version: 1.0 To: Gerard Jungman CC: GSL Discuss Mailing List Subject: Re: new double precision data structure? References: <48E25CA9.6080306@iki.fi> <490DE4BD.7070907@iki.fi> <497B00F6.2080400@iki.fi> <498727E5.6080407@iki.fi> <49AA9DB5.6030908@iki.fi> <49FB01D1.30000@iki.fi> <4A7ADFDC.9080408@iki.fi> <1251414774.23092.80.camel@manticore.lanl.gov> <1251414939.23092.82.camel@manticore.lanl.gov> <1253062179.23092.971.camel@manticore.lanl.gov> In-Reply-To: <1253062179.23092.971.camel@manticore.lanl.gov> Content-Type: multipart/mixed; boundary="------------000500050609060806010103" 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-q3/txt/msg00067.txt.bz2 This is a multi-part message in MIME format. --------------000500050609060806010103 Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Content-length: 674 On 09/16/2009 03:49 AM, Gerard Jungman wrote: > It smells like a view, and indeed it is. But it is a notion of > view which is more central to the design than the current one. > Everything should be a view. I think we really need to sort the data structure out now. I'm not sure what you have in mind, but to get some discussion going on, I've attached an (incomplete and untested) idea for a general n-dimensional double precision data structure. Is this anything like what you pictured? I'd like to hear people's comments on this. Would something like this be of use? Feel free to give your full critique. -- Tuomo.Keskitalo@iki.fi http://iki.fi/tuomo.keskitalo --------------000500050609060806010103 Content-Type: text/x-csrc; name="newblock.c" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="newblock.c" Content-length: 4355 /* newblock.c - An idea for n-dimensional double precision data structure. This is pseudocode to describe the idea, does not compile. */ /* Structure for n-dimensional double precision data storage */ typedef struct { size_t dim; /* data structure dimension */ size_t * sizes; /* dimensions of system, e.g. {3,4,5} for real 3*4*5 tensor */ double * data; /* continuous data storage */ size_t layout; /* data layout type - LAST_INDEX_CONTINUOUS = "column major" for matrices, or FIRST_INDEX_CONTINUOUS = "row major" (default) for matrices. Note: It might be possible to extend this to support other storage schemes, e.g. diagonal/banded storage BANDED_INDEX_CONTINUOUS. */ } base_block; /* Structure for a generic n-dimensional view object to access base block data. Access to data storage is done via this view object. There can be multiple views referring to the same data storage. */ typedef struct { size_t * firstcell; /* indices to first access cell (minimum indices) */ size_t * lastcell; /* indices to last access cell (maximum indices) */ size_t * stride; /* strides for each dimension for the view */ base_block * block; double * data; /* pointer to block->data */ size_t owner; /* indicator for view owning the data */ } base_view; /* vector object based on base block and view */ typedef struct { base_view * view; } vector; /* matrix object based on base block and view */ typedef struct { base_view * view; } matrix; matrix * matrix_alloc (size_t n1, size_t n2) { /* Allocates matrix object (both view and storage objects) in row-major layout. */ // 1. allocate base_block with dim=2, sizes={n1,n2} // 2. allocate base_view with firstcell={1,1}, // lastcell={n1,n2}, stride={0,0} (=whole matrix) // 3. set view owner=TRUE // 4. return base_view; } matrix * matrix_alloc_with_layout (size_t n1, size_t n2, size_t layout) { /* Allocates matrix object (both view and storage objects) in the specified layout. */ } int vector_alloc_view_of_matrix (vector * v, matrix * m, size_t * firstcell, size_t * lastcell, size_t * stride) { /* Creates a new vector view of the matrix m */ // 1. allocate base_view with given firstcell, lastcell and strides // 2. associate view with the storage // 3. set owner=FALSE // 4. return base_view; } /* Example code */ int main (void) { matrix * a; matrix * b; vector * v; int retval; /* allocate matrix with 3 rows and 4 columns */ a = matrix_alloc (3, 4); /* User can change the layout between row- and column-major at any point. After change to column-major, those GSL functions that require row-major storage would raise error if the matrix is passed to them. */ retval = matrix_set_layout (a, LAST_INDEX_CONTINUOUS); { double temp[] = { 0, 10, 20, 01, 11, 21, 02, 12, 22, 03, 13, 23 }; size_t templen = 12; retval = matrix_set_doublevec (a, temp, templen); /* This would initialize matrix a = [00 01 02 03 10 11 12 13 20 21 22 23] */ } /* Set up a vector view of the diagonal */ const size_t b0[] = { 1, 1 }; const size_t e0[] = { 3, 3 }; const size_t s0[] = { 1, 1 }; retval = vector_alloc_view_of_matrix (v, a, &b0, &e0, &s0); /* Access cells in matrix through the vector */ double val; val = vector_get(v, 1); // returns 00 from matrix a val = vector_get(v, 2); // returns 11 from matrix a val = vector_get(v, 3); // returns 22 from matrix a retval = vector_set(v, 2, 55); // sets value 55 in place of 11 to matrix a size_t l; l = vector_length(v); // returns 3; double diag[3]; retval = vector_get_doublevec (&diag, v); // copies diagonal element values to diag } /* complex number matrix (one possible implementation) based on base block and view */ typedef struct { base_view * view; } complex_matrix; complex_matrix * complex_matrix_alloc (size_t n1, size_t n2) { /* Allocates matrix object (both view and storage objects) */ // 1. allocate base_block with dim=3, sizes={n1,n2,2}, so that // the last dimension spans the real and imaginary parts // 2. allocate base_view with firstcell={1,1,1} and // lastcell={n1,n2,2} (=whole matrix) // 3. set view owner=TRUE // 4. return base_view; } --------------000500050609060806010103--