Voici ma fonction qui permet de calculer la pseudo inverse, le pointeur mat correspond a ma matrice de départ (elle est recu correctement par la fonction), le pointeur inv correspond a la matrice pseudo inversé, N et M sont la taille de la matrice (ou tableau). Le probleme c que je trouve pas les meme résultat que sous Matlab donc si vous voyez une erreur ds l'algo. Merci d'avance car je me doute que pas grand monde connait l'algo d'une pseudo-inverse
# include <math.h>
# include <stdio.h>
int Pseudo(double *Mat,double *Inv,int N,int M)
{
int i,j,i1,i2,k,k1 ;
double xmax,xnul,x,d2 ;
double B[4],C[4],D[4],xnorm[4] ;
double precision = 1e12;
/* Calcul de la norme des vecteurs colonne */
for(i=0 ; i<M ; i++)
{
xnorm[i]=0.0 ;
for(j=0 ; j<N ; j++) xnorm[i] += ( *(Mat+(j*M)+i) ) * ( *(Mat+(j*M)+i) ) ;
}
xmax = xnorm[0] ; /* Recherche de la colonne de norme maximale */
for(i=1 ; i<M ; i++)
if( xnorm[i]>xmax ) xmax = xnorm[i] ;
xnul = xmax/precision ;
/* Initialisation de la Pseudo-inverse */
/* La première ligne de Inv est la première colonne */
/* de Mat divise par sa norme (si non nul) */
if( xnorm[0]>xnul )
for(i=0 ; i<N ; i++) ( *(Inv+i) ) = ( *(Mat+(i*M)) )/xnorm[0] ;
else
for(i=0 ; i<N ; i++) ( *(Inv+i) ) = 0.0 ;
/* Calcul des autres lignes de la Pseudo-inverse */
for(k=1 ; k<M ; k++)
{
k1 = k-1 ; /* indice pour le partitionnement de la matrice */
for(i1=0 ; i1<=k1 ; i1++) /* Calcul de D */
{
D[i1] = 0.0 ;
for(i2=0 ; i2<M ; i2++)
D[i1] += ( *(Inv+(i1*N)+i2) ) * ( *(Mat+(i2*M)+k) ) ;
}
for(i2=0 ; i2<N ; i2++) /* Calcul de C */
{
C[i2] = 0.0 ;
for(i1=0 ; i1<=k1 ; i1++)
C[i2] += ( *(Mat+(i2*M)+i1) ) * D[i1] ;
C[i2] = ( *(Mat+(i2*M)+k) ) - C[i2] ;
}
/* Calcul de la norme de C */
d2 = 0.0 ;
for(i=0 ; i<N ; i++) d2 += C[i]*C[i] ;
if(d2>xnul) /* Si cette norme n'est pas nulle */
for(i=0 ; i<N ; i++)
{
B[i] = C[i]/d2 ;
}
else /* Si cette norme est nulle */
{
d2 = 0.0 ;
for(i=0 ; i<=k1 ; i++)
d2 += D[i]*D[i] ; /* Calcule de la norme de D */
x = d2 + 1.0 ; x= 1.0/x ;
for(i=0 ; i<N ; i++)
{
d2=0.0 ;
for(j=0 ; j<=k1 ; j++) d2 += D[j] * ( *(Inv+(j*N)+i) ) ;
B[i] = x*d2 ;
}
}
for(i=0 ; i<=k1 ; i++)
for(j=0 ; j<N ; j++)
( *(Inv+(i*N)+j) ) -= D[i]*B[j] ;
for(i=0 ; i<N ; i++) ( *(Inv+(k*N)+i) ) = B[i] ;
}
for (i=0;i<12;i++) printf ("%lf\n",*Inv++);
return(1) ; /* C'est fini */
}
