матрица, содержащая LU-декомпозицию

Question

Я делаю LU-декомпрессию, и я нашел этот код в 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 в документации.