我正在做 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