IDA_new/jerasure/Examples/raph_old.c

   1 /* *
   2  * Copyright (c) 2014, James S. Plank and Kevin Greenan
   3  * All rights reserved.
   4  *
   5  * Jerasure - A C/C++ Library for a Variety of Reed-Solomon and RAID-6 Erasure
   6  * Coding Techniques
   7  *
   8  * Revision 2.0: Galois Field backend now links to GF-Complete
   9  *
  10  * Redistribution and use in source and binary forms, with or without
  11  * modification, are permitted provided that the following conditions
  12  * are met:
  13  *
  14  *  - Redistributions of source code must retain the above copyright
  15  *    notice, this list of conditions and the following disclaimer.
  16  *
  17  *  - Redistributions in binary form must reproduce the above copyright
  18  *    notice, this list of conditions and the following disclaimer in
  19  *    the documentation and/or other materials provided with the
  20  *    distribution.
  21  *
  22  *  - Neither the name of the University of Tennessee nor the names of its
  23  *    contributors may be used to endorse or promote products derived
  24  *    from this software without specific prior written permission.
  25  *
  26  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  27  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  28  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  29  * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  30  * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  31  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  32  * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
  33  * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
  34  * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  35  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY
  36  * WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  37  * POSSIBILITY OF SUCH DAMAGE.
  38  */
  39
  40 /* Jerasure's authors:
  41
  42    Revision 2.x - 2014: James S. Plank and Kevin M. Greenan.
  43    Revision 1.2 - 2008: James S. Plank, Scott Simmerman and Catherine D. Schuman.
  44    Revision 1.0 - 2007: James S. Plank.
  45  */
  46
  47 #include <stdio.h>
  48 #include <stdlib.h>
  49 #include <string.h>
  50 #include "jerasure.h"
  51
  52 #define talloc(type, num) (type *) malloc(sizeof(type)*(num))
  53
  54 static void usage(char *s)
  55 {
  56   fprintf(stderr, "usage: jerasure_03 k w - Creates a kxk Cauchy matrix in GF(2^w). \n\n");
  57   fprintf(stderr, "       k must be < 2^w.  Element i,j is 1/(i+(2^w-j-1)).  (If that is\n");
  58   fprintf(stderr, "       If that is 1/0, then it sets it to zero).  \n");
  59   fprintf(stderr, "       It then tests whether that matrix is invertible.\n");
  60   fprintf(stderr, "       If it is invertible, then it prints out the inverse.\n");
  61   fprintf(stderr, "       Finally, it prints the product of the matrix and its inverse.\n");
  62   fprintf(stderr, "       \n");
  63   fprintf(stderr, "This demonstrates: jerasure_print_matrix()\n");
  64   fprintf(stderr, "                   jerasure_invertible_matrix()\n");
  65   fprintf(stderr, "                   jerasure_invert_matrix()\n");
  66   fprintf(stderr, "                   jerasure_matrix_multiply().\n");
  67   if (s != NULL) fprintf(stderr, "%s\n", s);
  68   exit(1);
  69 }
  70
  71
  72 typedef unsigned long mylong;
  73 #define LLUI (long long unsigned int)
  74
  75
  76 void display(mylong *mat, int r, int c) {
  77  for(int i=0;i<r;i++) {
  78     for(int j=0;j<c;j++) {
  79       printf("%016llu ",LLUI mat[i*c+j]);
  80     }
  81     printf("\n");
  82   }
  83  printf("\n");
  84 }
  85
  86 mylong *matrix_multiply(gf_t *gf, mylong *m1, mylong *m2, int r1, int c1, int r2, int c2, int w)
  87 {
  88   mylong *product;
  89   int i, j, k;
  90
  91   product = (mylong *) malloc(sizeof(mylong)*r1*c2);
  92   for (i = 0; i < r1*c2; i++) product[i] = 0;
  93
  94   for (i = 0; i < r1; i++) {
  95     for (j = 0; j < c2; j++) {
  96       for (k = 0; k < r2; k++) {
  97         product[i*c2+j] ^= gf->multiply.w64(gf,m1[i*c1+k], m2[k*c2+j]);
  98       }
  99     }
 100   }
 101   return product;
 102 }
 103
 104
 105
 106 int invert_matrix(gf_t *gf, mylong *mat, mylong *inv, int rows)
 107 {
 108   int cols, i, j, k, x, rs2;
 109   int row_start;
 110   mylong tmp, inverse;
 111
 112   cols = rows;
 113
 114   k = 0;
 115   for (i = 0; i < rows; i++) {
 116     for (j = 0; j < cols; j++) {
 117       inv[k] = (i == j) ? 1 : 0;
 118       k++;
 119     }
 120   }
 121 //  display(inv, rows, rows);
 122 //  printf("\n");
 123
 124   /* First -- convert into upper triangular  */
 125   for (i = 0; i < cols; i++) {
 126     row_start = cols*i;
 127
 128     /* Swap rows if we ave a zero i,i element.  If we can't swap, then the
 129        matrix was not invertible  */
 130
 131     if (mat[row_start+i] == 0) {
 132       for (j = i+1; j < rows && mat[cols*j+i] == 0; j++) ;
 133       if (j == rows) return -1;
 134       rs2 = j*cols;
 135       for (k = 0; k < cols; k++) {
 136         tmp = mat[row_start+k];
 137         mat[row_start+k] = mat[rs2+k];
 138         mat[rs2+k] = tmp;
 139         tmp = inv[row_start+k];
 140         inv[row_start+k] = inv[rs2+k];
 141         inv[rs2+k] = tmp;
 142       }
 143     }
 144
 145     /* Multiply the row by 1/element i,i  */
 146     tmp = mat[row_start+i];
 147     if (tmp != 1) {
 148       inverse = gf->divide.w64(gf,1, tmp);
 149       for (j = 0; j < cols; j++) {
 150         mat[row_start+j] = gf->multiply.w64(gf,mat[row_start+j], inverse);
 151         inv[row_start+j] = gf->multiply.w64(gf,inv[row_start+j], inverse);
 152       }
 153     }
 154
 155     /* Now for each j>i, add A_ji*Ai to Aj  */
 156     k = row_start+i;
 157     for (j = i+1; j != cols; j++) {
 158       k += cols;
 159       if (mat[k] != 0) {
 160         if (mat[k] == 1) {
 161           rs2 = cols*j;
 162           for (x = 0; x < cols; x++) {
 163             mat[rs2+x] ^= mat[row_start+x];
 164             inv[rs2+x] ^= inv[row_start+x];
 165           }
 166         } else {
 167           tmp = mat[k];
 168           rs2 = cols*j;
 169           for (x = 0; x < cols; x++) {
 170             mat[rs2+x] ^= gf->multiply.w64(gf,tmp, mat[row_start+x]);
 171             inv[rs2+x] ^= gf->multiply.w64(gf,tmp, inv[row_start+x]);
 172           }
 173         }
 174       }
 175     }
 176   }
 177
 178   /* Now the matrix is upper triangular.  Start at the top and multiply down  */
 179
 180   for (i = rows-1; i >= 0; i--) {
 181     row_start = i*cols;
 182     for (j = 0; j < i; j++) {
 183       rs2 = j*cols;
 184       if (mat[rs2+i] != 0) {
 185         tmp = mat[rs2+i];
 186         mat[rs2+i] = 0;
 187         for (k = 0; k < cols; k++) {
 188           inv[rs2+k] ^= gf->multiply.w64(gf,tmp, inv[row_start+k]);
 189         }
 190       }
 191     }
 192   }
 193
 194 /*  printf("mat\n");
 195   display(mat, rows, rows);
 196   printf("\n");
 197    printf("inv\n");
 198   display(inv, rows, rows);
 199   printf("\n");
 200 */
 201   return 0;
 202 }
 203
 204
 205
 206
 207 int invertible_matrix(gf_t *gf, int *mat, int rows, int w)
 208 {
 209   int cols, i, j, k, x, rs2;
 210   int row_start;
 211   ulong tmp, inverse;
 212
 213   cols = rows;
 214
 215   /* First -- convert into upper triangular  */
 216   for (i = 0; i < cols; i++) {
 217     row_start = cols*i;
 218
 219     /* Swap rows if we ave a zero i,i element.  If we can't swap, then the
 220        matrix was not invertible  */
 221
 222     if (mat[row_start+i] == 0) {
 223       for (j = i+1; j < rows && mat[cols*j+i] == 0; j++) ;
 224       if (j == rows) return 0;
 225       rs2 = j*cols;
 226       for (k = 0; k < cols; k++) {
 227         tmp = mat[row_start+k];
 228         mat[row_start+k] = mat[rs2+k];
 229         mat[rs2+k] = tmp;
 230       }
 231     }
 232
 233     /* Multiply the row by 1/element i,i  */
 234     tmp = mat[row_start+i];
 235     if (tmp != 1) {
 236       inverse =  gf->divide.w64(gf,1, tmp);
 237       for (j = 0; j < cols; j++) {
 238         mat[row_start+j] =  gf->multiply.w64(gf,mat[row_start+j], inverse);
 239       }
 240     }
 241
 242     /* Now for each j>i, add A_ji*Ai to Aj  */
 243     k = row_start+i;
 244     for (j = i+1; j != cols; j++) {
 245       k += cols;
 246       if (mat[k] != 0) {
 247         if (mat[k] == 1) {
 248           rs2 = cols*j;
 249           for (x = 0; x < cols; x++) {
 250             mat[rs2+x] ^= mat[row_start+x];
 251           }
 252         } else {
 253           tmp = mat[k];
 254           rs2 = cols*j;
 255           for (x = 0; x < cols; x++) {
 256             mat[rs2+x] ^=  gf->multiply.w64(gf,tmp,mat[row_start+x]);
 257           }
 258         }
 259       }
 260     }
 261   }
 262   return 1;
 263 }
 264
 265
 266
 267 int main(int argc, char **argv)
 268 {
 269   unsigned int k, w, i, j, n;
 270   mylong *matrix;
 271   mylong *matrix_copy;
 272   mylong *inverse;
 273   mylong *identity;
 274
 275
 276   k=8;
 277   // if (argc != 2) usage(NULL);
 278   //if (sscanf(argv[1], "%d", &k) == 0 || k <= 0) usage("Bad k");
 279   /* if (sscanf(argv[2], "%d", &w) == 0 || w <= 0 || w > 31) usage("Bad w");
 280   if (k >= (1 << w)) usage("K too big");
 281   */
 282   w=64;
 283
 284
 285   matrix = malloc(sizeof(mylong)*k*k);
 286   matrix_copy = malloc(sizeof(mylong)* k*k);
 287   inverse = malloc(sizeof(mylong)* k*k);
 288
 289   /*for (i = 0; i < k; i++) {
 290     for (j = 0; j < k; j++) {
 291       n = i ^ ((1 << w)-1-j);
 292       matrix[i*k+j] = (n == 0) ? 0 : galois_single_divide(1, n, w);
 293     }
 294     }*/
 295
 296   srand48(time(0));
 297   for(int i=0;i<k;i++)
 298     for(int j=0;j<k;j++) {
 299       //mat1[i*c+j]=i*c+j;
 300       matrix[i*k+j]=lrand48()<<32|lrand48();//k*k-(j*k+i);
 301     }
 302 /*   matrix[3]+=34;
 303   matrix[7]+=134;
 304   matrix[13]+=234;
 305   matrix[5]+=67;
 306   matrix[9]+=24;
 307   matrix[2]+=490;
 308 */
 309   gf_t gf;
 310   gf_init_easy(&gf, w);
 311
 312
 313   /* printf("<HTML><TITLE>jerasure_03");
 314   for (i = 1; i < argc; i++) printf(" %s", argv[i]);
 315   printf("</TITLE>\n");
 316   printf("<h3>jerasure_03");
 317   for (i = 1; i < argc; i++) printf(" %s", argv[i]);
 318   printf("</h3>\n");
 319   printf("<pre>\n");
 320   */
 321   printf("Matrix:\n");
 322   //jerasure_print_matrix(matrix, k, k, w);
 323   display(matrix, k, k);
 324
 325   memcpy(matrix_copy, matrix, sizeof(mylong)*k*k);
 326   i = 1;
 327   //jerasure_invertible_matrix(matrix_copy, k, w);
 328   invertible_matrix(&gf,matrix_copy, k, w);
 329   printf("\nInvertible: %s\n", (i == 1) ? "Yes" : "No");
 330   if (i == 1) {
 331     printf("\nInverse:\n");
 332     memcpy(matrix_copy, matrix, sizeof(mylong)*k*k);
 333     //   i = jerasure_invert_matrix(matrix_copy, inverse, k, w);
 334
 335     invert_matrix(&gf, matrix_copy, inverse, k);
 336
 337     //jerasure_print_matrix(inverse, k, k, w);
 338     display(inverse, k, k);
 339
 340     //identity = jerasure_matrix_multiply(inverse, matrix, k, k, k, k, w);
 341     identity=matrix_multiply(&gf, inverse, matrix, k, k, k, k, w);
 342     printf("\nInverse times matrix (should be identity):\n");
 343     display(identity, k, k);
 344     //jerasure_print_matrix(identity, k, k, w);
 345   }
 346   return 0;
 347 }
 348
 349