HP Forums
LUP matrix decomposition - Printable Version

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: LUP matrix decomposition (/thread-22885.html)



LUP matrix decomposition - Pekis - 12-21-2024 08:58 PM

Hello,

In advanced functions handbook of the HP-15C, I read that it uses the LU matrix decomposition, with pivoting (P).

I understand the advantages of pivoting and I often saw on Internet the method of selecting the row (in the remaining rows) with the biggest absolute element in the current column and swap current row with that one. Like in:

Code:
for_i  {
        int max_j = i;
        foreach(j, i, n)
            if (fabs(a[j][i]) > fabs(a[max_j][i])) max_j = j;

        if (max_j != i)
            for_k { _swap(p[i][k], p[max_j][k]); }
    }
...

But I also found sources proposing another method for selecting pivot and I don't understand quite well why it's done that way, with a dot product sum very similar to getting the "u"s and the "L"s of the LU matrix:

Code:

  for (i=0; i<N; i++) { 

    //start pivot section
    Umax = 0;
    for (r=i; r<N; r++) {
      Uii=A[r][i];
      q = 0;
      while (q<i) {
        Uii -= A[r][q]*A[q][r];
        q++;
      } 
      if (Math.abs(Uii)>Umax) {
        Umax = Math.abs(Uii);
        row = r;
      }
    }
    if (i!=row) {//swap rows
      exchanges++;
      for (q=0; q<N; q++) {
        tmp = P[i][q];
        P[i][q]=P[row][q];
        P[row][q]=tmp;
        tmp = A[i][q];
        A[i][q]=A[row][q];
        A[row][q]=tmp;
      } 
    }
   //end pivot section

    j = i;
    while (j<N) { //Determine U across row i
      q = 0;
      while (q<i) {
        A[i][j] -= A[i][q]*A[q][j];
        q++;
      } 
      j++;
    }
    j = i+1;
    while (j<N) { //Determine L down column i
      q = 0;
      while (q<i) {
        A[j][i] -= A[j][q]*A[q][i];
        q++;
      } 
      A[j][i] = A[j][i]/A[i][i];
      j++;
    } 
  }
  return {LU:A,P:P,exchanges:exchanges};
...

I hope you see what I mean ... Did you know this method ?

Thanks