Chapter 10  Algèbre linéaire

10.1  Produit matrice-vecteur sur anneau

for (i=1;i<=n;i++)
  for (j=1;j<=n;j++)
    y[i] = y[i]+a[i,j]*x[j];

float a[r][n];
float x[r],y[r];

matrice-vecteur(A,x,y) {
  q = my_num();
  p = tot_proc_num();
  for (step=0;step<p;step++) {
    send(x,r);
    for (i=0;i<r;i++) 
      for (j=0;j<r;j++)
        y[i] = y[i]+a[i,(q-step mod p)r+j]*x[j];
    receive(temp,r); 
    x = temp;
  }
}

10.2  Factorisation LU

for (k=0;k<n-1;k++) {
  prep(k): for (i=k+1;i<n;i++)
    a[i,k]=a[i,k]/a[k,k];
  for (j=k+1;j<n;j++)
    update(k,j): for (i=k+1;i<n;i++)
      a[i,j]=a[i,j]-a[i,k]*a[k,j];
}

q = my_num();
p = tot_proc_num();
for (k=0;k<n-1;k++) {
  if (k == q) { 
    prep(k): for (i=k+1;i<n;i++)
      buffer[i-k-1] = a[i,k]/a[k,k];
    broadcast(k,buffer,n-k); 
  }
  else { 
      receive(buffer,n-k);
      update(k,q): for (i=k+1;k<n;k++)
        a[i,q] = a[i,q]-buffer[i-k-1]*a[k,q]; }
}

q = my_num();
p = tot_proc_num();
l = 0;
for (k=0;k<n-1;k++) {
  if (alloc(k) == q) {
    prep(k): for (i=k+1;i<n;i++)
      buffer[i-k-1] = a[i,l]/a[k,l];
    l++; }
  broadcast(alloc(k),buffer,n-k);
  for (j=l;j<r;j++)
    update(k,j): for (i=k+1;k<n;k++)
      a[i,j] = a[i,j]-buffer[i-k]*a[k,j]; }

10.2.1  Cas de l'allocation cyclique par lignes

10.2.2  Recouvrement communication/calcul

q = my_num();
p = tot_proc_num();
l = 0;
for (k=0;k<n-1;k++) {
  if (k == q mod p) {
    prep(k): for (i=k+1;i<n;i++)
      buffer[i-k-1] = -a[i,l]/a[k,l];
    l++; send(buffer,n-k); }
  else { receive(buffer,n-k);
    if (q != k-1 mod p) send(buffer,n-k); }
  for (j=l;j<r;j++)
    update(k,j): for (i=k+1;k<n;k++)
      a[i,j] = a[i,j]+buffer[i-k-1]*a[k,j]; }

• A l'étape k=0: P₁ reçoit la colonne pivot 0 de P₀
• P₁ l'envoit à P₂
• Fait update(0,j) pour toutes les colonnes j qui lui appartiennent, c'est-à-dire j=1 mod p
• A l'étape k=1: fait prep(1)
• Envoie la colonne pivot 1 à P₂
• Fait update(1,j) pour toutes les colonnes j qui lui appartiennent, c'est-à-dire j=1 mod p
• etc.

• update(0,1)
• prep(1)
• Envoi vers P₂
• update(0,j) pour j=1 mod p et j>1
• etc.

10.3  Algorithmique sur grille 2D

10.3.1  Principe de l'algorithme de Cannon

/* diag(A) sur col 0, diag(B) sur ligne 0 */
Rotations(A,B); /* preskewing */

/* calcul du produit de matrice */
forall (k=1; k<=sqrt(P)) {
  LocalProdMat(A,B,C);
  VerticalRotation(B,downwards);
  HorizontalRotation(A,leftwards); }

/* mouvements des donnees apres les calculs */
Rotations(A,B); /* postskewing */

10.3.2  Principe de l'algorithme de Fox

/* pas de mouvements de donnees avant les calculs */

/* calcul du produit de matrices */
broadcast(A(x,x));
forall (k=1; k<sqrt(P)) {
  LocalProdMat(A,B,C);
  VerticalRotation(B,upwards);
  broadcast(A(k+x,k+x)); }
forall () {
  LocalProdMat(A,B,C);
  VerticalRotation(B,upwards); }

/* pas de mouvements de donnees apres les calculs */

10.3.3  Principe de l'algorithme de Snyder

/* mouvements des donnees avant les calculs */
Transpose(B);
/* calcul du produit de matrices */
forall () {
  LocalProdMat(A,B,C);
  VerticalRotation(B,upwards); }
forall (k=1;k<sqrt(P)) {
  GlobalSum(C(i,(i+k-1) mod sqrt(P)));
  LocalProdMat(A,B,C);
  VerticalRotation(B,upwards); }
GlobalSum(C(i,(i+sqrt(P)-1) mod sqrt(P)));
/* mouvements des donnees apres les calculs */
Transpose(B);

10.4  Algorithmique hétérogène

Distribute(B,t1,t2,...,tn)
/* initialisation: calcule ci */
for (i=1;i<=p;i++)
  c[i]=...
/* incrementer iterativement les ci minimisant le temps */
while (sum(c[])<B) {
  find k in {1,...,p} st t[k]*(c[k]+1)=min{t[i]*(c[i]+1)};
  c[k] = c[k]+1;
return(c[]);

Distribute(B,t1,t2,...tp)
/* Initialisation: aucune tache a distribuer m=0 */
for (i=1;i<=p;i++) c[i]=0;
/* construit iterativement l'allocation sigma */
for (m=1;m<=B;m++)
  find(k in {1,...p} st t[k]*(c[k]+1)=min{t[i]*(c[i]+1)});
  c[k]=c[k]+1;
  sigma[m]=k;
return(sigma[],c[]);

10.4.1   LU hétérogène (1D)