123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297 |
- #include<stdio.h>
- #include<stdlib.h>
- #include<math.h>
- // Необязательный флаг для дебага (можно удалить)
- #define DEBUG 0
- // Фиксированная размерность матрицы (обрабатывается препроцессором, просто alias)
- #define N 10
- // Тип вектора
- typedef double vec[N];
- // Тип матрицы
- typedef double mat[N][N];
- // Скалярное произведение
- double dot(vec v, vec u) {
- double res = 0;
- for(int i = 0; i < N; i++) res += v[i]*u[i];
- return res;
- }
- // Прибавление вектора
- void add(vec v, vec u) {
- for(int i = 0; i < N; i++) {
- v[i] += u[i];
- }
- }
- // Вычитание вектора
- void sub(vec v, vec u) {
- for(int i = 0; i < N; i++) {
- v[i] -= u[i];
- }
- }
- // Умножение на скаляр
- void mul(vec v, double a) {
- for(int i = 0; i < N; i++) {
- v[i] *= a;
- }
- }
- // Умножение матрицы на вектор (res = m*v)
- void apply(mat m, vec v, vec res) {
- for(int i = 0; i < N; i++) {
- res[i] = 0;
- for(int j = 0; j < N; j++) {
- res[i] += m[i][j] * v[j];
- }
- }
- }
- // Умножение матриц (res = m*n)
- void multiply(mat m, mat n, mat res) {
- for(int i = 0; i < N; i++) {
- for(int j = 0; j < N; j++) {
- res[i][j] = 0;
- for(int k = 0; k < N; k++) {
- res[i][j] += m[i][k] * n[k][j];
- }
- }
- }
- }
- // Копирование матрицы (copy = m)
- void copy_mat(mat m, mat copy) {
- for(int i = 0; i < N; i++) {
- for(int j = 0; j < N; j++) {
- copy[i][j] = m[i][j];
- }
- }
- }
- // Получение j-ого столбца матрицы (res = m[*][j])
- void get_column(mat m, int j, vec res) {
- for(int i = 0; i < N; i++) {
- res[i] = m[i][j];
- }
- }
- // Получение i-ой строки матрицы (res = m[i][*])
- void get_row(mat m, int i, vec res) {
- for(int j = 0; j < N; j++) {
- res[j] = m[i][j];
- }
- }
- // Задание i-ого столбца матрицы (m[*][i] = col)
- void set_column(mat m, int i, vec col) {
- for(int j = 0; j < N; j++) {
- m[j][i] = col[j];
- }
- }
- // Нормализация столбцов матрицы
- void normalize(mat m) {
- vec v;
- for(int i = 0; i < N; i++) {
- get_column(m, i, v);
- double norm = sqrt(dot(v,v));
- mul(v, 1/norm);
- set_column(m, i, v);
- }
- }
- // Ортогонализация Грамма-Шмидта
- void ortogonalize(mat m) {
- vec v;
- vec u;
- for(int i = 0; i < N; i++) {
- get_column(m, i, v);
- for(int j = 0; j < i; j++) {
- get_column(m, j, u);
- double c = dot(v, u)/dot(u, u);
- mul(u, c);
- sub(v, u);
- }
- set_column(m, i, v);
- }
- normalize(m);
- }
- // QR-разложение (m = q*r)
- void qr_decomposition(mat m, mat q, mat r) {
- copy_mat(m, q);
- ortogonalize(q);
- vec v;
- vec u;
- for(int j = 0; j < N; j++) {
- get_column(m, j, v);
- for(int i = 0; i <= j; i++) {
- get_column(q, i, u);
- r[i][j] = dot(v, u);
- }
- for(int i = j+1; i < N; i++) r[i][j] = 0;
- }
- }
- // Транспонирование матрицы
- void transpose(mat m) {
- for(int i = 0; i < N; i++) {
- for(int j = 0; j < i; j++) {
- double d = m[i][j];
- m[i][j] = m[j][i];
- m[j][i] = d;
- }
- }
- }
- // Обращение верхне-треугольной матрицы (inv = r^{-1})
- void inverse_r(mat r, mat inv) {
- for(int i = 0; i < N; i++) {
- for(int j = 0; j < N; j++) {
- inv[i][j] = i == j ? 1./r[i][i] : 0;
- }
- }
- vec v;
- vec u;
- for(int j = 0; j < N; j++) {
- for(int i = j-1; i >= 0; i--) {
- get_row(r, i, v);
- get_column(inv, j, u);
- inv[i][j] = -dot(v, u)/r[i][i];
- }
- }
- }
- // Обращение матрицы (inv = m^{-1})
- void inverse(mat m, mat inv) {
- mat q;
- mat r;
- qr_decomposition(m, q, r);
- mat inv_r;
- inverse_r(r, inv_r);
- transpose(q);
- multiply(inv_r, q, inv);
- }
- // Решение СЛАУ методом QR-разложения
- void solve(mat m, vec v, vec sol) {
- mat inv;
- inverse(m, inv);
- apply(inv, v, sol);
- }
- // Чтение матрицы из файла
- void mat_from_file(char* filename, mat m) {
- double d;
- FILE* f;
- if ((f = fopen(filename, "r")) == NULL) return;
- for (int i = 0; i < N; i++) {
- for (int j = 0; j < N; j++) {
- fscanf(f, "%lf,", &d);
- m[i][j] = d;
- }
- }
- fclose(f);
- }
- // Красивая печать матрицы
- void print_mat(mat m) {
- for(int j = 0; j < N; j++) {
- printf("-----------|");
- }
- printf("\n");
- for(int i = 0; i < N; i++) {
- for(int j = 0; j < N; j++) {
- printf("%10.4lf |", m[i][j]);
- }
- printf("\n");
- for(int j = 0; j < N; j++) {
- printf("-----------|");
- }
- printf("\n");
- }
- }
- // Красивая печать вектора
- void print_vec(vec v) {
- for(int j = 0; j < N; j++) {
- printf("%10.4lf |", v[j]);
- }
- printf("\n");
- }
- int main(){
- // Матрица системы
- mat m;
- mat_from_file("lab_3.txt", m);
- printf("Coefficients: \n");
- print_mat(m);
- #if 1
- // Решение по умолчанию (x[i] = 2*i + 1)
- vec x;
- for(int i = 0; i < N; i++) {
- x[i] = 2*i+1;
- }
- // Свободный столбец уравнения (b = m*x)
- vec b;
- apply(m, x, b);
- printf("Constant terms: \n");
- print_vec(b);
- // Решение с помощью QR-разложения
- vec sol;
- solve(m, b, sol);
- printf("Solution: \n");
- print_vec(sol);
- #endif
- // Тесты
- #if 0
- ortogonalize(m);
- print_mat(m);
- for(int i = 0; i < N; i++){
- for(int j = i; j < N; j++){
- vec v;
- vec u;
- get_column(m, i, v);
- get_column(m, j, u);
- printf("%d:%d - %lf\n", i, j, dot(u,v));
- }
- }
- #endif
- #if 0
- mat q;
- mat r;
- qr_decomposition(m, q, r);
- print_mat(m);
- print_mat(q);
- print_mat(r);
- mat inv;
- inverse_r(r, inv);
- print_mat(inv);
- #endif
- #if 0
- transpose(m);
- print_mat(m);
- #endif
- #if 0
- mat inv;
- inverse(m, inv);
- print_mat(inv);
- #endif
- }
|