define _CRT_SECURE_NO_WARNI NGS include iostream include iomanip inclu

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
#define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#define vd vector<double>
#define vvd vector<vector<double> >
using namespace std;
double eps = 1e-4;
int n;
bool zeros_above_diag(vvd a) { //проверяет, есть ли нули выше главной диагонали
for (int i = 0; i < n; ++i) {
for (int j = i + 1; j < n; ++j) {
if (fabs(a[i][j]) > eps) {
return true;
}
}
}
return false;
}
vvd operator *(vvd x, vvd y) {
vvd ret(n, vd(n, 0.0));
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
for (int k = 0; k < n; k++) {
ret[i][j] += x[i][k] * y[k][j];
}
}
}
return ret;
}
int main() {
freopen("input.txt", "r", stdin);
cin >> n;
vvd a(n, vd(n, 0.0));
vvd t(n, vd(n, 0.0));
for (int i = 0; i < n; ++i) {
t[i][i] = 1.0;
}
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
cin >> a[i][j];
}
}
while (zeros_above_diag(a)) {
int maxi = 0; //1) находим максимальное по модулю число в исходной матрице А
int maxj = 1;
double max_element = 0.0;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < i; ++j) {
if (max_element < fabs(a[i][j])) {
maxi = i;
maxj = j;
max_element = fabs(a[i][j]);
}
}
}
double p = 2.0 * a[maxi][maxj]; //2) вычичляем коэффициенты с и s
double q = a[maxi][maxi] - a[maxj][maxj];
double d = sqrt(p * p + q * q);
double c;
double s;
double tiny_eps = 1e-15;
if (fabs(q) < tiny_eps) {
c = 1.0 / sqrt(2.0);
s = 1.0 / sqrt(2.0);
}
else {
double r = fabs(q) / (2.0 * d);
c = sqrt(0.5 + r);
if (fabs(p) * 1000.0 < fabs(q)) {
s = fabs(p) * (p * q > 0.0 ? 1.0 : -1.0) / (2.0 * c * d);
}
else {
s = sqrt(0.5 - r) * (p * q > 0.0 ? 1.0 : -1.0);
}
}
vvd cur_t(n, vd(n, 0.0)); //3) составляем временную матрицу поворота
for (int i = 0; i < n; ++i) {
cur_t[i][i] = 1.0;
}
cur_t[maxi][maxi] = c;
cur_t[maxj][maxj] = c;
cur_t[maxi][maxj] = -s;
cur_t[maxj][maxi] = s;
t = t * cur_t;
double a1 = a[maxi][maxi]; //4) изменяем исходную матрицу а
double a2 = a[maxj][maxj];
a[maxi][maxi] = c * c * a1 + s * s * a2 + 2.0 * c * s * a[maxi][maxj];
a[maxj][maxj] = s * s * a1 + c * c * a2 - 2.0 * c * s * a[maxi][maxj];
a[maxi][maxj] = 0.0;
a[maxj][maxi] = 0.0;
for (int i = 0; i < n; ++i) {
if (i == maxi || i == maxj) {
continue;
}
a1 = a[i][maxi];
a2 = a[i][maxj];
a[i][maxi] = c * a1 + s * a2;
a[maxi][i] = c * a1 + s * a2;
a[i][maxj] = -s * a1 + c * a2;
a[maxj][i] = -s * a1 + c * a2;
}
}
cout << "solution\n egienvalues:" << endl;
cout << fixed << setprecision(5);
for (int i = 0; i < n; ++i) {
cout << fixed << setw(9) << setprecision(6) << a[i][i] << endl;
}
cout << endl << "eigenvectors:" << endl;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
cout << fixed << setw(9) << setprecision(6) << t[i][j] << " ";
}
cout << endl;
}
return 0;
}