none
Повреждение кучи RRS feed

  • Общие обсуждения

  • Имеется программа (недоработанная):

    main.cpp

    #include<stdio.h>
    #include<stdlib.h>
    #include<conio.h>
    #include<math.h>
    #include"matrix.h"
    
    const double d = 0.0001;
    const double e = 2.7182818284;
    const double eps = 0.0000001;
    const double h = 0.01;
    typedef double (*function)(double x1, double x2);
    
    double f1(double x1, double x2)
    {
    	return(pow(e, sin(x2)) - 4 * x1);
    }
    
    double f2(double x1, double x2)
    {
    	return(sqrt(pow(x1, 2) + 3 * pow(x2, 2)) * sin (x1));
    }
    
    double s(double x1, double x2)
    {
    	return(5 * x1 - x2);
    }
    
    void Runge_Kutta(function f1, function f2, function s, double x, double y, int j, Matrix*M)
    {
    	int c, i = 0;
    	double k1, m1, k2, m2, k3, m3, k4, m4;
    
    	for(c = 0; c < 110; c++)
    	{
    		k1 = f1(x, y);
    		m1 = f2(x, y);
    		k2 = f1((x + h * k1 / 2), (y + h * m1 / 2));
    		m2 = f2((x + h * k1 / 2), (y + h * m1 / 2));
    		k3 = f1((x + h * k2 / 2), (y + h * m2 / 2));
    		m3 = f2((x + h * k2 / 2), (y + h * m2 / 2));
    		k4 = f1((x + h * k3), (y + h * m3));
    		m4 = f2((x + h * k3), (y + h * m3));
    		x = x + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
    		y = y + h * (m1 + 2 * m2 + 2 * m3 + m4) / 6;
    		
    		if((c == 10) || (c == 20) || (c == 30) || (c == 40) || (c == 45) ||
    		   (c == 50) || (c == 60) || (c == 66) || (c == 90) || (c == 110))
    		   M->Matrix[i++][j] = s(x, y);
    	}
    }
    
    void main()
    {
    	double x = 2, y = -3;
    	Matrix dR, dO, K, K_1, M, Mt, R, S, S1, S2, T, T1;
    
    	K = new_Matrix(0, 10, 10);
    	R = new_Matrix(0, 10, 1);
    	S = new_Matrix(0, 10, 1);
    	S1 = new_Matrix(0, 10, 2);
    	S2 = new_Matrix(0, 10, 2);
    
    	
    	K.Matrix[0][0] = 1.4380981508E-06;
    	K.Matrix[1][1] = 9.1018497960E-07;
    	K.Matrix[2][2] = 5.6740138858E-07;
    	K.Matrix[3][3] = 3.6134436446E-07;
    	K.Matrix[4][4] = 2.9321213273E-07;
    	K.Matrix[5][5] = 2.4121969626E-07;
    	K.Matrix[6][6] = 1.7072159697E-07;
    	K.Matrix[7][7] = 1.4283725936E-07;
    	K.Matrix[8][8] = 8.5110641850E-08;
    	K.Matrix[9][9] = 6.6127236730E-08;
    
    	R.Matrix[0][0] = 1.1990453016e+1;
    	R.Matrix[1][0] = 9.5385595392e+0;
    	R.Matrix[2][0] = 7.5325142241e+0;
    	R.Matrix[3][0] = 6.0110233129e+0;
    	R.Matrix[4][0] = 5.4145243611e+0;
    	R.Matrix[5][0] = 4.9115614762e+0;
    	R.Matrix[6][0] = 4.1320620808e+0;
    	R.Matrix[7][0] = 3.7796009645e+0;
    	R.Matrix[8][0] = 2.9173796705e+0;
    	R.Matrix[9][0] = 2.5725721232e+0;	
    
    	K_1 = inv_Matrix(K);
    	
    	{
    		Runge_Kutta(f1, f2, s, 2, (-3), 0, (&S));
    		Runge_Kutta(f1, f2, s, (2 + d), (-3), 0, (&S1));
    		Runge_Kutta(f1, f2, s, 2, (-3 + d), 1, (&S1));
    		Runge_Kutta(f1, f2, s, (2 - d), (-3), 0, (&S2));
    		Runge_Kutta(f1, f2, s, 2, (-3 - d), 1, (&S2));
    
    		M = sub_Matrix(S1, S2);
    		div_a_Matrix(&M, (2 * d));
    
    		dR = sub_Matrix(R, S);
    
    		Mt = cpy_Matrix(M);
    		trans_Matrix(&Mt);
    
    
    		T = mult_Matrix(Mt, K_1);
    		T1 = mult_Matrix(T, M);
    		T1 = inv_Matrix(T1);
    
    		/*
    		free_Matrix(&M);
    		free_Matrix(&dR);
    		free_Matrix(&Mt);
    		free_Matrix(&dO);
    		*/
    	}
    
    	_getch();
    }

    в которой на шаге T1 = inv_Matrix(T1) происходит повреждение кучи.
    Все функции, которые тут используются описаны в моём модуле для работы с матрицами:

    matrix.cpp

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include "matrix.h"
    
    //Выделяет память для матрицы и заносит её размеры
    Matrix new_Matrix(int ext, int m, int n)
    {
    	Matrix M;
    	
    	if((m > 0) && (n > 0))
    	{
    		M.ext = ext;
    		M.m = m;
    		M.n = n;
    		M.Matrix = (double**)calloc(M.m, sizeof(double*));
    		if(M.Matrix != NULL)
    		{
    			for(int i = 0; i < M.m; i++)
    				M.Matrix[i] = (double*)calloc(M.n, sizeof(double));
    			if(M.ext == 1)
    				M.Ext = (double*)calloc(M.m, sizeof(double));
    			return M;
    		}
    		free(M.Matrix);
    	}
    	M.ext = M.m = M.n = 0;
    	M.Matrix = NULL;
    	M.Ext = NULL;
    	return M;
    }
    
    //Производит удаление матрицы из памяти
    void free_Matrix(Matrix*M)
    {
    	M->n = 0;
    	for(int i = 0; i < M->m; i++)
    		free(M->Matrix[i]);
    	M->m = 0;
    	free(M->Matrix);
    	if(M->ext != 0)
    		free(M->Ext);
    	M->ext = 0;
    }
    
    //Производит ввод матрицы с клавиатуры
    void rd_Matrix(Matrix*M)
    {
    	char ext;
    	int i, m, n;
    	
    	printf("The extended matrix (y/n): ");
    	scanf("%c", &ext);
    	printf("Enter the dimensions of the main matrix (M, N): ");
    	scanf("%i%i", &m, &n);
    	if((m > 0) && (n > 0))
    	{
    		(*M) = new_Matrix(ext == 'y' ? 1 : 0, m, n);
    		if(M->Matrix != NULL)
    		{
    			puts("Enter the matrix:");
    			for(i = 0; i < M->m; i++)
    			{
    				for(int j = 0; j < M->n; j++)
    					scanf("%lf", &(M->Matrix[i][j]));
    				if(ext == 'y')
    					scanf("%lf", &(M->Ext[i]));
    			}
    		}
    		else
    			puts("Error! No memory!");
    	}
    	else
    		puts("Error! The size of the matrix is zero or negative");
    }
    
    //Производит вывод матрицы на экран
    void wr_Matrix(Matrix M)
    {
    	if((M.m > 0) && (M.n > 0))
    		for(int i = 0; i < M.m; i++)
    		{
    			for(int j = 0; j < M.n; j++)
    				printf("%lg ", M.Matrix[i][j]);
    			if(M.ext == 1)
    				printf("| %lg", M.Ext[i]);
    			printf("\n");
    		}
    	else
    		puts("Error! The size of the matrix is zero or negative");
    }
    
    //Возвращает копию матрицы
    Matrix cpy_Matrix(Matrix M)
    {
    	Matrix T;
    
    	T = new_Matrix(M.ext, M.m, M.n);
    	if(T.Matrix != NULL)
    		for(int i = 0; i < M.m; i++)
    		{
    			for(int j = 0; j < M.n; j++)
    				T.Matrix[i][j] = M.Matrix[i][j];
    			if(T.ext == 1)
    				T.Ext[i] = M.Ext[i];
    		}
    	return T;
    }
    
    void swp_rows_Matrix(Matrix*M, int i1, int i2)
    {
    	double t;
    
    	if(M->ext == 1)
    	{
    		t = M->Ext[i2];
    		M->Ext[i2] = M->Ext[i1];
    		M->Ext[i1] = t;
    	}
    	for(int j = 0; j < M->n; j++)
    	{	
    		t = M->Matrix[i2][j];
    		M->Matrix[i2][j] = M->Matrix[i1][j];
    		M->Matrix[i1][j] = t;
    	}
    }
    
    //Возвращает 1 если матрицы эквивалентны, иначе 0
    int eq_Matrix(Matrix A, Matrix B)
    {
    	if((A.m == B.m) && (A.n == B.n) && (A.m > 0) && (A.n > 0))
    	{
    		for(int i = 0; i < A.m; i++)
    			for(int j = 0; j < A.n; j++)
    				if(A.Matrix[i][j] != B.Matrix[i][j])
    					return 0;
    		return 1;
    	}
    	else
    		return 0;
    }
    
    //Умножает матрицу на число
    void mult_a_Matrix(Matrix*M, double a)
    {
    	if((M->m > 0) && (M->n > 0))
    		for(int i = 0; i < M->m; i++)
    			for(int j = 0; j < M->n; j++)
    				M->Matrix[i][j] *= a;
    }
    
    //Делит матрицу на число
    void div_a_Matrix(Matrix*M, double a)
    {
    	if((M->m > 0) && (M->n > 0))
    		for(int i = 0; i < M->m; i++)
    			for(int j = 0; j < M->n; j++)
    				M->Matrix[i][j] /= a;
    }
    
    //Возвращает сумму матриц A и B
    Matrix add_Matrix(Matrix A, Matrix B)
    {
    	Matrix M;
    
    	if((A.m == B.m) && (A.n == B.n))
    		M = new_Matrix(0, A.m, A.n);
    	else
    		M = new_Matrix(0, 0, 0);
    	if(M.Matrix != NULL)
    	{
    		for(int i = 0; i < A.m; i++)
    			for(int j = 0; j < A.n; j++)
    				M.Matrix[i][j] = A.Matrix[i][j] + B.Matrix[i][j];
    	}
    	return M;
    }
    
    //Возвращает разность матриц A и B
    Matrix sub_Matrix(Matrix A, Matrix B)
    {
    	Matrix M;
    
    	if((A.m == B.m) && (A.n == B.n))
    		M = new_Matrix(0, A.m, A.n);
    	else
    		M = new_Matrix(0, 0, 0);
    	if(M.Matrix != NULL)
    	{
    			for(int i = 0; i < A.m; i++)
    				for(int j = 0; j < A.n; j++)
    					M.Matrix[i][j] = A.Matrix[i][j] - B.Matrix[i][j];
    	}
    	return M;
    }
    
    //Возвращает произведение матриц A и B
    Matrix mult_Matrix(Matrix A, Matrix B)
    {
    	Matrix M;
    
    	if(A.n == B.m)
    		M = new_Matrix(0, A.m, B.n);
    	else
    		M = new_Matrix(0, 0, 0);
    	if(M.Matrix != NULL)
    	{
    		for(int i = 0; i < A.m; i++)
    			for(int k = 0; k < B.m; k++)
    			{
    				M.Matrix[i][k] = 0;
    				for(int j = 0; j < B.n; j++)
    					M.Matrix[i][j] += A.Matrix[i][k] * B.Matrix[k][j];
    			}
    	}
    	return M;
    }
    
    //Транспонирует матрицу
    void trans_Matrix(Matrix*M)
    {
    	Matrix T;
    
    	T = new_Matrix(0, M->n, M->m);
    	if(T.Matrix != NULL)
    	{
    		for(int i = 0; i < M->m; i++)
    			for(int j = 0; j < M->n; j++)
    				T.Matrix[j][i] = M->Matrix[i][j];
    		free_Matrix(M);
    		(*M) = T;
    	}
    }
    
    //Возводит матрицу в степень
    void deg_Matrix(Matrix*M, double deg)
    {
    	if((M->m == M->n) && (M->m > 0))
    		for(int i = 1; i < deg; i++)
    			(*M) = mult_Matrix((*M), (*M));
    }
    
    //Возвращает единичную матрицу
    Matrix E_Matrix(int m, int n)
    {
    	Matrix M;
    
    	if(m == n)
    		M = new_Matrix(0, m, n);
    	else
    		M = new_Matrix(0, 0, 0);
    	if(M.Matrix != NULL)
    	{
    		if(M.Matrix != NULL)
    			for(int i = 0; i < m; i++)
    				M.Matrix[i][i] = 1;
    	}
    	return M;
    }
    
    //Возвращает определитель матрицы
    double det_Matrix(Matrix M)
    {
    	double r = 1, t;
    	int c = 0, i;
    	Matrix T;
    	
    	T = cpy_Matrix(M);
    	if((T.m == T.n) && (T.m > 0) && (T.Matrix != NULL))
    	{
    		int k = 0;
    		while(k != T.m - 1)
    		{
    			if(T.Matrix[k][k] == 0)
    				for(i = 0; i < T.m; i++)
    					if((T.Matrix[i][k] != 0) && (i > k))
    					{
    						swp_rows_Matrix((&T), i, k);
    						c++;
    						i = T.m;
    					}
    			r *= T.Matrix[k][k];
    			for(i = k + 1; i < T.m; i++)
    			{
    				t = (-T.Matrix[i][k] / T.Matrix[k][k]);
    				for(int j = k; j < T.n; j++)
    					T.Matrix[i][j] += t * T.Matrix[k][j];
    			}
    			k++;
    		}
    		r *= T.Matrix[k][k];
    		free_Matrix(&T);
    		return r * pow((float)(-1), c);
    	}
    	return 9999999;
    }
    
    //Возвращает обратную матрицу
    Matrix inv_Matrix(Matrix M)
    {
    	double t;
    	int i, j;
    	Matrix T, T1;
    
    	T = cpy_Matrix(M);
    	if((T.m == T.n) && (T.Matrix != NULL) && (det_Matrix(T) != 0))
    	{
    		T1 = E_Matrix(T.m, T.n);
    		if((T.Matrix != NULL) && (T1.Matrix != NULL))
    			for(int k = 0; k < T.m; k++)
    			{
    				if(T.Matrix[k][k] == 0)
    				for(i = 0; i < T.m; i++)
    					if((T.Matrix[i][k] != 0) && (i > k))
    					{
    						swp_rows_Matrix((&T), i, k);
    						swp_rows_Matrix((&T1), i, k);
    						i = T.m;
    					}
    				t = T.Matrix[k][k];
    				for(j = 0; j < T.n; j++)
    				{
    					T.Matrix[k][j] /= t;
    					T1.Matrix[k][j] /= t;
    				}
    				for(i = 0; i < T.m; i++)
    					if(i != k)
    					{
    						t = T.Matrix[i][k];
    						for(j = 0; j < T.n; j++)
    						{
    							(T.Matrix[i][j]) -= (T.Matrix[k][j]) * t;
    							(T1.Matrix[i][j]) -= (T1.Matrix[k][j]) * t;
    						}
    					}
    			}
    		else
    			free_Matrix(&T);
    	}
    	else
    		T1 = new_Matrix(1, 0, 0);
    	return T1;
    }
    
    /// Возвращает решение СЛАУ
    Matrix sltn_SLAU(Matrix M)
    {
    	double t;
    	int i, j;
    	Matrix T;
    
    	T = cpy_Matrix(M);
    	if((T.ext == 1) && (T.Matrix != NULL))
    	{
    		if((T.Matrix != NULL) && (T.Ext != NULL))
    			for(int k = 0; k < T.m; k++)
    			{
    				if(T.Matrix[k][k] == 0)
    				for(i = 0; i < T.m; i++)
    					if((T.Matrix[i][k] != 0) && (i > k))
    					{
    						swp_rows_Matrix((&T), i, k);
    						i = T.m;
    					}
    				t = T.Matrix[k][k];
    				T.Ext[k] /= t;
    				for(j = 0; j < T.n; j++)
    					T.Matrix[k][j] /= t;
    				for(i = 0; i < T.m; i++)
    					if(i != k)
    					{
    						t = T.Matrix[i][k];
    						T.Ext[i] -= T.Ext[k] * t;
    						for(j = 0; j < T.n; j++)
    							T.Matrix[i][j] -= T.Matrix[k][j] * t;
    					}
    			}
    		else
    			free_Matrix(&T);
    	}
    	return T;
    }

    matrix.h

    #ifndef _MATRIX_H
    
    	typedef struct{
    		int ext; //индикатор расширенной матрицы
    		int m; //количество строк
    		int n; //количество столбцов (без учёта расширения матрицы)
    		double**Matrix; //основная матрица
    		double*Ext; //расширение основной матрицы
    	}Matrix;
    
    	Matrix new_Matrix(int ext, int m, int n);
    	void free_Matrix(Matrix*M);
    	void rd_Matrix(Matrix*M);
    	void wr_Matrix(Matrix M);
    	Matrix cpy_Matrix(Matrix M);
    	int eq_Matrix(Matrix A, Matrix B);
    	void mult_a_Matrix(Matrix*M, double a);
    	void div_a_Matrix(Matrix*M, double a);
    	Matrix add_Matrix(Matrix A, Matrix B);
    	Matrix sub_Matrix(Matrix A, Matrix B);
    	Matrix mult_Matrix(Matrix A, Matrix B);
    	void trans_Matrix(Matrix*M);
    	void deg_Matrix(Matrix*M, double deg);
    	Matrix E_Matrix(int m, int n);
    	double det_Matrix(Matrix M);
    	Matrix inv_Matrix(Matrix M);
    	Matrix sltn_SLAU(Matrix M);
    
    #endif
    Программа реализует "Метод максимального правдоподобия" для оценки параметров системы. После инициализации переменной K_1 идёт один шаг этого метода. Пожалуйста, подскажите, где и в чём ошибка и почему повреждается куча?
    (Если не хватает информации о коде, пишите, дам всю возможную информацию)
    22 января 2015 г. 23:29

Все ответы