source/kernel/NormalForms.c

   1 /*
   2     This file is part of PolyLib.
   3
   4     PolyLib is free software: you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation, either version 3 of the License, or
   7     (at your option) any later version.
   8
   9     PolyLib is distributed in the hope that it will be useful,
  10     but WITHOUT ANY WARRANTY; without even the implied warranty of
  11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License
  15     along with PolyLib.  If not, see <http://www.gnu.org/licenses/>.
  16 */
  17
  18 #include <stdlib.h>
  19 #include <polylib/polylib.h>
  20
  21 /* nota bene: on stocke les matrices par lignes */
  22 /* nota bene: matrices are stored in row major order */
  23
  24 /*------------------------------------------------------------------------
  25         change les signes de la ligne i de la matrice A
  26         change the sign of row i of matrix A
  27 ------------------------------------------------------------------------*/
  28
  29 static void moins_l(Value *a,int i,int n,int p) {
  30
  31   int k;
  32   Value *c;
  33
  34   c=a+(i-1)*p;
  35
  36   for(k=1; k<=p; k++) {
  37     value_oppose(*c,*c);
  38     c++;
  39   }
  40   return;
  41 } /* moins_l */
  42
  43 /*------------------------------------------------------------------------
  44         change les signes de la colonne i de la matrice A
  45         change the sign of column i of matrix A
  46 ------------------------------------------------------------------------*/
  47
  48 static void moins_c(Value *a,int i,int n,int p) {
  49
  50   int k;
  51   Value *c;
  52
  53   c=a+(i-1);
  54
  55   for(k=1;k<=n;k++) {
  56     value_oppose(*c,*c);
  57     c=c+p;
  58   }
  59   return;
  60 } /* moins_c */
  61
  62 /*------------------------------------------------------------------------
  63         echange les lignes i et j de la matrice A
  64         exchange row i and j of matrix A
  65 ------------------------------------------------------------------------*/
  66
  67 static void echange_l(Value *a,int i,int j,int n,int p) {
  68
  69   int k;
  70   Value s;
  71   Value *c1,*c2;
  72
  73   value_init(s);
  74   c1=a+(i-1)*p;
  75   c2=a+(j-1)*p;
  76
  77   for(k=1;k<=p;k++) {
  78     value_assign(s,*c1);
  79     value_assign(*c1,*c2);
  80     value_assign(*c2,s);
  81     c1++;
  82     c2++;
  83   }
  84   value_clear(s);
  85   return;
  86 } /* echange_l */
  87
  88 /*------------------------------------------------------------------------
  89         echange les colonnes i et j de la matrice A
  90         exchange columns i and j of matrix A
  91 ------------------------------------------------------------------------*/
  92
  93 static void echange_c(Value *a,int i,int j,int n,int p) {
  94
  95   int k;
  96   Value s;
  97   Value *c1,*c2;
  98
  99   value_init(s);
 100   c1=a+(i-1);
 101   c2=a+(j-1);
 102
 103   for(k=1;k<=n;k++) {
 104     value_assign(s,*c1);
 105     value_assign(*c1,*c2);
 106     value_assign(*c2,s);
 107     c1=c1+p;
 108     c2=c2+p;
 109   }
 110   value_clear(s);
 111   return;
 112 } /* echange_c */
 113
 114 /*------------------------------------------------------------------------
 115         A est une matrice de taille n*p
 116         on ajoute a la ligne i, x fois la ligne j
 117         A is a matrix of size n*p
 118         we add x times the row j to the row i
 119 ------------------------------------------------------------------------*/
 120
 121 static void ligne(Value *a,int i,int j,Value x,int n,int p) {
 122
 123   int k;
 124   Value *c1,*c2;
 125
 126   c1=a+(i-1)*p;
 127   c2=a+(j-1)*p;
 128
 129   for(k=1;k<=p;k++) {
 130
 131     value_addmul(*c1, x, *c2);
 132     c1++;
 133     c2++;
 134   }
 135   return;
 136 } /* ligne */
 137
 138 /*------------------------------------------------------------------------
 139         A est une matrice de taille n*p
 140         on ajoute a la colonne i, x fois la colonne j
 141         A is a matrix of size n*p
 142         we add x times the column j to the column i
 143
 144 ------------------------------------------------------------------------*/
 145
 146 static void colonne(Value *a,int i,int j,Value x,int n,int p) {
 147
 148   int k;
 149   Value *c1,*c2;
 150
 151   c1=a+(i-1);
 152   c2=a+(j-1);
 153
 154   for(k=1;k<=n;k++) {
 155     value_addmul(*c1, x, *c2);
 156     c1=c1+p;
 157     c2=c2+p;
 158   }
 159   return;
 160 } /* colonne */
 161
 162 /*----------------------------------------------------------------------
 163         trouve le numero de colonne du plus petit element non nul de
 164         la ligne q, d'indice superieur a q, de la matrice A,  mais
 165         renvoie zero si tous les elements sont nuls sauf le qieme.
 166
 167         find the column number (greater than q) of the smallest non
 168         null element of row q . it returns
 169         0 if all the last elements of row q are null except the qth
 170
 171 ----------------------------------------------------------------------*/
 172
 173 static int petit_l(Value *a,int n,int p,int q) {
 174
 175   int numero=0, k, tousnuls;
 176   Value minus, comp;
 177   Value *c;
 178
 179   value_init(minus); value_init(comp);
 180   c=a+(q-1)*p+(p-1);
 181   tousnuls=1;
 182   for(k=p;k>q;k--) {
 183     value_absolute(comp,*c);
 184     if (value_notzero_p(comp)) {
 185       if (tousnuls==1) {
 186         value_assign(minus,comp);
 187         numero=k;
 188         tousnuls=0;
 189       }
 190       else if (value_ge(minus,comp)) {
 191         value_assign(minus,comp);
 192         numero=k;
 193       }
 194     }
 195     c--;
 196   }
 197   if (tousnuls==1) {
 198     value_clear(minus); value_clear(comp);
 199     return(0);
 200   }
 201   else  {
 202     value_absolute(comp,*c);
 203     if ((value_notzero_p(comp))&&(value_ge(minus,comp))) {
 204       value_clear(minus); value_clear(comp);
 205       return(q);
 206     }
 207     else {
 208       value_clear(minus); value_clear(comp);
 209       return(numero);
 210     }
 211   }
 212 } /* petit_l */
 213
 214 /*----------------------------------------------------------------------
 215         trouve le numero de ligne du plus petit element non nul de
 216         la colonne q, d'indice superieur a q, de la matrice A,  mais
 217         renvoie zero si tous les elements sont nuls sauf le qieme.
 218
 219         find the row number (greater than q) of the smallest non
 220         null element of column q . it returns
 221         0 if all the last elements of column q are null except the qth
 222
 223 ----------------------------------------------------------------------*/
 224
 225 static int petit_c(Value *a,int n,int p,int q) {
 226
 227   int numero=0, k, tousnuls;
 228   Value minus, comp;
 229   Value *c;
 230
 231   value_init(minus); value_init(comp);
 232   c = a+q-1+p*(n-1);
 233   tousnuls=1;
 234   for(k=n;k>q;k--) {
 235     value_absolute(comp,*c);
 236     if (value_notzero_p(comp)) {
 237       if (tousnuls==1) {
 238         value_assign(minus,comp);
 239         numero=k;
 240         tousnuls=0;
 241       }
 242       else if (value_ge(minus,comp)) {
 243         value_assign(minus,comp);
 244         numero=k;
 245       }
 246     }
 247     c=c-p;
 248   }
 249
 250   if (tousnuls==1) {
 251     value_clear(minus); value_clear(comp);
 252     return(0);
 253   }
 254   else {
 255     value_absolute(comp,*c);
 256     if ((value_notzero_p(comp)) && (value_ge(minus,comp))) {
 257       value_clear(minus); value_clear(comp);
 258       return(q);
 259     }
 260     else {
 261       value_clear(minus); value_clear(comp);
 262       return(numero);
 263     }
 264   }
 265 } /* petit_c */
 266
 267 /*-----------------------------------------------------------------------
 268         A est initialisee a une matrice "identite" de taille n*p
 269         A is initialized to an "identity" matrix of size n*p
 270 -----------------------------------------------------------------------*/
 271
 272 static void identite(Value *a,int n,int p) {
 273
 274   int i,j;
 275   Value *b;
 276
 277   b = a;
 278   for(i=1;i<=n;i++) {
 279     for(j=1;j<=p;j++) {
 280       if (i==j)
 281         value_set_si(*b,1);
 282       else
 283         value_set_si(*b,0);
 284       b++;
 285     }
 286   }
 287   return;
 288 } /* identite */
 289
 290 /*-----------------------------------------------------------------------
 291   transpose la sous-matrice de A dont les indices sont
 292   superieurs ou egaux a q
 293   transposition of the sub-matrix of A composed of the elements
 294   whose indices are greater or equal to q
 295 -----------------------------------------------------------------------*/
 296
 297 static void transpose(Value *a,int n,int q) {
 298
 299   int i;
 300   Value val;
 301   Value *b,*c;
 302
 303   value_init(val);
 304   if (q<n) {
 305     b=a+(q-1)+(q-1)*n;
 306     c=b;
 307     for(i=q+1;i<=n;i++) {
 308       b++;
 309       c=c+n;
 310       value_assign(val,*b);
 311       value_assign(*b,*c);
 312       value_assign(*c,val);
 313     }
 314     transpose(a,n,q+1);
 315   }
 316   value_clear(val);
 317   return;
 318 } /* transpose */
 319
 320 /*----------------------------------------------------------------------
 321         trouve le numero de colonne du premier element non divisible
 322         par val dans la sous-matrice d'indices > q mais renvoie 0 sinon
 323
 324         find the column number of the first non-divisible by val element
 325         in the sub-matrix of a composed of the elements of indices >q.
 326         Return 0 is there is no such element
 327 ----------------------------------------------------------------------*/
 328
 329 static int encore(Value *a,int n,int p,int q,Value val) {
 330
 331   int k,l;
 332   Value comp, tmp;
 333   Value *c;
 334
 335   if ((q>=n)||(q>=p)) return(0);
 336
 337   value_init(comp); value_init(tmp);
 338   c=a+q*p+q;
 339   k=q+1;
 340   l=q+1;
 341   while (k<=n) {
 342     value_absolute(comp,*c);
 343     if (value_zero_p(val)) {
 344       if (value_notzero_p(comp)) {
 345         value_clear(comp);
 346         value_clear(tmp);
 347         return(l);
 348       }
 349     }
 350     else {
 351       value_modulus(tmp,comp,val);
 352       if (value_notzero_p(tmp)) {
 353         value_clear(comp);
 354         value_clear(tmp);
 355         return(l);
 356       }
 357     }
 358     if (l==p) {
 359       k=k+1;
 360       l=q+1;
 361       c=c+q+1;
 362     }
 363     else {
 364       l++;
 365       c++;
 366     }
 367   }
 368   value_clear(comp);
 369   value_clear(tmp);
 370   return(0);
 371 } /* encore */
 372
 373 /*----------------------------------------------------------------------
 374         transforme la matrice A (de taille n*p), a partir de la position
 375         q*q en sa forme de smith, les matrices b et c subissant les memes
 376         transformations respectivement sur les lignes et colonnes
 377
 378         transform the matrix A (size n*p), from position (q,q) into
 379         its smith normal form. the same transformations are applied to
 380         matrices b and c, respectively on rows and on columns
 381 ----------------------------------------------------------------------*/
 382
 383 static void smith(Value *a,Value *b,Value *c,Value *b_inverse,Value *c_inverse,int n,int p,int q) {
 384
 385   int i,k,fini;
 386   Value x, pivot, tmp, x_inv;
 387   Value *f;
 388
 389   value_init(pivot); value_init(tmp);
 390   value_init(x); value_init(x_inv);
 391
 392   if ((q<=n)&&(q<=p)) {
 393     fini = 1;
 394
 395     while (fini!=0) {
 396       i=petit_c(a,n,p,q);
 397       while (i!=0) {
 398         if (i!=q) {
 399           echange_l(a,i,q,n,p);
 400           echange_l(b,i,q,n,n);
 401           echange_c(b_inverse,i,q,n,n);
 402         }
 403         f=a+(q-1)+(q-1)*p;
 404         value_assign(pivot,*f);
 405         if (value_neg_p(pivot)) {
 406           moins_l(a,q,n,p);
 407           moins_l(b,q,n,n);
 408           moins_c(b_inverse,q,n,n);
 409           value_oppose(pivot,pivot);
 410         }
 411         for(k=q+1;k<=n;k++) {
 412           f=f+p;
 413           if (value_notzero_p(*f)) {
 414             value_division(x,*f,pivot);
 415             value_modulus(tmp,*f,pivot);
 416             if (value_neg_p(tmp))
 417               value_decrement(x,x);
 418             value_oppose(x_inv,x);
 419             ligne(a,k,q,x_inv,n,p);
 420             ligne(b,k,q,x_inv,n,n);
 421             colonne(b_inverse,q,k,x,n,n);
 422           }
 423         }
 424         i=petit_c(a,n,p,q);
 425       }
 426       fini=0;
 427       i=petit_l(a,n,p,q);
 428       while (i!=0) {
 429         if (i!=q) {
 430           echange_c(a,i,q,n,p);
 431           echange_c(c,i,q,p,p);
 432           echange_l(c_inverse,i,q,p,p);
 433           fini=1;
 434         }
 435         f=a+(q-1)+(q-1)*p;
 436         value_assign(pivot,*f);
 437         if (value_neg_p(pivot)) {
 438           moins_c(a,q,n,p);
 439           moins_c(c,q,p,p);
 440           moins_l(c_inverse,q,p,p);
 441           value_oppose(pivot,pivot);
 442         }
 443         for(k=q+1;k<=p;k++) {
 444           f++;
 445           if (value_notzero_p(*f)) {
 446             value_division(x,*f,pivot);
 447             value_modulus(tmp,*f,pivot);
 448             if (value_neg_p(tmp))
 449               value_decrement(x,x);
 450             value_oppose(x_inv,x);
 451             colonne(a,k,q,x_inv,n,p);
 452             colonne(c,k,q,x_inv,p,p);
 453             ligne(c_inverse,q,k,x,p,p);
 454           }
 455         }
 456         i=petit_l(a,n,p,q);
 457       }
 458     }
 459     value_assign(pivot,*(a+(q-1)+(q-1)*p));
 460     if (value_neg_p(pivot)) {
 461       moins_l(a,q,n,p);
 462       moins_l(b,q,n,n);
 463       moins_c(b_inverse,q,n,n);
 464       value_oppose(pivot,pivot);
 465     }
 466
 467     i=encore(a,n,p,q,pivot);
 468     if (i!=0) {
 469       value_set_si(x,1);
 470       value_set_si(x_inv,-1);
 471       colonne(a,q,i,x,n,p);
 472       colonne(c,q,i,x,p,p);
 473       ligne(c_inverse,i,q,x_inv,p,p);
 474       smith(a,b,c,b_inverse,c_inverse,n,p,q);
 475     }
 476     else smith(a,b,c,b_inverse,c_inverse,n,p,q+1);
 477   }
 478   value_clear(pivot); value_clear(tmp);
 479   value_clear(x); value_clear(x_inv);
 480
 481   return;
 482 } /* smith */
 483
 484 /*----------------------------------------------------------------------
 485         decompose la matrice A en le produit d'une matrice B
 486         unimodulaire et d'une matrice triangulaire superieure
 487         D est l'inverse de B
 488         Decompose the matrix A in the product of a matrix B unimodular
 489         and a matrix upper triangular, D is the inverse of B
 490 ----------------------------------------------------------------------*/
 491
 492 static void hermite(Value *a,Value *b,Value *d,int n,int p,int q) {
 493
 494   int i,k;
 495   Value x, pivot, x_inv, tmp;
 496   Value *c1;
 497
 498   value_init(pivot);  value_init(tmp);
 499   value_init(x); value_init(x_inv);
 500
 501   if ((q<=p)&&(q<=n)) {
 502     i=petit_c(a,n,p,q);
 503     while (i!=0) {
 504       if (i!=q) {
 505         echange_l(a,i,q,n,p);
 506         echange_c(b,i,q,n,n);
 507         echange_l(d,i,q,n,n);
 508       }
 509
 510       c1=a+(q-1)+(q-1)*p;
 511       value_assign(pivot,*c1);
 512       if (value_neg_p(pivot)) {
 513         moins_l(a,q,n,p);
 514         moins_l(d,q,n,n);
 515         moins_c(b,q,n,n);
 516         value_oppose(pivot,pivot);
 517       }
 518       for(k=q+1;k<=n;k++) {
 519         c1=c1+p;
 520         if (value_notzero_p(*c1)) {
 521           value_division(x,*c1,pivot);
 522           value_modulus(tmp,*c1,pivot);
 523           if (value_neg_p(tmp))
 524             value_decrement(x,x);
 525           value_oppose(x_inv,x);
 526           ligne(a,k,q,x_inv,n,p);
 527           colonne(b,q,k,x,n,n);
 528           ligne(d,k,q,x_inv,n,n);
 529         }
 530       }
 531       i=petit_c(a,n,p,q);
 532     }
 533     c1=a+(q-1);
 534     value_assign(pivot,*(c1+(q-1)*p));
 535     if (value_neg_p(pivot)) {
 536       moins_l(a,q,n,p);
 537       moins_l(d,q,n,n);
 538       moins_c(b,q,n,n);
 539       value_oppose(pivot,pivot);
 540     }
 541     if (value_notzero_p(pivot)) {
 542       for(k=1;k<q;k++) {
 543         if (value_notzero_p(*c1)) {
 544           value_division(x,*c1,pivot);
 545           value_modulus(tmp,*c1,pivot);
 546           if (value_neg_p(tmp))
 547             value_decrement(x,x);
 548           value_oppose(x_inv,x);
 549           ligne(a,k,q,x_inv,n,p);
 550           colonne(b,q,k,x,n,n);
 551           ligne(d,k,q,x_inv,n,n);
 552         }
 553         c1=c1+p;
 554       }
 555     }
 556     hermite(a,b,d,n,p,q+1);
 557   }
 558   value_clear(pivot); value_clear(tmp);
 559   value_clear(x); value_clear(x_inv);
 560   return;
 561 } /* hermite */
 562
 563 /*----------------------------------------------------------------------
 564   B est une g_base de dimension n dont les p premiers vecteurs
 565   colonnes sont colineaires aux vecteurs colonnes de A
 566   C est l'invers de B
 567
 568   B is an integral (g_base ?) of dimension n whose p first columns
 569   vectors are proportionnal to the column vectors of A. C is the
 570   inverse of B.
 571   ----------------------------------------------------------------------*/
 572
 573 /** Convert PolmattoDarmat :
 574 ***  This function converts the matrix of a Polylib to a int * as necessary
 575 ***    for the functions in Darte's implementation.
 576 ***
 577 ***  Input  : A Polylib Matrix ( a pointer to it );
 578 ***  Output : An int * with the matrix copied into it
 579 **/
 580 static Value *ConvertPolMattoDarMat(Matrix *A ) {
 581
 582   int i,j;
 583   Value *result;
 584
 585   result = (Value *)malloc(sizeof(Value) * A->NbRows * A->NbColumns);
 586   for(i=0;i<A->NbRows * A->NbColumns;i++)
 587     value_init(result[i]);
 588
 589   for (i = 0; i < A->NbRows; i++)
 590     for (j = 0 ;  j < A->NbColumns; j++)
 591       value_assign(result[i*A->NbColumns + j],A->p[i][j]);
 592   return result;
 593 } /* ConvertPolMattoDarMat */
 594
 595 /** Convert DarmattoPolmat
 596 ***   This function converts the matrix from Darte representation to a matrix
 597 ***    in PolyLib.
 598 ***
 599 ***   Input  : The matrix (a pointer to it),  number of Rows, number of columns
 600 ***   Output : The matrix of the PolyLib
 601 ***
 602 **/
 603
 604 static Matrix *ConvertDarMattoPolMat (Value *A, int NbRows, int NbCols) {
 605
 606   int i,j;
 607   Matrix *result;
 608
 609   result = Matrix_Alloc(NbRows, NbCols);
 610
 611   for (i = 0; i < NbRows; i ++)
 612     for (j = 0; j < NbCols; j ++)
 613       value_assign(result->p[i][j],A[i*NbCols + j]);
 614   return result;
 615 } /* ConvertDarMattoPolMat */
 616
 617 /**
 618 *** Smith : This function takes a Matrix A of dim n * l as its input
 619 ***         and returns the three matrices U, V and Product such that
 620 ***         A = U * Product * V, where U is an unimodular matrix of dimension
 621 ***         n * n and V is an unimodular matrix of dimension l * l.
 622 ***         Product is a diagonal matrix of dimension n * l.
 623 ***
 624 ***         We use Alan Darte's implementation of Smith for computing
 625 ***         the Smith Normal Form
 626 **/
 627 void Smith(Matrix *A, Matrix **U, Matrix **V, Matrix **Product)
 628 {
 629   int i;
 630   Matrix *u, *v;
 631
 632   u = Identity(A->NbRows);
 633   v = Identity(A->NbColumns);
 634
 635   *U = Identity(A->NbRows);
 636   *V = Identity(A->NbColumns);
 637
 638   *Product = Matrix_Copy(A);
 639   smith((*Product)->p_Init, u->p_Init, v->p_Init, (*U)->p_Init, (*V)->p_Init,
 640         A->NbRows, A->NbColumns, 1);
 641
 642   Matrix_Free(u);
 643   Matrix_Free(v);
 644 } /* Smith */
 645
 646 /** Hermite :
 647 ***   This function takes a Matrix as its input and finds its HNF
 648 ***    ( Left form )
 649 ***
 650 ***   Input  : A Matrix A (The Matrix A is not necessarily a Polylib matrix.
 651 ***             It is just a matrix as far as Hermite is concerned. It
 652 ***             does not even check if the matrix is a Polylib matrix or not)
 653 ***   Output : The Hnf matrix H and the Unimodular matrix U such that
 654 ***               A = H * U.
 655 ***
 656 ***   We use Alan Darte's implementation of Hermite to compute the HNF.
 657 ***   Alan Darte's implementation computes the Upper Triangular HNF.
 658 ***   So We work on the fact that if A = H * U then
 659 ***       A (transpose) = U(transpose) * H (transpose)
 660 ***       There are a set of interface functions written in Interface.c
 661 ***        which convert a matrix from Polylib representationt to that of
 662 ***        Alan Darte's and vice versa.
 663 ***
 664 ***   This Function Does the Following
 665 ***    Step 1 : Given the matrix A it finds its Transpose.
 666 ***    Step 2 : Finds the HNF (Right Form) using Alan Darte's Algorithm.
 667 ***    Step 3 : The H1 and U1 obtained in Step2 are both Transposed to get
 668 ***              the actual H and U such that A = HU.
 669 **/
 670 void Hermite (Matrix *A, Matrix **H, Matrix **U) {
 671
 672   int i;
 673   Matrix *transpose, *tempH, *tempU ;
 674   Value *darte_matA, *darte_identite, *darte_id_inv;
 675
 676   transpose = Transpose(A);
 677   darte_matA = ConvertPolMattoDarMat(transpose);
 678
 679   darte_identite = (Value *)malloc(sizeof(Value) * A->NbColumns* A->NbColumns);
 680   darte_id_inv = (Value *) malloc(sizeof(Value) * A->NbColumns*A->NbColumns);
 681   for (i=0; i< A->NbColumns * A->NbColumns; i++)
 682     value_init(darte_identite[i]);
 683   for (i=0; i< A->NbColumns * A->NbColumns; i++)
 684     value_init(darte_id_inv[i]);
 685
 686   identite(darte_identite, transpose->NbRows, transpose->NbRows);
 687   identite(darte_id_inv, transpose->NbRows, transpose->NbRows);
 688   hermite(darte_matA, darte_identite, darte_id_inv,A->NbColumns,A->NbRows, 1);
 689   Matrix_Free (transpose);
 690   transpose = ConvertDarMattoPolMat(darte_matA, A->NbColumns, A->NbRows);
 691   tempU = ConvertDarMattoPolMat(darte_identite, A->NbColumns, A->NbColumns);
 692
 693   /* Finding H Transpose */
 694   tempH = Transpose(transpose);
 695   Matrix_Free(transpose);
 696
 697   /* Finding U Transpose */
 698   transpose = Transpose(tempU);
 699
 700   H[0] = Matrix_Copy(tempH);
 701   U[0] = Matrix_Copy(transpose);
 702
 703   Matrix_Free (tempH);
 704   Matrix_Free (transpose);
 705   Matrix_Free (tempU);
 706
 707   for (i=0; i< A->NbRows * A->NbColumns; i++)
 708     value_clear(darte_matA[i]);
 709   for (i=0; i< A->NbColumns * A->NbColumns; i++)
 710     value_clear(darte_identite[i]);
 711   for (i=0; i< A->NbColumns * A->NbColumns; i++)
 712     value_clear(darte_id_inv[i]);
 713   free (darte_matA);
 714   free (darte_identite);
 715   free (darte_id_inv);
 716   return;
 717 } /* Hermite */
 718
 719