Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]
Groups > comp.lang.c > #61593 > unrolled thread
| Started by | luser droog <luser.droog@gmail.com> |
|---|---|
| First post | 2015-05-05 23:18 -0700 |
| Last post | 2015-05-25 00:13 -0700 |
| Articles | 20 on this page of 50 — 10 participants |
Back to article view | Back to comp.lang.c
Column slices luser droog <luser.droog@gmail.com> - 2015-05-05 23:18 -0700
Column slices asetofsymbols@gmail.com - 2015-05-06 00:26 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-05-06 02:33 -0700
Re: Column slices Bartc <bc@freeuk.com> - 2015-05-06 10:46 +0100
Re: Column slices asetofsymbols@gmail.com - 2015-05-07 11:02 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-07 14:05 -0700
Re: Column slices Nobody <nobody@nowhere.invalid> - 2015-05-06 08:29 +0100
Re: Column slices Bartc <bc@freeuk.com> - 2015-05-06 10:30 +0100
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-05-06 03:11 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-08 09:05 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-05-11 16:50 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-11 17:24 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-11 19:01 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-05-19 09:56 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-19 10:43 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-19 11:16 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-19 11:27 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-05-19 21:28 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-20 08:05 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-05-19 21:20 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-20 08:55 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-20 08:58 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-21 01:22 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-21 01:34 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-26 19:17 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-31 13:56 -0700
Re: Column slices Keith Thompson <kst-u@mib.org> - 2015-05-31 14:23 -0700
Re: Column slices Ian Collins <ian-news@hotmail.com> - 2015-06-01 09:37 +1200
Re: Column slices asetofsymbols@gmail.com - 2015-06-01 01:23 -0700
Re: Column slices Ian Collins <ian-news@hotmail.com> - 2015-06-01 20:35 +1200
Re: Column slices mark.bluemel@gmail.com - 2015-06-01 03:50 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-06-01 18:04 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-06-01 18:16 -0700
Re: Column slices "Chris M. Thomasson" <nospam@nospam.nospam> - 2015-06-01 19:57 -0700
Re: Column slices Rosario19 <Ros@invalid.invalid> - 2015-06-02 08:30 +0200
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-06-02 00:22 -0700
Re: Column slices Ike Naar <ike@iceland.freeshell.org> - 2015-05-12 06:06 +0000
Re: Column slices Ike Naar <ike@iceland.freeshell.org> - 2015-05-12 21:15 +0000
Re: Column slices asetofsymbols@gmail.com - 2015-05-13 20:49 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-05-13 23:43 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-05-13 23:44 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-14 01:02 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-15 23:06 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-16 03:46 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-05-23 00:17 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-23 07:59 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-23 21:25 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-05-24 18:29 -0700
Re: Column slices luser droog <luser.droog@gmail.com> - 2015-05-24 18:08 -0700
Re: Column slices asetofsymbols@gmail.com - 2015-05-25 00:13 -0700
Page 1 of 3 [1] 2 3 Next page →
| From | luser droog <luser.droog@gmail.com> |
|---|---|
| Date | 2015-05-05 23:18 -0700 |
| Subject | Column slices |
| Message-ID | <3f240d60-031e-4605-a07d-a9d20bed581b@googlegroups.com> |
I've been thinking further about the arbitrary-dimensional array
structure (and am still stunned by the lackluster response so far),
and I think a suitable coup-de-grace for demonstrating the value
of indirection is to tease-out yet another data member for this
structure. This ties everything back to the 1962 APL book, the
chapter on Representation, which has kept me up so many, many nights.
A dynamic array object is defined by its rank, list of dimension sizes,
a pointer to the start of the data, AND a weighting vector which for
simple uses merely caches the multiplications from the array-access
formula but it affords us much greater flexibility as (I hope) we
shall see...
typedef struct arr {
int rank;
int *dims;
int *weight;
int *data;
} *parr;
All of the pointer members can be assigned locations within the
same allocated block of memory as the struct itself (which we will
call the header). But by replacing the earlier use of offsets
and struct-hackery, the implementation of algorithms can be made
independent of that actual memory layout within (or without) the
block.
The basic memory layout for a self-contained array object is
rank dims weight data
dims[0] dims[1] ... dims[rank-1]
weight[0] weight[1] ... weight[rank-1]
data[0] data[1] ... data[ product(dims)-1 ]
An indirect array sharing data (whole array or 1 or more row-slices)
will have the following memory layout
rank dims weight data
dims[0] dims[1] ... dims[rank-1]
weight[0] weight[1] ... weight[rank-1]
//no data! it's somewhere else!
And an indirect array containing an orthogonal slice along
another axis will have the same layout as a basic indirect array,
but with dims and weight suitably modified.
The access formula for an element with indices (i0 i1 ... iN)
is
pa->data[ i0*pa->weight[0] + i1*pa->weight[1] + ...
+ iN*pa->weight[N] ]
, assuming each index i[j] is between [ 0 ... dims[j] ).
Nothing implemented yet.
[toc] | [next] | [standalone]
| From | asetofsymbols@gmail.com |
|---|---|
| Date | 2015-05-06 00:26 -0700 |
| Message-ID | <6dfd7905-9c07-4f0c-92fa-b3d988a55164@googlegroups.com> |
| In reply to | #61593 |
Is weight[i] == dim[i]*sizeof(matrix element)?
[toc] | [prev] | [next] | [standalone]
| From | luser droog <luser.droog@gmail.com> |
|---|---|
| Date | 2015-05-06 02:33 -0700 |
| Message-ID | <3e43ec2b-e875-4da5-812d-1d8db80552a9@googlegroups.com> |
| In reply to | #61595 |
On Wednesday, May 6, 2015 at 2:26:12 AM UTC-5, asetof...@gmail.com wrote: > Is weight[i] == dim[i]*sizeof(matrix element)? It's the product of all "lower" dimensions. weight[i] = product(dims[i+1 .. rank-1] So for dims = (3, 4, 5) weight = (20, 5, 1) = (4*5, 5, 1) or for dims = (7, 6, 5, 4, 3, 2) weight = (720, 120, 24, 6, 2, 1)
[toc] | [prev] | [next] | [standalone]
| From | Bartc <bc@freeuk.com> |
|---|---|
| Date | 2015-05-06 10:46 +0100 |
| Message-ID | <micnrf$mdk$1@dont-email.me> |
| In reply to | #61603 |
On 06/05/2015 10:33, luser droog wrote: > On Wednesday, May 6, 2015 at 2:26:12 AM UTC-5, asetof...@gmail.com wrote: >> Is weight[i] == dim[i]*sizeof(matrix element)? > > It's the product of all "lower" dimensions. > weight[i] = product(dims[i+1 .. rank-1] > > So for dims = (3, 4, 5) > weight = (20, 5, 1) = (4*5, 5, 1) > > or for dims = (7, 6, 5, 4, 3, 2) > weight = (720, 120, 24, 6, 2, 1) > On 06/05/2015 10:30, Bartc wrote: > So for a 3x4x5 array it would be: > Presumably the W vector might be W0=12 W1=4 W2=1, to help locate Or W0=20 as you demonstrate above. (The boxes in the diagram I created were too badly drawn to see properly what size they were!) -- Bartc
[toc] | [prev] | [next] | [standalone]
| From | asetofsymbols@gmail.com |
|---|---|
| Date | 2015-05-07 11:02 -0700 |
| Message-ID | <687c2ccd-29e0-46a8-bc79-2dc23faf82b9@googlegroups.com> |
| In reply to | #61603 |
Thank you, clear
[toc] | [prev] | [next] | [standalone]
| From | asetofsymbols@gmail.com |
|---|---|
| Date | 2015-05-07 14:05 -0700 |
| Message-ID | <9f17ce7b-ddc2-4eb1-8e96-9cd77074af84@googlegroups.com> |
| In reply to | #61702 |
On Thursday, May 7, 2015 at 8:02:35 PM UTC+2, asetof...@gmail.com wrote:
> Thank you, clear
Here my start in implementing matrices ...
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <iostream.h>
#define G(a,b) if(a)goto b
#define R return
#define P printf
#define M malloc
#define Fr free
#define F for
#define S sizeof
#define u8 unsigned char
#define u16 unsigned short
#define u32 unsigned
#define i32 int
#define d64 double
/*
I[0],I[1],I[2],I[3],I[4],sizeElement, sz, led, [matrix data]
0 1 2 3 4 5 6 7 8
1 2 3 4 5 6 7 8 9
I[I[0]+1]=I[3]
I[I[0]+I[0]]=I[4]
I[0],I[1],..,I[p[0]], sizeElement, sz, led, [matrix data]
0 1 I[0] 2*I[0]+1 2*I[0]+2 2*I[0]+3 2*I[0]+4 indirizzo
1 2 I[0]+1 2*I[0]+2 2*I[0]+3 2*I[0]+4 2*I[0]+5 nelement
I[0]: is the number of index as 32 bit
I[1]: is the first max index as 32 bit
....
I[I[0]]: is the last max index
I[I[0]+1]: is w[0]
...
I[2*I[0]]: is w[I[0]]
I[2*I[0]+1]: is the size of element as 32 bit 5
I[2*I[0]+2]: is sz size of "data" area 6
I[2*I[0]+3]: is led: 0 matrix is ok, 1 overflow or errors in matrix 7
sono 2*I[0]+4 elementi [8 elementi]
I[2*I[0]+4]: is matrix "data" area 8
*/
// access to led as u32*
#define Led(a) ((a)+*(a)*2+3)
// access to size of matrix data as u32*
#define SZ(a) ((a)+*(a)*2+2)
// access to size element of matrix as u32*
#define SZe(a) ((a)+*(a)*2+1)
// access to data as u8*
#define Data(a) (u8*) ((a)+*(a)*2+4)
/*
reserve mem space to matrix of element msize for dimension
i0xi1x...xin taken as input
where i0, i1 etc are max+1 index and n<=255
the last max+1 index is -1 that end the input
Return -1 for error
Return 0 all ok
NB
**The first function to call on pm [construct matrix] **
**the last function to call on pm [destruct matrix] **
*/
i32 Mcon(u32** pm, u32 msize, ...)
{va_list ap;
u32 v, i, s, j;
u32 arr[256];
G(pm, l1);
l0: R -1; // for error
n0: va_end(ap); *pm=0; G(1,l0);
l1: G(msize, l3);Fr(*pm); *pm=0;/* free the Matrix */
l2: R 0; // for all ok
y0: va_end(ap); G(1,l2);
l3: va_start(ap, msize); G(msize>0xFFF, n0); G(msize==0, n0); i=0;
l4: v=va_arg(ap, u32); //v the first index;if i==0 means no element
G(v>0xFFFF && v!= -1, n0); G(v!=-1, la); G(i==0, n0);
s=1; j=1;
l5: s*=arr[j-1]; G(s>0xFFFFF, n0); ++j; G( j<=i, l5);
s*=msize; G(s>0xFFFFF, n0);
*pm=(u32*)M((2*i+4)*sizeof(u32)+s); G(*pm==0, n0);
(*pm)[0]=i; *SZe(*pm)=msize; *SZ(*pm)=s; *Led(*pm)=0;
j=1; s=1;
l6: (*pm)[j]=arr[j-1]; (*pm)[2*i+1-j]=s;
s*=arr[i-j]; ++j; G(j<=i, l6);
G(1,y0);
la: G(i>255, n0); arr[i++]=v; G(1,l4);
}
// return one pointer to element
// or 0 if find some error
/*
2x2 y*n +x (y, x)e{0..m -1}x{0..n -1}
2 elm(i0,i1)=i0*M1+i1 (i0,i1)e{0..M0-1}x{0..M1-1} Max=(M0-1)*M1+M1-1=M0*M1-1
se M0=7 M1=1 7x1 (1,1)
elm(i0,i1)=i0+i1 (i0,i1)e{0..6}x{0..0} Max=(7-1)*1+1-1=6
se M0=1 M1=7 1x7 (7,1)
elm(i0,i1)=i0*7+i1 (i0,i1)e{0..0}x{0..6} Max=(0)*7+7-1=6
3 elm(i0,i1,i2)=i0*M1*M2+i1*M1+i2 (i0,i1,i2)e{0..M0-1}x{0..M1-1}x{0..M2-1}
Max=(M0-1)*M1*M2+(M1-1)*M2+(M2-1)=M0*M1*M2-M1*M2+M1*M2-M2+M2-1=M0*M1*M2-1
*/
u8* Aelm(u32* pm, ...)
{va_list ap;
u8 *r;
u32 maxn, size, v, s, i;
G(pm,l1);
l0: R 0;
n0: va_end(ap); G(1,l0);
l1: maxn=pm[0]; size= *SZe(pm); G(size>0xFFFF, l0);
// 0 1
G(maxn>=256||maxn==0,l0); //for 2x2 y*n+x (y,x)e[0..m-1]x[0..n-1]
va_start(ap, pm); // gets all index i0 i1 maxn=2
i=1; s=0; /* I[0]+1..2*I[0] */
l3: v=va_arg(ap,u32); G(v>=pm[i],n0); s+=v*pm[maxn+i]; ++i; G(i<=maxn, l3);
v=s*size; r=Data(pm); r+=v;
va_end(ap);
R r;
}
// for now only nxm matrix of u8 u16 u32 d64[double]
void PrntMtrx(u32* m)
{va_list ap;
d64 dv;
u8 *p, *pmax;
u32 i,j,maxn, max, size, v;
u32 Index[256];
G(m, l1);
l0: R; // error
la: P("[Me]"); G(1,l0);
l1: G(m[0]>2, l0);
size= *SZe(m);
G(size>8||size==0, l0);
G( *Led(m) , la);
i=0; p=Data(m); pmax=p+(size*m[1]*m[2]);
m0: G(size==1,l2);G(size==2,l3);G(size==4,l4);G(size==8,l11);G(1,l0);
l11:dv=*(d64*)p; P("%10.2f",dv); G(1,l5);
l2: v =*(u8*) p; P("%3u ", v); G(1,l5);
l3: v =*(u16*)p; P("%5u ", v); G(1,l5);
l4: v =*(u32*)p; P("%10u ", v);
l5: G( (i+1)%m[2]!=0, l6); P("\n");
l6: p+=size; ++i; G(p<pmax, m0);
R;
}
template<class T> class Matrix{
public:
Matrix(u32 a, u32 b){Mcon(&mtr, S(T), a, b, -1);}
Matrix(u32 a, u32 b, u32 c){Mcon(&mtr, S(T), a, b, c, -1);}
Matrix(u32 a, u32 b, u32 c, u32 d){Mcon(&mtr, S(T), a, b, c, d, -1);}
~Matrix(){Mcon(&mtr, 0);}
T& elm(u32 i, u32 j){R *(T*) Aelm(mtr, i,j);}
T& elm(u32 i, u32 j, u32 k){R *(T*) Aelm(mtr, i,j,k); }
T& elm(u32 i, u32 j, u32 k, u32 v){R *(T*) Aelm(mtr, i,j,k, v); }
Matrix<T>& operator=(Matrix<T>& a)
{u32 i;
if(a.mtr== 0)
{
l0: *Led(mtr)=1; G(1,z); /* errore trovato e riportato */
}
if(a.mtr!=mtr)
{// deve essere delle stesse dimensioni e indici massimi
G(a.mtr[0]!=mtr[0], l0);
F(i=1;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
G(*Led(a.mtr)!=0, l0); G(*SZe(mtr)!=*SZe(a.mtr), l0);
memcpy(mtr, a.mtr, *SZ(a.mtr)+(2*a.mtr[0]+4)*sizeof(u32));
}
z:
R *this;
}
Matrix<T>& operator+=(Matrix<T>& a)
{u32 i;
T *p1,*p2;
u8 *pmax;
if(a.mtr== 0)
{
l0: *Led(mtr)=1; G(1,z); /* errore trovato e riportato */
}
// deve essere delle stesse dimensioni e indici massimi
G(a.mtr[0]!=mtr[0], l0);
F(i=1;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
G(*Led(a.mtr)!=0, l0); G(*SZe(mtr)!=*SZe(a.mtr), l0);
p1=(T*) Data(mtr); p2=(T*) Data(a.mtr);
i=*SZ(mtr); pmax=(u8*) p1; pmax+=i;
F(; p1<pmax ;++p1,++p2) *p1=*p1+*p2;
z:
R *this;
}
Matrix<T>& operator-=(Matrix<T>& a)
{u32 i;
T *p1,*p2;
u8 *pmax;
if(a.mtr== 0)
{
l0: *Led(mtr)=1; G(1,z); /* errore trovato e riportato */
}
// deve essere delle stesse dimensioni e indici massimi
G(a.mtr[0]!=mtr[0], l0);
F(i=1;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
G(*Led(a.mtr)!=0, l0); G(*SZe(mtr)!=*SZe(a.mtr), l0);
p1=(T*) Data(mtr); p2=(T*) Data(a.mtr);
i=*SZ(mtr); pmax=(u8*) p1; pmax+=i;
F(; p1<pmax ;++p1,++p2) *p1=*p1-*p2;
z:
R *this;
}
Matrix<T>& operator*=(T a)
{u32 i;
T *p1;
u8 *pmax;
p1=(T*) Data(mtr);
i=*SZ(mtr); pmax=(u8*) p1; pmax+=i;
F(; p1<pmax ;++p1) *p1 = *p1 * a;
R *this;
}
/* ritorna il sup dell'indice consentito */
/* indici = 0, 1, 2 ,3 ,4 */
u32 operator[](u32 i){if(mtr==0||mtr[0]<=i) R 0; R mtr[i+1];}
int print(void){PrntMtrx(mtr);}
u32 *mtr;
};
int main(void)
{Matrix<u32> m1(2,2), m2(2,3), m7(2,2);
Matrix<d64> m3(2,2);
Matrix<u8> m4(1,7), m5(7,1);
u32 i,j;
F(i=0; i<m2[0]; ++i)
F(j=0; j<m2[1]; ++j)
m2.elm(i,j)=10*i+j;
F(i=0; i<m4[0]; ++i)
F(j=0; j<m4[1]; ++j)
m4.elm(i,j)=10*i+j;
F(i=0; i<m5[0]; ++i)
F(j=0; j<m5[1]; ++j)
m5.elm(i,j)=10*i+j;
m1.elm(0,0)=123; m1.elm(0,1)=124;
m1.elm(1,0)=125; m1.elm(1,1)=126;
F(i=0; i<m3[0]; ++i)
F(j=0; j<m3[1]; ++j)
m3.elm(i,j)=10.23*i+j;
P("M1:\n"); m1.print(); P("\n\n");
P("M2:\n"); m2.print(); P("\n\n");
P("M3:\n"); m3.print(); P("\n\n");
P("M4:\n"); m4.print(); P("\n\n");
P("M5:\n"); m5.print(); P("\n\n");
m7=m1;
P("M7:\n"); m7.print(); P("\n\n");
m1*=10;
P("M1=M1*10:\n"); m1.print(); P("\n\n");
m1+=m1;
P("M1+=M1:\n"); m1.print(); P("\n\n");
R 0;
}
----------------------
M1:
123 124
125 126
M2:
0 1 2
10 11 12
M3:
0.00 1.00
10.23 11.23
M4:
0 1 2 3 4 5 6
M5:
0
10
20
30
40
50
60
M7:
123 124
125 126
M1=M1*10:
1230 1240
1250 1260
M1+=M1:
2460 2480
2500 2520
[toc] | [prev] | [next] | [standalone]
| From | Nobody <nobody@nowhere.invalid> |
|---|---|
| Date | 2015-05-06 08:29 +0100 |
| Message-ID | <pan.2015.05.06.07.28.58.999000@nowhere.invalid> |
| In reply to | #61593 |
On Tue, 05 May 2015 23:18:12 -0700, luser droog wrote: > A dynamic array object is defined by its rank, list of dimension sizes, a > pointer to the start of the data, AND a weighting vector which for simple > uses merely caches the multiplications from the array-access formula but > it affords us much greater flexibility as (I hope) we shall see... This is essentially NumPy's ndarray (N-dimensional array) class: http://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html That may be worth looking at if you're going to reinvent that particular wheel (the concept certainly pre-dates NumPy, but it's well documented, maintained, and doesn't require a fancy keyboard).
[toc] | [prev] | [next] | [standalone]
| From | Bartc <bc@freeuk.com> |
|---|---|
| Date | 2015-05-06 10:30 +0100 |
| Message-ID | <micmuc$jde$1@dont-email.me> |
| In reply to | #61593 |
On 06/05/2015 07:18, luser droog wrote:
> A dynamic array object is defined by its rank, list of dimension sizes,
> a pointer to the start of the data, AND a weighting vector which for
> simple uses merely caches the multiplications from the array-access
> formula but it affords us much greater flexibility as (I hope) we
> shall see...
>
> typedef struct arr {
> int rank;
> int *dims;
> int *weight;
> int *data;
> } *parr;
>
> All of the pointer members can be assigned locations within the
> same allocated block of memory as the struct itself (which we will
> call the header). But by replacing the earlier use of offsets
> and struct-hackery, the implementation of algorithms can be made
> independent of that actual memory layout within (or without) the
> block.
>
> The basic memory layout for a self-contained array object is
>
> rank dims weight data
> dims[0] dims[1] ... dims[rank-1]
> weight[0] weight[1] ... weight[rank-1]
> data[0] data[1] ... data[ product(dims)-1 ]
So for a 3x4x5 array it would be:
3 4 5
W0 W1 W2
D0 D1 ..... D59 ?
If the last index (the 0 to 4 one) varies more quickly, the data
consists of 3x4 rows of 5, or 3 blocks of (4 rows of 5). And the data
will be more like:
D(0,0,0) D(0,0,1),... D(0,0,4) D(0,1,0) .... D(2,3,3) D(2,3,4)
Presumably the W vector might be W0=12 W1=4 W2=1, to help locate
D[i][j][k], with the last W not really needed.
But where does the column slice come into it? To access any slice other
than a row, would need the ability to deal with an 'array' that has a
stride, or offset between elements, the stride being W0, W1 (or W2 for a
row).
>
> An indirect array sharing data (whole array or 1 or more row-slices)
> will have the following memory layout
>
> rank dims weight data
> dims[0] dims[1] ... dims[rank-1]
> weight[0] weight[1] ... weight[rank-1]
> //no data! it's somewhere else!
>
> And an indirect array containing an orthogonal slice along
> another axis will have the same layout as a basic indirect array,
> but with dims and weight suitably modified.
Don't forget 2D, 3D and N-dimensional slices, where the slice (for a 3D
array) might be rectangular 2D array (in, I think, 3 possible
orientations), or a cuboid, which is a 3D subarray (and might also come
in different orientations, although you can't tell apart from the
indexing order).
> Nothing implemented yet.
> typedef struct arr {
> int rank;
> int *dims;
> int *weight;
> int *data;
> } *parr;
Does it need a different struct depending on the array element type?
This also seems inefficient for small dimensions, for 2D for example
which perhaps is as high as is usually needed.
And for 1D, although I guess a different method is used for those,
unless it is possible for the number of dimensions to be variable? Ie.
be itself a (1D) vector, requiring a similarly sized vector of indices.
(I suppose an N-dimensional array of dimensions is a possibility but I
can't get my head around that!)
(I guess these arrays form part of a different language (APL?), and C is
used to implement it. For such purposes, I don't bother with
N-dimensional arrays at all, just have arrays of arrays to any degree.
It's probably not quite as efficient, but it's simpler. And there are no
multiplications involved.)
--
Bartc
[toc] | [prev] | [next] | [standalone]
| From | luser droog <luser.droog@gmail.com> |
|---|---|
| Date | 2015-05-06 03:11 -0700 |
| Message-ID | <39c32351-d2b8-4f2d-adc0-03fd98b01426@googlegroups.com> |
| In reply to | #61602 |
On Wednesday, May 6, 2015 at 4:31:02 AM UTC-5, Bart wrote:
> On 06/05/2015 07:18, luser droog wrote:
>
> > A dynamic array object is defined by its rank, list of dimension sizes,
> > a pointer to the start of the data, AND a weighting vector which for
> > simple uses merely caches the multiplications from the array-access
> > formula but it affords us much greater flexibility as (I hope) we
> > shall see...
> >
> > typedef struct arr {
> > int rank;
> > int *dims;
> > int *weight;
> > int *data;
> > } *parr;
> >
> > All of the pointer members can be assigned locations within the
> > same allocated block of memory as the struct itself (which we will
> > call the header). But by replacing the earlier use of offsets
> > and struct-hackery, the implementation of algorithms can be made
> > independent of that actual memory layout within (or without) the
> > block.
> >
> > The basic memory layout for a self-contained array object is
> >
> > rank dims weight data
> > dims[0] dims[1] ... dims[rank-1]
> > weight[0] weight[1] ... weight[rank-1]
> > data[0] data[1] ... data[ product(dims)-1 ]
>
>
> So for a 3x4x5 array it would be:
>
> 3 4 5
> W0 W1 W2
> D0 D1 ..... D59 ?
>
> If the last index (the 0 to 4 one) varies more quickly, the data
> consists of 3x4 rows of 5, or 3 blocks of (4 rows of 5). And the data
> will be more like:
>
> D(0,0,0) D(0,0,1),... D(0,0,4) D(0,1,0) .... D(2,3,3) D(2,3,4)
>
> Presumably the W vector might be W0=12 W1=4 W2=1, to help locate
> D[i][j][k], with the last W not really needed.
>
> But where does the column slice come into it? To access any slice other
> than a row, would need the ability to deal with an 'array' that has a
> stride, or offset between elements, the stride being W0, W1 (or W2 for a
> row).
I get confused when I try to do it directly, but you can get column
slices by taking a transpose and then taking row slices.
So if you have (d,w) = ((3,4,5),(20,5,1), then (d,w)' = ((5,3,4),(1,20,5))
and (d,w)'[0;*;*] = ((1,3,4)(1,20,5)). The transpose of this resulting
array is the [*;*;0] column slice of the original (d,w).
The weighting vector AFAIK is the sine qua non of taking column slices.
And yes, the last weight will usually be 1, unless you're working with
a transpose.
> >
> > An indirect array sharing data (whole array or 1 or more row-slices)
> > will have the following memory layout
> >
> > rank dims weight data
> > dims[0] dims[1] ... dims[rank-1]
> > weight[0] weight[1] ... weight[rank-1]
> > //no data! it's somewhere else!
> >
> > And an indirect array containing an orthogonal slice along
> > another axis will have the same layout as a basic indirect array,
> > but with dims and weight suitably modified.
>
> Don't forget 2D, 3D and N-dimensional slices, where the slice (for a 3D
> array) might be rectangular 2D array (in, I think, 3 possible
> orientations), or a cuboid, which is a 3D subarray (and might also come
> in different orientations, although you can't tell apart from the
> indexing order).
I think all those are possible, so long as the data can be accessed "linearly"
in the algebraic sense.
> > Nothing implemented yet.
>
> > typedef struct arr {
> > int rank;
> > int *dims;
> > int *weight;
> > int *data;
> > } *parr;
>
> Does it need a different struct depending on the array element type?
> This also seems inefficient for small dimensions, for 2D for example
> which perhaps is as high as is usually needed.
That gave me a lot of grief with inca2 which has 3 basic datatypes:
char, int, and double. It was indeed a royal pain to deal with the
sizes all over the place. For the new rewrite, I had a different
idea to encode everything into int-size handles. This is described
in this post in comp.lang.apl:
https://groups.google.com/d/topic/comp.lang.apl/KSwo40-eLKI/discussion
> And for 1D, although I guess a different method is used for those,
> unless it is possible for the number of dimensions to be variable? Ie.
> be itself a (1D) vector, requiring a similarly sized vector of indices.
> (I suppose an N-dimensional array of dimensions is a possibility but I
> can't get my head around that!)
Yes, the number of dimensions is variable and is stored in the .rank
struct member. The dimension list is 1D, a vector.
> (I guess these arrays form part of a different language (APL?), and C is
> used to implement it. For such purposes, I don't bother with
> N-dimensional arrays at all, just have arrays of arrays to any degree.
> It's probably not quite as efficient, but it's simpler. And there are no
> multiplications involved.)
>
Yes, this all comes directly from my experience trying to implement a
(better) APL interpreter in C. But I hope the ideas have more general
applicability. And perhaps I can build a library or at least a demo
of some linear algebra something or other that may help somebody
somewhere someday.
[toc] | [prev] | [next] | [standalone]
| From | asetofsymbols@gmail.com |
|---|---|
| Date | 2015-05-08 09:05 -0700 |
| Message-ID | <5f4c0f5f-1c50-4ab1-8f34-09e45d20796e@googlegroups.com> |
| In reply to | #61605 |
"The dimension list is 1D, a vector" What about 1D not exist in matrix? I in my code not allow that it does not exist. Exist nx1 As column matrix for example for 3x1 1 2 3 Or 1xn matrix example 1x3 1 2 3 For both program use 2 index not only one so is what you call 'rank' 2? I forget all math on matrix and follow the computer...
[toc] | [prev] | [next] | [standalone]
| From | luser droog <luser.droog@gmail.com> |
|---|---|
| Date | 2015-05-11 16:50 -0700 |
| Message-ID | <8a54a6ec-0827-42e9-8f9b-b5c4d0084e8e@googlegroups.com> |
| In reply to | #61593 |
On Wednesday, May 6, 2015 at 1:18:25 AM UTC-5, luser droog wrote:
> I've been thinking further about the arbitrary-dimensional array
> structure (and am still stunned by the lackluster response so far),
> and I think a suitable coup-de-grace for demonstrating the value
> of indirection is to tease-out yet another data member for this
> structure. This ties everything back to the 1962 APL book, the
> chapter on Representation, which has kept me up so many, many nights.
>
> A dynamic array object is defined by its rank, list of dimension sizes,
> a pointer to the start of the data, AND a weighting vector which for
> simple uses merely caches the multiplications from the array-access
> formula but it affords us much greater flexibility as (I hope) we
> shall see...
>
> typedef struct arr {
> int rank;
> int *dims;
> int *weight;
> int *data;
> } *parr;
>
> All of the pointer members can be assigned locations within the
> same allocated block of memory as the struct itself (which we will
> call the header). But by replacing the earlier use of offsets
> and struct-hackery, the implementation of algorithms can be made
> independent of that actual memory layout within (or without) the
> block.
>
> The basic memory layout for a self-contained array object is
>
> rank dims weight data
> dims[0] dims[1] ... dims[rank-1]
> weight[0] weight[1] ... weight[rank-1]
> data[0] data[1] ... data[ product(dims)-1 ]
>
> An indirect array sharing data (whole array or 1 or more row-slices)
> will have the following memory layout
>
> rank dims weight data
> dims[0] dims[1] ... dims[rank-1]
> weight[0] weight[1] ... weight[rank-1]
> //no data! it's somewhere else!
>
> And an indirect array containing an orthogonal slice along
> another axis will have the same layout as a basic indirect array,
> but with dims and weight suitably modified.
>
> The access formula for an element with indices (i0 i1 ... iN)
> is
>
> pa->data[ i0*pa->weight[0] + i1*pa->weight[1] + ...
> + iN*pa->weight[N] ]
>
> , assuming each index i[j] is between [ 0 ... dims[j] ).
>
> Nothing implemented yet.
A basic implementation. Insufficient error checking.
//arrind.c - Arrays, indirect.
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "ppnarg.h"
typedef struct arr {
int rank;
int *dims;
int *weight;
int *data;
} *arr;
/* multiply together rank integers in dims array
*/
int productdims(int rank, int *dims){
int i,z=1;
for(i=0; i<rank; i++)
z *= dims[i];
return z;
}
/* create a new array with specified dimensions */
#define array(...) (makearr)(PP_NARG(__VA_ARGS__),__VA_ARGS__)
/* create a new array with specified rank and dimensions
*/
arr (makearr)(int rank, ...){
va_list ap;
//int *dims=calloc(rank,sizeof(int));
int dims[rank];
int datasz;
int i;
int x;
arr z;
va_start(ap,rank);
for (i=0; i<rank; i++){
dims[i]=va_arg(ap,int);
}
va_end(ap);
datasz=productdims(rank,dims);
z=malloc(sizeof(struct arr)
+ (rank+rank+datasz)*sizeof(int));
z->rank = rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + rank;
z->data = z->weight + rank;
memmove(z->dims,dims,rank*sizeof(int));
for(x=1, i=rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}
//free(dims);
return z;
}
/* create a new array which shares the data of an existing array
*/
arr clone(arr a){
arr z=malloc(sizeof(struct arr)
+ (a->rank+a->rank)*sizeof(int));
z->rank = a->rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + z->rank;
z->data = a->data;
memmove(z->dims,a->dims,z->rank*sizeof(int));
memmove(z->weight,a->weight,z->rank*sizeof(int));
return z;
}
/* exchange the leftmost two dimensions (only two in 2D)
*/
void transpose2(arr a){
int t;
//if (a->rank != 2) error();
t = a->dims[0]; a->dims[0] = a->dims[1]; a->dims[1] = t;
t = a->weight[0]; a->weight[0] = a->weight[1]; a->weight[1] = t;
}
/* take a (row) slice (in 2D) */
arr slice(arr a, int i){
int rank = a->rank-1;
arr z=malloc(sizeof(struct arr)
+ (rank+rank)*sizeof(int));
z->rank = rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + z->rank;
memmove(z->dims,a->dims+1,z->rank*sizeof(int));
memmove(z->weight,a->weight+1,z->rank*sizeof(int));
z->data = a->data + i*a->weight[0];
return z;
}
/* access element of a indexed by int[] */
int *elema(arr a, int *ind){
int idx = 0;
int i;
for (i=0; i<a->rank; i++){
idx += ind[i] * a->weight[i];
}
return a->data + idx;
}
/* access element of a indexed by va_list */
int *elemv(arr a, va_list ap){
int idx = 0;
int i;
for(i=0; i<a->rank; i++){
int ind;
ind = va_arg(ap,int);
//if (ind > a->dims[i]) error();
idx += ind * a->weight[i];
}
return a->data + idx;
}
/* access element of a indexed by integer arguments */
int *elem(arr a, ...){
va_list ap;
int *z;
va_start(ap,a);
z = elemv(a,ap);
va_end(ap);
return z;
}
/* create a (contiguous) copy of a (not necessarily contiguous)
existing array
*/
arr copy(arr a){
int datasz = productdims(a->rank,a->dims);
arr z=malloc(sizeof(struct arr)
+ (a->rank+a->rank+datasz)*sizeof(int));
int i;
int x;
int ind[a->rank];
z->rank = a->rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + z->rank;
z->data = z->weight + z->rank;
memmove(z->dims,a->dims,z->rank*sizeof(int));
for(x=1, i=z->rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}
for (i=0;i<datasz;i++){
int j;
int idx = i;
for (j=z->rank-1; j>=0; j--) { //generate index
ind[j] = idx % z->dims[j];
idx /= z->dims[j];
}
z->data[i] = *elema(a,ind);
}
return z;
}
#define OPERATORS(_) \
/* f F id */ \
_('+',+,0) \
_('*',*,1) \
_('=',==,1) \
/**/
/* perform binary operation F upon corresponding elements of vectors X and Y */
#define binop(X,F,Y) (binop)(X,*#F,Y)
#define BINOP(f,F,id) case f: *elem(z,i) = *elem(x,i) F *elem(y,i); break;
arr (binop)(arr x, char f, arr y){
arr z=copy(x);
int n=x->dims[0];
int i;
for (i=0; i<n; i++){
switch(f){
OPERATORS(BINOP)
}
}
return z;
}
#undef BINOP
/* perform binary operation F upon adjacent elements of vector X, right to left,
reducing vector to a single value */
#define reduce(F,X) (reduce)(*#F,X)
#define REDID(f,F,id) case f: x = id; break;
#define REDOP(f,F,id) case f: x = *elem(a,i) F x; break;
int (reduce)(char f, arr a){
int n = a->dims[0];
int x;
int i;
switch(f){
OPERATORS(REDID)
}
if (n) {
x=*elem(a,n-1);
for (i=n-2;i>=0;i--){
switch(f){
OPERATORS(REDOP)
}
}
}
return x;
}
#undef REDID
#undef REDOP
/* perform a (2D) matrix multiplication upon x and y
using operations F and G.
Z[i,j] = F/ X[i,*] G Y'[j,*]
*/
#define matmul(X,F,G,Y) (matmul)(X,*#F,*#G,Y)
arr (matmul)(arr x, char f, char g, arr y){
int i,j;
arr z=array(x->dims[0],y->dims[1]);
y=clone(y);
transpose2(y);
for (i=0; i<x->dims[0]; i++){
for (j=0; j<y->dims[0]; j++){
arr xs = slice(x,i);
arr ys = slice(y,j);
*elem(z,i,j)=(reduce)(f,(binop)(xs,g,ys));
free(xs);
free(ys);
}
}
free(y);
return z;
}
int main(){
#if 0
/* testing basic functionality and copies, transposes, and slices */
int loop;
for (loop = 0;
loop < 1/*00000*/;
loop++)
{
{
int i,n=12;
arr a=array(n);
for (i=0;i<n;i++)
*elem(a,i) = n-i;
for (i=0;i<n;i++,printf(" "))
printf("%2d",*elem(a,i));
printf("\n\n");
free(a);
}
{
int i,j,n=6;
arr a=array(n,n);
arr b;
for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(a,i,j) = n*n - (i*n+j);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%2d",*elem(a,i,j));
printf("\n");
b=clone(a);
transpose2(b);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%2d",*elem(b,i,j));
printf("\n");
free(b);
free(a);
}
{
int i,j,k,n=3;
arr a=array(n,n,n);
arr b;
arr c;
arr d;
for (i=0;i<n;i++)
for (j=0;j<n;j++)
for (k=0;k<n;k++)
*elem(a,i,j,k) = n*n*n - ((i*n+j)*n+k);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d",*elem(a,i,j,k));
printf("\n");
/*
b=slice(a,0);
transpose2(b);
for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(b,i,j) = n*n - (i*n+j);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d",*elem(a,i,j,k));
printf("\n");
free(b);
*/
b=clone(a);
transpose2(b);
c=slice(b,1);
transpose2(c);
for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(c,i,j) = n*n - (i*n+j);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d",*elem(b,i,j,k));
printf("\n");
//d = copy(a);
d = copy(b);
for (i=0; i<n*n*n; i++){
printf("%2d", d->data[i]);
printf(" ");
if (((i+1)%n)==0) printf("\n");
if (((i+1)%(n*n))==0) printf("\n");
}
free(a);
free(b);
free(c);
free(d);
}
}
#endif
{ /* testing reduce() */
int i,n=3;
arr a=array(n);
arr b;
for (i=0; i<n; i++)
*elem(a,i) = i+1;
printf("%2d\n", reduce(*,a));
free(a);
n=6;
a=array(n);
for (i=0; i<n; i++)
*elem(a,i) = 5;
b=binop(a,=,a);
for (i=0; i<n; i++,printf(" "))
printf("%2d", *elem(b,i));
printf("\n%2d\n", reduce(=,b));
free(a);
free(b);
}
{ /* testing matmul() */
int i,j,n=3;
arr a=array(n,n);
arr b;
for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(a,i,j) = ((i*n)+j)+1;
b = matmul(a,+,*,a);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%2d",*elem(a,i,j));
printf("\n");
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%2d",*elem(b,i,j));
printf("\n");
free(a);
free(b);
}
return 0;
}
[toc] | [prev] | [next] | [standalone]
| From | asetofsymbols@gmail.com |
|---|---|
| Date | 2015-05-11 17:24 -0700 |
| Message-ID | <f68a20e5-279f-4cc4-b6b7-883af097eca2@googlegroups.com> |
| In reply to | #61867 |
luser droog wrote:"
perform a (2D) matrix multiplication upon x and y
using operations F and G.
Z[i,j] = F/ X[i,*] G Y'[j,*]
*/
#define matmul(X,F,G,Y) (matmul)(X,*#F,*#G,Y)
arr (matmul)(arr x, char f, char g, arr y){
int i,j;
arr z=array(x->dims[0],y->dims[1]);
y=clone(y);
transpose2(y);
for (i=0; i<x->dims[0]; i++){
for (j=0; j<y->dims[0]; j++){
arr xs = slice(x,i);
arr ys = slice(y,j);
*elem(z,i,j)=(reduce)(f,(binop)(xs,g,ys));
free(xs);
free(ys);
}
}
free(y);
return z;
}"
A mxn nxk multiplication,
can be defined without the use
of new matrix to free o to traspose
etc
In C++ Something as
i32 mul(matrix<T>& a, matrix<T>& b)
3 matrices a, b and *this as result
Gaming with their index
the loop result something as
for(i=0; i<i0; ++i)
{for(j=0; j<i1; ++j)
for(k=0,s=0u; k<?; ++k)
{s=s+a.e(i,k)*b.e(k,j);}
*this.e(i,j)=s;
}
[toc] | [prev] | [next] | [standalone]
| From | asetofsymbols@gmail.com |
|---|---|
| Date | 2015-05-11 19:01 -0700 |
| Message-ID | <75e5349d-50c3-4aa3-b262-a5b7a95470dd@googlegroups.com> |
| In reply to | #61870 |
I wrote:"
A mxn nxk multiplication,
can be defined without the use
of new matrix to free o to traspose
etc
In C++ Something as
i32 mul(matrix<T>& a, matrix<T>& b)
3 matrices a, b and *this as result
Gaming with their index
the loop result something as
for(i=0; i<i0; ++i)
{for(j=0; j<i1; ++j)
for(k=0,s=0u; k<?; ++k)
{s=s+a.e(i,k)*b.e(k,j);}
*this.e(i,j)=s;
}"
for(i=0; i<i0; ++i)
{for(j=0; j<i1; ++j)
{for(k=0,s=0u; k<?; ++k)
{s=s+a.e(i,k)*b.e(k,j);}
*this.e(i,j)=s;
}
}
[toc] | [prev] | [next] | [standalone]
| From | luser droog <luser.droog@gmail.com> |
|---|---|
| Date | 2015-05-19 09:56 -0700 |
| Message-ID | <d42b2604-878d-437a-bf4e-d32c66bd88d8@googlegroups.com> |
| In reply to | #61874 |
On Monday, May 11, 2015 at 9:01:47 PM UTC-5, asetof...@gmail.com wrote:
> I wrote:"
> A mxn nxk multiplication,
> can be defined without the use
> of new matrix to free o to traspose
> etc
> In C++ Something as
> i32 mul(matrix<T>& a, matrix<T>& b)
> 3 matrices a, b and *this as result
> Gaming with their index
> the loop result something as
> for(i=0; i<i0; ++i)
> {for(j=0; j<i1; ++j)
> {for(k=0,s=0u; k<?; ++k)
> {s=s+a.e(i,k)*b.e(k,j);}
> *this.e(i,j)=s;
> }
> }
Took a little more time and energy than anticipated,
but finally here's the cake. matmul() now implements
the generalized Inner Product for arguments of
arbitrary, compatible dimensions.
Getting a little long for this format. Probably the
last post I'll make here. This file is also available
at https://github.com/luser-dr00g/inca/blob/master/arrind.c
//arrind.c - Arrays, indirect.
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "ppnarg.h" //https://github.com/luser-dr00g/inca/ppnarg.h
typedef struct arr {
int rank;
int *dims;
int *weight;
int *data;
} *arr;
/* multiply together rank integers in dims array
*/
int productdims(int rank, int *dims){
int i,z=1;
for(i=0; i<rank; i++)
z *= dims[i];
return z;
}
arr arraya(int rank, int *dims){
int datasz;
int i;
int x;
arr z;
datasz=productdims(rank,dims);
z=malloc(sizeof(struct arr)
+ (rank+rank+datasz)*sizeof(int));
z->rank = rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + rank;
z->data = z->weight + rank;
memmove(z->dims,dims,rank*sizeof(int));
for(x=1, i=rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}
return z;
}
void loaddimsv(int rank, int dims[], va_list ap){
int i;
for (i=0; i<rank; i++){
dims[i]=va_arg(ap,int);
}
}
/* create a new array with specified dimensions */
#define array(...) (array)(PP_NARG(__VA_ARGS__),__VA_ARGS__)
/* create a new array with specified rank and dimensions
*/
arr (array)(int rank, ...){
va_list ap;
//int *dims=calloc(rank,sizeof(int));
int dims[rank];
int i;
int x;
arr z;
va_start(ap,rank);
loaddimsv(rank,dims,ap);
va_end(ap);
z = arraya(rank,dims);
//free(dims);
return z;
}
/* create an array header to access existing data in multidimensional layout */
arr cast(int *data, int rank, ...){
va_list ap;
int dims[rank];
int i;
int x;
arr z;
va_start(ap,rank);
loaddimsv(rank,dims,ap);
va_end(ap);
z=malloc(sizeof(struct arr)
+ (rank+rank)*sizeof(int));
z->rank = rank;
z->dims = (int *)(((char *)z) + sizeof(struct arr));
z->weight = z->dims + rank;
z->data = data;
memmove(z->dims,dims,rank*sizeof(int));
for(x=1, i=rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}
return z;
}
/* create a new array which shares the data of an existing array
*/
arr clone(arr a){
arr z=malloc(sizeof(struct arr)
+ (a->rank+a->rank)*sizeof(int));
z->rank = a->rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + z->rank;
z->data = a->data;
memmove(z->dims,a->dims,z->rank*sizeof(int));
memmove(z->weight,a->weight,z->rank*sizeof(int));
return z;
}
/* exchange the leftmost two dimensions (only two in 2D)
*/
void transpose2(arr a){
int t;
//if (a->rank != 2) error();
t = a->dims[0]; a->dims[0] = a->dims[1]; a->dims[1] = t;
t = a->weight[0]; a->weight[0] = a->weight[1]; a->weight[1] = t;
}
/* rotate dims and weights according to sign and magnitude of shift
transpose(1,a)==transpose(-1,a)==transpose2(a) for 2D
*/
void transpose(int shift, arr a){
int i;
int t;
while(shift){
if (shift>0){
t=a->dims[0];
for (i=1; i<a->rank; i++)
a->dims[i-1]=a->dims[i];
a->dims[a->rank-1]=t;
t=a->weight[0];
for (i=1; i<a->rank; i++)
a->weight[i-1]=a->weight[i];
a->weight[a->rank-1]=t;
--shift;
} else {
t=a->dims[a->rank-1];
for (i=a->rank-2; i>=0; i--)
a->dims[i+1]=a->dims[i];
a->dims[0]=t;
t=a->weight[a->rank-1];
for (i=a->rank-2; i>=0; i--)
a->weight[i+1]=a->weight[i];
a->weight[0]=t;
++shift;
}
}
}
/* take a (row) slice (in 2D) */
arr slice(arr a, int i){
int rank = a->rank-1;
arr z=malloc(sizeof(struct arr)
+ (rank+rank)*sizeof(int));
z->rank = rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + z->rank;
memmove(z->dims,a->dims+1,z->rank*sizeof(int));
memmove(z->weight,a->weight+1,z->rank*sizeof(int));
z->data = a->data + i*a->weight[0];
return z;
}
/* access element of a indexed by int[] */
int *elema(arr a, int *ind){
int idx = 0;
int i;
for (i=0; i<a->rank; i++){
idx += ind[i] * a->weight[i];
}
return a->data + idx;
}
/* access element of a indexed by va_list */
int *elemv(arr a, va_list ap){
int idx = 0;
int i;
for(i=0; i<a->rank; i++){
int ind;
ind = va_arg(ap,int);
//if (ind > a->dims[i]) error();
idx += ind * a->weight[i];
}
return a->data + idx;
}
/* access element of a indexed by integer arguments */
int *elem(arr a, ...){
va_list ap;
int *z;
va_start(ap,a);
z = elemv(a,ap);
va_end(ap);
return z;
}
/* compute vector index list for ravel index ind */
int *vector_index(int ind, int *dims, int n, int *vec){
int i,t=ind, *z=vec;
for (i=0; i<n; i++){
z[n-1-i] = t % dims[n-1-i];
t /= dims[n-1-i];
}
return z;
}
/* compute ravel index for vector index list */
int ravel_index(int *vec, int *dims, int n){
int i,z=*vec;
for (i=0; i<n-1;i++){
z*=dims[i+1];
z+=vec[i+1];
}
return z;
}
/* create a vector of all elements of x followed by all elements of y */
arr catv(arr x, arr y){
int xsz = productdims(x->rank,x->dims);
int ysz = productdims(y->rank,y->dims);
int datasz = xsz + ysz;
arr z=array(datasz);
int scratch[x->rank+y->rank];
int i;
for (i=0; i<xsz; i++)
*elem(z,i) = *elema(x,vector_index(i,x->dims,x->rank,scratch));
for (i=0; i<ysz; i++)
*elem(z,xsz+i) = *elema(y,vector_index(i,y->dims,y->rank,scratch));
return z;
}
/* create a (contiguous) copy of a (not necessarily contiguous)
existing array
*/
arr copy(arr a){
int datasz = productdims(a->rank,a->dims);
arr z=malloc(sizeof(struct arr)
+ (a->rank+a->rank+datasz)*sizeof(int));
int i;
int x;
int ind[a->rank];
z->rank = a->rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + z->rank;
z->data = z->weight + z->rank;
memmove(z->dims,a->dims,z->rank*sizeof(int));
for(x=1, i=z->rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}
for (i=0;i<datasz;i++){
int j;
#if 0
int idx = i;
for (j=z->rank-1; j>=0; j--) { //generate index
ind[j] = idx % z->dims[j];
idx /= z->dims[j];
}
#endif
vector_index(i,z->dims,z->rank,ind);
z->data[i] = *elema(a,ind);
}
return z;
}
#define OPERATORS(_) \
/* f F id */ \
_('+',+,0) \
_('*',*,1) \
_('=',==,1) \
/**/
/* perform binary operation F upon corresponding elements of vectors X and Y */
#define binop(X,F,Y) (binop)(X,*#F,Y)
#define BINOP(f,F,id) case f: *elem(z,i) = *elem(x,i) F *elem(y,i); break;
arr (binop)(arr x, char f, arr y){
arr z=copy(x);
int n=x->dims[0];
int i;
for (i=0; i<n; i++){
switch(f){
OPERATORS(BINOP)
}
}
return z;
}
#undef BINOP
/* perform binary operation F upon adjacent elements of vector X, right to left,
reducing vector to a single value */
#define reduce(F,X) (reduce)(*#F,X)
#define REDID(f,F,id) case f: x = id; break;
#define REDOP(f,F,id) case f: x = *elem(a,i) F x; break;
int (reduce)(char f, arr a){
int n = a->dims[0];
int x;
int i;
switch(f){
OPERATORS(REDID)
}
if (n) {
x=*elem(a,n-1);
for (i=n-2;i>=0;i--){
switch(f){
OPERATORS(REDOP)
}
}
}
return x;
}
#undef REDID
#undef REDOP
/* perform a (2D) matrix multiplication upon rows of x and columns of y
using operations F and G.
Z = X F.G Y
Z[i,j] = F/ X[i,*] G Y'[j,*]
more generally,
perform an inner product on arguments of compatible dimension.
Z = X[A;B;C;D;E;F] +.* Y[G;H;I;J;K] |(F = G)
Z[A;B;C;D;E;H;I;J;K] = +/ X[A;B;C;D;E;*] * Y[*;H;I;J;K]
*/
#define matmul(X,F,G,Y) (matmul)(X,*#F,*#G,Y)
arr (matmul)(arr x, char f, char g, arr y){
int i,j;
//arr z=array(x->dims[0],y->dims[y->rank-1]);
arr xdims = cast(x->dims,1,x->rank);
arr ydims = cast(y->dims,1,y->rank);
xdims->dims[0]--;
ydims->dims[0]--;
ydims->data++;
arr z=arraya(x->rank+y->rank-2,catv(xdims,ydims)->data);
int datasz = productdims(z->rank,z->dims);
int k=y->dims[0];
arr xs = array(k);
arr ys = array(k);
y=clone(y);
transpose(1,y);
for (i=0; i<datasz; i++){
int idx[z->rank];
vector_index(i,z->dims,z->rank,idx);
int *xdex=idx;
int *ydex=idx+x->rank;
memmove(ydex,ydex-1,y->rank);
for (j=0; j<k; j++){
xdex[x->rank-1]=j;
xs->data[j] = *elema(x,xdex);
}
for (j=0; j<k; j++){
ydex[y->rank-1]=j;
ys->data[j] = *elema(y,ydex);
}
z->data[i] = (reduce)(f,(binop)(xs,g,ys));
}
free(y);
return z;
}
arr iota(int n){
arr z = array(n);
int i;
for (i=0;i<n;i++)
*elem(z,i)=i;
return z;
}
int main(){
#ifdef TEST_BASIC
/* testing basic functionality and copies, transposes, and slices */
int loop;
for (loop = 0;
loop < 1/*00000*/;
loop++)
{
{
int i,n=12;
arr a=array(n);
for (i=0;i<n;i++)
*elem(a,i) = n-i;
for (i=0;i<n;i++,printf(" "))
printf("%2d",*elem(a,i));
printf("\n\n");
free(a);
}
{
int i,j,n=6;
arr a=array(n,n);
arr b;
for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(a,i,j) = n*n - (i*n+j);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%2d",*elem(a,i,j));
printf("\n");
b=clone(a);
transpose2(b);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%2d",*elem(b,i,j));
printf("\n");
free(b);
free(a);
}
{
int i,j,k,n=3;
arr a=array(n,n,n);
arr b;
arr c;
arr d;
for (i=0;i<n;i++)
for (j=0;j<n;j++)
for (k=0;k<n;k++)
*elem(a,i,j,k) = n*n*n - ((i*n+j)*n+k);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d",*elem(a,i,j,k));
printf("\n");
/*
b=slice(a,0);
transpose2(b);
for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(b,i,j) = n*n - (i*n+j);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d",*elem(a,i,j,k));
printf("\n");
free(b);
*/
b=clone(a);
transpose2(b);
c=slice(b,1);
transpose2(c);
for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(c,i,j) = n*n - (i*n+j);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d",*elem(b,i,j,k));
printf("\n");
//d = copy(a);
d = copy(b);
for (i=0; i<n*n*n; i++){
printf("%2d", d->data[i]);
printf(" ");
if (((i+1)%n)==0) printf("\n");
if (((i+1)%(n*n))==0) printf("\n");
}
free(a);
free(b);
free(c);
free(d);
}
}
#endif
#ifdef TEST_MATMUL
{ /* testing reduce() */
int i,n=3;
arr a=array(n);
arr b;
for (i=0; i<n; i++)
*elem(a,i) = i+1;
printf("%2d\n", reduce(*,a));
free(a);
n=6;
a=array(n);
for (i=0; i<n; i++)
*elem(a,i) = 5;
b=binop(a,=,a);
for (i=0; i<n; i++,printf(" "))
printf("%2d", *elem(b,i));
printf("\n%2d\n", reduce(=,b));
free(a);
free(b);
}
{ /* testing matmul() */
int i,j,n=3;
arr a=array(n,n);
arr b;
for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(a,i,j) = ((i*n)+j)+1;
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%3d",*elem(a,i,j));
printf("\n");
b = matmul(a,+,*,a);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%3d",*elem(b,i,j));
printf("\n");
free(a);
free(b);
}
{
int i,j,k,l,n=2;
arr a=iota(n*n*n);
arr b=cast(a->data,3,n,n,n);
arr c=matmul(b,+,*,b);
printf("%d\n",c->rank);
for (i=0; i<c->rank; i++,printf(" "))
printf("%d",c->dims[i]);
printf("\n");
for (i=0; i<c->dims[0]; i++,printf("\n"))
for (j=0; j<c->dims[1]; j++,printf("\n"))
for (k=0; k<c->dims[2]; k++,printf("\n"))
for (l=0; l<c->dims[3]; l++,printf(" "))
printf("%3d", *elem(c,i,j,k,l));
printf("\n");
free(a);
free(b);
free(c);
}
#endif
#ifdef TEST_CAST
{ /* testing cast() */
int i,j,k,n=3;
int q[n][n][n];
arr a=cast((int *)q,3,n,n,n);
for (i=0;i<n;i++)
for (j=0;j<n;j++)
for (k=0;k<n;k++)
*elem(a,i,j,k) = n*n*n - ((i*n+j)*n+k);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d",*elem(a,i,j,k));
printf("\n");
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d", q[i][j][k]);
printf("\n");
free(a);
}
#endif
#ifdef TEST_IOTA
{
arr a=iota(12);
arr b=catv(a,a);
int i;
for (i=0;i<b->dims[0];i++,printf(" "))
printf("%2d", *elem(b,i));
arr c=cast(b->data,3,2,3,4);
int j,k;
for (i=0; i<2; i++,printf("\n"))
for (j=0; j<3; j++,printf("\n"))
for (k=0; k<4; k++,printf(" "))
printf("%2d", *elem(c,i,j,k));
free(a);
free(b);
free(c);
}
#endif
return 0;
}
--
eof
[toc] | [prev] | [next] | [standalone]
| From | asetofsymbols@gmail.com |
|---|---|
| Date | 2015-05-19 10:43 -0700 |
| Message-ID | <03478eef-3187-4d9e-94ba-0c9d6a9cfef5@googlegroups.com> |
| In reply to | #62303 |
here the traspose of a=0 1 2 3 4 5 6 would be b= 0 1 2 3 4 5 6 for a 'rank'=2 i0=1 i1=7 w0=7 w1=1 for b 'rank'=2 i0=7 i1=1 w0=1 w1=1 So b here can not be obtained from a just swap i0<->i1 and w0<->w1
[toc] | [prev] | [next] | [standalone]
| From | asetofsymbols@gmail.com |
|---|---|
| Date | 2015-05-19 11:16 -0700 |
| Message-ID | <6e7a9fc1-50ed-484c-9242-511cbe98ed5b@googlegroups.com> |
| In reply to | #62311 |
Here this 1 2 3 4 5 6 7 8 has rank=2 i0=2 i1=4 w0=4 w1=1 5 1 6 2 7 3 8 4 has rank=2 i0=4 i1=2 w0=2 w1=1 And would be one the transpose of the other... I think for traspose swap i0 i1 and from them recalculate w0 and w1. This would be not sufficient because for me one has to swap elements in array too...
[toc] | [prev] | [next] | [standalone]
| From | asetofsymbols@gmail.com |
|---|---|
| Date | 2015-05-19 11:27 -0700 |
| Message-ID | <298bb251-4448-4cfd-893c-b747cb3af400@googlegroups.com> |
| In reply to | #62313 |
Here this 1 2 3 4 5 6 7 8 has rank=2 i0=2 i1=4 w0=4 w1=1 1 5 2 6 3 7 4 8 has rank=2 i0=4 i1=2 w0=2 w1=1 And would be one the transpose of the other... I think for traspose swap i0 i1 and from them recalculate w0 and w1. This would be not sufficient because for me one has to swap elements in array too... The linear elements of 1 are 1 2 3 4 5 6 7 8 The linear elements of 2 are 1 5 2 6 3 7 4 8 So in this case for to have the right matrix I would change linear array (or data) for the transpose matrix too not only i0 etc
[toc] | [prev] | [next] | [standalone]
| From | luser droog <luser.droog@gmail.com> |
|---|---|
| Date | 2015-05-19 21:28 -0700 |
| Message-ID | <02fba50f-73ed-4ad1-af91-39369115b455@googlegroups.com> |
| In reply to | #62313 |
On Tuesday, May 19, 2015 at 1:16:30 PM UTC-5, asetof...@gmail.com wrote: > Here this > 1 2 3 4 > 5 6 7 8 > has rank=2 i0=2 i1=4 w0=4 w1=1 > > 5 1 > 6 2 > 7 3 > 8 4 > has rank=2 i0=4 i1=2 w0=2 w1=1 > And would be one the transpose > of the other... Here you've done a different kind of transpose, about a different axis. Or you could consider it a transpose + reversal. A "simple" transpose would yield 1 5 2 6 3 7 4 8 > I think for traspose swap i0 i1 > and from them recalculate w0 and w1. > This would be not sufficient because > for me one has to swap elements in > array too... I think this kind of transpose can also be accommodated by using a negative weight, but I haven't fully thought it through. You might also need to adjust the pointer to account for the "new" first element.
[toc] | [prev] | [next] | [standalone]
| From | asetofsymbols@gmail.com |
|---|---|
| Date | 2015-05-20 08:05 -0700 |
| Message-ID | <ddaf96db-4dcf-485c-88cd-4b1623b51d87@googlegroups.com> |
| In reply to | #62325 |
Transpose here is ok... It seems code write itself... 800 lines in C++ like Matrix on numeric type operator defined =,+,-,*, determinant, transpose print base type build with overflow that matrix class can be used detecting overflow To do: Inverse, <<, >>, == 'complementi algebrici'
[toc] | [prev] | [next] | [standalone]
| From | luser droog <luser.droog@gmail.com> |
|---|---|
| Date | 2015-05-19 21:20 -0700 |
| Message-ID | <775dbce5-2d73-42d9-ad27-a72eac224834@googlegroups.com> |
| In reply to | #62311 |
On Tuesday, May 19, 2015 at 12:43:19 PM UTC-5, asetof...@gmail.com wrote:
> here the traspose of
> a=0 1 2 3 4 5 6
> would be b=
> 0
> 1
> 2
> 3
> 4
> 5
> 6
> for a 'rank'=2 i0=1 i1=7 w0=7 w1=1
> for b 'rank'=2 i0=7 i1=1 w0=1 w1=1
>
> So b here can not be obtained from
> a just swap i0<->i1 and w0<->w1
I think it still works for this case. Remember that
this is not a normal 'transpose' like everybody else
does. We don't rearrange the data at all. This is more
properly called a 'dope-vector pseudo-transpose'.
Instead of changing the data, we just change the
constants in the access formula, rearranging the
coefficients of the polynomial. In a real sense, this
is just an application of the commutativity and
associativity of addition.
So for concreteness, assume the data is arranged
sequentially starting at hypothetical address 500.
500: 0
501: 1
502: 2
503: 3
504: 4
505: 5
506: 6
if a is rank 2, dims {1, 7), weight (7, 1), then the
sum of the indices multiplied by the associated weights
added to the initial pointer (500) yield the appropriate
addresses for each element
a[0][0] == *(500+0*7+0*1)
a[0][1] == *(500+0*7+1*1)
a[0][2] == *(500+0*7+2*1)
a[0][3] == *(500+0*7+3*1)
a[0][4] == *(500+0*7+4*1)
a[0][5] == *(500+0*7+5*1)
a[0][6] == *(500+0*7+6*1)
So the dope-vector pseudo-transpose rearranges the
weights and dimensions to match the new ordering of
indices, but the sum remains the same. The initial
pointer remains the same. The data does not move.
b[0][0] == *(500+0*1+0*7)
b[1][0] == *(500+1*1+0*7)
b[2][0] == *(500+2*1+0*7)
b[3][0] == *(500+3*1+0*7)
b[4][0] == *(500+4*1+0*7)
b[5][0] == *(500+5*1+0*7)
b[6][0] == *(500+6*1+0*7)
[toc] | [prev] | [next] | [standalone]
Page 1 of 3 [1] 2 3 Next page →
Back to top | Article view | comp.lang.c
csiph-web