0
私はLUの分解を行っていますが、googelでこのコードを見つけましたが、出力 'pvt'と 'a'によってそれを理解したいのですが、私のpvtは正しくないので、いずれもここ.. おかげLUの分解を含む行列
私を修正する可能性がここに私のコード
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
これは[Code Review]（http：//codereview.stackexchange.com/） –