/*
$Id$
Copyright 1989-2016 MINES ParisTech
This file is part of Linear/C3 Library.
Linear/C3 Library is free software: you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
any later version.
Linear/C3 Library is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.
See the GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with Linear/C3 Library. If not, see .
*/
/* package matrix */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include
#include
#include "linear_assert.h"
#include "boolean.h"
#include "arithmetique.h"
#include "matrix.h"
/* void matrix_unimodular_triangular_inversion(Pmatrix u ,Pmatrix inv_u,
* * bool infer)
* u soit le matrice unimodulaire triangulaire.
* si infer = true (triangulaire inferieure),
* infer = false (triangulaire superieure).
* calcul de l'inversion de matrice u telle que I = U x INV_U .
* Les parametres de la fonction :
*
* Pmatrix u : matrice unimodulaire triangulaire -- inpout
* int n : dimension de la matrice caree -- inpout
* bool infer : type de triangulaire -- input
* matrice inv_u : l'inversion de matrice u -- output
*/
void matrix_unimodular_triangular_inversion(u,inv_u,infer)
Pmatrix u;
Pmatrix inv_u;
bool infer;
{
int i, j;
Value x;
int n = MATRIX_NB_LINES(u);
/* test de l'unimodularite et de la trangularite de u */
assert(matrix_triangular_unimodular_p(u,infer));
matrix_identity(inv_u,0);
if (infer){
for (i=n; i>=1;i--)
for (j=i-1; j>=1; j--){
x = MATRIX_ELEM(u,i,j);
if (value_notzero_p(x))
matrix_subtraction_column(inv_u,j,i,x);
}
}
else{
for (i=1; i<=n; i++)
for(j=i+1; j<=n; j++){
x = MATRIX_ELEM(u,i,j);
if (value_notzero_p(x))
matrix_subtraction_column(inv_u,j,i,x);
}
}
}
/* void matrix_diagonal_inversion(Pmatrix s,Pmatrix inv_s)
* calcul de l'inversion du matrice en forme reduite smith.
* s est un matrice de la forme reduite smith, inv_s est l'inversion
* de s ; telles que s * inv_s = I.
*
* les parametres de la fonction :
* matrice s : matrice en forme reduite smith -- input
* int n : dimension de la matrice caree -- input
* matrice inv_s : l'inversion de s -- output
*/
void matrix_diagonal_inversion(s,inv_s)
Pmatrix s;
Pmatrix inv_s;
{
Value d, d1;
Value gcd, lcm;
int i;
int n = MATRIX_NB_LINES(s);
/* tests des preconditions */
assert(matrix_diagonal_p(s));
assert(matrix_hermite_rank(s)==n);
assert(value_pos_p(MATRIX_DENOMINATOR(s)));
matrix_nulle(inv_s);
/* calcul de la ppcm(s[1,1],s[2,2],...s[n,n]) */
lcm = MATRIX_ELEM(s,1,1);
for (i=2; i<=n; i++)
lcm = ppcm(lcm,MATRIX_ELEM(s,i,i));
d = MATRIX_DENOMINATOR(s);
gcd = pgcd(lcm,d);
d1 = value_div(d,gcd);
for (i=1; i<=n; i++) {
Value tmp = value_div(lcm,MATRIX_ELEM(s,i,i));
MATRIX_ELEM(inv_s,i,i) = value_mult(d1,tmp);
}
MATRIX_DENOMINATOR(inv_s) = value_div(lcm,gcd);
}
/* void matrix_triangular_inversion(Pmatrix h, Pmatrix inv_h,bool infer)
* calcul de l'inversion du matrice en forme triangulaire.
* soit h matrice de la reduite triangulaire; inv_h est l'inversion de
* h ; telle que : h * inv_h = I.
* selon les proprietes de la matrice triangulaire:
* Aii = a11* ...aii-1*aii+1...*ann;
* Aij = 0 i>j pour la matrice triangulaire inferieure (infer==true)
* i