c++ - 包含 LU 分解的矩阵

Question

我正在做 LU decom，我在 googel 上找到了这段代码，但想通过输出“pvt”和“a”来理解它，但它似乎是我的 pvt 不正确，所以我得到了一些不同的东西所以请问谁能纠正我 .. 谢谢

这是我的代码

    int* LUfactor ( double **a, int n, int ps )


/*PURPOSE:   compute an LU decomposition for the coefficient matrix a

CALLING SEQUENCE:
       pvt = LUfactor ( a, n, ps );

INPUTS:
       a        coefficient matrix
                type:  **doble
       n        number of equations in system
                type:  int
       ps       flag indicating which pivoting strategy to use
                     ps == 0:   no pivoting
                     ps == 1;   partial pivoting
                     ps == 2;   scaled partial pivoting
                type:  int

  OUTPUT:
       pvt      vector which indicates the permutation of the rows
                performed during the decomposition process
                type:  *int
       a        matrix containing LU decomposition of the input    coefficient
                matrix - the L matrix in the decomposition consists of 1's
                along the main diagonal together with the strictly lower
                triangular portion of the output matrix a; the U matrix
                in the decomposition is theupper triangular portion of the
                output matrix a
                type:  **double
   */                 



 {
        int pass, row, col, *pvt, j, temp; 
    double *s,rmax,ftmp, mult, sum;

       /*initialize row pointer array*/


  pvt = new int [n];
  for ( row = 0; row < n; row++ )
      pvt[row] = row;


   /* if scaled partial pivoting option was selected,
   initialize scale vector*/


 if ( ps == 2 ) {
    s = new double [n];
    for ( row = 0; row < n; row++ ) {
        s[row] = fabs( a[row][0] );
        for ( col = 1; col < n; col++ )
            if ( fabs( a[row][col] ) > s[row] ) 
               s[row] = fabs( a[row][col] );
    }
 }


   /*elimination phase*/


 for ( pass = 0; pass < n; pass++ ) {


 /*  perform requested pivoting strategy

  even if no pivoting option is requested, still must check for
  zero pivot*/


     if ( ps != 0 ) {
        rmax = ( ps == 1 ? fabs( a[pvt[pass]][pass] ) : 
                           fabs( a[pvt[pass]][pass] ) / s[pvt[pass]] );
        j = pass;
        for ( row = pass+1; row < n; row++ ) { 
            ftmp = ( ps == 1 ? fabs( a[pvt[row]][pass] ) : 
                               fabs( a[pvt[row]][pass] ) / s[pvt[row]] );
            if ( ftmp > rmax ) {
               rmax = ftmp;
               j = row;
            }
        }

        if ( j != pass ) {
           temp = pvt[j];
           pvt[j] = pvt[pass];
           pvt[pass] = temp;
        }
     }
     else {
        if ( a[pvt[pass]][pass] == 0.0 ) {
           for ( row = pass+1; row < n; row++ )
               if ( a[pvt[row]][pass] != 0.0 ) break;
           temp = pvt[row];
           pvt[row] = pvt[pass];
           pvt[pass] = temp;  
        }
     }

     for ( row = pass + 1; row < n; row++ ) {
         mult = - a[pvt[row]][pass] / a[pvt[pass]][pass];
         a[pvt[row]][pass] = -mult;
         for ( col = pass+1; col < n; col++ )
             a[pvt[row]][col] += mult * a[pvt[pass]][col];
     }
 }

 if ( ps == 2 ) delete [] s;
 return ( pvt );
}

这是我的主要内容

     double **af;
 int *pvt;

     int i, j, n;

/*
  allocate space for coefficient matrix
           */

 n = 4;

     af = new double* [n];
 pvt = new int [n];

 for ( i = 0; i < n; i++ )
     af[i] = new double [n];

      af[0][0] = 2.00;    af[0][1] = 1.00;   af[0][2] = 1.00; af[0][3] = -2.00;
 af[1][0] = 4.00;    af[1][1] = 0.00;   af[1][2] = 2.00; af[1][3] = 1.00;
 af[2][0] = 3.00;    af[2][1] = 2.00;   af[2][2] = 2.00; af[2][3] = 0.00;
 af[3][0] = 1.00;    af[3][1] = 3.00;   af[3][2] = 2.00; af[3][3] = 0.00;

 pvt =LUfactor ( af, n, 0 );

    cout << "pvt" << endl;
      for ( i = 0; i < n; i++ )
          cout << pvt[i] << endl;
          cout << endl << endl <<  endl;


      cout << "a" << endl;
      for ( i = 0; i < n; i++ )
          cout << af[i][i] << endl;
          cout << endl << endl <<  endl;

/////// 
        out put

  pvt
0
3
1
2

 LU matrix is 
 2  1  1  -2  0  
 2  -0.8  1.2  5.8  0  
 1.5  0.2  0.166667  1.83333  0  
 0.5  2.5  1.5  1  0  
 Segmentation fault

////////////////////////////////////

The out put I'm looking for is 
  Matrix A
     0         2         0         1
     2         2         3         2
     4        -3         0         1
     6         1        -6        -5

determinant: -234
pivot vector:   3  2  1  0

 Lower triangular matrix
      6         0         0         0
      4    -3.667         0         0
      2     1.667     6.818         0
      0         2     2.182      1.56

Upper triangular matrix
     1    0.1667        -1   -0.8333
     0         1    -1.091    -1.182
     0         0         1    0.8267
     0         0         0         1


Product of L U

     6         1        -6        -5
     4        -3         0         1
     2         2         3         2
     0         2         0         1

Right-hand-side number 1

0.0000   -2.0000   -7.0000    6.0000

 Solution vector

-0.5000    1.0000    0.3333   -2.0000

Answer 1

您没有阅读精美的文档。上面写的很清楚

CALLING SEQUENCE:
   pvt = LUfactor ( a, n, ps );

您错误地使用了函数。您分配并填充了 pvt，然后您忽略了 LUfactor 的返回值。您不分配 pvt； LUfactor 函数执行的操作。您需要根据文档调用 LUfactor。