#include <string>
#include <fstream>
#include <iostream>
#include <iomanip>
#include "nr.h"
using namespace std;
// Driver for routine ludcmp
int main(void)
{
int j,k,l,m,n,dum;
string txt;
DP d;
ifstream fp("matrx1.dat");
if (fp.fail())
NR::nrerror("Data file matrx1.dat not found");
cout << fixed << setprecision(6);
getline(fp,txt);
while (!fp.eof()) {
getline(fp,txt);
fp >> n >> m;
getline(fp,txt);
Vec_INT indx(n),jndx(n);
Mat_DP a(n,n),xl(n,n),xu(n,n),x(n,n);
getline(fp,txt);
for (k=0;k<n;k++)
for (l=0;l<n;l++) fp >> a[k][l];
getline(fp,txt);
getline(fp,txt);
for (l=0;l<m;l++)
for (k=0;k<n;k++) fp >> x[k][l];
getline(fp,txt);
getline(fp,txt);
// Print out a-matrix for comparison with product of
// lower and upper decomposition matrices
cout << "original matrix:" << endl;
for (k=0;k<n;k++) {
for (l=0;l<n;l++) cout << setw(12) << a[k][l];
cout << endl;
}
// Perform the decomposition
NR::ludcmp(a,indx,d);
// Compose separately the lower and upper matrices
for (k=0;k<n;k++) {
for (l=0;l<n;l++) {
if (l > k) {
xu[k][l]=a[k][l];
xl[k][l]=0.0;
} else if (l < k) {
xu[k][l]=0.0;
xl[k][l]=a[k][l];
} else {
xu[k][l]=a[k][l];
xl[k][l]=1.0;
}
}
}
// Compute product of lower and upper matrices for
// comparison with original matrix
for (k=0;k<n;k++) {
jndx[k]=k;
for (l=0;l<n;l++) {
x[k][l]=0.0;
for (j=0;j<n;j++)
x[k][l] += (xl[k][j]*xu[j][l]);
}
}
cout << endl << "product of lower and upper ";
cout << "matrices (rows unscrambled):" << endl;
for (k=0;k<n;k++) {
dum=jndx[indx[k]];
jndx[indx[k]]=jndx[k];
jndx[k]=dum;
}
for (k=0;k<n;k++)
for (j=0;j<n;j++)
if (jndx[j] == k) {
for (l=0;l<n;l++)
cout << setw(12) << x[j][l];
cout << endl;
}
cout << endl << "lower matrix of the decomposition:" << endl;
for (k=0;k<n;k++) {
for (l=0;l<n;l++)
cout << setw(12) << xl[k][l];
cout << endl;
}
cout << endl << "upper matrix of the decomposition:" << endl;
for (k=0;k<n;k++) {
for (l=0;l<n;l++) cout << setw(12) << xu[k][l];
cout << endl;
}
cout << endl << "***********************************" << endl;
cout << "press return for next problem:" << endl;
cin.get();
}
fp.close();
return 0;
}
LU-Decomposition.zip_The First
版权申诉
166 浏览量
2022-09-23
11:51:54
上传
评论
收藏 7KB ZIP 举报
林当时
- 粉丝: 98
- 资源: 1万+