高斯消元法求逆矩阵

有多组测试数据。每组测试数据先输入一个整数n，表示方阵的阶。然后下面输入n阶方阵。输出其逆矩阵。若无逆矩阵，则输出No inverse matrix。

 

#include <iostream>
#include <cmath>
#include <algorithm>

using namespace std;

const double eps = 1e-6;

bool is_zero( const double num )
{
    return fabs(num) < eps;
}

void create( double ** & matrix, const int n )
{
    matrix = new double* [n];
    for ( int i = 0; i < n; ++i )
        matrix[i] = new double[n];
}

void input ( double ** matrix, const int n )
{
    for ( int i = 0; i < n; ++i )
    {
        for ( int j  = 0; j < n; ++ j )
            cin >> matrix[i][j];
    }
}

bool inverse ( double ** matrix1, double ** matrix2, const int n )
{
    int i, j;
    for ( i = 0; i < n; ++ i )
    {
        for ( j = 0; j < n; ++ j )
        {
            if ( i == j )
                matrix2[i][j] = 1;
            else
                matrix2[i][j] = 0;
        }
    }
    for ( i = 0; i < n; ++i )
    {
        int rowmaxpos = i;
        for ( j = i + 1; j < n; ++j )
        {
            if ( matrix1[i][j] > matrix1[i][rowmaxpos] )
                rowmaxpos = j;
        }
        for ( j = i; j < n; ++ j )
        {
            swap( matrix1[j][rowmaxpos], matrix1[j][i]);
            swap( matrix2[j][rowmaxpos], matrix2[j][i]);
        }
        if ( !is_zero(matrix1[i][i]) )
        {
            int divisor = matrix1[i][i];
            for ( j = i; j < n; ++ j )
            {
                matrix1[i][j] /= divisor;
                matrix2[i][j] /= divisor;
            }
            for ( j = i + 1; j < n; ++ j )
            {
                int multiple = matrix1[j][i];
                for ( int k = i; k < n; ++ k )
                {
                    matrix1[i][j] -= matrix1[i][k] * multiple;
                    matrix2[i][j] -= matrix2[i][k] * multiple;
                }
            }
        }
        else
            return false;
    }
    return true;
}

void output( double ** matrix, const int n )
{
    for ( int i = 0; i < n; ++i )
    {
        for ( int j = 0; j < n; ++ j )
            cout << matrix[i][j] << ' ';
        cout<<endl;
    }
}

void destroy( double ** matrix, const int n )
{
    for ( int i = 0; i < n; ++ i )
        delete [] matrix[i];
    delete [] matrix;
}

int main()
{
    int n;
    double ** matrix1;
    double ** matrix2;
    while ( cin >> n )
    {
        create( matrix1, n );
        create( matrix2, n );
        input( matrix1, n);
        if ( inverse(matrix1, matrix2, n) )
            output( matrix2, n );
        else
            cout << "No inverse matrix" << endl;
        destroy( matrix1, n );
        destroy( matrix2, n );
    }
    return 0;
}