using System;

namespace RealMatrixProcessor
{
	/// <summary>
	/// Operations with real matrix
	/// </summary>
	public class RealMatrix
	{
		public RealMatrix()
		{
			//
			// TODO: Add constructor logic here
			//
		}

		/// <summary>
		/// Inverts matrix
		/// </summary>
		/// <param name="a">Matrix to invert</param>
		/// <param name="a1">Result of inversion</param>
		public static void Invert(double[,] a, double[,] a1)
		{
			double e, y, det, w, d = 0, d1;
			int i, j, k, ir, ip, n = a.GetLength(0);
			int[] jz = new int[n];
			double[] c = new double[n];
			double[] ab = new double[n];
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					a1[i, j] = a[i, j];
				}
			}
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					e = Math.Abs(a1[i, j]);
					if (d < e) d = e;
				}
			}
			d1 = 1 / d;
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					a1[i, j] *= d1;
				}
			}
			e = 1.0e-26;
			det=1;
			for (i = 0; i < n;  i++)
			{
				jz[i]=i;
			}
			for (i = 0; i < n; i++)
			{
				k = i;
				y = a1[i, i];
				ir = i - 1;
				ip = i + 1;
				if (ip < n)
				{
					for (j = ip; j < n; j++)
					{
						w = a1[i, j];
						if (Math.Abs(w) > Math.Abs(y))
						{
							k = j;
							y = w;
						}
					}
				}
				det *= y;
				y = 1.0 / y;
				for (j = 0; j < n; j++)
				{
					c[j] = a1[j, k];
					a1[j, k] = a1[j, i];
					a1[j, i] = -c[j] * y;
					ab[j] = a1[i, j] * y;
					a1[i, j] = ab[j];
				}
				a1[i, i] = y;
				j = jz[i];
				jz[i] = jz[k];
				jz[k] = j;
				k = 0;
				do
				{
					if ((k <= ir) | (k >= ip))
					{
						j = 0;
						do
						{
							if ((j <= ir) | (j >= ip))
							{
								a1[k, j] -= ab[j] * c[k];
							}
							j++;
						}
						while (j < n);
					}
					k++;
				}
				while(k < n);
			}
			i = 0;
			do
			{
				k = jz[i];
				if (k != i)
				{
					for(j = 0; j < n; j++)
					{
						w = a1[i, j];
						a1[i, j] = a1[k, j];
						a1[k, j] = w;
					}
					ip = jz[i];
					jz[i] = jz[k];
					jz[k] = ip;
					det = -det;
				}
				else
				{
					i++;
				}
			}
			while (i < n);
			d1= 1.0 / d;
			for (i = 0; i < n; i++)
			{
				for (j = 0; j < n; j++)
				{
					a1[i, j] *= d1;
				}
			}
		}

		static public double Det(double[,] a)
		{
			double[,] A = new double[a.GetLength(0),a.GetLength(1)];
			for (int ii = 0; ii < a.GetLength(0); ii++)
			{
				for (int jj = 0; jj < a.GetLength(1); jj++)
				{
					A[ii,jj] = a[ii,jj];
				}
			}
			double MAX = 0;
			double D = 1;
			double T = 0;
			int k, i, j = 0;
			double z = 0;
			int n = a.GetLength(0);
			for (k = 0; k < n; k++) 
			{
				MAX = 0;
				for (i = k; i < n; i++) 
				{
					T = A[i, k];
					if (!(T == 0)) 
					{
						MAX = T;
						j = i;
						goto calc;
					}
				}
			calc:
				if (MAX == 0)
				{
					return z;
				}
				if (j != k) 
				{
					D =- D;
					for (i = k; i < n; i++) 
					{
						T = A[j, i];
						A[j, i] = A[k, i];
						A[k, i] = T;						
					}
				}
				for (i = k + 1; i < n; i++) 
				{ 
					T = A[i, k] / MAX;
					for (j = k + 1; j < n; j++)
					{
						A[i, j] = A[i, j] - T * A[k, j];
					}
				}
				D = D * A[k, k];
			}
			return  D;
		}

		static public void Multiply(double[,] a, double[,] b, double [,] c)
		{
			if ((a.GetLength(1) != b.GetLength(0)) | (a.GetLength(0) != c.GetLength(0)) | 
				(b.GetLength(1) != c.GetLength(1)))
			{
				throw new Exception("Illegal matrix product dimension");
			}
			int i, j, k = 0;
			for (i = 0; i < c.GetLength(0); i++)
			{
				for (j = 0; j < c.GetLength(1); j++)
				{
					c[i,j] = 0;
				}
			}

			for (i = 0; i < a.GetLength(0); i++)
			{
				for (j = 0; j < b.GetLength(1); j++)
				{
					for (k = 0; k < a.GetLength(1); k++)
					{
						c[i,j] += a[i, k] * b[k, j];
					}
				}
			}
		}

		static public void Multiply(double[,] a, double[] b, double[] c)
		{
			if ((a.GetLength(1) != b.Length) | (a.GetLength(0) != c.Length))
			{
				throw new Exception("Illegal dimension of vector or matrix product");
			}
			int i, j = 0;
			for (i = 0; i < c.Length; i++)
			{
                c[i] = 0;
			}
			for (i = 0; i < a.GetLength(0); i++)
			{
				for (j = 0; j < b.Length; j++)
				{
					c[i] += a[i, j] * b[j];
				}
			}
		}

		static public void Multiply(double[] b, double[,] a, double[] c)
		{
			if ((a.GetLength(0) != b.Length) | (a.GetLength(1) != c.Length))
			{
				throw new Exception("Illegal dimension of vector or matrix product");
			}
			int i, j = 0;
			for (i = 0; i < c.Length; i++)
			{
				c[i] = 0;
			}
			for (i = 0; i < a.GetLength(1); i++)
			{
				for (j = 0; j < b.Length; j++)
				{
					c[i] += a[i, j] * b[j];
				}
			}
		}

	}
}
 

