rt/snapraid
1
0
Fork 0
mirror of https://github.com/amadvance/snapraid.git synced 2026-10-26 19:01:57 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions
snapraid/raid/raid.c

496 lines
16 KiB
C
Raw Blame History

// SPDX-License-Identifier: GPL-2.0-or-later
// Copyright (C) 2013 Andrea Mazzoleni
#include "internal.h"
#include "gf.h"
/*
* Extended Cauchy Matrix Construction for MDS RAID in GF(2^8)
*
* The Cauchy RAID modes operate in GF(2^8) and support up to six parity
* levels. RAID_MODE_CAUCHY_RAID uses polynomial 0x11d, while
* RAID_MODE_CAUCHY_AES uses polynomial 0x11b.
*
* For RAID5 and RAID6, it follows the method described in H. Peter Anvin's
* paper "The mathematics of RAID-6" [1]. Please refer to this paper for a
* complete explanation.
*
* Triple parity was first evaluated using an extension of the same approach,
* with additional parity coefficients set as powers of 2^-1, defined as:
*
* P = sum(Di)
* Q = sum(2^i * Di)
* R = sum(2^-i * Di) for 0 <= i < N
*
* This approach works efficiently for triple parity because multiplications
* and divisions by 2 in GF(2^8) are very fast. It is similar to ZFS RAIDZ3,
* which uses powers of 4 instead of 2^-1.
*
* Unfortunately, this method does not extend beyond triple parity because
* no choice of power coefficients guarantees solvable equations for all
* combinations of missing disks. This is expected: a Vandermonde matrix,
* used in this method, does not ensure that all submatrices are nonsingular
* [2, Chap. 11, Problem 7], which is required for an MDS (Maximum Distance
* Separable) code [2, Chap. 11, Theorem 8].
*
* To overcome this limitation, a Cauchy matrix [3][4] is used for parity
* computation. A Cauchy matrix ensures that all square submatrices are
* nonsingular, so the linear system is always solvable for any combination
* of missing disks. This guarantees an MDS code.
*
* RAID_MODE_CAUCHY_RAID uses polynomial 0x11d and the existing primitive
* generator g=2 construction described below.
*
* RAID_MODE_CAUCHY_AES instead uses polynomial 0x11b and a G23 Extended
* Cauchy construction. In the AES field, element 2 has multiplicative order
* 51 while element 3 is primitive. Therefore H=<2> has five cosets in
* GF(256)*. The Q coefficients q[0..254] enumerate those cosets, starting
* with q[0]=1. Within a coset q advances by multiplication by 2. The four
* transitions between cosets advance by multiplication by 3, after the
* coefficients indexed 50, 101, 152, and 203.
*
* The AES data X set consumes inverses of the G23 sequence from the beginning:
*
* X[i] = inv(q[i])
*
* The Cauchy Y set starts with zero and consumes inverses backward from the
* end of the sequence:
*
* Y[0] = 0
* Y[j] = inv(q[255-j]) for j > 0
*
* Thus the Y values associated with R, S, T, and U are inv(q[254]),
* inv(q[253]), inv(q[252]), and inv(q[251]), respectively. For PARITY=n and
* DISK=257-n, increasing parity removes one X element from the tail and
* appends it to Y, leaving the existing rows unchanged on the remaining
* columns. X and Y are disjoint and together contain all 256 field elements, so
* C[j,i]=1/(X[i]+Y[j]) defines valid Cauchy rows. Adding P as an all-ones row
* forms an Extended Cauchy matrix. Q is the row for Y[0]=0 and therefore
* Q=sum(q[i]*Di). Rows R through U are scaled so their first coefficient is 1.
* This nonzero row scaling preserves the MDS property, and the resulting 6x251
* matrix is MDS.
*
* The numerical examples below refer to RAID_MODE_CAUCHY_RAID,
* using polynomial 0x11d and generator g = 2.
*
* How the RAID/g=2 matrix is obtained:
*
* The matrix is constructed to match Linux RAID coefficients for the
* first two rows, while ensuring that all square submatrices are nonsingular.
*
* 1. Start by forming a Cauchy matrix with elements defined as 1/(x_i + y_j),
* where all x_i and y_j are distinct elements (the textbook definition
* of a Cauchy matrix).
*
* 2. For the first row (j=0), set x_i = g^-i and y_0 = 0, resulting in:
*
* row j=0 -> 1/(x_i + y_0) = 1/(g^-i + 0) = g^i
*
* which reproduces the RAID-6 coefficients.
*
* 3. For subsequent rows (j>0), set y_j = g^j, yielding:
*
* rows j>0 -> 1/(x_i + y_j) = 1/(g^-i + g^j)
*
* Since g is primitive and has multiplicative order 255,
* g^-i != g^j for the supported indexes where i + j < 255.
*
* 4. Place a row filled with 1 at the top of the matrix, transforming it
* into an Extended Cauchy Matrix. This preserves the nonsingularity
* of all square submatrices.
*
* This first row reproduces the RAID-5 coefficients.
*
* 5. Adjust each row with a factor so that the first column contains all 1s.
* This transformation also preserves the nonsingularity property,
* maintaining the MDS code property.
*
* Concrete example in GF(256) for k=6, m=4:
*
* First, create a 3x6 Cauchy matrix using x_i = g^-i and y_0 = 0, y_j = g^j for j>0.
* For RAID_MODE_CAUCHY_RAID, where g=2, this gives:
*
* x = { 1, 142, 71, 173, 216, 108 }
* y = { 0, 2, 4 }
*
* The Cauchy matrix is:
*
* 1 2 4 8 16 32
* 244 83 78 183 118 47
* 167 39 213 59 153 82
*
* Divide row 2 by 244 and row 3 by 167. Then extend it with a row of ones
* on top, and it remains MDS. This forms the code for m=4, with RAID-6 as
* a subset:
*
* 1 1 1 1 1 1
* 1 2 4 8 16 32
* 1 245 210 196 154 113
* 1 187 166 215 199 7
*
* Extending the same computation to k=251 and m=6 gives:
*
* This results in the matrix A[row, col] defined as:
*
* 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01...
* 01 02 04 08 10 20 40 80 1d 3a 74 e8 cd 87 13 26 4c 98 2d 5a b4 75...
* 01 f5 d2 c4 9a 71 f1 7f fc 87 c1 c6 19 2f 40 55 3d ba 53 04 9c 61...
* 01 bb a6 d7 c7 07 ce 82 4a 2f a5 9b b6 60 f1 ad e7 f4 06 d2 df 2e...
* 01 97 7f 9c 7c 18 bd a2 58 1a da 74 70 a3 e5 47 29 07 f5 80 23 e9...
* 01 2b 3f cf 73 2c d6 ed cb 74 15 78 8a c1 17 c9 89 68 21 ab 76 3b...
*
* This matrix supports six levels of parity, one for each row, for up to 251
* data disks, one for each column. All the 377,342,351,231 square submatrices
* are nonsingular, verified also by brute-force testing.
*
* The RAID/g=2 construction can be extended beyond six parity levels by
* adding rows and removing one data column for each additional row, while
* preserving the field-size constraint i + j < 255. The current
* implementation remains limited to RAID_PARITY_MAX = 6.
* (See mktables.c for more details on how the matrix is generated.)
*
* The downside is that generic multiplications are required to compute the
* parity, not only fast multiplications by 2^n or 2^-n, which negatively
* affect performance. However, optimized parallel multiplications using SSSE3
* or AVX2 instructions [1][5] make this approach competitive with the triple
* parity computation using power coefficients.
*
* Another advantage of the RAID/g=2 Cauchy matrix is that the first two rows can
* replicate the RAID5 and RAID6 approach, resulting in a compatible extension.
* SSSE3 or AVX2 instructions are only needed for triple parity or beyond.
*
* Parity computation is as follows:
*
* P = sum(Di)
* Q = sum(A[1,i] * Di)
* R = sum(A[2,i] * Di)
* S = sum(A[3,i] * Di)
* T = sum(A[4,i] * Di)
* U = sum(A[5,i] * Di) for 0 <= i < N
*
* For RAID_MODE_CAUCHY_RAID, A[1,i]=2^i. For RAID_MODE_CAUCHY_AES,
* A[1,i]=q[i] from the G23 sequence.
*
* Recovery from six disk failures at indices x, y, z, h, v, w (0 <= x < y < z < h < v < w < N)
* involves computing parity of the remaining N-6 disks:
*
* Pa = sum(Di)
* Qa = sum(A[1,i] * Di)
* Ra = sum(A[2,i] * Di)
* Sa = sum(A[3,i] * Di)
* Ta = sum(A[4,i] * Di)
* Ua = sum(A[5,i] * Di) for i != x,y,z,h,v,w
*
* Defining:
*
* Pd = Pa + P
* Qd = Qa + Q
* Rd = Ra + R
* Sd = Sa + S
* Td = Ta + T
* Ud = Ua + U
*
* yields:
*
* Pd = Dx + Dy + Dz + Dh + Dv + Dw
* Qd = A[1,x] * Dx + A[1,y] * Dy + A[1,z] * Dz + A[1,h] * Dh + A[1,v] * Dv + A[1,w] * Dw
* Rd = A[2,x] * Dx + A[2,y] * Dy + A[2,z] * Dz + A[2,h] * Dh + A[2,v] * Dv + A[2,w] * Dw
* Sd = A[3,x] * Dx + A[3,y] * Dy + A[3,z] * Dz + A[3,h] * Dh + A[3,v] * Dv + A[3,w] * Dw
* Td = A[4,x] * Dx + A[4,y] * Dy + A[4,z] * Dz + A[4,h] * Dh + A[4,v] * Dv + A[4,w] * Dw
* Ud = A[5,x] * Dx + A[5,y] * Dy + A[5,z] * Dz + A[5,h] * Dh + A[5,v] * Dv + A[5,w] * Dw
*
* This linear system is always solvable since the coefficient matrix is
* nonsingular due to the properties of A[].
*
* Example performance on a Intel Core i7-10700 CPU @ 2.90GHz, stripe 256 KiB, 8 data disks:
*
* int8 int32 int64 sse2 sse2e ssse3 ssse3e avx2 avx2e
* gen1 24995 46562 62683 77069
* gen2 8243 15879 25788 27292 44907
* genz 5509 10588 14485 13920 25175
* gen3 1190 13201 14307 27268
* gen4 914 10061 11006 21328
* gen5 746 7970 8828 17040
* gen6 630 6627 7516 14584
*
* Values are in MiB/s of data processed by a single thread.
* Results can be reproduced using "raid/test/speedtest.c".
*
* Triple parity with power coefficients "1, 2, 2^-1" is slightly faster than
* the Cauchy matrix computation if SSSE3 or AVX2 are available:
*
* int8 int32 int64 sse2 ssse3 avx2
* genz 2337 2874 10920 18944
*
* In conclusion, the use of power coefficients, specifically powers
* of 1, 2, and 2^-1, is the best option to implement triple parity on CPUs
* without SSSE3 and AVX2.
* On modern CPUs with SSSE3 or AVX2, the Cauchy matrix is the best
* option because it provides a fast and general approach that works for any
* number of parities.
*
* References:
* [1] Anvin, "The mathematics of RAID-6", 2004
* [2] MacWilliams, Sloane, "The Theory of Error-Correcting Codes", 1977
* [3] Blomer, "An XOR-Based Erasure-Resilient Coding Scheme", 1995
* [4] Roth, "Introduction to Coding Theory", 2006
* [5] Plank, "Screaming Fast Galois Field Arithmetic Using Intel SIMD Instructions", 2013
*/
/**
* Buffer filled with 0 used in recovering.
*/
static void *raid_zero_block;
void raid_zero(void *zero)
{
raid_zero_block = zero;
}
/**
* Inverts the square matrix M of size nxn into V.
*
* This is not a general matrix inversion because we assume the matrix M
* to have all the square submatrix not singular.
* We use Gauss elimination to invert.
*
* @M Matrix to invert with @n rows and @n columns.
* @V Destination matrix where the result is put.
* @n Number of rows and columns of the matrix.
*/
void raid_invert(uint8_t *M, uint8_t *V, int n)
{
int i, j, k;
/* set the identity matrix in V */
for (i = 0; i < n; ++i)
for (j = 0; j < n; ++j)
V[i * n + j] = i == j;
/* for each element in the diagonal */
for (k = 0; k < n; ++k) {
uint8_t f;
/*
* The diagonal element cannot be 0 because
* we are inverting matrices with all the square
* submatrices not singular
*/
BUG_ON(M[k * n + k] == 0);
/* make the diagonal element to be 1 */
f = inv(M[k * n + k]);
for (j = 0; j < n; ++j) {
M[k * n + j] = mul(f, M[k * n + j]);
V[k * n + j] = mul(f, V[k * n + j]);
}
/*
* Make all the elements over and under the diagonal
* to be zero
*/
for (i = 0; i < n; ++i) {
if (i == k)
continue;
f = M[i * n + k];
for (j = 0; j < n; ++j) {
M[i * n + j] ^= mul(f, M[k * n + j]);
V[i * n + j] ^= mul(f, V[k * n + j]);
}
}
}
}
/**
* Computes the parity without the missing data blocks
* and store it in the buffers of such data blocks.
*
* This is the parity expressed as Pa,Qa,Ra,Sa,Ta,Ua in the equations.
*/
void raid_delta_gen(int nr, int *id, int *ip, int nd, size_t size, void **v)
{
void *p[RAID_PARITY_MAX];
void *pa[RAID_PARITY_MAX];
int i, j;
int np;
void *latest;
/* zero buffer must be initialized prior to delta recovery */
BUG_ON(raid_zero_block == 0);
/*
* Total number of parities we are going to process
* they are both the used and the unused ones
*/
np = ip[nr - 1] + 1;
/* latest missing data block */
latest = v[id[nr - 1]];
/* setup pointers for delta computation */
for (i = 0, j = 0; i < np; ++i) {
/* keep a copy of the original parity vector */
p[i] = v[nd + i];
if (ip[j] == i) {
/*
* Set used parities to point to the missing
* data blocks.
*
* The related data blocks are instead set
* to point to the "zero" buffer.
*/
/*
* The latest parity to use ends the for loop and
* then it cannot happen to process more of them
*/
BUG_ON(j >= nr);
/* buffer for missing data blocks */
pa[j] = v[id[j]];
/* set at zero the missing data blocks */
v[id[j]] = raid_zero_block;
/* compute the parity over the missing data blocks */
v[nd + i] = pa[j];
/* check for the next used entry */
++j;
} else {
/*
* Unused parities are going to be rewritten with
* not significant data, because we don't have
* functions able to compute only a subset of
* parities.
*
* To avoid this, we reuse parity buffers,
* assuming that all the parity functions write
* parities in order.
*
* We assign the unused parity block to the same
* block of the latest used parity that we know it
* will be written.
*
* This means that this block will be written
* multiple times and only the latest write will
* contain the correct data.
*/
v[nd + i] = latest;
}
}
/* all the parities have to be processed */
BUG_ON(j != nr);
/*
* Recompute the parity, note that np may be smaller than the
* total number of parities available
*/
raid_gen(nd, np, size, v);
/* restore data buffers as before */
for (j = 0; j < nr; ++j)
v[id[j]] = pa[j];
/* restore parity buffers as before */
for (i = 0; i < np; ++i)
v[nd + i] = p[i];
}
/**
* Recover failure of one data block for PAR1.
*
* Starting from the equation:
*
* Pd = Dx
*
* and solving we get:
*
* Dx = Pd
*/
void raid_rec1of1(int *id, int nd, size_t size, void **v)
{
void *p;
void *pa;
/*
* For PAR1 we can directly compute the missing block
* and we don't need to use the zero buffer
*/
p = v[nd];
pa = v[id[0]];
/* use the parity as missing data block */
v[id[0]] = p;
/* compute the parity over the missing data block */
v[nd] = pa;
/* compute */
raid_gen(nd, 1, size, v);
/* restore as before */
v[id[0]] = pa;
v[nd] = p;
}
/**
* Recover failure of two data blocks for PAR2.
*
* Starting from the equations:
*
* Pd = Dx + Dy
* Qd = A[1,id[0]] * Dx + A[1,id[1]] * Dy
*
* and solving we get:
*
* Qd + A[1,id[0]] * Pd
* Dy = ------------------------
* A[1,id[0]] + A[1,id[1]]
*
* Dx = Dy + Pd
*
* with conditions:
*
* A[1,id[0]] + A[1,id[1]] != 0
*
* That are always satisfied for any 0<=id[0]<id[1]<255.
*/
void raid_rec2of2_int8(int *id, int *ip, int nd, size_t size, void **vv)
{
uint8_t **v = (uint8_t **)vv;
size_t i;
uint8_t *p;
uint8_t *pa;
uint8_t *q;
uint8_t *qa;
const uint8_t *T[2];
/* get multiplication tables from the actual active Q coefficients */
T[0] = table(mul(A(1, id[0]), inv(A(1, id[0]) ^ A(1, id[1]))));
T[1] = table(inv(A(1, id[0]) ^ A(1, id[1])));
/* compute delta parity */
raid_delta_gen(2, id, ip, nd, size, vv);
p = v[nd];
q = v[nd + 1];
pa = v[id[0]];
qa = v[id[1]];
for (i = 0; i < size; ++i) {
/* delta */
uint8_t Pd = p[i] ^ pa[i];
uint8_t Qd = q[i] ^ qa[i];
/* reconstruct */
uint8_t Dy = T[0][Pd] ^ T[1][Qd];
uint8_t Dx = Pd ^ Dy;
/* set */
pa[i] = Dx;
qa[i] = Dy;
}
}
Reference in a new issue View git blame Copy permalink