fork(55) download
copy 
  1. #include <cstdio>
  2. #include <algorithm>
  3. #include <cmath>
  4. //очистить выделенную память
  5. void clear(double **arr, int n)
  6. {
  7.     for(int i = 0; i < n; i++)
  8.         delete[] arr[i];
  9.     delete[] arr;
  10. }
  11. //создать копию массива
  12. double** clone(double **arr, int n)
  13. {
  14.     double **newArr = new double*[n];
  15.     for(int row = 0; row < n; row++)
  16.     {
  17.         newArr[row] = new double[n];
  18.         for(int col = 0; col < n; col++)
  19.             newArr[row][col] = arr[row][col];
  20.     }
  21.     return newArr;
  22. }
  23. //print the matrix
  24. void show(double **matrix, int n)
  25. {
  26.     for(int row = 0; row < n; row++){
  27.             for(int col = 0; col < n; col++)
  28.                 printf("%lf\t", matrix[row][col]);
  29.             printf("\n");
  30.     }
  31.     printf("\n");
  32. }
  33. //матричное умножение матриц
  34. double** matrix_multi(double **A, double **B, int n)
  35. {
  36.     double **result = new double*[n];
  37.     //заполнение нулями
  38.     for(int row = 0; row < n; row++){
  39.         result[row] = new double[n];
  40.         for(int col = 0; col < n; col++){
  41.             result[row][col] = 0;
  42.         }
  43.     }
  44.  
  45.     for(int row = 0; row < n; row++){
  46.         for(int col = 0; col < n; col++){
  47.             for(int j = 0; j < n; j++){
  48.                 result[row][col] += A[row][j] * B[j][col];
  49.             }
  50.         }
  51.     }
  52.     return result;
  53. }
  54. //умножение матрицы на число
  55. void scalar_multi(double **m, int n, double a){
  56.     for(int row = 0; row < n; row++)
  57.         for(int col = 0; col < n; col++){
  58.             m[row][col] *= a;
  59.         }
  60. }
  61. //вычисление суммы двух квадратных матриц
  62. void sum(double **A, double **B, int n)
  63. {
  64.     for(int row = 0; row < n; row++)
  65.         for(int col = 0; col < n; col++)
  66.             A[row][col] += B[row][col];
  67. }
  68.  
  69. //вычисление определителя
  70. double det(double **matrix, int n) //квадратная матрица размера n*n
  71. {
  72.     double **B = clone(matrix, n);
  73.     //приведение матрицы к верхнетреугольному виду
  74.     for(int step = 0; step < n - 1; step++)
  75.         for(int row = step + 1; row < n; row++)
  76.         {
  77.             double coeff = -B[row][step] / B[step][step]; //метод Гаусса
  78.             for(int col = step; col < n; col++)
  79.                 B[row][col] += B[step][col] * coeff;
  80.         }
  81.     //Рассчитать определитель как произведение элементов главной диагонали
  82.     double Det = 1;
  83.     for(int i = 0; i < n; i++)
  84.         Det *= B[i][i];
  85.     //Очистить память
  86.     clear(B, n);
  87.     return Det;
  88. }
  89.  
  90. int main()
  91. {
  92.     //Исходная матрица, динамический двухмерный массив
  93.     int n; scanf("%d", &n);
  94.     double **A = new double*[n];
  95.     for(int row = 0; row < n; row++)
  96.     {
  97.         A[row] = new double[n];
  98.             for(int col = 0; col < n; col++)
  99.                 scanf("%lf", &A[row][col]);
  100.     }
  101.  
  102.     /* Численное вычисление обратной матрицы по методу Ньютона-Шульца
  103.         1. Записать начальное приближение [Pan, Schreiber]:
  104.             1) Транспонировать данную матрицу
  105.             2) Нормировать по столбцам и строкам
  106.         2. Повторять процесс до достижения заданной точности.
  107.     */
  108.  
  109.     double N1 = 0, Ninf = 0; //норма матрицы по столбцам и по строкам
  110.     double **A0 = clone(A, n);       //инициализация начального приближения
  111.     for(size_t row = 0; row < n; row++){
  112.         double colsum = 0, rowsum = 0;
  113.         for(size_t col = 0; col < n; col++){
  114.             rowsum += fabs(A0[row][col]);
  115.             colsum += fabs(A0[col][row]);
  116.         }
  117.         N1 = std::max(colsum, N1);
  118.         Ninf = std::max(rowsum, Ninf);
  119.     }
  120.     //транспонирование
  121.     for(size_t row = 0; row < n - 1; row++){
  122.         for(size_t col = row + 1; col < n; col++)
  123.             std::swap(A0[col][row], A0[row][col]);
  124.     }
  125.     scalar_multi(A0, n, (1 / (N1 * Ninf))); //нормирование матрицы
  126.     //инициализация удвоенной единичной матрицы нужного размера
  127.     double **E2 = new double*[n];
  128.     for(int row = 0; row < n; row++)
  129.     {
  130.         E2[row] = new double[n];
  131.             for(int col = 0; col < n; col++){
  132.                 if(row == col)
  133.                     E2[row][col] = 2;
  134.                 else
  135.                     E2[row][col] = 0;
  136.             }
  137.     }
  138.     double **inv = clone(A0, n); //A_{0}
  139.     double EPS = 0.001;   //погрешность
  140.     if(det(A, n) != 0){ //если матрица не вырождена
  141.         while(fabs(det(matrix_multi(A, inv, n), n) - 1) >= EPS) //пока |det(A * A[k](^-1)) - 1| >= EPS
  142.         {
  143.             double **prev = clone(inv, n); //A[k-1]
  144.             inv = matrix_multi(A, prev, n);   //A.(A[k-1]^(-1))
  145.             scalar_multi(inv, n, -1);         //-A.(A[k-1]^(-1))
  146.             sum(inv, E2, n);                   //2E - A.(A[k-1]^(-1))
  147.             inv = matrix_multi(prev, inv, n); //(A[k-1]^(-1)).(2E - A.(A[k-1]^(-1)))
  148.             clear(prev, n);
  149.         }
  150.         //вывод матрицы на экран
  151.         show(inv, n);
  152.     }
  153.     else
  154.         printf("Impossible\n");
  155.     clear(A, n);
  156.     clear(E2, n);
  157.     return 0;
  158. }
Success #stdin #stdout 0s 3436KB
comments (?)
stdin
copy 
2
1 2
2 1
stdout
copy 
-0.333067	0.666400	
0.666400	-0.333067
https://ideone.com/7YphoX
Вычисление обратной матрицы методом Ньютона-Шульца
language:
C++14 (gcc 8.3)
created:
9 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!