Algorithme résolution système linéaire

Question

Bonjour à tous,

Je recherche un algorithme de résolution direct me permettant de résoudre un système linéaire AX=B avec A symétrique, de largeur de bande connue.
A est n'est pas stocké dans un tableau de taille N*N mais je stocke uniquement les termes du triangle supérieur de la matrice dans un tableau de taille N*LB. LB est la largeur de bande de ma matrice. Ceci me permet de gagner énormément de mémoire vive.

Merci à celui ou celle qui pourrait m'aider.

DD

ricky78 · Answer

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(...)

rad_issam · Answer

Algorithme résolution système linéaire

2 réponses

Votre réponse

Discussions similaires