ricky78
|
bonjour DD05
Tu peux essayer de t'inspirer du code suivant en l'adaptant pour des matrices reelle ou l'utiliser directement en mettant la partie complexe à 0. Par contre il faut que tu modifie la forme de ta matrice qui ne peut rester triangulaire en utilisant directement ce code ou si tu veut la laisser triangulaire tu n'auras qu'à modifier le point de départ des boucles pour travailler uniquement avec la triangulaire haute.
cordialement
ricky78
//=========================================
// Division en complexe
// (qr + J*qi) = (ar + J*ai) / (br + J*bi)
// J = sqrt(-1)
//=========================================
void div_complex(float ar, float ai, float br, float bi, float *qr, float *qi)
// qr+J.qi = (ar+J.ai) / (br+J.bi) = (ar+J.ai).(br-J.bi) / (br+J.bi).(br-J.bi)
// = ((ar.br+ai.bi + J.(ai.br-ar.bi)) / (br² + bi²)
{
float m;
m = br*br + bi*bi;
*qr = (ar*br + ai*bi) / m;
*qi = (ai*br - ar*bi) / m;
}
//=========================================
// Multiplication en complexe
// (pr + J*pi) = (ar + J*ai) * (br + J*bi)
// J = sqrt(-1)
//=========================================
void mul_complex(float ar, float ai, float br, float bi, float *pr, float *pi)
// pr+J.pi = (ar+J.ai).(br+J.bi) = (ar.br - ai.bi) + J.(ar.bi + ai.br)
{
*pr = ar*br - ai*bi;
*pi = ar*bi + ai*br;
}
//===============================================
// Resolution d'un systeme d'equations lineaires
// en complexe
// |Y| = |A|.|X| par la methode de Gauss
//===============================================
void resol_equa_lin_complex(int n, float *Ma, float *Mx, float *My)
{
//-----------------------------------------------------
// Ma: matrice carree n*n dans un tableau (WW)*(MM)
// Ma[i, j] (i et j de 1 a n) a pour adresse (i*WW + 2.j)
// Mx[i] (i de 1 a n) a pour adresse (2.i)
// Exemple: MM=6 WW=2*MM=12 n=3
//
// j: 0 1 2 3 4 5
// i=0 - - - - - - - - - - - - i=0 - -
// i=1 - - A a A a A a - - - - i=1 X x
// i=2 - - A a A a A a - - - - i=2 X x
// i=3 - - A a A a A a - - - - i=3 X x
// i=4 - - - - - - - - - - - - i=4 - -
// i=5 - - - - - - - - - - - - i=5 - -
//
// Partie reelle: A X
// Partie imaginaire: a x
//-----------------------------------------------------
int i, j, k, m;
float sr, si, dpr, dpi, dr, di, pr, pi;
//......................................................
// Nota - Pour les commentaires:
// MA(y, x) = Ma[y*WW+2*x] + J*Ma[y*WW+2*x+1]
// MX(x) = Mx[2*x] + J*Mx[2*x+1]
// MY(x) = My[2*x] + J*My[2*x+1] J = sqrt(-1)
//......................................................
for (k = 1; k <= (n-1); k ++)
{
// Calcul de DP = 1 / MA(k, k)
div_complex(1.0, 0.0, Ma[k*WW+2*k], Ma[k*WW+2*k+1], &dpr, &dpi);
for (j = k; j <= n; j ++)
{
// MA(k,j) = MA(k,j).DP
mul_complex(Ma[k*WW+2*j], Ma[k*WW+2*j+1], dpr, dpi,
&Ma[k*WW+2*j], &Ma[k*WW+2*j+1]);
}
// MY(k) = MY(k).DP
mul_complex(My[2*k], My[2*k+1], dpr, dpi, &My[2*k], &My[2*k+1]);
for (j = k+1; j <= n; j ++)
{
// D = MA(j,k)
dr = Ma[j*WW+2*k];
di = Ma[j*WW+2*k+1];
for (m = k; m <= n; m ++)
{
// P = MA(k,m).D
mul_complex(Ma[k*WW+2*m], Ma[k*WW+2*m+1], dr, di, &pr, &pi);
// MA(j,m) = MA(j,m)-P
Ma[j*WW+2*m] -= pr;
Ma[j*WW+2*m+1] -= pi;
}
// P = MY(k).D
mul_complex(My[2*k], My[2*k+1], dr, di, &pr, &pi);
// MY(j) = MY(j)-P
My[2*j] -= pr;
My[2*j+1] -= pi;
}
}
// Solution
// MX(n) = MY(n)/MA(n,n)
div_complex(My[2*n], My[2*n+1], Ma[n*WW+2*n], Ma[n*WW+2*n+1],
&Mx[2*n], &Mx[2*n+1]);
for (i = 1; i <= (n-1); i ++)
{
// S = 0
sr = 0.0;
si = 0.0;
for (j = (n-i+1); j <= n; j ++)
{
// P = MA((n-i),j).MX(j)
mul_complex(Ma[(n-i)*WW+2*j], Ma[(n-i)*WW+2*j+1], Mx[2*j], Mx[2*j+1],
&pr, &pi);
// S = S+P
sr += pr;
si += pi;
}
// MX(n-i) = MY(n-i)-S
Mx[2*(n-i)] = My[2*(n-i)] - sr;
Mx[2*(n-i)+1] = My[2*(n-i)+1] - si;
}
}
// fin de resol_equa_lin_complex(...)
|