include cstdio include algorithm include cmath очистить выделенную пам

  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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#include <cstdio>
#include <algorithm>
#include <cmath>
//очистить выделенную память
void clear(double **arr, int n)
{
for(int i = 0; i < n; i++)
delete[] arr[i];
delete[] arr;
}
//создать копию массива
double** clone(double **arr, int n)
{
double **newArr = new double*[n];
for(int row = 0; row < n; row++)
{
newArr[row] = new double[n];
for(int col = 0; col < n; col++)
newArr[row][col] = arr[row][col];
}
return newArr;
}
//print the matrix
void show(double **matrix, int n)
{
for(int row = 0; row < n; row++){
for(int col = 0; col < n; col++)
printf("%lf\t", matrix[row][col]);
printf("\n");
}
printf("\n");
}
//матричное умножение матриц
double** matrix_multi(double **A, double **B, int n)
{
double **result = new double*[n];
//заполнение нулями
for(int row = 0; row < n; row++){
result[row] = new double[n];
for(int col = 0; col < n; col++){
result[row][col] = 0;
}
}
for(int row = 0; row < n; row++){
for(int col = 0; col < n; col++){
for(int j = 0; j < n; j++){
result[row][col] += A[row][j] * B[j][col];
}
}
}
return result;
}
//умножение матрицы на число
void scalar_multi(double **m, int n, double a){
for(int row = 0; row < n; row++)
for(int col = 0; col < n; col++){
m[row][col] *= a;
}
}
//вычисление суммы двух квадратных матриц
void sum(double **A, double **B, int n)
{
for(int row = 0; row < n; row++)
for(int col = 0; col < n; col++)
A[row][col] += B[row][col];
}
//вычисление определителя
double det(double **matrix, int n) //квадратная матрица размера n*n
{
double **B = clone(matrix, n);
//приведение матрицы к верхнетреугольному виду
for(int step = 0; step < n - 1; step++)
for(int row = step + 1; row < n; row++)
{
double coeff = -B[row][step] / B[step][step]; //метод Гаусса
for(int col = step; col < n; col++)
B[row][col] += B[step][col] * coeff;
}
//Рассчитать определитель как произведение элементов главной диагонали
double Det = 1;
for(int i = 0; i < n; i++)
Det *= B[i][i];
//Очистить память
clear(B, n);
return Det;
}
int main()
{
//Исходная матрица, динамический двухмерный массив
int n; scanf("%d", &n);
double **A = new double*[n];
for(int row = 0; row < n; row++)
{
A[row] = new double[n];
for(int col = 0; col < n; col++)
scanf("%lf", &A[row][col]);
}
/* Численное вычисление обратной матрицы по методу Ньютона-Шульца
1. Записать начальное приближение [Pan, Schreiber]:
1) Транспонировать данную матрицу
2) Нормировать по столбцам и строкам
2. Повторять процесс до достижения заданной точности.
*/
double N1 = 0, Ninf = 0; //норма матрицы по столбцам и по строкам
double **A0 = clone(A, n); //инициализация начального приближения
for(size_t row = 0; row < n; row++){
double colsum = 0, rowsum = 0;
for(size_t col = 0; col < n; col++){
rowsum += fabs(A0[row][col]);
colsum += fabs(A0[col][row]);
}
N1 = std::max(colsum, N1);
Ninf = std::max(rowsum, Ninf);
}
//транспонирование
for(size_t row = 0; row < n - 1; row++){
for(size_t col = row + 1; col < n; col++)
std::swap(A0[col][row], A0[row][col]);
}
scalar_multi(A0, n, (1 / (N1 * Ninf))); //нормирование матрицы
//инициализация удвоенной единичной матрицы нужного размера
double **E2 = new double*[n];
for(int row = 0; row < n; row++)
{
E2[row] = new double[n];
for(int col = 0; col < n; col++){
if(row == col)
E2[row][col] = 2;
else
E2[row][col] = 0;
}
}
double **inv = clone(A0, n); //A_{0}
double EPS = 0.001; //погрешность
if(det(A, n) != 0){ //если матрица не вырождена
while(fabs(det(matrix_multi(A, inv, n), n) - 1) >= EPS) //пока |det(A * A[k](^-1)) - 1| >= EPS
{
double **prev = clone(inv, n); //A[k-1]
inv = matrix_multi(A, prev, n); //A.(A[k-1]^(-1))
scalar_multi(inv, n, -1); //-A.(A[k-1]^(-1))
sum(inv, E2, n); //2E - A.(A[k-1]^(-1))
inv = matrix_multi(prev, inv, n); //(A[k-1]^(-1)).(2E - A.(A[k-1]^(-1)))
clear(prev, n);
}
//вывод матрицы на экран
show(inv, n);
}
else
printf("Impossible\n");
clear(A, n);
clear(E2, n);
return 0;
}