2 This file is part of PolyLib.
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.
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.
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/>.
20 Both this software and its documentation are
22 Copyright 1993 by IRISA /Universite de Rennes I -
23 France, Copyright 1995,1996 by BYU, Provo, Utah
32 #include <polylib/polylib.h>
39 * Allocate space for matrix dimensioned by 'NbRows X NbColumns'.
41 Matrix
*Matrix_Alloc(unsigned NbRows
,unsigned NbColumns
) {
47 Mat
=(Matrix
*)malloc(sizeof(Matrix
));
49 errormsg1("Matrix_Alloc", "outofmem", "out of memory space");
53 Mat
->NbColumns
=NbColumns
;
54 if (NbRows
==0 || NbColumns
==0) {
56 Mat
->p_Init
= (Value
*)0;
59 q
= (Value
**)malloc(NbRows
* sizeof(*q
));
62 errormsg1("Matrix_Alloc", "outofmem", "out of memory space");
65 p
= value_alloc(NbRows
* NbColumns
, &Mat
->p_Init_size
);
69 errormsg1("Matrix_Alloc", "outofmem", "out of memory space");
74 for (i
=0;i
<NbRows
;i
++) {
86 * Free the memory space occupied by Matrix 'Mat'
88 void Matrix_Free(Matrix
*Mat
)
91 value_free(Mat
->p_Init
, Mat
->p_Init_size
);
99 void Matrix_Extend(Matrix
*Mat
, unsigned NbRows
)
104 q
= (Value
**)realloc(Mat
->p
, NbRows
* sizeof(*q
));
106 errormsg1("Matrix_Extend", "outofmem", "out of memory space");
110 if (Mat
->p_Init_size
< NbRows
* Mat
->NbColumns
) {
111 p
= (Value
*)realloc(Mat
->p_Init
, NbRows
* Mat
->NbColumns
* sizeof(Value
));
113 errormsg1("Matrix_Extend", "outofmem", "out of memory space");
117 Vector_Set(Mat
->p_Init
+ Mat
->NbRows
*Mat
->NbColumns
, 0,
118 Mat
->p_Init_size
- Mat
->NbRows
*Mat
->NbColumns
);
119 for (i
= Mat
->p_Init_size
; i
< Mat
->NbColumns
*NbRows
; ++i
)
120 value_init(Mat
->p_Init
[i
]);
121 Mat
->p_Init_size
= Mat
->NbColumns
*NbRows
;
123 Vector_Set(Mat
->p_Init
+ Mat
->NbRows
*Mat
->NbColumns
, 0,
124 (NbRows
- Mat
->NbRows
) * Mat
->NbColumns
);
125 for (i
=0;i
<NbRows
;i
++) {
126 Mat
->p
[i
] = Mat
->p_Init
+ (i
* Mat
->NbColumns
);
128 Mat
->NbRows
= NbRows
;
132 * Print the contents of the Matrix 'Mat'
134 void Matrix_Print(FILE *Dst
, const char *Format
, Matrix
*Mat
)
138 unsigned NbRows
, NbColumns
;
140 fprintf(Dst
,"%d %d\n", NbRows
=Mat
->NbRows
, NbColumns
=Mat
->NbColumns
);
145 for (i
=0;i
<NbRows
;i
++) {
147 for (j
=0;j
<NbColumns
;j
++) {
149 value_print(Dst
," "P_VALUE_FMT
" ",*p
++);
152 value_print(Dst
,Format
,*p
++);
160 * Read the contents of the Matrix 'Mat'
162 Matrix
*Matrix_Read_Input(Matrix
*Mat
) {
169 for (i
=0;i
<Mat
->NbRows
;i
++) {
172 c
= fgets(s
, 1024, stdin
);
175 /* jump to the first non space char */
176 while(isspace(*c
) && *c
!='\n' && *c
)
179 while(*c
=='#' || *c
== '\n'); /* continue if the line is empty or is a comment */
181 errormsg1( "Matrix_Read", "baddim", "not enough rows" );
187 /* skip to EOL and ignore the rest of the line if longer than 1024 */
188 errormsg1( "Matrix_Read", "warning", "line longer than 1024 chars (ignored remaining chars till EOL)" );
191 while( n
!=EOF
&& n
!='\n' );
195 for (j
=0;j
<Mat
->NbColumns
;j
++)
198 if(*c
=='\n' || *c
=='#' || *c
=='\0') {
199 errormsg1("Matrix_Read", "baddim", "not enough columns");
202 /* go the the next space or end */
203 for( z
=c
; *z
; z
++ )
205 if( *z
=='\n' || *z
=='#' || isspace(*z
) )
211 z
--; /* hit EOL :/ go back one char */
212 value_read(*(p
++),c
);
213 /* go point to the next non space char */
220 } /* Matrix_Read_Input */
224 * Read the contents of the matrix 'Mat' from standard input.
225 * a '#' is a comment (till EOL)
227 Matrix
*Matrix_Read(void) {
230 unsigned NbRows
, NbColumns
;
233 if (fgets(s
, 1024, stdin
) == NULL
)
235 while ((*s
=='#' || *s
=='\n') ||
236 (sscanf(s
, "%d %d", &NbRows
, &NbColumns
)<2)) {
237 if (fgets(s
, 1024, stdin
) == NULL
)
240 Mat
= Matrix_Alloc(NbRows
,NbColumns
);
242 errormsg1("Matrix_Read", "outofmem", "out of memory space");
245 if( !Matrix_Read_Input(Mat
) )
254 * Basic hermite engine
256 static int hermite(Matrix
*H
,Matrix
*U
,Matrix
*Q
) {
258 int nc
, nr
, i
, j
, k
, rank
, reduced
, pivotrow
;
260 Value
*temp1
, *temp2
;
263 /* Computes form: A = Q H and U A = H and U = Q */
266 errormsg1("Domlib", "nullH", "hermite: ? Null H");
271 temp1
= (Value
*) malloc(nc
* sizeof(Value
));
272 temp2
= (Value
*) malloc(nr
* sizeof(Value
));
273 if (!temp1
||!temp2
) {
274 errormsg1("Domlib", "outofmem", "out of memory space");
278 /* Initialize all the 'Value' variables */
279 value_init(pivot
); value_init(x
);
282 value_init(temp1
[i
]);
284 value_init(temp2
[i
]);
287 fprintf(stderr
,"Start -----------\n");
288 Matrix_Print(stderr
,0,H
);
290 for (k
=0, rank
=0; k
<nc
&& rank
<nr
; k
=k
+1) {
291 reduced
= 1; /* go through loop the first time */
293 fprintf(stderr
, "Working on col %d. Rank=%d ----------\n", k
+1, rank
+1);
298 /* 1. find pivot row */
299 value_absolute(pivot
,H
->p
[rank
][k
]);
301 /* the kth-diagonal element */
304 /* find the row i>rank with smallest nonzero element in col k */
305 for (i
=rank
+1; i
<nr
; i
++) {
306 value_absolute(x
,H
->p
[i
][k
]);
307 if (value_notzero_p(x
) &&
308 (value_lt(x
,pivot
) || value_zero_p(pivot
))) {
309 value_assign(pivot
,x
);
314 /* 2. Bring pivot to diagonal (exchange rows pivotrow and rank) */
315 if (pivotrow
!= rank
) {
316 Vector_Exchange(H
->p
[pivotrow
],H
->p
[rank
],nc
);
318 Vector_Exchange(U
->p
[pivotrow
],U
->p
[rank
],nr
);
320 Vector_Exchange(Q
->p
[pivotrow
],Q
->p
[rank
],nr
);
323 fprintf(stderr
,"Exchange rows %d and %d -----------\n", rank
+1, pivotrow
+1);
324 Matrix_Print(stderr
,0,H
);
327 value_assign(pivot
,H
->p
[rank
][k
]); /* actual ( no abs() ) pivot */
329 /* 3. Invert the row 'rank' if pivot is negative */
330 if (value_neg_p(pivot
)) {
331 value_oppose(pivot
,pivot
); /* pivot = -pivot */
333 value_oppose(H
->p
[rank
][j
],H
->p
[rank
][j
]);
335 /* H->p[rank][j] = -(H->p[rank][j]); */
338 value_oppose(U
->p
[rank
][j
],U
->p
[rank
][j
]);
340 /* U->p[rank][j] = -(U->p[rank][j]); */
343 value_oppose(Q
->p
[rank
][j
],Q
->p
[rank
][j
]);
345 /* Q->p[rank][j] = -(Q->p[rank][j]); */
347 fprintf(stderr
,"Negate row %d -----------\n", rank
+1);
348 Matrix_Print(stderr
,0,H
);
352 if (value_notzero_p(pivot
)) {
354 /* 4. Reduce the column modulo the pivot */
355 /* This eventually zeros out everything below the */
356 /* diagonal and produces an upper triangular matrix */
358 for (i
=rank
+1;i
<nr
;i
++) {
359 value_assign(x
,H
->p
[i
][k
]);
360 if (value_notzero_p(x
)) {
361 value_modulus(aux
,x
,pivot
);
363 /* floor[integer division] (corrected for neg x) */
364 if (value_neg_p(x
) && value_notzero_p(aux
)) {
367 value_division(x
,x
,pivot
);
368 value_decrement(x
,x
);
371 value_division(x
,x
,pivot
);
372 for (j
=0; j
<nc
; j
++) {
373 value_multiply(aux
,x
,H
->p
[rank
][j
]);
374 value_subtract(H
->p
[i
][j
],H
->p
[i
][j
],aux
);
377 /* U->p[i][j] -= (x * U->p[rank][j]); */
379 for (j
=0; j
<nr
; j
++) {
380 value_multiply(aux
,x
,U
->p
[rank
][j
]);
381 value_subtract(U
->p
[i
][j
],U
->p
[i
][j
],aux
);
384 /* Q->p[rank][j] += (x * Q->p[i][j]); */
387 value_addmul(Q
->p
[rank
][j
], x
, Q
->p
[i
][j
]);
393 "row %d = row %d - %d row %d -----------\n", i
+1, i
+1, x
, rank
+1);
394 Matrix_Print(stderr
,0,H
);
399 } /* if (pivot != 0) */
400 } /* while (reduced) */
402 /* Last finish up this column */
403 /* 5. Make pivot column positive (above pivot row) */
404 /* x should be zero for i>k */
406 if (value_notzero_p(pivot
)) {
407 for (i
=0; i
<rank
; i
++) {
408 value_assign(x
,H
->p
[i
][k
]);
409 if (value_notzero_p(x
)) {
410 value_modulus(aux
,x
,pivot
);
412 /* floor[integer division] (corrected for neg x) */
413 if (value_neg_p(x
) && value_notzero_p(aux
)) {
414 value_division(x
,x
,pivot
);
415 value_decrement(x
,x
);
420 value_division(x
,x
,pivot
);
422 /* H->p[i][j] -= x * H->p[rank][j]; */
423 for (j
=0; j
<nc
; j
++) {
424 value_multiply(aux
,x
,H
->p
[rank
][j
]);
425 value_subtract(H
->p
[i
][j
],H
->p
[i
][j
],aux
);
428 /* U->p[i][j] -= x * U->p[rank][j]; */
430 for (j
=0; j
<nr
; j
++) {
431 value_multiply(aux
,x
,U
->p
[rank
][j
]);
432 value_subtract(U
->p
[i
][j
],U
->p
[i
][j
],aux
);
435 /* Q->p[rank][j] += x * Q->p[i][j]; */
437 for (j
=0; j
<nr
; j
++) {
438 value_addmul(Q
->p
[rank
][j
], x
, Q
->p
[i
][j
]);
442 "row %d = row %d - %d row %d -----------\n", i
+1, i
+1, x
, rank
+1);
443 Matrix_Print(stderr
,0,H
);
448 } /* if (pivot!=0) */
451 /* Clear all the 'Value' variables */
452 value_clear(pivot
); value_clear(x
);
455 value_clear(temp1
[i
]);
457 value_clear(temp2
[i
]);
463 void right_hermite(Matrix
*A
,Matrix
**Hp
,Matrix
**Up
,Matrix
**Qp
) {
469 /* Computes form: A = QH , UA = H */
474 *Hp
= H
= Matrix_Alloc(nr
,nc
);
476 errormsg1("DomRightHermite", "outofmem", "out of memory space");
480 /* Initialize all the 'Value' variables */
483 Vector_Copy(A
->p_Init
,H
->p_Init
,nr
*nc
);
487 *Up
= U
= Matrix_Alloc(nr
, nr
);
489 errormsg1("DomRightHermite", "outofmem", "out of memory space");
493 Vector_Set(U
->p_Init
,0,nr
*nr
); /* zero's */
494 for(i
=0;i
<nr
;i
++) /* with diagonal of 1's */
495 value_set_si(U
->p
[i
][i
],1);
501 /* Actually I compute Q transpose... its easier */
503 *Qp
= Q
= Matrix_Alloc(nr
,nr
);
505 errormsg1("DomRightHermite", "outofmem", "out of memory space");
509 Vector_Set(Q
->p_Init
,0,nr
*nr
); /* zero's */
510 for (i
=0;i
<nr
;i
++) /* with diagonal of 1's */
511 value_set_si(Q
->p
[i
][i
],1);
518 /* Q is returned transposed */
521 for (i
=0; i
<nr
; i
++) {
522 for (j
=i
+1; j
<nr
; j
++) {
523 value_assign(tmp
,Q
->p
[i
][j
]);
524 value_assign(Q
->p
[i
][j
],Q
->p
[j
][i
] );
525 value_assign(Q
->p
[j
][i
],tmp
);
531 } /* right_hermite */
533 void left_hermite(Matrix
*A
,Matrix
**Hp
,Matrix
**Qp
,Matrix
**Up
) {
535 Matrix
*H
, *HT
, *Q
, *U
;
539 /* Computes left form: A = HQ , AU = H ,
541 using right form A = Q H , U A = H */
546 /* HT = A transpose */
547 HT
= Matrix_Alloc(nc
, nr
);
549 errormsg1("DomLeftHermite", "outofmem", "out of memory space");
555 value_assign(HT
->p
[j
][i
],A
->p
[i
][j
]);
559 *Up
= U
= Matrix_Alloc(nc
,nc
);
561 errormsg1("DomLeftHermite", "outofmem", "out of memory space");
565 Vector_Set(U
->p_Init
,0,nc
*nc
); /* zero's */
566 for (i
=0;i
<nc
;i
++) /* with diagonal of 1's */
567 value_set_si(U
->p
[i
][i
],1);
573 *Qp
= Q
= Matrix_Alloc(nc
, nc
);
575 errormsg1("DomLeftHermite", "outofmem", "out of memory space");
579 Vector_Set(Q
->p_Init
,0,nc
*nc
); /* zero's */
580 for (i
=0;i
<nc
;i
++) /* with diagonal of 1's */
581 value_set_si(Q
->p
[i
][i
],1);
586 /* H = HT transpose */
587 *Hp
= H
= Matrix_Alloc(nr
,nc
);
589 errormsg1("DomLeftHermite", "outofmem", "out of memory space");
595 value_assign(H
->p
[i
][j
],HT
->p
[j
][i
]);
600 for (i
=0; i
<nc
; i
++) {
601 for (j
=i
+1; j
<nc
; j
++) {
602 value_assign(tmp
,U
->p
[i
][j
]);
603 value_assign(U
->p
[i
][j
],U
->p
[j
][i
] );
604 value_assign(U
->p
[j
][i
],tmp
);
612 * Given a integer matrix 'Mat'(k x k), compute its inverse rational matrix
613 * 'MatInv' k x (k+1). The last column of each row in matrix MatInv is used
614 * to store the common denominator of the entries in a row. The output is 1,
615 * if 'Mat' is non-singular (invertible), otherwise the output is 0. Note::
616 * (1) Matrix 'Mat' is modified during the inverse operation.
617 * (2) Matrix 'MatInv' must be preallocated before passing into this function.
619 int MatInverse(Matrix
*Mat
,Matrix
*MatInv
) {
625 if(Mat
->NbRows
!= Mat
->NbColumns
) {
626 fprintf(stderr
,"Trying to invert a non-square matrix !\n");
630 /* Initialize all the 'Value' variables */
631 value_init(x
); value_init(gcd
); value_init(piv
);
632 value_init(m1
); value_init(m2
);
636 /* Initialise MatInv */
637 Vector_Set(MatInv
->p
[0],0,k
*(k
+1));
639 /* Initialize 'MatInv' to Identity matrix form. Each diagonal entry is set*/
640 /* to 1. Last column of each row (denominator of each entry in a row) is */
643 value_set_si(MatInv
->p
[i
][i
],1);
644 value_set_si(MatInv
->p
[i
][k
],1); /* denum */
646 /* Apply Gauss-Jordan elimination method on the two matrices 'Mat' and */
647 /* 'MatInv' in parallel. */
650 /* Check if the diagonal entry (new pivot) is non-zero or not */
651 if(value_zero_p(Mat
->p
[i
][i
])) {
653 /* Search for a non-zero pivot down the column(i) */
655 if(value_notzero_p(Mat
->p
[j
][i
]))
658 /* If no non-zero pivot is found, the matrix 'Mat' is non-invertible */
662 /* Clear all the 'Value' variables */
663 value_clear(x
); value_clear(gcd
); value_clear(piv
);
664 value_clear(m1
); value_clear(m2
);
668 /* Exchange the rows, row(i) and row(j) so that the diagonal element */
669 /* Mat->p[i][i] (pivot) is non-zero. Repeat the same operations on */
670 /* matrix 'MatInv'. */
673 /* Interchange rows, row(i) and row(j) of matrix 'Mat' */
674 value_assign(x
,Mat
->p
[j
][c
]);
675 value_assign(Mat
->p
[j
][c
],Mat
->p
[i
][c
]);
676 value_assign(Mat
->p
[i
][c
],x
);
678 /* Interchange rows, row(i) and row(j) of matrix 'MatInv' */
679 value_assign(x
,MatInv
->p
[j
][c
]);
680 value_assign(MatInv
->p
[j
][c
],MatInv
->p
[i
][c
]);
681 value_assign(MatInv
->p
[i
][c
],x
);
685 /* Make all the entries in column(i) of matrix 'Mat' zero except the */
686 /* diagonal entry. Repeat the same sequence of operations on matrix */
689 if (j
==i
) continue; /* Skip the pivot */
690 value_assign(x
,Mat
->p
[j
][i
]);
691 if(value_notzero_p(x
)) {
692 value_assign(piv
,Mat
->p
[i
][i
]);
693 value_gcd(gcd
, x
, piv
);
694 if (value_notone_p(gcd
) ) {
695 value_divexact(x
, x
, gcd
);
696 value_divexact(piv
, piv
, gcd
);
698 for(c
=((j
>i
)?i
:0);c
<k
;++c
) {
699 value_multiply(m1
,piv
,Mat
->p
[j
][c
]);
700 value_multiply(m2
,x
,Mat
->p
[i
][c
]);
701 value_subtract(Mat
->p
[j
][c
],m1
,m2
);
704 value_multiply(m1
,piv
,MatInv
->p
[j
][c
]);
705 value_multiply(m2
,x
,MatInv
->p
[i
][c
]);
706 value_subtract(MatInv
->p
[j
][c
],m1
,m2
);
709 /* Simplify row(j) of the two matrices 'Mat' and 'MatInv' by */
710 /* dividing the rows with the common GCD. */
711 Vector_Gcd(&MatInv
->p
[j
][0],k
,&m1
);
712 Vector_Gcd(&Mat
->p
[j
][0],k
,&m2
);
713 value_gcd(gcd
, m1
, m2
);
714 if(value_notone_p(gcd
)) {
716 value_divexact(Mat
->p
[j
][c
], Mat
->p
[j
][c
], gcd
);
717 value_divexact(MatInv
->p
[j
][c
], MatInv
->p
[j
][c
], gcd
);
724 /* Simplify every row so that 'Mat' reduces to Identity matrix. Perform */
725 /* the same sequence of operations on the matrix 'MatInv'. */
727 value_assign(MatInv
->p
[j
][k
],Mat
->p
[j
][j
]);
729 /* Make the last column (denominator of each entry) of every row greater */
731 Vector_Normalize_Positive(&MatInv
->p
[j
][0],k
+1,k
);
734 /* Clear all the 'Value' variables */
735 value_clear(x
); value_clear(gcd
); value_clear(piv
);
736 value_clear(m1
); value_clear(m2
);
742 * Given (m x n) integer matrix 'X' and n x (k+1) rational matrix 'P', compute
743 * the rational m x (k+1) rational matrix 'S'. The last column in each row of
744 * the rational matrices is used to store the common denominator of elements
747 void rat_prodmat(Matrix
*S
,Matrix
*X
,Matrix
*P
) {
750 int last_column_index
= P
->NbColumns
- 1;
751 Value lcm
, gcd
,last_column_entry
,s1
;
754 /* Initialize all the 'Value' variables */
755 value_init(lcm
); value_init(gcd
);
756 value_init(last_column_entry
); value_init(s1
);
757 value_init(m1
); value_init(m2
);
759 /* Compute the LCM of last column entries (denominators) of rows */
760 value_assign(lcm
,P
->p
[0][last_column_index
]);
761 for(k
=1;k
<P
->NbRows
;++k
) {
762 value_assign(last_column_entry
,P
->p
[k
][last_column_index
]);
763 value_gcd(gcd
, lcm
, last_column_entry
);
764 value_divexact(m1
, last_column_entry
, gcd
);
765 value_multiply(lcm
,lcm
,m1
);
768 /* S[i][j] = Sum(X[i][k] * P[k][j] where Sum is extended over k = 1..nbrows*/
769 for(i
=0;i
<X
->NbRows
;++i
)
770 for(j
=0;j
<P
->NbColumns
-1;++j
) {
772 /* Initialize s1 to zero. */
774 for(k
=0;k
<P
->NbRows
;++k
) {
776 /* If the LCM of last column entries is one, simply add the products */
777 if(value_one_p(lcm
)) {
778 value_addmul(s1
, X
->p
[i
][k
], P
->p
[k
][j
]);
781 /* Numerator (num) and denominator (denom) of S[i][j] is given by :- */
782 /* numerator = Sum(X[i][k]*P[k][j]*lcm/P[k][last_column_index]) and */
783 /* denominator= lcm where Sum is extended over k = 1..nbrows. */
785 value_multiply(m1
,X
->p
[i
][k
],P
->p
[k
][j
]);
786 value_division(m2
,lcm
,P
->p
[k
][last_column_index
]);
787 value_addmul(s1
, m1
, m2
);
790 value_assign(S
->p
[i
][j
],s1
);
793 for(i
=0;i
<S
->NbRows
;++i
) {
794 value_assign(S
->p
[i
][last_column_index
],lcm
);
796 /* Normalize the rows so that last element >=0 */
797 Vector_Normalize_Positive(&S
->p
[i
][0],S
->NbColumns
,S
->NbColumns
-1);
800 /* Clear all the 'Value' variables */
801 value_clear(lcm
); value_clear(gcd
);
802 value_clear(last_column_entry
); value_clear(s1
);
803 value_clear(m1
); value_clear(m2
);
809 * Given a matrix 'Mat' and vector 'p1', compute the matrix-vector product
810 * and store the result in vector 'p2'.
812 void Matrix_Vector_Product(Matrix
*Mat
,Value
*p1
,Value
*p2
) {
814 int NbRows
, NbColumns
, i
, j
;
815 Value
**cm
, *q
, *cp1
, *cp2
;
818 NbColumns
=Mat
->NbColumns
;
822 for(i
=0;i
<NbRows
;i
++) {
825 value_multiply(*cp2
,*q
,*cp1
);
829 /* *cp2 = *q++ * *cp1++ */
830 for(j
=1;j
<NbColumns
;j
++) {
831 value_addmul(*cp2
, *q
, *cp1
);
838 } /* Matrix_Vector_Product */
841 * Given a vector 'p1' and a matrix 'Mat', compute the vector-matrix product
842 * and store the result in vector 'p2'
844 void Vector_Matrix_Product(Value
*p1
,Matrix
*Mat
,Value
*p2
) {
846 int NbRows
, NbColumns
, i
, j
;
847 Value
**cm
, *cp1
, *cp2
;
850 NbColumns
=Mat
->NbColumns
;
853 for (j
=0;j
<NbColumns
;j
++) {
855 value_multiply(*cp2
,*(*cm
+j
),*cp1
);
858 /* *cp2= *(*cm+j) * *cp1++; */
859 for (i
=1;i
<NbRows
;i
++) {
860 value_addmul(*cp2
, *(*(cm
+i
)+j
), *cp1
);
866 } /* Vector_Matrix_Product */
869 * Given matrices 'Mat1' and 'Mat2', compute the matrix product and store in
872 void Matrix_Product(Matrix
*Mat1
,Matrix
*Mat2
,Matrix
*Mat3
) {
875 unsigned NbRows
, NbColumns
;
876 Value
**q1
, **q2
, *p1
, *p3
,sum
;
878 NbRows
= Mat1
->NbRows
;
879 NbColumns
= Mat2
->NbColumns
;
881 Size
= Mat1
->NbColumns
;
882 if(Mat2
->NbRows
!=Size
||Mat3
->NbRows
!=NbRows
||Mat3
->NbColumns
!=NbColumns
) {
883 fprintf(stderr
, "? Matrix_Product : incompatible matrix dimension\n");
891 /* Mat3[i][j] = Sum(Mat1[i][k]*Mat2[k][j] where sum is over k = 1..nbrows */
892 for (i
=0;i
<NbRows
;i
++) {
893 for (j
=0;j
<NbColumns
;j
++) {
896 for (k
=0;k
<Size
;k
++) {
897 value_addmul(sum
, *p1
, *(*(q2
+k
)+j
));
900 value_assign(*p3
,sum
);
906 } /* Matrix_Product */
909 * Given a rational matrix 'Mat'(k x k), compute its inverse rational matrix
912 * if 'Mat' is non-singular (invertible), otherwise the output is 0. Note::
913 * (1) Matrix 'Mat' is modified during the inverse operation.
914 * (2) Matrix 'MatInv' must be preallocated before passing into this function.
916 int Matrix_Inverse(Matrix
*Mat
,Matrix
*MatInv
) {
923 if(Mat
->NbRows
!= Mat
->NbColumns
) {
924 fprintf(stderr
,"Trying to invert a non-square matrix !\n");
928 /* Initialize all the 'Value' variables */
929 value_init(x
); value_init(gcd
); value_init(piv
);
930 value_init(m1
); value_init(m2
);
934 /* Initialise MatInv */
935 Vector_Set(MatInv
->p
[0],0,k
*k
);
937 /* Initialize 'MatInv' to Identity matrix form. Each diagonal entry is set*/
938 /* to 1. Last column of each row (denominator of each entry in a row) is */
941 value_set_si(MatInv
->p
[i
][i
],1);
942 /* value_set_si(MatInv->p[i][k],1); // denum */
944 /* Apply Gauss-Jordan elimination method on the two matrices 'Mat' and */
945 /* 'MatInv' in parallel. */
948 /* Check if the diagonal entry (new pivot) is non-zero or not */
949 if(value_zero_p(Mat
->p
[i
][i
])) {
951 /* Search for a non-zero pivot down the column(i) */
953 if(value_notzero_p(Mat
->p
[j
][i
]))
956 /* If no non-zero pivot is found, the matrix 'Mat' is non-invertible */
960 /* Clear all the 'Value' variables */
961 value_clear(x
); value_clear(gcd
); value_clear(piv
);
962 value_clear(m1
); value_clear(m2
);
966 /* Exchange the rows, row(i) and row(j) so that the diagonal element */
967 /* Mat->p[i][i] (pivot) is non-zero. Repeat the same operations on */
968 /* matrix 'MatInv'. */
971 /* Interchange rows, row(i) and row(j) of matrix 'Mat' */
972 value_assign(x
,Mat
->p
[j
][c
]);
973 value_assign(Mat
->p
[j
][c
],Mat
->p
[i
][c
]);
974 value_assign(Mat
->p
[i
][c
],x
);
976 /* Interchange rows, row(i) and row(j) of matrix 'MatInv' */
977 value_assign(x
,MatInv
->p
[j
][c
]);
978 value_assign(MatInv
->p
[j
][c
],MatInv
->p
[i
][c
]);
979 value_assign(MatInv
->p
[i
][c
],x
);
983 /* Make all the entries in column(i) of matrix 'Mat' zero except the */
984 /* diagonal entry. Repeat the same sequence of operations on matrix */
987 if (j
==i
) continue; /* Skip the pivot */
988 value_assign(x
,Mat
->p
[j
][i
]);
989 if(value_notzero_p(x
)) {
990 value_assign(piv
,Mat
->p
[i
][i
]);
991 value_gcd(gcd
, x
, piv
);
992 if (value_notone_p(gcd
) ) {
993 value_divexact(x
, x
, gcd
);
994 value_divexact(piv
, piv
, gcd
);
996 for(c
=((j
>i
)?i
:0);c
<k
;++c
) {
997 value_multiply(m1
,piv
,Mat
->p
[j
][c
]);
998 value_multiply(m2
,x
,Mat
->p
[i
][c
]);
999 value_subtract(Mat
->p
[j
][c
],m1
,m2
);
1002 value_multiply(m1
,piv
,MatInv
->p
[j
][c
]);
1003 value_multiply(m2
,x
,MatInv
->p
[i
][c
]);
1004 value_subtract(MatInv
->p
[j
][c
],m1
,m2
);
1007 /* Simplify row(j) of the two matrices 'Mat' and 'MatInv' by */
1008 /* dividing the rows with the common GCD. */
1009 Vector_Gcd(&MatInv
->p
[j
][0],k
,&m1
);
1010 Vector_Gcd(&Mat
->p
[j
][0],k
,&m2
);
1011 value_gcd(gcd
, m1
, m2
);
1012 if(value_notone_p(gcd
)) {
1014 value_divexact(Mat
->p
[j
][c
], Mat
->p
[j
][c
], gcd
);
1015 value_divexact(MatInv
->p
[j
][c
], MatInv
->p
[j
][c
], gcd
);
1022 /* Find common denom for each row */
1023 den
= (Value
*)malloc(k
*sizeof(Value
));
1025 for(j
=0 ; j
<k
; ++j
) {
1027 value_assign(den
[j
],Mat
->p
[j
][j
]);
1029 /* gcd is always positive */
1030 Vector_Gcd(&MatInv
->p
[j
][0],k
,&gcd
);
1031 value_gcd(gcd
, gcd
, den
[j
]);
1032 if (value_neg_p(den
[j
]))
1033 value_oppose(gcd
,gcd
); /* make denominator positive */
1034 if (value_notone_p(gcd
)) {
1036 value_divexact(MatInv
->p
[j
][c
], MatInv
->p
[j
][c
], gcd
); /* normalize */
1037 value_divexact(den
[j
], den
[j
], gcd
);
1039 value_gcd(gcd
, x
, den
[j
]);
1040 value_divexact(m1
, den
[j
], gcd
);
1041 value_multiply(x
,x
,m1
);
1043 if (value_notone_p(x
))
1044 for(j
=0 ; j
<k
; ++j
) {
1045 for (c
=0; c
<k
; c
++) {
1046 value_division(m1
,x
,den
[j
]);
1047 value_multiply(MatInv
->p
[j
][c
],MatInv
->p
[j
][c
],m1
); /* normalize */
1051 /* Clear all the 'Value' variables */
1052 for(j
=0 ; j
<k
; ++j
) {
1053 value_clear(den
[j
]);
1055 value_clear(x
); value_clear(gcd
); value_clear(piv
);
1056 value_clear(m1
); value_clear(m2
);
1060 } /* Matrix_Inverse */