63 lines
1.7 KiB
C++
63 lines
1.7 KiB
C++
// Problem 2: Scaled partial pivoting
|
|
#include <iostream>
|
|
#include <iomanip>
|
|
#include <cmath>
|
|
#include <vector>
|
|
using namespace std;
|
|
|
|
double roundSig(double val, int sig) {
|
|
if (val == 0) return 0;
|
|
double d = ceil(log10(fabs(val)));
|
|
double p = pow(10, sig - d);
|
|
return round(val * p) / p;
|
|
}
|
|
|
|
void solve(vector<vector<double>> A, int n, int sig) {
|
|
cout << "\n" << sig << " significant digits:\n";
|
|
|
|
vector<double> s(n);
|
|
for (int i = 0; i < n; i++)
|
|
for (int j = 0; j < n; j++) s[i] = max(s[i], fabs(A[i][j]));
|
|
|
|
vector<int> idx(n);
|
|
for (int i = 0; i < n; i++) idx[i] = i;
|
|
|
|
for (int k = 0; k < n - 1; k++) {
|
|
int maxRow = k;
|
|
double maxRatio = fabs(A[idx[k]][k]) / s[idx[k]];
|
|
for (int i = k + 1; i < n; i++) {
|
|
double ratio = fabs(A[idx[i]][k]) / s[idx[i]];
|
|
if (ratio > maxRatio) { maxRatio = ratio; maxRow = i; }
|
|
}
|
|
swap(idx[k], idx[maxRow]);
|
|
|
|
for (int i = k + 1; i < n; i++) {
|
|
double m = roundSig(A[idx[i]][k] / A[idx[k]][k], sig);
|
|
for (int j = k; j <= n; j++)
|
|
A[idx[i]][j] = roundSig(A[idx[i]][j] - m * A[idx[k]][j], sig);
|
|
}
|
|
}
|
|
|
|
vector<double> x(n);
|
|
for (int i = n - 1; i >= 0; i--) {
|
|
x[i] = A[idx[i]][n];
|
|
for (int j = i + 1; j < n; j++)
|
|
x[i] = roundSig(x[i] - A[idx[i]][j] * x[j], sig);
|
|
x[i] = roundSig(x[i] / A[idx[i]][i], sig);
|
|
}
|
|
|
|
for (int i = 0; i < n; i++)
|
|
cout << "x" << (i+1) << " = " << fixed << setprecision(6) << x[i] << "\n";
|
|
}
|
|
|
|
int main() {
|
|
vector<vector<double>> A = {
|
|
{4.13, -2.20, 0.95, 3.02},
|
|
{6.14, 4.45, -1.45, -4.02},
|
|
{1.03, 1.86, 0.44, 5.22}
|
|
};
|
|
solve(A, 3, 6);
|
|
solve(A, 3, 3);
|
|
return 0;
|
|
}
|