C++: СЛАУ методом Гаусса :)

Старая добрая проблема о решении СЛАУ методом Гаусса, всплыла в аспекте "а есть ли что-то на Си" и т.п.

Есть, прилагаю консольное решение C/C++ отдельной заметкой.

Данные берутся из файла текущей папки с именем gauss.dat. Его формат следующий:

• Строка описания системы (до 80 символов)
• Строка с размерностью системы, например, 4
• Строки матрицы
• Элементы вектора правой части в столбик.

Например для СЛАУ 4 на 4 может получиться файл:

My System
4
0 35 17.3 2
2 0 1 -1.5
3 17 12.1 14
-9 77.3 122.11 23
1
2
3
4

Сам код мне кажется довольно очевидным, просто приведу его. В конcоли Studio 2015 и т.п. работает.

Интересное в решении - как на C/C++ нумеровать элементы массивов с единицы, а не с нуля :)

#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define ICHAR 80  /*  Длина строки описания системы  */

#define IDEBUG 1  /* Печатать ли шаги декомпозиции */

int matrix_print_off(int nr, int nc, double **A) {
	int i, j;
	if (nr <= 0) return (-1);
	if (nc <= 0) return (-2);
	for (i = 1; i <= nr; i++) {
		for (j = 1; j <= nc; j++) {
			printf("%9.4f  ", A[i][j]);
		}
		printf("\n");
	}
	return (0);
}

int vector_print_off(int nr, double *x) {
	int i;

	if (nr <= 0) return (-1);

	for (i = 1; i <= nr; i++) {
		printf("%9.4f  \n", x[i]);
	}
	printf("\n");
	return (0);
}

double *matrix_on_vector_mult(int nr, int nc, double **A, double *x) {
	if (nr <= 0) return NULL;
	if (nc <= 0) return NULL;
	double *b = (double *)calloc(nr, sizeof(double)); b--;
	for (int i = 1; i <= nr; i++) {
		b[i] = 0;
		for (int j = 1; j <= nc; j++) b[i] += A[i][j] * x[j];
	}
	return b;
}

void gauss(double **a, double *b, double *x, int n) {
	int   i, j, k, m, rowx;
	double xfac, temp, temp1, amax;

	rowx = 0;   
	for (k = 1; k <= n - 1; ++k) {
		amax = (double)fabs(a[k][k]);
		m = k;
		for (i = k + 1; i <= n; i++) {  
			xfac = (double)fabs(a[i][k]);
			if (xfac > amax) { amax = xfac; m = i; }
		}
		if (m != k) { 
			rowx = rowx + 1;
			temp1 = b[k];
			b[k] = b[m];
			b[m] = temp1;
			for (j = k; j <= n; j++) {
				temp = a[k][j];
				a[k][j] = a[m][j];
				a[m][j] = temp;
			}
		}
		for (i = k + 1; i <= n; ++i) {
			xfac = a[i][k] / a[k][k];

			for (j = k + 1; j <= n; ++j) {
				a[i][j] = a[i][j] - xfac*a[k][j];
			}
			b[i] = b[i] - xfac*b[k];
		}

		if (IDEBUG == 1) {
			printf("\n Decomposition step %d\n\n", k);
			matrix_print_off(n, n, a);
		}
	}

	for (j = 1; j <= n; ++j) {
		k = n - j + 1;
		x[k] = b[k];
		for (i = k + 1; i <= n; ++i) {
			x[k] = x[k] - a[k][i] * x[i];
		}
		x[k] = x[k] / a[k][k];
	}

	if (IDEBUG == 1) printf("\nNumber of row exchanges = %d\n", rowx);
}

int main(void) {
	double  **a, **a0, *b, *x;
	float  aij, bi;
	char   desc[ICHAR];
	int    i, j, n;
	FILE   *finput;

    finput = fopen("gauss.dat", "r");
	if (finput == NULL) {
		printf("Data file gauss.dat not found\n");
		return(-1);
	}
		
	fgets(desc, ICHAR, finput);
	fscanf(finput, "%d", &n);
	printf("%s", desc);
	printf("\nDimension of matrix = %d\n\n", n);

	a = (double **)calloc(n, sizeof(double *)); --a; //Для нумерации с единицы :)
	a0 = (double **)calloc(n, sizeof(double *)); --a0; //Копия матрицы для проверки решения

	for (i = 1; i <= n; ++i) {
		a[i] = (double *)calloc(n, sizeof(double)); --a[i];
		a0[i] = (double *)calloc(n, sizeof(double)); --a0[i];
	}

	b = (double *)calloc(n, sizeof(double)); --b;
	x = (double *)calloc(n, sizeof(double)); --x;

	for (i = 1; i <= n; i++) {
		for (j = 1; j <= n; j++) {
			fscanf(finput, "%f ", &aij);
			a[i][j] = (double)aij;
			a0[i][j] = (double)aij;
		}
	}

	for (i = 1; i <= n; i++) {
		fscanf(finput, "%f ", &bi);
		b[i] = (double)bi;
	}

	fclose(finput); 

	printf("\nMatrices read from input file\n");

	printf("\nMatrix A\n\n");
	matrix_print_off(n, n, a);

	printf("\nVector b\n\n");
	vector_print_off(n, b);

	gauss(a, b, x, n);

	printf("\nSolution x\n\n");
	vector_print_off(n, x);

	if (IDEBUG == 1) {
		double *b0 = matrix_on_vector_mult(n, n, a0, x);
		printf("\nChecking\n\n");
		vector_print_off (n, b0);
	}

	getchar();
	return(0);
}

28.09.2017, 15:44 [13383 просмотра]