#1 10.11.10 22:23
Нахождение определителя матрицы через разложение по столбцам [C++]
Привет всем, вот такой у меня код, решающий эту задачу-
double Matrix_A(double **table1,int razm)
{
extern int nom;
int n_fakt=1;
for (int v=1;v<razm;v++)
n_fakt=n_fakt*v;
double ***table = new double** [n_fakt];
for(int c=0;c<razm;c++){table[c]=new double*[razm];for(int x=0;x<razm;x++){table[c][x]=new double[razm];};};
for(int c1=0;c1<razm;c1++){for(int c2=0;c2<razm;c2++){table[0][c1][c2]=table1[c1][c2];};};
double A=Minor(table, razm, 0, 0);
return(A);};
//=============
double Minor(double ***table, int razm, int i_minor,int j_minor)
{
int znak=1;
int vnut_nom;
vnut_nom=nom;
double A=0;
if(razm==0)return(0);
if(razm==1)return(table[nom][0][0]);
for(int q=0;q<(j_minor+i_minor);q++){znak=znak*(-1);};
for(int M=0;M<razm;M++)
{
vnut_nom++;
for(int b=1;b<razm;b++)
{
for(int c=0;c<razm;c++)
{
if(c==M)c++;
int n=0;
table[vnut_nom][b][n]=table[nom][b][c];
n++;
}
}
}
for(int z=0;z<razm;z++)
A=A+znak*table[0][i_minor+1][j_minor+1+z]*Minor(table[nom-razm+1+z][0][0],razm-1,i_minor+1,j_minor+1+z);
return(A);
};
Немного опишу её, а то вникнуть сложно сразу - здесь функция Matrix_A получает двойной массив и его размерность (предполагается что матрица квадратная)и создаёт тройной массив, где по адресу table[0][][] хранится исходная матрица, а в остальных все возможные её миноры (матрица разбивается на n! миноров 1х1) потом эта функция вызывает подфункцию - Minor которая разбивает получиную матрицу поочерёдно исключая столбцы и строки на все возможные миноры и опять сама себя вызывает передавая эти миноры как параметры и так пока не получится матрица 1х1, когда получилось считается формула. Так вот у меня возникла проблема с передачей этих массивов...даже боюсь как бы не была эта проблема только верхушкой айсберга - ХЕЛП ПЛИЗ! Помогите довести подпрограмму до ума, ну или,накрай, предложите другой способ нахождения определителя матрицы (а лучше минора, чтобы можно было и обратную потом написать).
Исправлено NoobsEnslaver (10.11.10 22:24)
Offline
#2 11.11.10 03:09
Re: Нахождение определителя матрицы через разложение по столбцам [C++]
NoobsEnslaver, есть два алгоритма для нахождения обратной матрицы, сам ими полбзуюсь:
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/* Функции для работы с матрицами и векторами из библиотеки uBLAS
* Поиск обратной матрицы первый способ
* Поиск обратной матрицы второй способ*/
#ifndef INVERT_MATRIX_H
#define INVERT_MATRIX_H
#include "boost/numeric/ublas/vector.hpp"
#include "boost/numeric/ublas/vector_proxy.hpp"
#include "boost/numeric/ublas/matrix.hpp"
#include "boost/numeric/ublas/triangular.hpp"
#include "boost/numeric/ublas/lu.hpp"
#include "boost/numeric/ublas/io.hpp"
namespace ublas = boost::numeric::ublas;
/* Matrix inversion routine.
Uses lu_factorize and lu_substitute in uBLAS to invert a matrix */
template<typename T1, typename T2>
bool InvertMatrix (const ublas::matrix<T1>& input, ublas::matrix<T2>& inverse)
{
using namespace boost::numeric::ublas;
typedef permutation_matrix<std::size_t> pmatrix;
// create a working copy of the input
matrix<T2> A(input);
// create a permutation matrix for the LU-factorization
pmatrix pm(A.size1());
// perform LU-factorization
int res = lu_factorize(A,pm);
if( res != 0 ) return false;
// create identity matrix of "inverse"
inverse.assign(ublas::identity_matrix<T2>(A.size1()));
// backsubstitute to get the inverse
lu_substitute(A, pm, inverse);
return true;
};
/* Matrix inversion routine.
Uses lu_factorize and lu_substitute in uBLAS to invert a matrix */
template<typename T>
bool InvertMatrix (ublas::matrix<T>& inverse)
{
using namespace boost::numeric::ublas;
typedef permutation_matrix<std::size_t> pmatrix;
// create a working copy of the input
matrix<T> A(inverse);
// create a permutation matrix for the LU-factorization
pmatrix pm(A.size1());
// perform LU-factorization
int res = lu_factorize(A,pm);
if( res != 0 ) return false;
// create identity matrix of "inverse"
inverse.assign(ublas::identity_matrix<T>(A.size1()));
// backsubstitute to get the inverse
lu_substitute(A, pm, inverse);
return true;
};
/* Поиск обратной матрицы
используется метод Фаддеева для массивов значений */
template<typename T>
bool faddeev_method(T** R, T** A, int n)
{
int i, j, k, m;
T p, s;
T** B = new T* [n];
T** U = new T* [n];
for(i = 0; i < n; ++i)
{
B[i] = new T[n];
U[i] = new T[n];
}
for(i = 0; i < n; ++i)
{
for(j = 0; j < n; ++j)
{
B[i][j] = 0.0;
}
B[i][i] = 1.0;
}
for(k = 1; k <= n; ++k) // общий цикл вычислений
{
for(i = 0; i < n; ++i) // U = A * B
{
for(j = 0; j < n; ++j)
{
s = 0.0;
for(m = 0; m < n; m++) s += A[i][m] * B[m][j];
U[i][j] = s;
}
}
s = 0;
for(i = 0; i < n; ++i)
{
s += U[i][i]; // след матрицы U
}
p = s / k; // характеристический коэффициент
if(k == n) break; // на последней итерации не вычислять
for(i = 0; i < n; ++i) // B = U - pE
{
for(j = 0; j < n; ++j)
{
i != j ? B[i][j] = U[i][j] : B[i][j] = U[i][j]-p;
}
}
}
bool z = true;// предположим, что обратная матрица существует
if(p == 0.0) z = false; // обратная матрица не существует
for(i = 0; i < n; ++i) // обратная матрица
{
for(j = 0; j < n; ++j)
{
R[i][j] = B[i][j] / p;
}
}
for(i = 0; i < n; ++i)
{
delete [] U[i]; // удалить строки U
}
delete [] U; // удалить указатели строк матрицы U
for(i = 0; i < n; ++i)
{
delete [] B[i]; // удалить строки B
}
delete [] B; // удалить указатели строк матрицы B
return z;
};
/* Поиск обратной матрицы
используется метод Фаддеева для класса типа матрица из uBLAS */
template<typename T>
bool faddeev_method (ublas::matrix<T>& result, const ublas::matrix<T>& matr/*, int n*/)
{
using namespace boost::numeric::ublas;
typedef permutation_matrix<std::size_t> pmatrix;
const int n = matr.size1();
matrix<T> B(n,n);
matrix<T> U(n,n);
int i, j, k;
T p, s;
for(i = n; --i >= 0;)
{
for(j = n; --j >= 0;)
{
B(i, j) = 0.0;
};
B(i, i) = 1.0;
};
for(k = 1; k <= n; ++k) // общий цикл вычислений
{
axpy_prod(matr, B, U, true); // U = matr * B
s = 0;
for(i = 0; i < n; ++i)
{
s += U(i, i);
};
p = s / k; // характеристический коэффициент
if(k == n) break; // на последней итерации не вычислять
B = U; // B = U - pE
for(i = 0; i < n; ++i)
{
B(i, i) = U(i, i) - p;
};
};
bool z = true;// предположим, что обратная матрица существует
if(p == 0.0) z = false; // обратная матрица не существует
result = B / p;
return z;
};
/* Поиск обратной матрицы
используется метод Фаддеева для класса типа матрица из uBLAS */
template<typename T>
bool faddeev_method (ublas::matrix<T>& matr)
{
using namespace boost::numeric::ublas;
typedef permutation_matrix<std::size_t> pmatrix;
const int n = matr.size1();
matrix<T> B(n,n);
matrix<T> U(n,n);
int i, j, k;
T p, s;
for(i = n; --i >= 0;)
{
for(j = n; --j >= 0;)
{
B(i, j) = 0.0;
};
B(i, i) = 1.0;
};
for(k = 1; k <= n; ++k) // общий цикл вычислений
{
axpy_prod(matr, B, U, true); // U = matr * B
s = 0;
for(i = 0; i < n; ++i)
{
s += U(i, i);
};
p = s / k; // характеристический коэффициент
if(k == n) break; // на последней итерации не вычислять
B = U; // B = U - pE
for(i = 0; i < n; ++i)
{
B(i, i) = U(i, i) - p;
};
};
bool z = true;// предположим, что обратная матрица существует
if(p == 0.0) z = false; // обратная матрица не существует
matr = B / p;
return z;
};
#endif //INVERT_MATRIX_H
первая вторая четвертая и пятая шаблоны для матриц из uBLAS.
третий шаблон для обычных массивов например double**. первый параметр - результат второй - исходная матрица третий - размерность. можешь изменить функцию так чтобы результат хранился на месте исходной матрицы.
все функции работают отлично, предпочитаю метод фаддеева - и работает вроде пошустрее и матрица обратная считается поточнее.
Исправлено Flinn (11.11.10 03:11)
Offline
