Метод LU-разложения. Показать, что если λ– собственное число матрицы A, то 1/λ– собственное число обратной матрицы A

//Задача. Написать программу, которая для заданной матрицы A = (aij) nxm (2<n<11) позволит определить собственные значения и проверить утверждение, сформулированное в варианте.
//Метод разложения LU. Доказать, что если λ является собственным значением матрицы A, то 1/λ является обратным собственным значением матрицы A.
 
#include <iostream>
#include <cmath>
#include <vector>
 
using namespace std;
 
//DECLARING FUNCTIONS 
//ОБЪЯВЛЕНИЕ ФУНКЦИЙ 
 
float exist(int n, int m); //check of n and m
//проверка n и m
 
vector< vector<float> > LU_decompos(vector< vector<float> > matr); //LU-decomposition
//LU-разложение
 
float det(vector< vector<float> > matr); //matrix determinant
//определитель матрицы
 
//MAIN 
//ОСНОВНАЯ ЧАСТЬ 
 
int main()
{
    int n, m;
    float x;
    bool real_neg;
 
    cout << "Enter matrix size (2<n<11)" << endl;
 
    cout << "n="; cin >> n;    
    cout << "m="; cin >> m;
 
    if (!exist(n, m)) //check of n and m
    {
        cout << "Error. Wrong enter" << endl;
        return 0;    
    }
 
    vector< vector<float> > matr(n, vector<float>(m));
 
    cout << "Enter matrix elements:" << endl;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++)
        {
            cin >> matr[i][j];
        }
 
    vector<float> y(n);
    float z;
 
    //LU_decomposition 
    //LU-разложение
    matr = LU_decompos(matr);
 
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
            x += -matr[i][j];
            
            y[i] = x + matr[i][i];
            x = 0;
    }
 
//determinant calculation 
//расчет определителя
 
float F = det(matr);
 
//determinant calculation complete 
//завершен расчет определителя
 
//eigenvalues and statement check 
//определение собственных значений и проверка утверждения
 
real_neg = false;
 
for (int i = 0; i < n; i++)
{
    if (abs(y[i]) < 0.0001)  
    {
        real_neg = true;
        break;
    }
 
    z = 1 / y[i];
    cout << "Eigenvalue: " <<y[i]<< ", inverse eigenvalue: " <<z<< endl;
}
 
if (real_neg)
{
    for (int i=0; i<n; i++)
    {
        if (abs(y[i]) < 0.00001)
        {
            cout << "Eigenvalue " << y[i] << " is equal to 0" << endl;
        }        
    }
} 
 
//proof result 
//результат доказательства 
 
if (F == 0)
{
    cout << "Statement is true" << endl;
} 
else 
{
    cout << "Statement is false" << endl;
}
 
return 0;
}
 
float exist(int n, int m)
{
    if ((n >= 2 && n <= 11) && (m >= 2 && m <= 11))
        return true;
    
    return false;
}
 
vector< vector<float> > LU_decompos(vector< vector<float> > matr)
{
    int n = matr.size();
    float mult;
 
    for (int k = 0; k < n; k++)
        for (int i = k + 1; i < n; i++)
        {
            mult = matr[i][k] / matr[k][k];
            matr[i][k] = mult;
 
        for (int j = k + 1; j < n; j++)
            matr[i][j] -= mult * matr[k][j];
        }    
    
    return matr;
}
 
float det(vector< vector<float> > matr)
{
    int n = matr.size();
    float det = 1;
 
    for (int i = 0; i < n; i++)
        det *= matr[i][i];
 
    return det;
}

for (int j = 0; j < m; j++)

#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;
 
double** Create( size_t n);               // Создание матрицы
void Free( double** M, size_t n );        // Удаление матрицы
void Input( double** M, size_t n);        // ввод матрицы
void Print( double** M, size_t n);        // Печать матрицы 
double** MatrixDuplicate(double**matrix,size_t n);// копирование матрицы
double** MatrixDecompose(double** matrix,int* perm,int &toggle,int n);// Разложение LUP Дулитла
double** MatrixInverse(double** matrix,int n);  // Обращение матрицы
double MatrixDeterminant(double** matrix,int n);// Определитель матрицы
double* SystemSolve(double** A,double* b,int n);// Решение систем линейных уравнений
double* HelperSolve(double** luMatrix,double* b,int n);//Вспомогательная функция для решения системы уравнений
 
 
//-------------------------------------
int main()
{
    int n;
    cout<<"n=";
    cin>>n;
    double** matr=Create(n);
    Input(matr,n);
    cout<<endl<<"Original matrix:"<<endl;
    Print(matr,n);
 
    int* perm=new int[n];
    int toggle;
    double** lum=MatrixDecompose(matr,perm,toggle,n);
    //cout<<endl<<"LU matrix:"<<endl;
    //Print(lum,n);
    cout<<endl<<"Eigen Values:"<<endl;
    for(int i=0;i<n;i++)
        cout<<lum[i][i]<<' ';
    double det=MatrixDeterminant(matr,n);
    cout<<endl<<"Determinant="<<setprecision(4)<<det<<endl;
    if(fabs(det)<1e-6)
        cout<<"Can't to find inverse matrix\n";
    else
    {
        double** minv=MatrixInverse(matr,n);
        cout<<endl<<"Inverse matrix:"<<endl;
        Print(minv,n);
        double** lumi=MatrixDecompose(minv,perm,toggle,n);
        cout<<endl<<"Eigen Values:"<<endl;
        for(int i=0;i<n;i++)
            cout<<lumi[i][i]<<' ';
        double deti=MatrixDeterminant(minv,n);
        cout<<endl<<"Inverse Determinant="<<setprecision(4)<<deti<<endl;
        Free(minv,n);
        Free(lumi,n);
    }
    delete[] perm;
    Free(matr,n);
    Free(lum,n);
    cin.get();
    return 0;
}
//-------- Печать матрицы ------------------------------------------------
void Print( double** M, size_t n) {
    for ( size_t i = 0; i < n; ++i ) {
        for ( size_t j = 0; j < n; ++j ) {
            cout<<setw(7)<<setprecision(3)<<M[i][j]<<' ';
        }
        cout<<endl;
    }
}
// Создание матрицы
double** Create( size_t n) {
    double** M = new double* [n];
    for ( size_t i = 0; i < n; ++i ) {
        M[i] = new double[n];
    }
    return M;
}
// Удаление матрицы
void Free( double** M, size_t n ) {
    for ( size_t i = 0; i < n; ++i ) {
        delete [] M[i];
    }
    delete [] M;
}
//---- ввод матрицы--------------------------------------
void Input( double** M, size_t n) {
    cout << "Enter matrix elements:\n";
    for ( size_t i = 0; i < n; ++i ) {
        for ( size_t j = 0; j < n; ++j ) {
            //cout << "M[" << i << "][" << j << "] = ";
            cin >> M[i][j];
        }
    }
    cin.get();
}
//-------копирование матрицы ------------------------------
double** MatrixDuplicate(double**matrix,size_t n)
{
    double** result=Create(n);
    for(size_t i=0;i<n;++i)
        for(size_t j=0;j<n;++j)
            result[i][j]=matrix[i][j];
    return result;
}
//----------------Разложение LUP Дулитла-----------------------
double** MatrixDecompose(double** matrix,int* perm,int &toggle,int n)
{
    double**result=MatrixDuplicate(matrix,n);
    for(int i=0;i<n;++i){perm[i]=i;}
    toggle=1;
    for(int j=0;j<n-1;++j)
    {
        double colMax=fabs(result[j][j]);
        int pRow=j;
        for(int i=j+1;i<n;++i)
        {
            if(result[i][j]>colMax)
            {
                colMax=result[i][j];
                pRow=i;
            }
        }
        if(pRow!=j)//перестановка строк
        {
            double* rowPtr=result[pRow];
            result[pRow]=result[j];
            result[j]=rowPtr;
 
            int tmp=perm[pRow];//Меняем информацию о перестановке
            perm[pRow]=perm[j];
            perm[j]=tmp;
 
            toggle=-toggle;//переключатель перестановки строк
        }
        if(fabs(result[j][j])<1.0E-20)
            return NULL;
        for(int i=j+1;i<n;++i)
        {
            result[i][j]/=result[j][j];
            for(int k=j+1;k<n;++k)
                result[i][k]-=result[i][j]*result[j][k];
        }
    }
    return result;
}
//--------------- Вспомогательная функция для решения системы уравнений   ------------------------
double* HelperSolve(double** luMatrix,double* b,int n)
{
    double* x=new double[n];
    for(int i=0;i<n;++i)
        x[i]=b[i];
    for(int i=1;i<n;++i)
    {
        double sum=x[i];
        for(int j=0;j<i;++j)
            sum-=luMatrix[i][j]*x[j];
        x[i]=sum;
    }
    x[n-1]/=luMatrix[n-1][n-1];
    for(int i=n-2;i>=0;--i)
    {
        double sum=x[i];
        for(int j=i+1;j<n;++j)
            sum-=luMatrix[i][j]*x[j];
        x[i]=sum/luMatrix[i][i];
    }
    return x;
}
//------------Обращение матрицы ---------------------
double** MatrixInverse(double** matrix,int n)
{
    double** result=MatrixDuplicate(matrix,n);
    int* perm=new int[n];
    int toggle;
    double** lum=MatrixDecompose(matrix,perm,toggle,n);
    if(lum==NULL)
    {
        cout<<"Unable to compute inverse\n";
        cin.get();
        return 0;
    }
    double* b=new double[n];
    for(int i=0;i<n;++i)
    {
        for(int j=0;j<n;++j)
        {
            if(i==perm[j])
                b[j]=1.0;
            else
                b[j]=0.0;
        }
        double* x=HelperSolve(lum,b,n);
        for(int j=0;j<n;++j)
            result[j][i]=x[j];
        delete[] x;
    }
    delete[] perm;
    delete[] b;
    Free(lum,n);
    return result;
}
//-------------Определитель матрицы -------------------------------
double MatrixDeterminant(double** matrix,int n)
{
    int* perm=new int[n];
    int toggle;
    double** lum=MatrixDecompose(matrix,perm,toggle,n);
    if(lum==NULL)
    {
        cout<<"Unable to compute Determinant\n";
        cin.get();
        return 0;
    }
    double result=toggle;
    for(int i=0;i<n;++i)
        result*=lum[i][i];
    delete[] perm;
    Free(lum,n);
    return result;
}
//-------------Решение систем линейных уравнений -------
double* SystemSolve(double** A,double* b,int n)
{
    //Решаем Ax=b
    int* perm=new int[n];
    int toggle;
    double** luMatrix=MatrixDecompose(A,perm,toggle,n);
    if(luMatrix==NULL)
    {
        cout<<"Unable to compute system\n";
        cin.get();
        return 0;
    }
    Print(luMatrix,n);
    double* bp=new double[n];
    for(int i=0;i<n;++i)
        bp[i]=b[perm[i]];
    double* x=HelperSolve(luMatrix,bp,n);
    delete[] perm;
    delete[] bp;
    Free(luMatrix,n);
    return x;
}