Implementing LU factorization with partial pivoting in C using only one matrix

Question

I have designed the following C function in order to compute the PA = LU factorization, using only one matrix to store and compute the data:

double plupmc(int n, double **c, int *p, double tol) {
  int i, j, k, pivot_ind = 0, temp_ind;
  int ii, jj;
  double pivot, *temp_row;

  for (j = 0; j < n-1; ++j) {

    pivot = 0.;
    for (i = j; i < n; ++i)
      if (fabs(c[i][j]) > fabs(pivot)) {
        pivot = c[i][j];
        pivot_ind = i;
      }

    temp_row = c[j];
    c[j] = c[pivot_ind];
    c[pivot_ind] = temp_row;

    temp_ind  = p[j];
    p[j] = p[pivot_ind];
    p[pivot_ind] = temp_ind;

    for (k = j+1; k < n; ++k) {
      c[k][j] /= c[j][j];
      c[k][k] -= c[k][j]*c[j][k];
    }

  }
  return 0.;
}

where n is the order of the matrix, c is a pointer to the matrix and p is a pointer to a vector storing the permutations done when partial pivoting the system. The variable tol is not relevant for now. The program works storing in c both the lower and upper triangular parts of the factorization, where U corresponds to the upper triangular part of c and L corresponds to the strictly lower triangular part of c, adding 1's in the diagonal. For what I have been able to test, the part of the program corresponding to partial pivoting is working properly, however, the algorithm used to compute the entries of the matrix is not giving the expected results, and I cannot see why. For instance, if I try to compute the LU factorization of the matrix

1. 2. 3.
4. 5. 6.
7. 8. 9.

I get

    1.    0.     0.         7. 8. 9.
l : 0.143 1.     0.     u : 0. 2. 1.714*
    0.571 0.214* 1.         0. 0. 5.663*

the product of which does not correspond to any permutation of the matrix c. In fact, the wrong entries seem to be the ones marked with a star.

I would appreciate any suggestion to fix this problem.

You need to swap the whole rows !

user1196549
– user1196549

2021-10-12 07:04:26 +00:00
Commented Oct 12, 2021 at 7:04 — user1196549
– user1196549, Commented Oct 12, 2021 at 7:04

marc_s · Accepted Answer · 2021-12-02 22:02:54Z

I found the problem with your code, there was a little conceptual error in the way that you were normalizing the row while computing the actual decomposition:

for (k = j+1; k < n; ++k) {
  c[k][j] /= c[j][j];
  c[k][k] -= c[k][j]*c[j][k];
}

became:

for (k = j+1; k < n; ++k) {
  temp=c[k][j]/=c[j][j];
  for(int q=j+1;q<n;q++){
        c[k][q] -= temp*c[j][q];
      }
}

which returns the result:

7.000000 8.000000 9.000000
0.142857 0.857143 1.714286
0.571429 0.500000 -0.000000

If you have any questions I am happy to help.

Full implementation here:

#include<stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

double plupmc(int n, double **c, int *p, double tol) {
  int i, j, k, pivot_ind = 0, temp_ind;
  int ii, jj;
  double *vv=calloc(n,sizeof(double));
  double pivot, *temp_row;
  double temp;

  for (j = 0; j < n; ++j) {
    pivot = 0;
    for (i = j; i < n; ++i)
      if (fabs(c[i][j]) > fabs(pivot)) {
        pivot = c[i][j];
        pivot_ind = i;
      }

    temp_row = c[j];
    c[j] = c[pivot_ind];
    c[pivot_ind] = temp_row;

    temp_ind  = p[j];
    p[j] = p[pivot_ind];
    p[pivot_ind] = temp_ind;

    for (k = j+1; k < n; ++k) {
      temp=c[k][j]/=c[j][j];
      for(int q=j+1;q<n;q++){
            c[k][q] -= temp*c[j][q];
          }
    }
    for(int q=0;q<n;q++){
      for(int l=0;l<n;l++){
        printf("%lf ",c[q][l]);
      }
      printf("\n");
    }

  }
  return 0.;
}

int main() {
  double **x;
  x=calloc(3,sizeof(double));
  for(int i=0;i<3;i++){
    x[i]=calloc(3,sizeof(double));
  }
  memcpy(x[0],(double[]){1,2,3},3*sizeof(double));
  memcpy(x[1],(double[]){4,5,6},3*sizeof(double));
  memcpy(x[2],(double[]){7,8,9},3*sizeof(double));

  int *p=calloc(3,sizeof(int));
  memcpy(p,(int[]){0,1,2},3*sizeof(int));
  plupmc(3,x,p,1);
  for(int i=0;i<3;i++){
    free(x[i]);
  }
  free(p);
  free(x);


}

Collectives™ on Stack Overflow

Implementing LU factorization with partial pivoting in C using only one matrix

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related