GaussJ_procedure_from_Numerical_Recipes

using System;
using System.Numerics;
using ExtendedNumerics;
namespace NamespaceName
{
    public class ClassName
    {
        public static void GaussJ(BigRational[,] a,out byte error)
        {
            /*
             * Matrix inversion by Gauss-Jordan elimination
             * The input matrix a[0..n-1,0..n-1] has n by n elements
             * On output, a is replaced by its matrix inverse
             */
            if(a.GetLength(0) != a.GetLength(1))
            {
                error = 2;
                return;
            }
            BigRational big,dum,pivinv;
            int n,icol,irow;
            int[] indxc,indxr,ipiv;
            /*The integer arrays ipiv[0..n-1],indxr[0..n-1] and indxc[0..n-1]
             * are used for bookkeeping on the pivoting
             */
            n = a.GetLength(0);
            error = 0;

            indxc = new int[n];
            indxr = new int[n];
            ipiv = new int[n];
            for(int j=0;j<n;j++) ipiv[j] = 0;
            for(int i=0;i<n;i++){   //This is main loop over columns to be reduced
                big = 0;
                icol = -1;
                irow = -1;
                for(int j=0;j<n;j++) //This is tho outer loop of the searchh for a povot element
                    if(ipiv[j] != 1)
                        for(int k=0;k<n;k++)
                            if(ipiv[k]==0)
                                if(BigRational.Abs(a[j,k]) >= big){
                                    big = BigRational.Abs(a[j,k]);
                                    irow = j;
                                    icol = k;
                                }
                                else if(ipiv[k] > 1){
                                    Console.WriteLine("Pause 1 in GAUSSJ - singular matrix");
                                    Console.ReadKey();
                                    error = 1;
                                    return;
                                }
                            ipiv[icol]++;
                            /*
                             * We now have the pivot element,so we interchange rows,if needed to put the pivot element on the
                             * diagonal The columns are not phisically interchanged, only relabeled: indxc[i] , the column of the
                             * ith pivot element, is the ith column that is reduced, while indxr[i]  is the row in which that pivot
                             * element was originally located. If indxr[i] != indxc[i] there is implied column interchange. With
                             * this form of bookkeeping the inverse matrix will be scrambled by columns
                             */
                            if(irow != icol){
                                for(int l = 0;l<n;l++){
                                    dum = a[irow,l];
                                    a[irow,l] = a[icol,l];
                                    a[icol,l] = dum;
                                }
                            }
                            indxr[i] = irow; //We are now ready to divide the pivot row  by the pivot element, located at irow and icol
                            indxc[i] = icol;
                            if(a[icol,icol] == 0){
                                Console.WriteLine("Pause 2 in GAUSSJ - singular matrix");
                                Console.ReadKey();
                                error = 1;
                                return;
                            }
                pivinv = 1/a[icol,icol];
                a[icol,icol] = 1;
                for(int l=0;l<n;l++)
                    a[icol,l] *= pivinv;
                for(int ll=0;ll<n;ll++) //Next we reduce rows...
                    if(ll != icol){     // ...except fotr the pivot one, of course
                        dum = a[ll,icol];
                        a[ll,icol] = 0;
                        for(int l=0;l<n;l++)
                            a[ll,l] -= a[icol,l]*dum;
                        }
                }
                /*
                 * This is the end of the main loop over columns of reduction. It only remains to unscramble the solution
                 * in view of the column interchanges. We do this by interchanging pairs of columns in reverse order
                 * that permutation was built up
                 */
                for(int l=n-1;l>=0;l--)
                    if(indxr[l] != indxc[l])
                        for(int k=0;k<n;k++){
                            dum = a[k,indxr[l]];
                            a[k,indxr[l]] = a[k,indxc[l]];
                            a[k,indxc[l]] = dum;
                        }

        }//And we are done

        public static void PrintMatrix(BigRational[,] a)
        {
            int m = a.GetLength(0),n = a.GetLength(1);
            for(int i = 0;i < m;i++)
            {
                for(int j = 0;j < m;j++)
                    Console.Write("{0} , ",a[i,j].GetImproperFraction());
                Console.WriteLine();
            }
            Console.WriteLine();

        }

        public static void Main(string[] args)
        {
            char esc;
            int n;
            BigRational[,] a;
            byte error;
            do
            {
                Console.WriteLine("Odwracanie macierzy metodą Gaussa-Jordana z pełnym wyborem elementu głównego");
                Console.WriteLine("Program korzysta z procedury zamieszczonej w Numerical Recipes in Pascal przetłumaczonej na C#");
                Console.WriteLine("Podaj stopien odwracanej macierzy");
                n = Convert.ToInt32(Console.ReadLine());
                a = new BigRational[n,n];
                Console.WriteLine("Wczytaj macierz do odwrócenia");
                for(int i=0;i<n;i++)
                {
                    Console.WriteLine("Wczytaj wiersz {0} ",i+1);
                    for(int j=0;j<n;j++)
                    {
                        Console.Write("a[{0},{1}] = ",i+1,j+1);
                        a[i,j] = BigRational.Parse(Console.ReadLine());
                    }
                }
                Console.WriteLine();
                Console.WriteLine("Wczytana macierz :");
                PrintMatrix(a);
                GaussJ(a,out error);
                if(error == 0)
                {
                    Console.WriteLine("Macierz odwrotna to:");
                    PrintMatrix(a);
                }
                esc = (char)Console.ReadKey(true).Key;
            }
            while(esc != (char)ConsoleKey.Escape);
        }
    }
}