omp/MatrixOps.cs

namespace omp
{
    public static class MatrixOps
    {
        public static void Print(double[,] a)
        {
            int rows = a.GetLength(0);
            int columns = a.GetLength(1);

            Console.WriteLine("Matrix of size {0}x{1}:", rows, columns);
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < columns; j++)
                {
                    Console.Write("\t{0}", a[i, j]);
                }
                Console.WriteLine("");
            }
        }

        public static double Adj(double[,] a, int i, int j)
        {
            double sgn = ((i + j) % 2 == 0) ? 1 : -1;

            return sgn * Det(Minor(a, j, i));
        }

        public static double[,] Minor(double[,] a, int i, int j)
        {
            int rows = a.GetLength(0);
            int columns = a.GetLength(1);

            double[,] minor = new double[rows - 1, columns - 1];

            int k = 0;
            int l;

            for (int r = 0; r < rows; ++r)
            {
                if (r == i)
                {
                    continue;
                }

                l = 0;

                for (int c = 0; c < columns; ++c)
                {
                    if (c == j)
                    {
                        continue;
                    }

                    minor[k, l] = a[r, c];

                    l++;
                }

                k++;
            }

            return minor;
        }

        public static double Det(double[,] a)
        {
            int n = a.GetLength(0);

            switch (n)
            {
                case 2:
                    return Det2(a);
                case 3:
                    return Det3(a);
                default:
                    double det = 0;

                    double sgn;
                    for (int j = 0; j < n; ++j)
                    {
                        sgn = (j % 2 == 0) ? 1 : -1;
                        det += sgn * a[0, j] * Det(Minor(a, 0, j));
                    }

                    return det;
            }
        }

        public static double Det2(double[,] a)
        {
            return a[0, 0] * a[1, 1] - a[0, 1] * a[1, 0];
        }

        public static double Det3(double[,] a)
        {
            if (a.GetLength(0) != 3 || a.GetLength(1) != 3)
                throw new ArgumentException("Matrix must be 3x3");

            return
                a[0, 0] * a[1, 1] * a[2, 2]
                + a[0, 1] * a[1, 2] * a[2, 0]
                + a[0, 2] * a[1, 0] * a[2, 1]
                - a[0, 2] * a[1, 1] * a[2, 0]
                - a[0, 0] * a[1, 2] * a[2, 1]
                - a[2, 2] * a[0, 1] * a[1, 0];
        }

        public static double[,] MultTrans(double[,] a)
        {
            int rows = a.GetLength(0);
            int columns = a.GetLength(1);

            double[,] r = new double[rows, rows];

            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < rows; j++)
                {
                    for (int k = 0; k < columns; k++)
                    {
                        r[i, j] += a[i, k] * a[j, k];
                    }
                }
            }

            return r;
        }
    }
}