使用涂料向量访问多维数组的任意轴向切片? [英] Use a dope vector to access arbitrary axial slices of a multidimensional array?

问题描述

我正在构建一套函数来处理 多维数组数据结构，我希望能够定义数组的任意切片，这样我就可以实现两个任意矩阵的广义内积(也称为张量或nd数组).

我读过的一篇 APL 论文(老实说，我找不到哪篇——我读过很多)定义了左矩阵 X 上的矩阵乘积，维度为 A;B;C;D;E;F 和右矩阵 Y 维度为 G;H;I;J;K 其中 F==G 为

Z <- X +.× YZ[A;B;C;D;E;H;I;J;K] <- +/X[A;B;C;D;E;*] × Y[*;H;I;J;K]

其中+/ 是和，× 将逐个元素应用于两个相同长度的向量.

所以我需要左侧的行"切片和右侧的列"切片.我当然可以先转置，然后再用行"切片来模拟列"切片，但我宁愿做得更优雅.

维基百科关于切片的文章引出了关于涂料载体 这似乎是我正在寻找的灵丹妙药，但没有太多可做的.

如何使用 dope vector 实现任意切片?

(很久以后我才注意到数组的步幅，其中有一些细节.)

解决方案

定义

<块引用>

一般的数组切片可以通过一个 dope 向量或描述符引用每个数组来实现(无论是否内置在语言中)——一个包含第一个数组元素地址的记录，然后是每个索引的范围和索引公式中的相应系数.这种技术还允许立即数组转置、索引反转、子采样等.对于像 C 这样的语言，索引总是从零开始，具有 d 个索引的数组的 dope 向量至少有 1 + 2d 个参数.
http://en.wikipedia.org/wiki/Array_slicing#Details

这是一个密集的段落，但实际上全部都在那里.所以我们需要这样的数据结构:

struct {类型 *数据；//第一个数组元素的地址国际排名；//维数int *昏暗；//每个维度的大小int *重量；//索引公式中对应的系数};

其中TYPE 是任何元素类型，矩阵的字段.为了简单和具体，我们将只使用 int.出于我自己的目的，我将各种类型的编码设计为 整数句柄 所以 int 为 me 做这项工作，YMMV.

typedef struct arr {国际排名；int *昏暗；int *重量；int *数据；} * arr;

所有的指针成员都可以在与结构本身相同的分配内存块(我们将调用标题).但是通过替换之前使用的偏移量和struct-hackery，可以进行算法的实现独立于内部(或外部)的实际内存布局堵塞.

自包含数组对象的基本内存布局是

rank 使权重数据变暗暗淡[0] 暗淡[1] ...暗淡[rank-1]weight[0] weight[1] ... weight[rank-1]数据[0] 数据[1] ... 数据[ 产品(dims)-1 ]

共享数据的间接数组(整个数组或 1 个或多个行切片)将有以下内存布局

rank 使权重数据变暗暗淡[0] 暗淡[1] ...暗淡[rank-1]weight[0] weight[1] ... weight[rank-1]//没有数据！它在别处！

以及一个包含正交切片的间接数组另一个轴将具有与基本间接数组相同的布局，但适当修改了暗淡和重量.

具有索引的元素的访问公式 (i0 i1 ... iN)是

a->data[ i0*a->weight[0] + i1*a->weight[1] + ...+ iN*a->weight[N]]

，假设每个索引 i[j] 在 [ 0 ... dims[j] 之间).

在正常布局的行主数组的weight向量中，每个元素都是所有较低维度的乘积.

for (int i=0; i

所以对于一个 3×4×5 的数组，元数据是

{ .rank=3, .dims=(int[]){3,4,5}, .weight=(int[]){4*5, 5, 1} }

或者对于 7×6×5×4×3×2 数组，元数据为

{ .rank=6, .dims={7,6,5,4,3,2}, .weight={720, 120, 24, 6, 2, 1} }

施工

因此，要创建其中之一，我们需要 上一个问题 从维度列表计算大小.

/* 将 dims 数组中的秩整数相乘 */int productdims(int rank, int *dims){int i,z=1;for(i=0; i

要分配，只需malloc足够的内存并将指针设置到适当的位置.

/* 创建数组给定 rank 和 int[] dims */arr arraya(int rank, int dims[]){整数数据；国际我;整数 x;arr z;datasz=productdims(rank,dims);z=malloc(sizeof(struct arr)+ (rank+rank+datasz)*sizeof(int));z->等级=等级；z->dims = z + 1;z->weight = z->dims + rank；z-> 数据 = z-> 权重 + 等级；memmove(z->dims,dims,rank*sizeof(int));for(x=1, i=rank-1; i>=0; i--){z->weight[i] = x;x *= z-> dims[i];}返回 z;}

并使用与上一个答案相同的技巧，我们可以制作一个可变参数接口以简化使用.

/* 将 va_list 中的秩整数加载到 int[] dims */void loaddimsv(int rank, int dims[], va_​​list ap){国际我;for (i=0; i

甚至通过使用 ppnargrank 参数一>.

#define array(...) (array)(PP_NARG(__VA_ARGS__),__VA_ARGS__)/* 创建一个指定维度的新数组 */

现在构建其中之一非常容易.

arr a = array(2,3,4);//创建一个动态的 [2][3][4] 数组

访问元素

通过调用elema 的函数来检索元素，该函数将每个索引乘以相应的权重，对它们求和，然后索引data 指针.我们返回一个指向元素的指针，以便调用者可以读取或修改它.

/* 访问由 int[] 索引的元素 */int *elema(arr a, int *ind){int idx = 0;国际我;对于 (i=0; i<a->rank; i++){idx += ind[i] * a->weight[i];}返回 a->data + idx;}

同样的可变参数技巧可以使界面更简单elem(a,i,j,k).

轴向切片

要进行切片，首先我们需要一种方法来指定要提取哪些维度以及要折叠哪些维度.如果我们只需要从一个维度中选择单个索引或所有元素，那么 slice 函数可以将 rank ints 作为参数并解释-1 表示选择整个维度或 0..dimsi-1 表示选择单个索引.

/* 获取以下 spec[] 指令的计算切片如果 spec[i] >= 0 并且 spec[i] rank，然后 spec[i] 选择来自维度 i 的索引.如果 spec[i] == -1，则 spec[i] 选择整个维度 i.*/arr slicea(arr a, int spec[]){int i,j;国际排名；for (i=0,rank=0; i<a->rank; i++)等级+=规格[i]==-1；int dims[等级];整数重量[等级];for (i=0,j=0; i=a->rank) 中断；dims[i] = a->dims[j];重量[i] = a->重量[j];}arr z = casta(a->data, rank, dims);memcpy(z->weight,weight,rank*sizeof(int));对于(j=0；j<a->rank；j++){如果(规格[j]！=-1)z->data += spec[j] * a->weight[j];}返回 z;}

参数数组中-1s对应的所有维度和权重都被收集起来，用于创建一个新的数组头.所有 >= 0 的参数乘以其相关的权重并添加到 data 指针，offset 指向正确元素的指针.

就数组访问公式而言，我们将其视为多项式.

offset = constant + sum_i=0,n( weight[i] * index[i] )

因此，对于我们从中选择单个元素的任何维度(+ 所有较低的维度)，我们分解出现在不变的索引并将其添加到公式中的常数项(在我们的 C 表示中是data 指针本身).

辅助函数casta 使用共享的data 创建新的数组头.slicea 当然会改变权重值，但是通过计算权重本身，casta 变成了一个更普遍可用的函数.它甚至可以用来构造一个动态数组结构，直接在常规的 C 风格多维数组上操作，从而强制转换.

/* 创建一个数组头来访问多维布局中的现有数据 */arr casta(int *data, int rank, int dims[]){int i,x;arr z=malloc(sizeof(struct arr)+ (rank+rank)*sizeof(int));z->等级=等级；z->dims = z + 1;z->weight = z->dims + rank；z->数据=数据；memmove(z->dims,dims,rank*sizeof(int));for(x=1, i=rank-1; i>=0; i--){z->weight[i] = x;x *= z-> dims[i];}返回 z;}

转置

涂料向量也可用于实现转置.维度(和索引)的顺序可以更改.

请记住，这不是像其他人一样的正常移调"做.我们根本不重新排列数据.这是更多正确地称为涂料向量伪转置".而不是改变数据，我们只是改变访问公式中的常量，重新排列多项式的系数.真正意义上的这只是交换性的一个应用，并且加法的结合性.

所以为了具体起见，假设数据是排列好的顺序从假设地址 500 开始.

500: 0501:1502:2503:3504:4505:5506:6

如果 a 是 2 级，使 {1, 7) 变暗，权重 (7, 1)，则指数之和乘以相关权重添加到初始指针 (500) 产生适当的每个元素的地址

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)

所以涂料向量伪转置重新排列重量和尺寸以匹配新的排序指数，但总和保持不变.最初的指针保持不变.数据不会移动.

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)

内积(又名矩阵乘法)

因此，通过使用通用切片或转置+行"-切片(更容易)，可以实现广义内积.

首先我们需要两个辅助函数，用于对两个向量应用二元运算产生向量结果，并使用二元运算减少向量产生标量结果.

如上一个问题 我们将传入运算符，因此相同的函数可以用于许多不同的运算符.对于这里的样式，我将运算符作为单个字符传递，因此已经存在从 C 运算符到我们函数的运算符的间接映射.这是一个x宏表.

#define OPERATORS(_) /* f F id */\_('+',+,0) \_('*',*,1) \_('=',==,1) /**/#define binop(X,F,Y) (binop)(X,*#F,Y)arr (binop)(arr x, char f, arr y);/* 对向量 X 和 Y 的对应元素执行二元运算 F */

表中的额外元素用于 reduce 函数，用于空向量参数的情况.在这种情况下，reduce 应该返回操作符的 identity 元素，0 代表 +，1 代表 *.

#define reduce(F,X) (reduce)(*#F,X)int (reduce)(char f, arr a);/* 对向量 X 的相邻元素进行二元运算 F，从右到左，将向量减少为单个值 */

所以 binop 在操作符字符上进行循环和切换.

/* 对向量 X 和 Y 的对应元素进行二元运算 F */#define BINOP(f,F,id) case f: *elem(z,i) = *elem(x,i) F *elem(y,i);休息;arr (binop)(arr x, char f, arr y){arr z=copy(x);int n=x->dims[0];国际我;对于 (i=0; i

如果有足够的元素，reduce 函数会向后循环，如果有，则将初始值设置为最后一个元素，并将初始值预设为运算符的标识元素.

/* 对向量 X 的相邻元素进行二元运算 F，从右到左，将向量减少为单个值 */#define REDID(f,F,id) case f: x = id;休息;#define REDOP(f,F,id) case f: x = *elem(a,i) F x;休息;int (reduce)(char f, arr a){int n = a->dims[0];整数 x;国际我;开关(f){运营商(REDID)}如果 (n) {x=*elem(a,n-1);对于 (i=n-2;i>=0;i--){开关(f){操作员(重做)}}}返回 x;}#undef REDID#undef REDOP

借助这些工具，可以以更高级别的方式实现内积.

/* 对 x 的行和 y 的列执行(2D)矩阵乘法使用操作 F 和 G.Z = X F.G YZ[i,j] = F/X[i,*] G Y'[j,*]更普遍，对兼容维度的参数执行内积.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]*/arr(matmul)(arr x，char f，char g，arr y){int i,j;arr xdims = cast(x->dims,1,x->rank);arr ydims = cast(y->dims,1,y->rank);xdims->dims[0]--;ydims->dims[0]--;ydims->数据++；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 = NULL;arr ys = NULL;for (i=0; irank+y->rank];vector_index(i,z->dims,z->rank,idx);int *xdex=idx;int *ydex=idx+x->rank-1;memmove(ydex+1,ydex,y->rank);xdex[x->rank-1] = -1;免费(xs)；免费(是)；xs = slicea(x,xdex);ys = slicea(y,ydex);z->data[i] = (reduce)(f,(binop)(xs,g,ys));}免费(xs)；免费(是)；免费(xdims)；免费(ydims)；返回 z;}

这个函数还使用了cast函数，它为casta提供了一个可变参数接口.

/* 创建一个数组头来访问多维布局中的现有数据 */arr cast(int *data, int rank, ...){va_list ap;int dims[等级];va_start(ap,rank);loaddimsv(rank,dims,ap);va_end(ap);返回卡斯塔(数据，排名，暗淡)；}

并且它还使用vector_index 将一维索引转换为索引的 nD 向量.

/* 计算 ravel index ind 的向量索引列表 */int *vector_index(int ind, int *dims, int n, int *vec){int i,t=ind, *z=vec;对于 (i=0; i

github 文件.额外的支持功能也在github文件中.

<小时>

这个问答对是在实施 inca 过程中出现的一系列相关问题的一部分a> 用 C 编写的 APL 语言解释器.其他:如何使用动态分配的任意维数组? 和 如何将 C 数学运算符 (+-*/%) 传递给函数 result=mathfunc(x,+,y);? .其中一些材料之前已发布到 comp.lang.c.comp.lang.apl 中的更多背景信息.

I'm building a suite of functions to work with a multidimensional-array data structure and I want to be able to define arbitrary slices of the arrays so I can implement a generalized inner product of two arbitrary matrices (aka Tensors or n-d arrays).

An APL paper I read (I honestly can't find which -- I've read so many) defines the matrix product on left-matrix X with dimensions A;B;C;D;E;F and right-matrix Y with dimensions G;H;I;J;K where F==G as

Z <- X +.× Y
Z[A;B;C;D;E;H;I;J;K] <- +/ X[A;B;C;D;E;*] × Y[*;H;I;J;K]

where +/ is sum of, and × applies element-by-element to two vectors of the same length.

So I need "row" slices of the left and "column" slices of the right. I can of course take a transpose and then a "row" slice to simulate a "column" slice, but I'd rather do it more elegantly.

Wikipedia's article on slicing leads to a stub about dope vectors which seem to be the miracle cure I'm looking for, but there's not a lot there to go on.

How do I use a dope vector to implement arbitrary slicing?

(Much later I did notice Stride of an array which has some details.)

解决方案

Definition

General array slicing can be implemented (whether or not built into the language) by referencing every array through a dope vector or descriptor — a record that contains the address of the first array element, and then the range of each index and the corresponding coefficient in the indexing formula. This technique also allows immediate array transposition, index reversal, subsampling, etc. For languages like C, where the indices always start at zero, the dope vector of an array with d indices has at least 1 + 2d parameters.
http://en.wikipedia.org/wiki/Array_slicing#Details

That's a dense paragraph, but it's actually all in there. So we need a data structure like this:

struct {
    TYPE *data;  //address of first array element
    int rank; //number of dimensions
    int *dims; //size of each dimension
    int *weight; //corresponding coefficient in the indexing formula
};

Where TYPE is whatever the element type is, the field of the matrices. For simplicity and concreteness, we'll just use int. For my own purposes, I've devised an encoding of various types into integer handles so int does the job for me, YMMV.

typedef struct arr { 
    int rank; 
    int *dims; 
    int *weight; 
    int *data; 
} *arr; 

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

a->data[ i0*a->weight[0] + i1*a->weight[1] + ... 
          + iN*a->weight[N] ] 

, assuming each index i[j] is between [ 0 ... dims[j] ).

In the weight vector for a normally laid-out row-major array, each element is the product of all lower dimensions.

for (int i=0; i<rank; i++)
    weight[i] = product(dims[i+1 .. rank-1]);

So for a 3×4×5 array, the metadata would be

{ .rank=3, .dims=(int[]){3,4,5}, .weight=(int[]){4*5, 5, 1} }

or for a 7×6×5×4×3×2 array, the metadata would be

{ .rank=6, .dims={7,6,5,4,3,2}, .weight={720, 120, 24, 6, 2, 1} }

Construction

So, to create one of these, we need the same helper function from the previous question to compute the size from a list of dimensions.

/* 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;
}

To allocate, simply malloc enough memory and set the pointers to the appropriate places.

/* create array given rank and int[] dims */
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 = z + 1;
    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;
}

And using the same trick from the previous answer, we can make a variable-argument interface to make usage easier.

/* load rank integers from va_list into int[] dims */
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 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;
}

And even automatically produce the rank argument by counting the other arguments using the awesome powers of ppnarg.

#define array(...) (array)(PP_NARG(__VA_ARGS__),__VA_ARGS__) /* create a new array with specified dimensions */

Now constructing one of these is very easy.

arr a = array(2,3,4);  // create a dynamic [2][3][4] array

Accessing elements

An element is retrieved with a function call to elema which multiplies each index by the corresponding weight, sums them, and indexes the data pointer. We return a pointer to the element so it can be read or modified by the caller.

/* 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;
}

The same varargs trick can make an easier interface elem(a,i,j,k).

Axial Slices

To take a slice, first we need a way of specifying which dimensions to extract and which to collapse. If we just need to select a single index or all elements from a dimension, then the slice function can take rank ints as arguments and interpret -1 as selecting the whole dimension or 0..dims_i-1 as selecting a single index.

/* take a computed slice of a following spec[] instructions
   if spec[i] >= 0 and spec[i] < a->rank, then spec[i] selects
      that index from dimension i.
   if spec[i] == -1, then spec[i] selects the entire dimension i.
 */
arr slicea(arr a, int spec[]){
    int i,j;
    int rank;
    for (i=0,rank=0; i<a->rank; i++)
        rank+=spec[i]==-1;
    int dims[rank];
    int weight[rank];
    for (i=0,j=0; i<rank; i++,j++){
        while (spec[j]!=-1) j++;
        if (j>=a->rank) break;
        dims[i] = a->dims[j];
        weight[i] = a->weight[j];
    }   
    arr z = casta(a->data, rank, dims);
    memcpy(z->weight,weight,rank*sizeof(int));
    for (j=0; j<a->rank; j++){
        if (spec[j]!=-1)
            z->data += spec[j] * a->weight[j];
    }   
    return z;
}

All the dimensions and weights corresponding to the -1s in the argument array are collected and used to create a new array header. All arguments >= 0 are multiplied by their associated weight and added to the data pointer, offsetting the pointer to the correct element.

In terms of the array access formula, we're treating it as a polynomial.

offset = constant + sum_i=0,n( weight[i] * index[i] )

So for any dimension from which we're selecting a single element (+ all lower dimensions), we factor-out the now-constant index and add it to the constant term in the formula (which in our C representation is the data pointer itself).

The helper function casta creates the new array header with shared data. slicea of course changes the weight values, but by calculating weights itself, casta becomes a more generally usable function. It can even be used to construct a dynamic array structure that operates directly on a regular C-style multidimensional array, thus casting.

/* create an array header to access existing data in multidimensional layout */
arr casta(int *data, int rank, int dims[]){
    int i,x;
    arr z=malloc(sizeof(struct arr)
            + (rank+rank)*sizeof(int));

    z->rank = rank;
    z->dims = z + 1;
    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;
}

Transposes

The dope vector can also be used to implement transposes. The order of the dimensions (and indices) can be changed.

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) 

Inner Product (aka Matrix Multiplication)

So, by using general slices or transpose+"row"-slices (which are easier), generalized inner product can be implemented.

First we need the two helper functions for applying a binary operation to two vectors producing a vector result, and reducing a vector with a binary operation producing a scalar result.

As in the previous question we'll pass in the operator, so the same function can be used with many different operators. For the style here, I'm passing operators as single characters, so there's already an indirect mapping from C operators to our function's operators. This is an x-macro table.

#define OPERATORS(_) 
    /* f  F id */ 
    _('+',+,0) 
    _('*',*,1) 
    _('=',==,1) 
    /**/

#define binop(X,F,Y) (binop)(X,*#F,Y)
arr (binop)(arr x, char f, arr y); /* perform binary operation F upon corresponding elements of vectors X and Y */

The extra element in the table is for the reduce function for the case of a null vector argument. In that case, reduce should return the operator's identity element, 0 for +, 1 for *.

#define reduce(F,X) (reduce)(*#F,X)
int (reduce)(char f, arr a); /* perform binary operation F upon adjacent elements of vector X, right to left,
                                   reducing vector to a single value */

So the binop does a loop and a switch on the operator character.

/* perform binary operation F upon corresponding elements of vectors X and 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

And the reduce function does a backwards loop if there are enough elements, having set the initial value to the last element if there was one, having preset the initial value to the operator's identity element.

/* perform binary operation F upon adjacent elements of vector X, right to left,
   reducing vector to a single value */
#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

And with these tools, inner product can be implemented in a higher-level manner.

/* 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]
 */
arr (matmul)(arr x, char f, char g, arr y){
    int i,j;
    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 = NULL;
    arr ys = NULL;

    for (i=0; i<datasz; i++){
        int idx[x->rank+y->rank];
        vector_index(i,z->dims,z->rank,idx);
        int *xdex=idx;
        int *ydex=idx+x->rank-1;
        memmove(ydex+1,ydex,y->rank);
        xdex[x->rank-1] = -1;
        free(xs);
        free(ys);
        xs = slicea(x,xdex);
        ys = slicea(y,ydex);
        z->data[i] = (reduce)(f,(binop)(xs,g,ys));
    }

    free(xs);
    free(ys);
    free(xdims);
    free(ydims);
    return z;
}

This function also uses the functions cast which presents a varargs interface to casta.

/* create an array header to access existing data in multidimensional layout */
arr cast(int *data, int rank, ...){
    va_list ap;
    int dims[rank];

    va_start(ap,rank);
    loaddimsv(rank,dims,ap);
    va_end(ap);

    return casta(data, rank, dims);
}

And it also uses vector_index to convert a 1D index into an nD vector of indices.

/* 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;
}

github file. Additional supporting functions are also in the github file.

---

This Q/A pair is part of a series of related issues which arose in implementing inca an interpreter for the APL language written in C. Others: How can I work with dynamically-allocated arbitrary-dimensional arrays? , and How to pass a C math operator (+-*/%) into a function result=mathfunc(x,+,y);? . Some of this material was previously posted to comp.lang.c. More background in comp.lang.apl.

这篇关于使用涂料向量访问多维数组的任意轴向切片?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！