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 
moins_l(Value * a,int i,int n,int p)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 
moins_c(Value * a,int i,int n,int p)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 
echange_l(Value * a,int i,int j,int n,int p)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 
echange_c(Value * a,int i,int j,int n,int p)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 
ligne(Value * a,int i,int j,Value x,int n,int p)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 
colonne(Value * a,int i,int j,Value x,int n,int p)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 
petit_l(Value * a,int n,int p,int q)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 
petit_c(Value * a,int n,int p,int q)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 
identite(Value * a,int n,int p)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 
transpose(Value * a,int n,int q)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 
encore(Value * a,int n,int p,int q,Value val)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 
smith(Value * a,Value * b,Value * c,Value * b_inverse,Value * c_inverse,int n,int p,int q)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 
hermite(Value * a,Value * b,Value * d,int n,int p,int q)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 **/
ConvertPolMattoDarMat(Matrix * A)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 
ConvertDarMattoPolMat(Value * A,int NbRows,int NbCols)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 **/
Smith(Matrix * A,Matrix ** U,Matrix ** V,Matrix ** Product)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 **/
Hermite(Matrix * A,Matrix ** H,Matrix ** U)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 
720