Category: Computer Science

using System;
using System.Diagnostics;
using System.Windows.Forms;

namespace MatrixInverseComparison
{
    public partial class MainForm : Form
    {
        public MainForm()
        {
            InitializeComponent();
        }
        static private string FormatNumber(double x)
        {
            string result = string.Empty;

            if (x > 0)
                result += x.ToString("F5").PadLeft(10);
            else
                result += x.ToString("F5").PadLeft(10);

            return result;
        }

        private void MultiplyPrintMatricies(
            double[,] A, double[,] B, int n)
        {
            double[,] I = new double[n, n];

            textBox1.Text += "Matrix Product:\r\n";

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    double sum = 0.0;

                    for (int k = 0; k < n; k++)
                        sum += A[i, k] * B[k, j];

                    I[i, j] = sum;
                }
            }

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    textBox1.Text += FormatNumber(I[i, j]) + " ";
                }

                textBox1.Text += "\r\n";
            }
        }

        private void PrintMatrix(string title, double[,] A, int n)
        {
            textBox1.Text += title + "\r\n";

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    textBox1.Text += FormatNumber(A[i, j]) + " ";
                }

                textBox1.Text += "\r\n";
            }
        }

        private void Compute(double[,] MI, int n)
        {
            double determinantGE = 1;
            double[] b = new double[n];
            double[] x = new double[n];
            double[,] MB1 = new double[n, n];
            double[,] MB2 = new double[n, n];
            double[,] MG = new double[n, n];
            double[,] MS = new double[n, n];
            double[,] IG = new double[n, n];
            double[,] IS = new double[n, n];
            int[] pivot = new int[n];
            Stopwatch sw = new();
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    MB1[i, j] = MB2[i, j] = MG[i, j] = MS[i, j] = MI[i, j];
            }
            PrintMatrix("Initial Matrix: ", MI, n);
            sw.Start();
            int flag = DirectMethods.Factor(MG, n, pivot);
            if (flag != 0)
            {
                for (int i = 0; i < n; i++)
                    determinantGE *= MG[i, i];
                determinantGE *= flag;
            }
            sw.Stop();
            PrintMatrix("Gaussian Elimination Final:", MG, n);
            for (int i = 0; i < n; i++)
                b[i] = 0;
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    b[j] = 1.0;
                    DirectMethods.Substitute(MG, b, x, n, pivot);
                    for (int k = 0; k < n; k++)
                        IG[k, j] = x[k];
                    b[j] = 0.0;
                }
            }
            PrintMatrix("Gaussian Elimination Inverse:", IG, n);
            textBox1.Text += "Determinant: " +
                FormatNumber(determinantGE) + "\r\n";
            MultiplyPrintMatricies(IG, MI, n);
            textBox1.Text += "Runtime (MS) = " +
                sw.ElapsedMilliseconds + "\r\n";
            sw.Start();
            double determinant1 = BareissAlgorithm.Determinant1(MB1, n);
            sw.Stop();
            PrintMatrix("Bareiss Algorithm Final 1:", MB1, n);
            textBox1.Text += "Determinant: " +
                FormatNumber(determinant1) + "\r\n";
            textBox1.Text += "Runtime (MS) = " +
                sw.ElapsedMilliseconds + "\r\n";
            sw.Start();
            double determinant2 = BareissAlgorithm.Determinant2(MB2, n);
            sw.Stop();
            PrintMatrix("Bareiss Algorithm Final 2:", MB2, n);
            textBox1.Text += "Determinant: " +
                FormatNumber(determinant2) + "\r\n";
            textBox1.Text += "Runtime (MS) = " +
                sw.ElapsedMilliseconds + "\r\n";
        }
        private void button1_Click(object sender, EventArgs e)
        {
            int n = (int)numericUpDown1.Value;
            int seed = (int)numericUpDown2.Value;
            double[,] A = new double[n, n];
            Random random = new Random(seed);

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    double q = n * random.NextDouble();

                    while (q == 0.0)
                        q = n * random.NextDouble();

                    A[i, j] = q;
                }
            }
            
            Compute(A, n);
        }

        private void button2_Click(object sender, EventArgs e)
        {
            textBox1.Text = string.Empty;
        }
    }
}

using System;

namespace MatrixInverseComparison
{
    public class DirectMethods
    {
        // Substitute and Factor translated from FORTRAN 77
        // source code found in "Elementary Numerical Analysis:
        // An Algorithmic Approach" by S. D. Conte and Carl de
        // Boor. Translator: James Pate Williams, Jr. (c)
        // August 14 - 17, 2023
        static public void Substitute(
            double[,] w,
            double[] b,
            double[] x,
            int n,
            int[] pivot)
        {
            double sum;
            int i, j, n1 = n - 1;

            if (n == 1)
            {
                x[0] = b[0] / w[0, 0];
                return;
            }

            // forward substitution
            x[0] = b[pivot[0]];

            for (i = 0; i < n; i++)
            {
                for (j = 0, sum = 0.0; j < i; j++)
                    sum += w[i, j] * x[j];
                x[i] = b[pivot[i]] - sum;
            }

            // backward substitution
            x[n1] /= w[n1, n1];

            for (i = n - 2; i >= 0; i--)
            {
                for (j = i + 1, sum = 0.0; j < n; j++)
                    sum += w[i, j] * x[j];
                x[i] = (x[i] - sum) / w[i, i];
            }
        }

        // Factor returns +1 if an even number of exchanges
        // Factor returns -1 if an odd  number of exchanges
        // Factor retrurn 0  if matrix is singular
        static public int Factor(
            double[,] w, int n, int[] pivot)
        // returns 0 if matrix is singular
        {
            double awikod, col_max, ratio, row_max, temp;
            double[] d = new double[n];
            int flag = 1, i, i_star, j, k;

            for (i = 0; i < n; i++)
            {
                pivot[i] = i;
                row_max = 0;
                for (j = 0; j < n; j++)
                    row_max = Math.Max(row_max, Math.Abs(w[i, j]));
                if (row_max == 0)
                {
                    flag = 0;
                    row_max = 1;
                }
                d[i] = row_max;
            }
            if (n <= 1)
                return flag;
            // factorization
            for (k = 0; k < n - 1; k++)
            {
                // determine pivot row the row i_star
                col_max = Math.Abs(w[k, k]) / d[k];
                i_star = k;
                for (i = k + 1; i < n; i++)
                {
                    awikod = Math.Abs(w[i, k]) / d[i];
                    if (awikod > col_max)
                    {
                        col_max = awikod;
                        i_star = i;
                    }
                }
                if (col_max == 0)
                    flag = 0;
                else
                {
                    if (i_star > k)
                    {
                        // make k the pivot row by
                        // interchanging with i_star
                        flag *= -1;
                        i = pivot[i_star];
                        pivot[i_star] = pivot[k];
                        pivot[k] = i;
                        temp = d[i_star];
                        d[i_star] = d[k];
                        d[k] = temp;
                        for (j = 0; j < n; j++)
                        {
                            temp = w[i_star, j];
                            w[i_star, j] = w[k, j];
                            w[k, j] = temp;
                        }
                    }
                    // eliminate x[k]
                    for (i = k + 1; i < n; i++)
                    {
                        w[i, k] /= w[k, k];
                        ratio = w[i, k];
                        for (j = k + 1; j < n; j++)
                            w[i, j] -= ratio * w[k, j];
                    }
                }
            }

            if (w[n - 1, n - 1] == 0)
                flag = 0;

            return flag;
        }
    }
}

namespace MatrixInverseComparison
{
    // One implementation is based on https://en.wikipedia.org/wiki/Bareiss_algorithm
    // Another perhaps better implementation is found on the following webpage
    // https://cs.stackexchange.com/questions/124759/determinant-calculation-bareiss-vs-gauss-algorithm
    class BareissAlgorithm
    {
        static public double Determinant1(double[,] M, int n)
        {
            double M00;

            for (int k = 0; k < n; k++)
            {
                if (k - 1 == -1)
                    M00 = 1;
                else
                    M00 = M[k - 1, k - 1];

                for (int i = k + 1; i < n; i++)
                {
                    for (int j = k + 1; j < n; j++)
                    {
                        M[i, j] = (
                            M[i, j] * M[k, k] -
                            M[i, k] * M[k, j]) / M00;
                    }
                }
            }

            return M[n - 1, n - 1];
        }

        static public double Determinant2(double[,] A, int dim)
        {
            if (dim <= 0)
            {
                return 0;
            }

            double sign = 1;

            for (int k = 0; k < dim - 1; k++)
            {
                //Pivot - row swap needed
                if (A[k, k] == 0)
                {
                    int m;
                    for (m = k + 1; m < dim; m++)
                    {
                        if (A[m, k] != 0)
                        {
                            double[] tempRow = new double[dim];

                            for (int i = 0; i < dim; i++)
                                tempRow[i] = A[m, i];
                            for (int i = 0; i < dim; i++)
                                A[m, i] = A[k, i];
                            for (int i = 0; i < dim; i++)
                                A[k, i] = tempRow[i];
                            sign = -sign;
                            break;
                        }
                    }
                    //No entries != 0 found in column k -> det = 0
                    if (m == dim)
                    {
                        return 0;
                    }
                }
                //Apply formula
                for (int i = k + 1; i < dim; i++)
                {
                    for (int j = k + 1; j < dim; j++)
                    {
                        A[i, j] = A[k, k] * A[i, j] - A[i, k] * A[k, j];

                        if (k >= 1)
                        {
                            A[i, j] /= A[k - 1, k - 1];
                        }
                    }
                }
            }

            return sign * A[dim - 1, dim - 1];
        }
    }
}