caculate matrix inverse by Gaussian-Jordan elimination

main.cpp

//
//  main.cpp
//  imat
//
//  Created by LARRYHOU on 2021/11/1.
//

#include <iostream>

using Float = float;

void print(Float *v, int n, int m)
{
    auto p = v;
    for (auto i = 0; i < n; i++)
    {
        for (auto j = 0; j < m; j++) { printf("%8.4f ", *p++); }
        printf("\n");
    }
    printf("\n");
}

// Gaussian-Jordan elimination
bool inverse(Float *v, Float *o, int n)
{
    const auto m = n << 1;
    
    Float t[n][m];
    for (auto i = 0; i < n; i++)
    {
        for (auto j = 0; j < n; j++) { t[i][j] = v[i*n+j]; }
        for (auto j = n; j < m; j++) { t[i][j] = j-n == i ? 1 : 0; }
    }
    print(t[0], n, m);
    
    // Gaussian elimination
    for (auto i = 0; i < n; i++)
    {
        int x = i;
        for (auto r = i+1; r < n; r++) { if (abs(t[r][i]) > abs(t[x][i])) { x = r; } }
        if (x != i) {
            for (auto j = i; j < m; j++) { std::swap(t[i][j], t[x][j]); }
        }
        
        for (auto r = i+1; r < n; r++) {
            if (t[r][i] == 0) {continue;}
            for (auto j = m-1; j >=i; j--) {
                t[r][j] = t[r][j] * t[i][i] / t[r][i] - t[i][j];
            }
            print(t[0], n, m);
        }
        if (t[i][i] == 0) {return false;}
    }
    
    // Jordan elimination
    for (auto i = n-1; i >= 0; i--)
    {
        for (auto j = i+1; j < n; j++)
        {
            for (auto c = m-1; c >= j; c--) { t[i][c] -= t[i][j] * t[j][c]; }
        }
        print(t[0], n, m);
        for (auto j = m-1; j >= i; j--) { t[i][j] /= t[i][i]; }
        print(t[0], n, m);
    }
    
    auto p = o;
    auto s = n * sizeof(Float);
    for (auto i = 0; i < n; i++,p+=n) { memcpy(p, t[i]+n, s); }
    print(o, n, n);
    return true;
}

int main(int argc, const char * argv[]) {
    {
        Float M[3*3] = {2,0,0,0,4,0,0,0,5};
        Float O[sizeof(M)/sizeof(Float)];
        inverse(M, O, 3);
    }
    
    {
        Float M[3*3] = {2,-1,0,-1,2,-1,0,-1,2};
        Float O[sizeof(M)/sizeof(Float)];
        inverse(M, O, 3);
    }
    
    {
        Float M[4*4] = {1,0,0,30,0,1,0,20,0,0,1,10,0,0,0,1};
        Float O[sizeof(M)/sizeof(Float)];
        inverse(M, O, 4);
    }
    
    {
        Float M[4*4] = {1,3,1,9,1,1,-1,20,0,0,1,10,0,0,0,1};
        Float O[sizeof(M)/sizeof(Float)];
        inverse(M, O, 4);
    }
    
    return 0;
}