Консультация # 186281: Здравствуйте! Прошу помощи в следующем вопросе: Реализовывал два алгоритм LU разложения матрицы. А именно метод исключения, и рекурсивный метод. Мне нужно сравнить работоспособность этих алгоритмов, и определить их сложность. Но у меня возникла проблема. Алгоритмы которые я делал с сайта: URL...
Здравствуйте! Прошу помощи в следующем вопросе: Реализовывал два алгоритм LU разложения матрицы. А именно метод исключения, и рекурсивный метод. Мне нужно сравнить работоспособность этих алгоритмов, и определить их сложность. Но у меня возникла проблема. Алгоритмы которые я делал с сайта: URL >> Теория (сейчас этой страницы нет на сайте, поэтому залил ее на файлообменник) раскладывают не все матрицы. Можно ли как-нибудь доработать мою
программу, чтобы она раскладывала на треугольные матрицы все квадратные матрицы? Подсчитать сложность алгоритма.
Код программы:
Код :
#include <iostream>
#include <stdlib.h>
#include <time.h>
using namespace std;
int k = 0;
void functionForTheMatrixMRekurs(double **M, int n, double **A) //recursive function of the matrix M
{
if (k != (n - 1)) // проверка на количесвто вызово рекурсии
{
for (int i = k + 1; i < n; i++) {
M[i][k] = A[i][k] / A[k][k]; //вычисляються элементы стоящие на строке равной номеру вызова рекурсии
for (int j = k + 1; j < n; j++)
M[i][j] = A[i][j] - A[i][k] * A[k][j]; //вычисление элементов стоящих на позициях i j используя элементы стоящие на позициях i j k
}
} else // если количество вызово рекурсии превышает количество строк в матрице значит всё вычислили и выходим из функции
return;
for (int l = 0; l < n; l++) // в любой млучае в матрице М 0 строка будет как в исходной матрице
M[0][l] = A[0][l];
k++; //счётчик рекрсии
functionForTheMatrixMRekurs(M, n, A); //рекурсивный вызов функции
}
void functionForTheMatrixL_and_U(double **L, int n, double **A, double **U) //function matrix for the calculation of M and L
{
for (int i = 0; i < n; i++)
for (int j = i; j < n; j++)
L[j][i] = A[j][i] / A[i][i];
for (int k = 1; k < n; k++) {
for (int i = k; i < n; i++)
for (int j = k - 1; j < n; j++)
U[i][j] = U[i][j] - L[i][k - 1] * U[k - 1][j]; //вычисление элементов проводиться простым способом
}
for (int i = 0; i < 2; i++) {
for (int j = i + 1; j < 3; j++) {
L[i][j] = 0; // в матриц Л в любой случаем в нижем левом углу будут нули
}
}
for (int k = 0; k < 3; k++)
L[k][k] = 1; // а по диагонали 1
}
void functionForTheMatrixLRekurs(double **L, int n, double **M)//recursive function of the matrix L
{
int l;
if (k == n)
return;
else {
for (int j = 0; j < n; j++) {
L[k][j] = M[k][j];
if (k < n - 1) {
for (int i = k + 1; i < n; i++)
L[k][i] = 0;
}
}
L[k][k] = 1;
k++;
functionForTheMatrixLRekurs(L, n, M);
}
}
void functionForTheMatrixURekurs(double **U, int n, double **M)//recursive function of the matrix U
{
if (k == n) //проверка количесва рекурсивных вызовов
return;
else {
for (int i = k; i < n; i++)
for (int j = k - 1; j < n; j++)
U[i][j] = U[i][j] - M[i][k - 1] * U[k - 1][j];
}
k++;
functionForTheMatrixURekurs(U, n, M);
}
void showMatrix(double **matrix, int n, int m) {
for (int i = 0; i < n; i++) {
cout << endl;
for (int j = 0; j < m; j++) {
cout << matrix[i][j] << " ";
}
}
cout << endl;
}
int main() {
double **A, **MR, **LR, **UR, **U, **L;
A = new double*[3];
MR = new double*[3];
LR = new double*[3];
UR = new double*[3];
L = new double*[3];
U = new double*[3];
for (int i = 0; i < 3; i++) {
A[i] = new double [3];
MR[i] = new double [3];
LR[i] = new double [3];
UR[i] = new double [3];
L[i] = new double [3];
U[i] = new double [3];
}
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) // здесь мы вводим наши элементы
{
cin >> A[i][j];
UR[i][j] = A[i][j]; // часть элементов У будет будет такоеже как в матрицу А
U[i][j] = A[i][j];
}
}
clock_t startRecurs, finishRecurs, startIs, finishIs;
startRecurs = clock();
{
functionForTheMatrixMRekurs(MR, 3, A);
k = 0;
functionForTheMatrixLRekurs(LR, 3, MR);
k = 1;
functionForTheMatrixURekurs(UR, 3, MR);
}
finishRecurs = clock();
startIs = clock();
{
functionForTheMatrixL_and_U(L, 3, A, U);
}
finishIs = clock();
for (int i = 1; i < 3; i++) {
for (int j = 0; j < i; j++) //нижний угол всегда будет нулевой
{
UR[i][j] = 0;
U[i][j] = 0;
}
}
//вывод
showMatrix(A, 3, 3);
showMatrix(MR, 3, 3);
showMatrix(LR, 3, 3);
showMatrix(UR, 3, 3);
showMatrix(L, 3, 3);
showMatrix(U, 3, 3);
double d0 = (double)(finishRecurs - startRecurs) /CLOCKS_PER_SEC;
double d1 = (double)(finishIs - startIs) /CLOCKS_PER_SEC;
cout << "Recursive method: " << d0 << endl;
cout << "Exception method: " << d1 << endl;
for (int i = 0; i < 3; i++)
delete [] A[i];
delete [] A;
for (int i = 0; i < 3; i++)
delete [] MR[i];
delete [] MR;
for (int i = 0; i < 3; i++)
delete [] LR[i];
delete [] LR;
for (int i = 0; i < 3; i++)
delete [] UR[i];
delete [] UR;
for (int i = 0; i < 3; i++)
delete [] U[i];
delete [] U;
for (int i = 0; i < 3; i++)
delete [] L[i];
delete [] L;
return 0;
}
Здравствуйте, Сушко Василий! Несколько переделал Ваш код: 1) В первом расчете сначала копирую матрицу в MR, затем рекурсивно ее же и обрабатываю В результате получаем в верхнем треугольнике с диагональю - матрицу U, в нижнем, без единичной диагонали, матрицу L Осталось только сформировать эти две матрицы (рекурсия здесь уже не нужна) 2) Второй расчет сделал на основе вот этой теории Добавил
только обнуление необходимых элементов и задание для L диагональных единичных элементов 3) Расчет длительности сделал на основе performance
Код :
#include <windows.h> ;для performance
#include <iostream>
#include <stdlib.h>
#include <time.h>
using namespace std;
int k = 0;
//переменные для расчета длительности
struct performance
{
_LARGE_INTEGER start;
_LARGE_INTEGER stop;
_LARGE_INTEGER frequency;
};
#define Start(perf) QueryPerformanceCounter(&perf.start)
#define Stop(perf) QueryPerformanceCounter(&perf.stop)
//запрос частоты
bool QueryPerfCounter(performance * perf)
{
return (0 != QueryPerformanceFrequency(&perf->frequency));
}
//длительность в мкс
double Duration(performance * perf)
{
return (((double)(perf->stop.QuadPart - perf->start.QuadPart) * 1.0e6 /
(double) perf->frequency.QuadPart));
}
void functionForTheMatrixMRekurs(double **M, int n) //recursive function of the matrix M
{
if (k != (n - 1)) // проверка на количесвто вызово рекурсии
{
for (int i = k + 1; i < n; i++) {
M[i][k] /= M[k][k]; //вычисляються элементы стоящие на строке равной номеру вызова рекурсии
for (int j = k + 1; j < n; j++)
M[i][j] -= M[i][k] * M[k][j]; //вычисление элементов стоящих на позициях i j используя элементы стоящие на позициях i j k
}
} else // если количество вызово рекурсии превышает количество строк в матрице значит всё вычислили и выходим из функции
return;
k++; //счётчик рекрсии
functionForTheMatrixMRekurs(M, n); //рекурсивный вызов функции
}
void functionForTheMatrixL(double **L, int n, double **M)
{
int i, j;
for(i=0; i<n; i++)
{
for (j=0; j < n; j++)
{
if (i == j)
L[i][j] = 1;
else if (i < j)
L[i][j] = 0;
else
L[i][j] = M[i][j];
}
}
}
void functionForTheMatrixU(double **U, int n, double **M)
{
int i, j;
for(i=0; i<n; i++)
{
for (j=0; j < n; j++)
{
if (i > j)
U[i][j] = 0;
else
U[i][j] = M[i][j];
}
}
}
void functionForTheMatrixL_and_U(double **A, int n, double **L, double **U) //function matrix for the calculation of M and L
{
int i, j, k;
for (j = 0; j < n; j++)
{
U[0][j] = A[0][j]; //первая строка U равна первой строки А
L[j][j] = 1; //единицы на диагональ L
}
for (i = 1; i < n; i++) //первый столбец L делится на первый элемент первой строки
L[i][0] = A[i][0] / U[0][0];
for (i = 1; i < n; i++) //для всех последующих надо учитывать предыдущий результат
{
for (j = 0; j < n; j++) //по столбцам
{
if (j < i)
U[i][j] = 0; //обнуляем нижний треугольник
else
{
U[i][j] = A[i][j]; //U(i,j) = A(i,j) - sum(L(i,k)*U(k,j)), k=0,...,i-1
for (k = 0; k < i; k++)
U[i][j] -= L[i][k] * U[k][j];
}
}
for (j = 0; j < i; j++) //элементы L выше диагонали = 0
L[j][i] = 0;
for (j = i+1; j < n; j++) //по строкам
{
L[j][i] = A[j][i]; //L(j,i) = [A(j,i) - sum(L(j,k)*U(k,i))] / U(i,i), k=0,...,i-1
for (k = 0; k < i; k++)
L[j][i] -= L[j][k] * U[k][i];
L[j][i] /= U[i][i];
}
}
}
void showMatrix(double **matrix, int n, int m)
{
int i, j;
cout.precision(4);
for (i = 0; i < n; i++) {
cout << endl;
for (j = 0; j < m; j++) {
cout << matrix[i][j] << "\t";
}
}
cout << endl;
}
int main()
{
performance perf;
int i, j;
double **A, **MR, **LR, **UR, **U, **L;
A = new double*[3];
MR = new double*[3];
LR = new double*[3];
UR = new double*[3];
L = new double*[3];
U = new double*[3];
for (i = 0; i < 3; i++)
{
A[i] = new double [3];
MR[i] = new double [3];
LR[i] = new double [3];
UR[i] = new double [3];
L[i] = new double [3];
U[i] = new double [3];
}
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++) // здесь мы вводим наши элементы
{
cin >> A[i][j];
MR[i][j] = A[i][j]; //сдублируем в MR
}
}
QueryPerfCounter(&perf);
Start(perf);
{
functionForTheMatrixMRekurs(MR, 3);
functionForTheMatrixL(LR, 3, MR);
functionForTheMatrixU(UR, 3, MR);
}
Stop(perf);
// Calculate time in mikroseconds
double d0 = Duration(&perf);
Start(perf);
{
functionForTheMatrixL_and_U(A, 3, L, U);
}
Stop(perf);
// Calculate time in mikroseconds
double d1 = Duration(&perf);
//вывод
showMatrix(A, 3, 3);
showMatrix(MR, 3, 3);
showMatrix(LR, 3, 3);
showMatrix(UR, 3, 3);
showMatrix(L, 3, 3);
showMatrix(U, 3, 3);
cout << "Recursive method: " << d0 << " mks" << endl;
cout << "Exception method: " << d1 << " mks" << endl;
for (i = 0; i < 3; i++)
delete [] A[i];
delete [] A;
for (i = 0; i < 3; i++)
delete [] MR[i];
delete [] MR;
for (i = 0; i < 3; i++)
delete [] LR[i];
delete [] LR;
for (i = 0; i < 3; i++)
delete [] UR[i];
delete [] UR;
for (i = 0; i < 3; i++)
delete [] U[i];
delete [] U;
for (i = 0; i < 3; i++)
delete [] L[i];
delete [] L;
return 0;
}
Консультировал: Лысков Игорь Витальевич (Старший модератор)
Дата отправки: 03.06.2012, 01:13
Команда портала RFPRO.RU благодарит Вас за то, что Вы пользуетесь нашими услугами. Вы только что прочли очередной выпуск рассылки. Мы старались.
Пожалуйста, оцените его. Если совет помог Вам, если Вам понравился ответ, Вы можете поблагодарить автора -
для этого в каждом ответе есть специальные ссылки. Вы можете оставить отзыв о работе портале. Нам очень важно знать Ваше мнение.
Вы можете поближе познакомиться с жизнью портала, посетив наш форум, почитав журнал,
который издают наши эксперты. Если у Вас есть желание помочь людям, поделиться своими знаниями, Вы можете зарегистрироваться экспертом.
Заходите - у нас интересно!