使用内情向量访问多维数组的任意轴位片？ [英] Use a dope vector to access arbitrary axial slices of a multidimensional array?
问题描述
我建立了一套功能以<一个工作href=\"http://stackoverflow.com/questions/30023867/how-can-i-work-with-dynamically-allocated-arbitrary-dimensional-arrays/30023868#30023868\">multidimensional-array数据结构，我希望能够定义任意的片阵列的的这样我就可以实现任意两个矩阵（指广义的内积又名张量或第二阵列的）。
这是APL论文我看了（老实说，我无法找到它 - 我已经看了这么多）定义的矩阵乘积左矩阵 X 与尺寸 A; B，C，D组;ê;˚F和右矩阵是，尺寸˚F==摹为
z，其中， - X +×大
Z [A; B; C; d，E，H，I;焦耳; K]＆LT; - + / X [A，B，C，D组;ê; *]×Y [*; H，I;焦耳; K]
，其中 + / 是的总和的，和×适用元件逐元件相同长度的两个向量
所以我需要在左边的行片和右侧的专栏片。我当然可以采取转，然后一个行片，以模拟列片，但我宁愿做更优雅。
维基百科上关于文章切片导致有关的预测矢量这似乎是灵丹妙药我在寻找，但不是很多那里去。
我如何使用的内情向量执行任意切片？
（很久以后我也注意到跨越其中有一些细节的阵列。）
定义
包含第一个数组元素的地址的记录，然后各指标和范围 -通用阵列切片可以通过引用每个数组通过的内情向量或描述来实现（不管是否内置到语言）索引公式中相应的系数。这种技术还允许直接阵列换位，索引反转，频带分割等。对于像C，语言，其中指数总是开始于零，阵列具有d指数具有至少1 + 2d的参数的内情向量。结果
http://en.wikipedia.org/wiki/Array_slicing#Details这是一个密集的段落，但实际上一切都在那里。因此，我们需要这样的数据结构：
结构{
TYPE *数据; //一个数组元素的地址
INT排名; //维数
INT *黯淡; //每一维的尺寸
INT *重; //索引公式中对应的系数
};其中，
键入是任何元素类型的矩阵中，的字段的。为了简化和具体化，我们就用INT。对于我自己而言，我已经设计了多种类型的编码成的整数处理，以便INT做这项工作了的我的，因人而异。typedef结构改编{
INT排名;
INT *黯淡;
INT *重;
为int *的数据;
} * ARR;所有指针构件可以是内所分配的位置
内存与结构本身（同一分配的块，我们将
调用标头）。但通过更换先前使用偏移
和struct-两轮牛车，算法的实现可以进行
独立之内（或没有）的实际内存布局
块。对于一个自包含的数组对象的基本内存布局是
排名变暗重量数据
DIMS [0]变暗[1] ... DIMS [秩1]
重量[0]重量[1] ...重[秩1]
数据[0]数据[1] ...数据[产品（DIMS）-1]这是间接阵列共享数据（整个阵列或1以上的行分片）
将有以下的内存布局排名变暗重量数据
DIMS [0]变暗[1] ... DIMS [秩1]
重量[0]重量[1] ...重[秩1]
//没有数据！这是别的地方！和沿包含正交切片间接阵列
另一轴将具有布局作为基本间接阵列相同，
但变暗，重量适当修改。对于指数的元素取数公式（I0 I1 ... IN）
是A-＆GT;数据[I0 * A-＆GT;重[0] + I1 * A-＆GT;重[1] + ...
+ IN * A-＆GT;重[N]，假设每个指标i [J]之间[0 ...变暗[J]）。
在
重量的正常下岗了的行主的阵列向量，每个元素都是低维度的产品。的for（int i = 0; I＆LT;职;我++）
重量[I] =产物（变暗第[i + 1 ..秩1]）;因此，对于一个3×4×5阵列中，元数据将是
{.rank = 3，.dims =（INT []）{} 3,4,5，。重量=（INT []）{4 * 5,5，1 }}或一个7×6×5×4×3×2阵列，所述元数据将是
{.rank = 6，.dims = {7,6,5,4,3,2}。重量= {720，120，24，6，2，1 }}建筑
所以，要创建其中之一，我们需要从<一个同样的辅助函数href=\"http://stackoverflow.com/questions/30023867/how-can-i-work-with-dynamically-allocated-arbitrary-dimensional-arrays\">$p$pvious问题从尺寸列表计算的大小。
/ *相乘居变暗整数数组* /
INT productdims（INT排名，为int *变暗）{
INT I，Z = 1;
对于（i = 0; I＆LT;职;我++）
Z * =变暗[I]
返回Z者除外;
}要分配，只需
的malloc足够的内存和指针设置到适当的地方。/ *创建数组给定秩和int []变暗* /
ARR arraya（INT排名，诠释变暗[]）{
INT datasz;
INT I;
INT X;
ARR Z者除外;
datasz = productdims（秩，变暗）;
Z =的malloc（sizeof的（结构ARR）
+（+等级+级datasz）* sizeof的（INT））; Z-GT＆;秩秩=;
Z-GT＆; DIMS = Z + 1;
Z-GT＆;体重= Z-GT＆;变暗+等级;
Z-GT＆;数据= Z-GT＆;重量+等级;
memmove与（Z-GT＆;变暗，变暗，等级*的sizeof（INT））;
为（X = 1，I =秩1; I＆GT; = 0; I - ）{
Z-GT＆;重量[I] = X;
X * = Z-GT＆;变暗[I]
} 返回Z者除外;
}和使用从previous回答同样的伎俩，我们可以做一个可变参数的接口，让使用更方便。
/ *负荷排名整数从va_list的转换成int []变暗* /
无效loaddimsv（INT排名，诠释变暗[]，va_list的AP）{
INT I;
对于（i = 0; I＆LT;职;我++）{
DIMS [I] =在va_arg（AP，INT）;
}
}/ *创建一个新的阵列指定等级和尺寸* /
ARR（阵列）（INT级，...）{
va_list的AP;
//为int *黯淡=释放calloc（军衔，sizeof的（INT））;
INT变暗[排名]
INT I;
INT X;
ARR Z者除外; 的va_start（AP，等级）;
loaddimsv（等级，变暗，AP）;
va_end用来（AP）; Z = arraya（秩，变暗）;
//自由（DIMS）;
返回Z者除外;
}和甚至自动生成的级别的使用的 ppnarg
的#define阵列（...）（阵列）（PP_NARG（__ __ VA_ARGS）__ VA_ARGS__）/ *创建一个指定尺寸的新数组* /现在构建，其中之一是很容易的。
ARR A =阵列（2,3,4）; //创建动态[2] [3] [4]的数组访问元素
元素是用一个函数调用中检索到
elema它乘以相应的权重各指标，对其求和，并索引数据指针。我们返回一个指向元件，因此它可以被读取或由呼叫者改性/ *的为int索引访问元素[] * /
为int * elema（ARR一，为int * IND）{
INT IDX = 0;
INT I;
对于（i = 0; I＆LT; A-＆GT;排名;我++）{
IDX + = IND [I] * A-＆GT;重[I]
}
返回A-＆GT;数据+ IDX;
}同样的可变参数的技巧可以让一个更简单的接口
ELEM（A，I，J，K）。轴位片
要分得一杯羹，首先我们需要指定的方式尺寸以提取和崩溃。如果我们只需要选择一个单一的指标或维度的所有元素，那么
片函数可以接受的级别的INT取值作为参数，跨preT -1作为选择全尺寸或0 ..的级别的-1作为选择单个指标。/ *采取以下规格的计算切片[]指示
如果规范[Ⅰ]≥= 0和规范[1] - ; A-＆GT;排名，然后SPEC [I]选择
该指数从维我。
如果规范[I] == -1，则SPEC [I]选择整个维我。
* /
ARR slicea（ARR一，诠释规范[]）{
INT I，J;
INT排名;
对于（i = 0，排名= 0; I＆LT; A-＆GT;排名;我++）
排名+ =规范[I] == - 1;
INT变暗[排名]
诠释重量[排名]
对于（i = 0，J = 0; I＆LT;职;我++，J ++）{
而（SPEC [J] = - 1！）J ++;
如果（J＆GT; = A-＆GT;排名）打破;
DIMS [I] = A-＆GT;变暗[J]。
重量[I] = A-＆GT;重[J]。
}
ARR Z =卡斯塔（A-＆GT;的数据，排名，变暗）;
的memcpy（Z-GT＆;重，重量，等级*的sizeof（INT））;
为（J = 0; J＆LT; A-＆GT;排名; J ++）{
如果（规范[J]！= - 1）
Z-GT＆;数据+ =规范[J] * A-＆GT;重[J]。
}
返回Z者除外;
}所有的尺寸和对应于自变量阵列中的-1s权重被收集并用于产生一个新的数组首部。所有参数> = 0通过其关联权重相乘，并添加到
数据指针的偏移的指针，正确的元素。在数组访问公式而言，我们把它当作一个多项式。
偏移=常数+ sum_i = 0，N（重量[I] *指数[I]）因此，对于任何尺寸从中我们选择单个元件（+所有较低维），我们因子出现在常数索引和在式（它添加到常数项这在我们的C重新presentation是
数据指针本身）。辅助函数
卡斯塔创建了共享数据的新数组头。slicea当然，改变权重值，而是通过计算权重本身，卡斯塔成为一个更加普遍使用的功能。它甚至可以用来构建直接操作一个常规的C风格的多维阵列上，从而动态数组结构的铸造的/ *创建一个数组头多维布局*访问现有的数据/
ARR卡斯塔（INT *数据，诠释排名，诠释变暗[]）{
INT I，X;
ARR Z =的malloc（sizeof的（结构ARR）
+（排名+等级）* sizeof的（INT））; Z-GT＆;秩秩=;
Z-GT＆; DIMS = Z + 1;
Z-GT＆;体重= Z-GT＆;变暗+等级;
Z-GT＆;数据=数据;
memmove与（Z-GT＆;变暗，变暗，等级*的sizeof（INT））;
为（X = 1，I =秩1; I＆GT; = 0; I - ）{
Z-GT＆;重量[I] = X;
X * = Z-GT＆;变暗[I]
} 返回Z者除外;
}转置
的内情向量还可以用于实现转置。尺寸（和索引）的顺序是可以改变的。
请记住，这不是一个正常的转和其他人一样
确实。我们并不在所有重新排列的数据。这是更
正确地称为涂料载体伪转置。
相反，改变数据，我们只是改变
取数公式常量，重新排列
的多项式的系数。在真正的意义上，这
仅仅是可交换的应用
另外的关联性。因此，对于具体性，假定布置在数据
依次开始假想地址500500：0
501：1
502：2
503：3
504：4
505：5
506：6如果一个秩2，变暗{1，7），重量（7,1），则
指数乘以相应的权重的总和
加到初始指针（500），得到相应的
为每个元素的地址A [0] [0] == *（500 + 0 * 7 + 0 * 1）
一个[0] [1] == *（500 + 0 * 7 + 1 * 1）
一个[0] [2] == *（500 + 0 * 7 + 2 * 1）
一个[0] [3] == *（500 + 0 * 7 + 3 * 1）
一个[0] [4] == *（500 + 0 * 7 + 4 * 1）
一个[0] [5] == *（500 + 0 * 7 + 5 * 1）
一个[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）内积（又名矩阵乘法）
因此，通过利用一般片或转+行-slices（其是更容易），广义内积可以实施。
首先，我们需要用于施加一个二进制操作以两个向量产生的向量结果，并减少与二进制运算的矢量产生一个标量结果的两个辅助功能。
如在<一href=\"http://stackoverflow.com/questions/29960273/how-to-pass-a-c-math-operator-into-a-function-result-mathfuncx-y\">$p$pvious问题的，我们将通过在操作员，因此同样的功能可与许多不同的运营商使用。对于这里的风格，我路过运营商为单个字符，所以已经有从C运营商的我们的函数的的运营商的间接映射。这是一个的x宏表。
的#define运算符（_）\\
/ *器f F ID * / \\
_（+，+，0）\\
_（'*'，*，1）\\
_（=，==，1）\\
/ ** /＃定义binop（X，F，Y）（binop）（X，*＃F，Y）
ARR（binop）（ARR X，焦楼ARR Y）; / *根据相应的向量X和Y *的元素执行二进制运算F /在表中的额外元素是一个零矢量参数的情况下，
减少功能。在这种情况下，减少应为+，1返回运营商的标识元素的，0*。的#define减少（F，X）（减少）（*＃F，X）
INT（减少）（焦楼ARR一）; / *对向量X的相邻的元素进行二进制运算男，从右到左，
还原矢量为单个值* /所以
binop做一个循环，并在操作符的开关。/ *根据相应的向量X和Y *的元素执行二进制运算F /
＃定义BINOP（F，F，id）的情况下，F：* ELEM（Z，I）= * ELEM（X，I）F * ELEM（Y，I）;打破;
ARR（binop）（ARR X，焦楼ARR Y）{
ARR Z =副本（X）;
INT N = X-GT＆;变暗[0];
INT I;
对于（i = 0; I＆LT; N;我++）{
开关（F）{
运算符（BINOP）
}
}
返回Z者除外;
}
和#undef BINOP和的降低函数做一个向后循环，如果有足够的元素，具有初始值设置为最后一个元素，如果有一个，具有preset中的初始值，以操作者的身份元素
/ *在执行向量X的相邻的元素二元运算男，从右到左，
还原矢量为单个值* /
＃定义重做（F，F，id）的情况下，F：X = ID;打破;
＃定义REDOP（F，F，id）的情况下，F：X = * ELEM（A，I）的F X;打破;
INT（减少）（焦楼ARR一）{
INT N = A-＆GT;变暗[0];
INT X;
INT I;
开关（F）{
运算符（重做）
}
如果（n）的{
X = * ELEM（一个中，n-1）;
为（ⅰ=正2; I＆GT; = 0; I - ）{
开关（F）{
运算符（REDOP）
}
}
}
返回X;
}
和#undef重做
和#undef REDOP和使用这些工具，内积可在更高级别的方式实现的。
/ *对x的行和y列执行（2D）矩阵乘法
使用操作F和G
ž= X F.Gÿ
Z [I，J] = F / X [I，*; G Y'[J，*] 更普遍，
上执行兼容尺寸的参数的内积。
Z = X [A; B; C; d，E，F] + * Y [G，H，I;焦耳; K] |（F = G）
Z [A; B; C; d，E，H，I;焦耳; K] = + / X [A，B，C，D组;ê; *] * Y [*; H，I;焦耳; K]
* /
ARR（MATMUL）（ARR X，焦楼炭克，ARR Y）{
INT I，J;
ARR xdims = CAST（X-＆GT;变暗，1，X轴和GT;等级）;
ARR ydims = CAST（Y-GT;变暗，1，Y＆GT;等级）;
xdims-＆GT;变暗[0] - ;
ydims-＆GT;变暗[0] - ;
ydims-＆GT;数据++;
ARR Z = arraya（X-GT＆;排名+ Y-GT;秩2，CATV（xdims，ydims） - GT;数据）;
INT datasz = productdims（Z-GT＆;秩，Z-＆GT;变暗）;
INT K = Y-＆GT;变暗[0];
ARR XS = NULL;
ARR YS = NULL; 对于（i = 0; I＆LT; datasz;我++）{
INT IDX [X-＆GT;排名+ Y-GT;等级]
vector_index（I，Z-GT＆;变暗，Z-＆GT;排名，IDX）;
为int * xdex = IDX;
为int * ydex = IDX + X＆GT;秩1;
memmove与（ydex + 1，ydex，Y＆GT;等级）;
xdex [X-＆GT;秩1] = -1;
免费（XS）;
免费（YS）;
XS = slicea（X，xdex）;
YS = slicea（Y，ydex）;
Z-GT＆;数据[I] =（减少）（F，（binop）（XS，G，YS））;
} 免费（XS）;
免费（YS）;
免费（xdims）;
免费（ydims）;
返回Z者除外;
}此功能还使用了
投其中presents一个可变参数的接口卡斯塔的功能。/ *创建一个数组头多维布局*访问现有的数据/
ARR投（INT *数据，INT级，...）{
va_list的AP;
INT变暗[排名] 的va_start（AP，等级）;
loaddimsv（等级，变暗，AP）;
va_end用来（AP）; 返回卡斯塔（数据，排名，变暗）;
}和它也使用
vector_index来一维指数转换为指数的第二矢量。/ *计算矢量指数拉威尔指数IND *列表/
为int * vector_index（INT IND，为int *黯淡，整数N，为int * VEC）{
INT I，T = IND，* Z = VEC;
对于（i = 0; I＆LT; N;我++）{
Z [n-1个-1] = T％变暗[N-1-i的];
T / DIMS = [N-1-I]
}
返回Z者除外;
}GitHub的文件。其它配套功能也都在GitHub的文件。
此Q / A对是一系列相关问题的一部分，出现在实施印加对于用C编写的其他的APL语言的跨preTER：<一href=\"http://stackoverflow.com/questions/30023867/how-can-i-work-with-dynamically-allocated-arbitrary-dimensional-arrays\">How我可以动态分配的任意维数组工作？和<一个href=\"http://stackoverflow.com/questions/29960273/how-to-pass-a-c-math-operator-into-a-function-result-mathfuncx-y\">How通过一个C数学运算符（+ - * /％）到一个函数的结果= mathfunc（X +，Y）;？。一些这方面的材料是previously发布到 comp.lang。 ç。在 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
Xwith dimensionsA;B;C;D;E;Fand right-matrixYwith dimensionsG;H;I;J;KwhereF==GasZ <- 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#DetailsThat'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
TYPEis whatever the element type is, the field of the matrices. For simplicity and concreteness, we'll just useint. For my own purposes, I've devised an encoding of various types into integer handles sointdoes 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
weightvector 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
mallocenough 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
elemawhich multiplies each index by the corresponding weight, sums them, and indexes thedatapointer. 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
slicefunction can take rankints as arguments and interpret -1 as selecting the whole dimension or 0..rank-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
datapointer, 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
datapointer itself).The helper function
castacreates the new array header with shareddata.sliceaof course changes the weight values, but by calculating weights itself,castabecomes 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: 6if 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
reducefunction 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
binopdoes 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 BINOPAnd 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 REDOPAnd 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
castwhich presents a varargs interface tocasta./* 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_indexto 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屋！
