SHARE
TWEET

Untitled

a guest Oct 13th, 2019 85 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. #define data_t double
  6. #define MODIFIER_A "%4.01lf"
  7. #define MODIFIER_B "%9.05lf"
  8.  
  9.  
  10. void swap(data_t* x, data_t* y) {
  11.     data_t buffer = *x;
  12.     *x = *y;
  13.     *y = buffer;
  14. }
  15.  
  16. // Добавить к одной строке другую, применяется в случае нуля на диагонали
  17. void swap_rows_inplace(data_t* A, size_t matrix_size, size_t row_first_index, size_t row_second_index) {
  18.     data_t* row_first = A + matrix_size * row_first_index;
  19.     data_t* row_second = A + matrix_size * row_second_index;
  20.  
  21.     // cache friendly optimization for sse2 set of instructions
  22.     for (size_t index = 0; index < matrix_size; ++index) {
  23.         swap(row_first + index, row_second + index);
  24.     }
  25. }
  26.  
  27.  
  28. void gauss_forward(data_t *A, data_t *B, size_t matrix_size) {
  29.     for (size_t row_index = 0; row_index < matrix_size; ++row_index) {
  30.         if (A[matrix_size * row_index + row_index] == 0) {
  31.  
  32.             for (size_t next_row_index = row_index + 1; next_row_index < matrix_size; ++next_row_index) {
  33.                 if (A[matrix_size * next_row_index + row_index] != 0) {
  34.                     swap_rows_inplace(A, matrix_size, row_index, next_row_index);
  35.                     swap(B + row_index, B + next_row_index);
  36.                     break;
  37.                 }
  38.             }
  39.  
  40.             if (A[matrix_size * row_index + row_index] == 0) {
  41.                 printf("Oops, det A = 0\n");
  42.                 exit(EXIT_FAILURE);
  43.             }
  44.  
  45.         }
  46.  
  47.         for (size_t next_row_index = row_index + 1; next_row_index < matrix_size; ++next_row_index) {
  48.             data_t multiplier = A[next_row_index * matrix_size + row_index] / A[row_index * matrix_size + row_index];
  49.  
  50.             for (size_t k = row_index; k < matrix_size; ++k) {
  51.                 A[next_row_index * matrix_size + k] -= A[row_index * matrix_size + k] * multiplier;
  52.             }
  53.  
  54.             B[next_row_index] -= multiplier * B[row_index];
  55.         }
  56.  
  57.     }
  58. }
  59.  
  60. void gauss_backward(data_t* A, data_t* B, size_t matrix_size) {
  61.     for (int row_index = (int)matrix_size - 1; row_index >= 0; --row_index) {
  62.         data_t diag_elem = A[row_index * matrix_size + row_index];
  63.  
  64.         A[row_index * matrix_size + row_index] /= diag_elem;
  65.         B[row_index] /= diag_elem;
  66.  
  67.         for (int prev_row_index = row_index - 1; prev_row_index >= 0; --prev_row_index) {
  68.             data_t local_mul =  A[prev_row_index * matrix_size + row_index];
  69.             A[prev_row_index * matrix_size + row_index] = 0;
  70.             B[prev_row_index] -= B[row_index] * local_mul;
  71.         }
  72.     }
  73. }
  74.  
  75. void print_system(data_t* A, data_t* B, size_t matrix_size) {
  76.     for (size_t i = 0; i < matrix_size; ++i) {
  77.         for (size_t j = 0; j < matrix_size; ++j) {
  78.             printf(MODIFIER_A " ", A[matrix_size * i + j]);
  79.         }
  80.         printf(" |" MODIFIER_B "\n", B[i]);
  81.     }
  82.     printf("\n");
  83. }
  84.  
  85. int main(void) {
  86.     size_t matrix_size = 20;
  87.     data_t* A = (data_t*) malloc (sizeof(data_t) * matrix_size * matrix_size);
  88.     data_t* B = (data_t*) malloc (sizeof(data_t) * matrix_size);
  89.  
  90.     if (A == NULL || B == NULL) {
  91.         printf("Not enough memory\n");
  92.         free(A);
  93.         free(B);
  94.         return EXIT_FAILURE;
  95.     }
  96.  
  97.     // Заполняем А
  98.     A[0 * matrix_size + 0] = 1;
  99.     for (size_t i = 1; i + 1 < matrix_size; ++i) {
  100.         A[i * matrix_size + i] = -2; // Главная диагональ
  101.         A[i * matrix_size + i + 1] = 1; // Над главной диагональю
  102.         A[i * matrix_size + i - 1] = 1; // под главной диагональю
  103.     }
  104.  
  105.     for (size_t i = 0; i < matrix_size; ++i) {
  106.         A[(matrix_size - 1) * matrix_size + i] = 1;
  107.     }
  108.  
  109.     // Заполняем В
  110.     B[0] = 1;
  111.     for (size_t i = 1; i + 1 < matrix_size; ++i) {
  112.         B[i] = 2.0 / (i + 1) / (i + 1);
  113.     }
  114.     B[matrix_size - 1] = -19.0 / 3;
  115.  
  116.     print_system(A, B, matrix_size);
  117.     gauss_forward(A, B, matrix_size);
  118.  
  119.     // print_system(A, B, matrix_size); // диагональная после прямого хода гаусса
  120.     gauss_backward(A, B, matrix_size);
  121.  
  122.     print_system(A, B, matrix_size);
  123.  
  124.     // В содержит решение методом гаусса
  125.  
  126.  
  127. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top