Files
UltimateFishing/Assets/Scripts/Assembly-CSharp/CSML/Matrix.cs
2026-02-21 16:45:37 +08:00

2628 lines
54 KiB
C#

using System;
using System.Collections;
namespace CSML
{
public class Matrix
{
public enum DefinitenessType
{
PositiveDefinite = 0,
PositiveSemidefinite = 1,
NegativeDefinite = 2,
NegativeSemidefinite = 3,
Indefinite = 4
}
private ArrayList Values;
private int rowCount;
private int columnCount;
public int RowCount
{
get
{
return rowCount;
}
}
public int ColumnCount
{
get
{
return columnCount;
}
}
public virtual Complex this[int i, int j]
{
get
{
if (i > 0 && i <= rowCount && j > 0 && j <= columnCount)
{
return (Complex)((ArrayList)Values[i - 1])[j - 1];
}
throw new ArgumentOutOfRangeException("Indices must not exceed size of matrix.");
}
set
{
if (i <= 0 || j <= 0)
{
throw new ArgumentOutOfRangeException("Indices must be real positive.");
}
if (i > rowCount)
{
for (int k = 0; k < i - rowCount; k++)
{
Values.Add(new ArrayList(columnCount));
for (int l = 0; l < columnCount; l++)
{
((ArrayList)Values[rowCount + k]).Add(Complex.Zero);
}
}
rowCount = i;
}
if (j > columnCount)
{
for (int m = 0; m < rowCount; m++)
{
for (int n = 0; n < j - columnCount; n++)
{
((ArrayList)Values[m]).Add(Complex.Zero);
}
}
columnCount = j;
}
((ArrayList)Values[i - 1])[j - 1] = value;
}
}
public virtual Complex this[int i]
{
get
{
if (RowCount == 1)
{
return (Complex)((ArrayList)Values[0])[i - 1];
}
if (ColumnCount == 1)
{
return (Complex)((ArrayList)Values[i - 1])[0];
}
throw new InvalidOperationException("General matrix acces requires double indexing.");
}
set
{
if (rowCount == 1)
{
if (i > columnCount)
{
for (int j = 0; j < i - columnCount; j++)
{
((ArrayList)Values[0]).Add(Complex.Zero);
}
columnCount = i;
}
((ArrayList)Values[0])[i - 1] = value;
return;
}
if (columnCount == 1)
{
if (i > rowCount)
{
for (int k = 0; k < i - rowCount; k++)
{
Values.Add(new ArrayList(columnCount));
((ArrayList)Values[rowCount + k]).Add(Complex.Zero);
}
rowCount = i;
}
((ArrayList)Values[i - 1])[0] = value;
return;
}
throw new InvalidOperationException("Cannot access multidimensional matrix via single index.");
}
}
public Matrix()
{
Values = new ArrayList();
rowCount = 0;
columnCount = 0;
}
public Matrix(int m, int n)
{
rowCount = m;
columnCount = n;
Values = new ArrayList(m);
for (int i = 0; i < m; i++)
{
Values.Add(new ArrayList(n));
for (int j = 0; j < n; j++)
{
((ArrayList)Values[i]).Add(Complex.Zero);
}
}
}
public Matrix(int n)
{
rowCount = n;
columnCount = n;
Values = new ArrayList(n);
for (int i = 0; i < n; i++)
{
Values.Add(new ArrayList(n));
for (int j = 0; j < n; j++)
{
((ArrayList)Values[i]).Add(Complex.Zero);
}
}
}
public Matrix(Complex x)
{
rowCount = 1;
columnCount = 1;
Values = new ArrayList(1);
Values.Add(new ArrayList(1));
((ArrayList)Values[0]).Add(x);
}
public Matrix(Complex[,] values)
{
if (values == null)
{
Values = new ArrayList();
columnCount = 0;
rowCount = 0;
}
rowCount = (int)values.GetLongLength(0);
columnCount = (int)values.GetLongLength(1);
Values = new ArrayList(rowCount);
for (int i = 0; i < rowCount; i++)
{
Values.Add(new ArrayList(columnCount));
for (int j = 0; j < columnCount; j++)
{
((ArrayList)Values[i]).Add(values[i, j]);
}
}
}
public Matrix(Complex[] values)
{
if (values == null)
{
Values = new ArrayList();
columnCount = 0;
rowCount = 0;
}
rowCount = values.Length;
columnCount = 1;
Values = new ArrayList(rowCount);
for (int i = 0; i < rowCount; i++)
{
Values.Add(new ArrayList(1));
((ArrayList)Values[i]).Add(values[i]);
}
}
public Matrix(double x)
{
rowCount = 1;
columnCount = 1;
Values = new ArrayList(1);
Values.Add(new ArrayList(1));
((ArrayList)Values[0]).Add(new Complex(x));
}
public Matrix(double[,] values)
{
if (values == null)
{
Values = new ArrayList();
columnCount = 0;
rowCount = 0;
}
rowCount = (int)values.GetLongLength(0);
columnCount = (int)values.GetLongLength(1);
Values = new ArrayList(rowCount);
for (int i = 0; i < rowCount; i++)
{
Values.Add(new ArrayList(columnCount));
for (int j = 0; j < columnCount; j++)
{
((ArrayList)Values[i]).Add(new Complex(values[i, j]));
}
}
}
public Matrix(double[] values)
{
if (values == null)
{
Values = new ArrayList();
columnCount = 0;
rowCount = 0;
}
rowCount = values.Length;
columnCount = 1;
Values = new ArrayList(rowCount);
for (int i = 0; i < rowCount; i++)
{
Values.Add(new ArrayList(1));
((ArrayList)Values[i]).Add(new Complex(values[i]));
}
}
public Matrix(string matrix_string)
{
matrix_string = matrix_string.Replace(" ", string.Empty);
string[] array = matrix_string.Split(';');
rowCount = array.Length;
Values = new ArrayList(rowCount);
columnCount = 0;
for (int i = 0; i < rowCount; i++)
{
Values.Add(new ArrayList());
}
for (int j = 1; j <= rowCount; j++)
{
string[] array2 = array[j - 1].Split(',');
for (int k = 1; k <= array2.Length; k++)
{
this[j, k] = new Complex(Convert.ToDouble(array2[k - 1]));
}
}
}
public static Matrix E(int n, int j)
{
Matrix matrix = Zeros(n, 1);
matrix[j] = Complex.One;
return matrix;
}
public static Complex KroneckerDelta(int i, int j)
{
return new Complex(Math.Min(Math.Abs(i - j), 1));
}
public static Matrix ChessboardMatrix(int m, int n, bool even)
{
Matrix matrix = new Matrix(m, n);
if (even)
{
for (int i = 1; i <= m; i++)
{
for (int j = 1; j <= n; j++)
{
matrix[i, j] = KroneckerDelta((i + j) % 2, 0);
}
}
}
else
{
for (int k = 1; k <= m; k++)
{
for (int l = 1; l <= n; l++)
{
matrix[k, l] = KroneckerDelta((k + l) % 2, 1);
}
}
}
return matrix;
}
public static Matrix ChessboardMatrix(int n, bool even)
{
Matrix matrix = new Matrix(n);
if (even)
{
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
matrix[i, j] = KroneckerDelta((i + j) % 2, 0);
}
}
}
else
{
for (int k = 1; k <= n; k++)
{
for (int l = 1; l <= n; l++)
{
matrix[k, l] = KroneckerDelta((k + l) % 2, 1);
}
}
}
return matrix;
}
public static Matrix Zeros(int m, int n)
{
return new Matrix(m, n);
}
public static Matrix Zeros(int n)
{
return new Matrix(n);
}
public static Matrix Ones(int m, int n)
{
Matrix matrix = new Matrix(m, n);
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
((ArrayList)matrix.Values[i])[j] = Complex.One;
}
}
return matrix;
}
public static Matrix Ones(int n)
{
Matrix matrix = new Matrix(n);
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
((ArrayList)matrix.Values[i])[j] = Complex.One;
}
}
return matrix;
}
public static Matrix Identity(int n)
{
return Diag(Ones(n, 1));
}
public static Matrix Eye(int n)
{
return Identity(n);
}
public static Matrix VerticalConcat(Matrix A, Matrix B)
{
Matrix matrix = A.Column(1);
for (int i = 2; i <= A.ColumnCount; i++)
{
matrix.InsertColumn(A.Column(i), i);
}
for (int j = 1; j <= B.ColumnCount; j++)
{
matrix.InsertColumn(B.Column(j), matrix.ColumnCount + 1);
}
return matrix;
}
public static Matrix VerticalConcat(Matrix[] A)
{
if (A == null)
{
throw new ArgumentNullException();
}
if (A.Length == 1)
{
return A[0];
}
Matrix matrix = VerticalConcat(A[0], A[1]);
for (int i = 2; i < A.Length; i++)
{
matrix = VerticalConcat(matrix, A[i]);
}
return matrix;
}
public static Matrix HorizontalConcat(Matrix A, Matrix B)
{
Matrix matrix = A.Row(1);
for (int i = 2; i <= A.RowCount; i++)
{
matrix.InsertRow(A.Row(i), i);
}
for (int j = 1; j <= B.RowCount; j++)
{
matrix.InsertRow(B.Row(j), matrix.RowCount + 1);
}
return matrix;
}
public static Matrix HorizontalConcat(Matrix[] A)
{
if (A == null)
{
throw new ArgumentNullException();
}
if (A.Length == 1)
{
return A[0];
}
Matrix matrix = HorizontalConcat(A[0], A[1]);
for (int i = 2; i < A.Length; i++)
{
matrix = HorizontalConcat(matrix, A[i]);
}
return matrix;
}
public static Matrix Diag(Matrix diag_vector)
{
int num = diag_vector.VectorLength();
if (num == 0)
{
throw new ArgumentException("diag_vector must be 1xN or Nx1");
}
Matrix matrix = new Matrix(num, num);
for (int i = 1; i <= num; i++)
{
matrix[i, i] = diag_vector[i];
}
return matrix;
}
public static Matrix Diag(Matrix diag_vector, int offset)
{
int num = diag_vector.VectorLength();
if (num == 0)
{
throw new ArgumentException("diag_vector must be 1xN or Nx1.");
}
Matrix matrix = new Matrix(num + Math.Abs(offset), num + Math.Abs(offset));
num = matrix.RowCount;
if (offset >= 0)
{
for (int i = 1; i <= num - offset; i++)
{
matrix[i, i + offset] = diag_vector[i];
}
}
else
{
for (int j = 1; j <= num + offset; j++)
{
matrix[j - offset, j] = diag_vector[j];
}
}
return matrix;
}
public static Matrix TriDiag(Complex l, Complex d, Complex u, int n)
{
if (n <= 1)
{
throw new ArgumentException("Matrix dimension must greater than one.");
}
return Diag(l * Ones(n - 1, 1), -1) + Diag(d * Ones(n, 1)) + Diag(u * Ones(n - 1, 1), 1);
}
public static Matrix TriDiag(Matrix l, Matrix d, Matrix u)
{
int num = l.VectorLength();
int num2 = d.VectorLength();
int num3 = u.VectorLength();
if (num * num2 * num3 == 0)
{
throw new ArgumentException("At least one of the paramter matrices is not a vector.");
}
if (num != num3)
{
throw new ArgumentException("Lower and upper secondary diagonal must have the same length.");
}
if (num + 1 != num2)
{
throw new ArgumentException("Main diagonal must have exactly one element more than the secondary diagonals.");
}
return Diag(l, -1) + Diag(d) + Diag(u, 1);
}
public static Complex Dot(Matrix v, Matrix w)
{
int num = v.VectorLength();
int num2 = w.VectorLength();
if (num == 0 || num2 == 0)
{
throw new ArgumentException("Arguments need to be vectors.");
}
if (num != num2)
{
throw new ArgumentException("Vectors must be of the same length.");
}
Complex zero = Complex.Zero;
for (int i = 1; i <= num; i++)
{
zero += v[i] * w[i];
}
return zero;
}
public static Complex Fib(int n)
{
Matrix matrix = Ones(2, 2);
matrix[2, 2] = Complex.Zero;
return (matrix ^ (n - 1))[1, 1];
}
public static Matrix RandomGraph(int n)
{
Matrix matrix = Random(n, n);
return matrix - Diag(matrix.DiagVector());
}
public static Matrix RandomGraph(int n, double p)
{
Matrix matrix = new Matrix(n);
Random random = new Random();
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
if (i == j)
{
matrix[i, j] = Complex.Zero;
}
else if (random.NextDouble() < p)
{
matrix[i, j] = new Complex(random.NextDouble());
}
else
{
matrix[i, j] = new Complex(double.PositiveInfinity);
}
}
}
return matrix;
}
public static Matrix Random(int m, int n)
{
Matrix matrix = new Matrix(m, n);
Random random = new Random();
for (int i = 1; i <= m; i++)
{
for (int j = 1; j <= n; j++)
{
matrix[i, j] = new Complex(random.NextDouble());
}
}
return matrix;
}
public static Matrix Random(int n)
{
Matrix matrix = new Matrix(n);
Random random = new Random();
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
matrix[i, j] = new Complex(random.NextDouble());
}
}
return matrix;
}
public static Matrix Random(int n, int lo, int hi)
{
Matrix matrix = new Matrix(n);
Random random = new Random();
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
matrix[i, j] = new Complex(random.Next(lo, hi));
}
}
return matrix;
}
public static Matrix RandomZeroOne(int m, int n, double p)
{
Matrix matrix = new Matrix(m, n);
Random random = new Random();
for (int i = 1; i <= m; i++)
{
for (int j = 1; j <= n; j++)
{
if (random.NextDouble() <= p)
{
matrix[i, j] = Complex.One;
}
}
}
return matrix;
}
public static Matrix RandomZeroOne(int n, double p)
{
Matrix matrix = new Matrix(n, n);
Random random = new Random();
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
if (random.NextDouble() <= p)
{
matrix[i, j] = Complex.One;
}
}
}
return matrix;
}
public static Matrix Random(int m, int n, int lo, int hi)
{
Matrix matrix = new Matrix(m, n);
Random random = new Random();
for (int i = 1; i <= m; i++)
{
for (int j = 1; j <= n; j++)
{
matrix[i, j] = new Complex(random.Next(lo, hi));
}
}
return matrix;
}
public static Matrix Vandermonde(Complex[] x)
{
if (x == null || x.Length < 1)
{
throw new ArgumentNullException();
}
int num = x.Length - 1;
Matrix matrix = new Matrix(num + 1);
for (int i = 0; i <= num; i++)
{
for (int j = 0; j <= num; j++)
{
matrix[i + 1, j + 1] = Complex.Pow(x[i], j);
}
}
return matrix;
}
public static Matrix[] Floyd(Matrix adjacence_matrix)
{
if (!adjacence_matrix.IsSquare())
{
throw new ArgumentException("Expected square matrix.");
}
if (!adjacence_matrix.IsReal())
{
throw new ArgumentException("Adjacence matrices are expected to be real.");
}
int num = adjacence_matrix.RowCount;
Matrix matrix = adjacence_matrix.Clone();
Matrix matrix2 = new Matrix(num);
for (int i = 1; i <= num; i++)
{
for (int j = 1; j <= num; j++)
{
for (int k = 1; k <= num; k++)
{
double num2 = matrix[j, i].Re + matrix[i, k].Re;
if (num2 < matrix[j, k].Re)
{
matrix[j, k].Re = num2;
matrix2[j, k].Re = i;
}
}
}
}
return new Matrix[2] { matrix, matrix2 };
}
public static ArrayList FloydPath(Matrix P, int i, int j)
{
if (!P.IsSquare())
{
throw new ArgumentException("Path matrix must be square.");
}
if (!P.IsReal())
{
throw new ArgumentException("Adjacence matrices are expected to be real.");
}
ArrayList arrayList = new ArrayList();
arrayList.Add(i);
while (P[i, j] != 0.0)
{
i = Convert.ToInt32(P[i, j]);
arrayList.Add(i);
}
arrayList.Add(j);
return arrayList;
}
public static Matrix DFS(Matrix adjacence_matrix, int root)
{
if (!adjacence_matrix.IsSquare())
{
throw new ArgumentException("Adjacence matrices are expected to be square.");
}
if (!adjacence_matrix.IsReal())
{
throw new ArgumentException("Adjacence matrices are expected to be real.");
}
int num = adjacence_matrix.RowCount;
if (root < 1 || root > num)
{
throw new ArgumentException("Root must be a vertex of the graph, e.i. in {1, ..., n}.");
}
Matrix matrix = new Matrix(num);
bool[] array = new bool[num + 1];
Stack stack = new Stack();
stack.Push(root);
array[root] = true;
ArrayList[] array2 = new ArrayList[num + 1];
for (int i = 1; i <= num; i++)
{
array2[i] = new ArrayList();
for (int j = 1; j <= num; j++)
{
if (adjacence_matrix[i, j].Re != 0.0 && adjacence_matrix[i, j].Im != double.PositiveInfinity)
{
array2[i].Add(j);
}
}
}
while (stack.Count > 0)
{
int num2 = (int)stack.Peek();
if (array2[num2].Count > 0)
{
int num3 = (int)array2[num2][0];
if (!array[num3])
{
array[num3] = true;
matrix[num2, num3].Re = 1.0;
stack.Push(num3);
}
array2[num2].RemoveAt(0);
}
else
{
stack.Pop();
}
}
return matrix;
}
public static Matrix BFS(Matrix adjacence_matrix, int root)
{
if (!adjacence_matrix.IsSquare())
{
throw new ArgumentException("Adjacence matrices are expected to be square.");
}
if (!adjacence_matrix.IsReal())
{
throw new ArgumentException("Adjacence matrices are expected to be real.");
}
int num = adjacence_matrix.RowCount;
if (root < 1 || root > num)
{
throw new ArgumentException("Root must be a vertex of the graph, e.i. in {1, ..., n}.");
}
Matrix matrix = new Matrix(num);
bool[] array = new bool[num + 1];
Queue queue = new Queue();
queue.Enqueue(root);
array[root] = true;
ArrayList[] array2 = new ArrayList[num + 1];
for (int i = 1; i <= num; i++)
{
array2[i] = new ArrayList();
for (int j = 1; j <= num; j++)
{
if (adjacence_matrix[i, j].Re != 0.0 && adjacence_matrix[i, j].Re != double.PositiveInfinity)
{
array2[i].Add(j);
}
}
}
while (queue.Count > 0)
{
int num2 = (int)queue.Peek();
if (array2[num2].Count > 0)
{
int num3 = (int)array2[num2][0];
if (!array[num3])
{
array[num3] = true;
matrix[num2, num3].Re = 1.0;
queue.Enqueue(num3);
}
array2[num2].RemoveAt(0);
}
else
{
queue.Dequeue();
}
}
return matrix;
}
public static Matrix ZeroOneRandom(int m, int n, double p)
{
Random random = new Random();
Matrix matrix = Zeros(m, n);
for (int i = 1; i <= m; i++)
{
for (int j = 1; j <= n; j++)
{
if (random.NextDouble() <= p)
{
matrix[i, j] = Complex.One;
}
}
}
return matrix;
}
public static Matrix ZeroOneRandom(int n, double p)
{
Random random = new Random();
Matrix matrix = Zeros(n);
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
if (random.NextDouble() <= p)
{
matrix[i, j] = Complex.One;
}
}
}
return matrix;
}
private static Matrix[] HouseholderVector(Matrix x)
{
int num = x.VectorLength();
if (num == 0)
{
throw new InvalidOperationException("Expected vector as argument.");
}
Matrix matrix = x / x.Norm();
Matrix matrix2 = matrix.Extract(2, num, 1, 1);
Complex complex = Dot(matrix2, matrix2);
Matrix matrix3 = Zeros(num, 1);
matrix3[1] = Complex.One;
matrix3.Insert(2, 1, matrix2);
double x2 = 0.0;
if (complex != 0.0)
{
Complex complex2 = Complex.Sqrt(matrix[1] * matrix[1] + complex);
if (matrix[1].Re <= 0.0)
{
matrix3[1] = matrix[1] - complex2;
}
else
{
matrix3[1] = -complex / (matrix[1] + complex2);
}
x2 = 2.0 * matrix3[1].Re * matrix3[1].Re / (complex.Re + matrix3[1].Re * matrix3[1].Re);
matrix3 /= matrix3[1];
}
return new Matrix[2]
{
matrix3,
new Matrix(x2)
};
}
public static Matrix BlockMatrix(Matrix A, Matrix B, Matrix C, Matrix D)
{
if (A.RowCount != B.RowCount || C.RowCount != D.RowCount || A.ColumnCount != C.ColumnCount || B.ColumnCount != D.ColumnCount)
{
throw new ArgumentException("Matrix dimensions must agree.");
}
Matrix matrix = new Matrix(A.RowCount + C.RowCount, A.ColumnCount + B.ColumnCount);
for (int i = 1; i <= matrix.rowCount; i++)
{
for (int j = 1; j <= matrix.columnCount; j++)
{
if (i <= A.RowCount)
{
if (j <= A.ColumnCount)
{
matrix[i, j] = A[i, j];
}
else
{
matrix[i, j] = B[i, j - A.ColumnCount];
}
}
else if (j <= C.ColumnCount)
{
matrix[i, j] = C[i - A.RowCount, j];
}
else
{
matrix[i, j] = D[i - A.RowCount, j - C.ColumnCount];
}
}
}
return matrix;
}
public static Matrix Solve(Matrix A, Matrix b)
{
Matrix matrix = A.Clone();
Matrix matrix2 = b.Clone();
if (!matrix.IsSquare())
{
throw new InvalidOperationException("Cannot uniquely solve non-square equation system.");
}
int n = matrix.RowCount;
Matrix matrix3 = matrix.LUSafe();
matrix2 = matrix3 * matrix2;
(matrix.ExtractLowerTrapeze() - Diag(matrix.DiagVector()) + Identity(n)).ForwardInsertion(matrix2);
matrix.ExtractUpperTrapeze().BackwardInsertion(matrix2);
return matrix2;
}
public Matrix Re()
{
Matrix matrix = new Matrix(rowCount, columnCount);
for (int i = 1; i <= rowCount; i++)
{
for (int j = 1; j <= columnCount; j++)
{
matrix[i, j] = new Complex(this[i, j].Re);
}
}
return matrix;
}
public Matrix Im()
{
Matrix matrix = new Matrix(rowCount, columnCount);
for (int i = 1; i <= rowCount; i++)
{
for (int j = 1; j <= columnCount; j++)
{
matrix[i, j] = new Complex(this[i, j].Im);
}
}
return matrix;
}
public Matrix[] HessenbergHouseholder()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot perform Hessenberg Householder decomposition of non-square matrix.");
}
int num = rowCount;
Matrix matrix = Identity(num);
Matrix matrix2 = Clone();
Matrix[] array = new Matrix[2];
for (int i = 1; i <= num - 2; i++)
{
array = HouseholderVector(matrix2.Extract(i + 1, num, i, i));
Matrix a = Identity(i);
Matrix matrix3 = Zeros(i, num - i);
int n = array[0].VectorLength();
Matrix matrix4 = Identity(n) - array[1][1, 1] * array[0] * array[0].Transpose();
matrix2.Insert(i + 1, i, matrix4 * matrix2.Extract(i + 1, num, i, num));
matrix2.Insert(1, i + 1, matrix2.Extract(1, num, i + 1, num) * matrix4);
Matrix matrix5 = BlockMatrix(a, matrix3, matrix3.Transpose(), matrix4);
matrix *= matrix5;
}
return new Matrix[2] { matrix2, matrix };
}
public Matrix Extract(int i1, int i2, int j1, int j2)
{
if (i2 < i1 || j2 < j1 || i1 <= 0 || j2 <= 0 || i2 > rowCount || j2 > columnCount)
{
throw new ArgumentException("Index exceeds matrix dimension.");
}
Matrix matrix = new Matrix(i2 - i1 + 1, j2 - j1 + 1);
for (int k = i1; k <= i2; k++)
{
for (int l = j1; l <= j2; l++)
{
matrix[k - i1 + 1, l - j1 + 1] = this[k, l];
}
}
return matrix;
}
public Matrix ExtractLowerTrapeze()
{
Matrix matrix = new Matrix(rowCount, columnCount);
for (int i = 1; i <= rowCount; i++)
{
for (int j = 1; j <= i; j++)
{
matrix[i, j] = this[i, j];
}
}
return matrix;
}
public Matrix ExtractUpperTrapeze()
{
Matrix matrix = new Matrix(rowCount, columnCount);
for (int i = 1; i <= rowCount; i++)
{
for (int j = i; j <= columnCount; j++)
{
matrix[i, j] = this[i, j];
}
}
return matrix;
}
public Matrix[] ColumnVectorize()
{
Matrix[] array = new Matrix[columnCount];
for (int i = 1; i <= array.Length; i++)
{
array[i] = Column(i);
}
return array;
}
public Matrix[] RowVectorize()
{
Matrix[] array = new Matrix[rowCount];
for (int i = 1; i <= array.Length; i++)
{
array[i] = Row(i);
}
return array;
}
public void VerticalFlip()
{
Values.Reverse();
}
public void HorizontalFlip()
{
for (int i = 0; i < rowCount; i++)
{
((ArrayList)Values[i]).Reverse();
}
}
public void SwapColumns(int j1, int j2)
{
if (j1 <= 0 || j1 > columnCount || j2 <= 0 || j2 > columnCount)
{
throw new ArgumentException("Indices must be positive and <= number of cols.");
}
if (j1 != j2)
{
j1--;
j2--;
for (int i = 0; i < rowCount; i++)
{
object value = ((ArrayList)Values[i])[j1];
((ArrayList)Values[i])[j1] = ((ArrayList)Values[i])[j2];
((ArrayList)Values[i])[j2] = value;
}
}
}
public void SwapRows(int i1, int i2)
{
if (i1 <= 0 || i1 > rowCount || i2 <= 0 || i2 > rowCount)
{
throw new ArgumentException("Indices must be positive and <= number of rows.");
}
if (i1 != i2)
{
ArrayList value = (ArrayList)Values[--i1];
Values[i1] = Values[--i2];
Values[i2] = value;
}
}
public void DeleteRow(int i)
{
if (i <= 0 || i > rowCount)
{
throw new ArgumentException("Index must be positive and <= number of rows.");
}
Values.RemoveAt(i - 1);
rowCount--;
}
public void DeleteColumn(int j)
{
if (j <= 0 || j > columnCount)
{
throw new ArgumentException("Index must be positive and <= number of cols.");
}
for (int i = 0; i < rowCount; i++)
{
((ArrayList)Values[i]).RemoveAt(j - 1);
}
columnCount--;
}
public Matrix ExtractRow(int i)
{
Matrix result = Row(i);
DeleteRow(i);
return result;
}
public Matrix ExtractColumn(int j)
{
if (j <= 0 || j > columnCount)
{
throw new ArgumentException("Index must be positive and <= number of cols.");
}
Matrix result = Column(j);
DeleteColumn(j);
return result;
}
public void InsertRow(Matrix row, int i)
{
int num = row.VectorLength();
if (num == 0)
{
throw new InvalidOperationException("Row must be a vector of length > 0.");
}
if (i <= 0)
{
throw new ArgumentException("Row index must be positive.");
}
if (i > rowCount)
{
this[i, num] = Complex.Zero;
}
else if (num > columnCount)
{
this[i, num] = Complex.Zero;
rowCount++;
}
else
{
rowCount++;
}
Values.Insert(--i, new ArrayList(num));
for (int j = 1; j <= num; j++)
{
((ArrayList)Values[i]).Add(row[j]);
}
for (int k = num; k < columnCount; k++)
{
((ArrayList)Values[i]).Add(Complex.Zero);
}
}
public void Insert(int i, int j, Matrix M)
{
for (int k = 1; k <= M.rowCount; k++)
{
for (int l = 1; l <= M.columnCount; l++)
{
this[i + k - 1, j + l - 1] = M[k, l];
}
}
}
public void InsertColumn(Matrix col, int j)
{
int num = col.VectorLength();
if (num == 0)
{
throw new InvalidOperationException("Row must be a vector of length > 0.");
}
if (j <= 0)
{
throw new ArgumentException("Row index must be positive.");
}
if (j > columnCount)
{
this[num, j] = Complex.Zero;
}
else
{
columnCount++;
}
if (num > rowCount)
{
this[num, j] = Complex.Zero;
}
j--;
for (int i = 0; i < num; i++)
{
((ArrayList)Values[i]).Insert(j, col[i + 1]);
}
for (int k = num; k < rowCount; k++)
{
((ArrayList)Values[k]).Insert(j, 0);
}
}
public Matrix Inverse()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot invert non-square matrix.");
}
Complex complex = Determinant();
if (complex == Complex.Zero)
{
throw new InvalidOperationException("Cannot invert (nearly) singular matrix.");
}
int num = columnCount;
if (num == 1)
{
return new Matrix(1.0 / complex);
}
if (IsReal() && IsOrthogonal())
{
return Transpose();
}
if (IsUnitary())
{
return ConjTranspose();
}
if (IsDiagonal())
{
Matrix matrix = DiagVector();
for (int i = 1; i <= num; i++)
{
matrix[i] = 1.0 / matrix[i];
}
return Diag(matrix);
}
Complex[,] array = new Complex[num, num];
for (int j = 0; j < num; j++)
{
for (int k = 0; k < num; k++)
{
array[j, k] = Math.Pow(-1.0, j + k) * Minor(k + 1, j + 1).Determinant();
}
}
return new Matrix(array) / complex;
}
public Matrix InverseLeverrier()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot invert non-square matrix.");
}
int num = rowCount;
Matrix matrix = Identity(num);
Matrix matrix2 = matrix;
Complex complex;
for (int i = 1; i < num; i++)
{
Matrix matrix3 = this * matrix2;
complex = 1.0 / (double)i * matrix3.Trace();
matrix2 = complex * matrix - matrix3;
}
complex = (this * matrix2).Trace() / num;
if (complex != Complex.Zero)
{
return matrix2 / complex;
}
throw new InvalidOperationException("WARNING: Matrix nearly singular or badly scaled.");
}
public Matrix Minor(int i, int j)
{
Matrix matrix = Clone();
matrix.DeleteRow(i);
matrix.DeleteColumn(j);
return matrix;
}
public Matrix Clone()
{
Matrix matrix = new Matrix();
matrix.rowCount = rowCount;
matrix.columnCount = columnCount;
for (int i = 0; i < rowCount; i++)
{
matrix.Values.Add(((ArrayList)Values[i]).Clone());
}
return matrix;
}
public Matrix DiagVector()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot get diagonal of non-square matrix.");
}
Matrix matrix = new Matrix(columnCount, 1);
for (int i = 1; i <= columnCount; i++)
{
matrix[i] = this[i, i];
}
return matrix;
}
public Matrix Column(int j)
{
Matrix matrix = new Matrix(rowCount, 1);
for (int i = 1; i <= rowCount; i++)
{
matrix[i] = this[i, j];
}
return matrix;
}
public Matrix Row(int i)
{
if (i <= 0 || i > rowCount)
{
throw new ArgumentException("Index exceed matrix dimension.");
}
Matrix matrix = new Matrix(columnCount, 1);
for (int j = 1; j <= columnCount; j++)
{
matrix[j] = this[i, j];
}
return matrix;
}
public Matrix Transpose()
{
Matrix matrix = new Matrix(columnCount, rowCount);
for (int i = 1; i <= columnCount; i++)
{
for (int j = 1; j <= rowCount; j++)
{
matrix[i, j] = this[j, i];
}
}
return matrix;
}
public Matrix Conjugate()
{
Matrix matrix = new Matrix(rowCount, columnCount);
for (int i = 1; i <= rowCount; i++)
{
for (int j = 1; j <= columnCount; j++)
{
matrix[i, j] = Complex.Conj(this[i, j]);
}
}
return matrix;
}
public Matrix ConjTranspose()
{
return Transpose().Conjugate();
}
public void LU()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot perform LU-decomposition of non-square matrix.");
}
int num = columnCount;
for (int i = 1; i <= num; i++)
{
if (this[i, i] == 0.0)
{
throw new DivideByZeroException("Warning: Matrix badly scaled or close to singular. Try LUSafe() instead. Check if det != 0.");
}
for (int j = 1; j < i; j++)
{
for (int k = j + 1; k <= num; k++)
{
this[k, i] -= this[k, j] * this[j, i];
}
}
for (int l = i + 1; l <= num; l++)
{
this[l, i] /= this[i, i];
}
}
}
public Matrix LUSafe()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot perform LU-decomposition of non-square matrix.");
}
int num = columnCount;
Matrix matrix = Identity(num);
for (int i = 1; i <= num; i++)
{
if (i < num)
{
int num2 = i;
for (int j = i + 1; j <= num; j++)
{
if (Complex.Abs(this[j, i]) > Complex.Abs(this[num2, i]))
{
num2 = j;
}
}
if (num2 > i)
{
matrix.SwapRows(i, num2);
SwapRows(i, num2);
}
if (this[i, i] == 0.0)
{
throw new DivideByZeroException("Warning: Matrix close to singular.");
}
}
for (int k = 1; k < i; k++)
{
for (int l = k + 1; l <= num; l++)
{
this[l, i] -= this[l, k] * this[k, i];
}
}
for (int m = i + 1; m <= num; m++)
{
this[m, i] /= this[i, i];
}
}
return matrix;
}
public void Cholesky()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot perform Cholesky decomposition of non-square matrix.");
}
if (!IsSPD())
{
throw new InvalidOperationException("Cannot perform Cholesky decomposition of matrix not being symmetric positive definite.");
}
int num = rowCount;
for (int i = 1; i < num; i++)
{
this[i, i] = Complex.Sqrt(this[i, i]);
for (int j = 1; j <= num - i; j++)
{
this[i + j, i] /= this[i, i];
}
for (int k = i + 1; k <= num; k++)
{
for (int l = 0; l <= num - k; l++)
{
this[k + l, k] -= this[k + l, i] * this[k, i];
}
}
}
this[num, num] = Complex.Sqrt(this[num, num]);
}
public void CholeskyUndo()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot undo cholesky decomposition on non-square matrix.");
}
this[1, 1] = Sqr(this[1, 1]);
for (int i = 2; i <= rowCount; i++)
{
Complex zero = Complex.Zero;
for (int j = 1; j <= i - 1; j++)
{
zero += Sqr(this[i, j]);
}
this[i, i] = Sqr(this[i, i]) + zero;
}
SymmetrizeDown();
}
public void ForwardInsertion(Matrix b)
{
if (!IsLowerTriangular())
{
throw new InvalidOperationException("Cannot perform forward insertion for matrix not being lower triangular.");
}
if (DiagProd() == 0.0)
{
throw new InvalidOperationException("Warning: Matrix is nearly singular.");
}
int num = rowCount;
if (b.VectorLength() != num)
{
throw new ArgumentException("Parameter must vector of the same height as matrix.");
}
for (int i = 1; i <= num - 1; i++)
{
b[i] /= this[i, i];
for (int j = 1; j <= num - i; j++)
{
b[i + j] -= b[i] * this[i + j, i];
}
}
b[num] /= this[num, num];
}
public void BackwardInsertion(Matrix b)
{
if (!IsUpperTriangular())
{
throw new InvalidOperationException("Cannot perform backward insertion for matrix not being upper triangular.");
}
if (DiagProd() == 0.0)
{
throw new InvalidOperationException("Warning: Matrix is nearly singular.");
}
int num = rowCount;
if (b.VectorLength() != num)
{
throw new ArgumentException("Parameter must vector of the same height as matrix.");
}
for (int num2 = num; num2 >= 2; num2--)
{
b[num2] /= this[num2, num2];
for (int i = 1; i <= num2 - 1; i++)
{
b[i] -= b[num2] * this[i, num2];
}
}
b[1] /= this[1, 1];
}
public void SymmetrizeDown()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot symmetrize non-square matrix.");
}
for (int i = 1; i <= columnCount; i++)
{
for (int j = i + 1; j <= columnCount; j++)
{
this[j, i] = this[i, j];
}
}
}
public void SymmetrizeUp()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot symmetrize non-square matrix.");
}
for (int i = 1; i <= rowCount; i++)
{
for (int j = i + 1; j <= columnCount; j++)
{
this[i, j] = this[j, i];
}
}
}
public Matrix[] QRGramSchmidt()
{
int num = rowCount;
int num2 = columnCount;
Matrix matrix = Clone();
Matrix matrix2 = new Matrix(num, num2);
Matrix matrix3 = new Matrix(num2, num2);
for (int i = 1; i <= num; i++)
{
matrix2[i, 1] = matrix[i, 1];
}
matrix3[1, 1] = Complex.One;
for (int j = 1; j <= num2; j++)
{
matrix3[j, j] = new Complex(matrix.Column(j).Norm());
for (int k = 1; k <= num; k++)
{
matrix2[k, j] = matrix[k, j] / matrix3[j, j];
}
for (int l = j + 1; l <= num2; l++)
{
matrix3[j, l] = Dot(matrix2.Column(j), matrix.Column(l));
for (int m = 1; m <= num; m++)
{
matrix[m, l] -= matrix2[m, j] * matrix3[j, l];
}
}
}
return new Matrix[2] { matrix2, matrix3 };
}
public Matrix Eigenvalues()
{
return QRIterationBasic(40).DiagVector();
}
public Matrix Eigenvector(Complex eigenvalue)
{
throw new NotImplementedException();
}
public Matrix SolveCG(Matrix b)
{
if (!IsSPD())
{
throw new InvalidOperationException("CG method only works for spd matrices.");
}
if (!IsReal())
{
throw new InvalidOperationException("CG method only works for real matrices.");
}
int m = rowCount;
int num = 150;
double num2 = 1E-06;
Matrix matrix = Ones(m, 1);
Matrix matrix2 = b - this * matrix;
Matrix matrix3 = matrix2;
double num3 = matrix2.Norm();
num3 *= num3;
num2 *= num2;
Matrix matrix4 = Zeros(m, 1);
if (num3 <= num2)
{
return matrix;
}
for (int i = 0; i < num; i++)
{
matrix4 = this * matrix3;
double re = Dot(matrix4, matrix3).Re;
if (Math.Abs(re) <= num2)
{
return matrix;
}
double num4 = num3 / re;
matrix += num4 * matrix3;
matrix2 -= num4 * matrix4;
double num5 = num3;
num3 = matrix2.Norm();
num3 *= num3;
if (num3 <= num2)
{
return matrix;
}
matrix3 = matrix2 + num3 / num5 * matrix3;
}
return matrix;
}
public Matrix QRIterationBasic(int max_iterations)
{
if (!IsReal())
{
throw new InvalidOperationException("Basic QR iteration is possible only for real matrices.");
}
Matrix matrix = Clone();
Matrix[] array = new Matrix[2];
for (int i = 0; i < max_iterations; i++)
{
array = matrix.QRGramSchmidt();
matrix = array[1] * array[0];
}
return matrix;
}
public Matrix QRIterationHessenberg(int max_iterations)
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot perform QR iteration of non-square matrix.");
}
int num = RowCount;
Matrix[] array = HessenbergHouseholder();
Matrix matrix = array[0];
for (int i = 1; i <= max_iterations; i++)
{
Matrix[] array2 = matrix.QRGivens();
matrix = array2[1];
for (int j = 1; j <= num - 1; j++)
{
matrix.Gacol(array2[2][j], array2[3][j], 1, j + 1, j, j + 1);
}
}
return matrix;
}
public Matrix[] QRGivens()
{
Matrix matrix = Clone();
int num = matrix.ColumnCount;
Matrix matrix2 = Zeros(num - 1, 1);
Matrix matrix3 = Zeros(num - 1, 1);
for (int i = 1; i <= num - 1; i++)
{
Complex[] array = GivensCS(matrix[i, i], matrix[i + 1, i]);
matrix2[i] = array[0];
matrix3[i] = array[1];
Garow(matrix2[i], matrix3[i], 1, i + 1, i, i + 1);
}
return new Matrix[4]
{
GivProd(matrix2, matrix3, num),
matrix,
matrix2,
matrix3
};
}
private Matrix GivProd(Matrix c, Matrix s, int n)
{
int num = n - 1;
int num2 = n - 2;
Matrix matrix = Eye(n);
matrix[num, num] = c[num];
matrix[n, n] = c[num];
matrix[num, n] = s[num];
matrix[n, num] = -s[num];
for (int num3 = num2; num3 >= 1; num3--)
{
int num4 = num3 + 1;
matrix[num3, num3] = c[num3];
matrix[num4, num3] = -s[num3];
Matrix matrix2 = matrix.Extract(num4, num4, num4, n);
matrix.Insert(num3, num4, s[num3] * matrix2);
matrix.Insert(num4, num4, c[num3] * matrix2);
}
return matrix;
}
private void Garow(Complex c, Complex s, int i, int k, int j1, int j2)
{
for (int l = j1; l <= j2; l++)
{
Complex complex = this[i, l];
Complex complex2 = this[k, l];
this[i, l] = c * complex - s * complex2;
this[k, l] = s * complex + c * complex2;
}
}
public void Gacol(Complex c, Complex s, int j1, int j2, int i, int k)
{
for (int l = j1; l <= j2; l++)
{
Complex complex = this[l, i];
Complex complex2 = this[l, k];
this[l, i] = c * complex - s * complex2;
this[l, k] = s * complex + c * complex2;
}
}
private Complex[] GivensCS(Complex xi, Complex xk)
{
Complex zero = Complex.Zero;
Complex complex = Complex.Zero;
if (xk == 0.0)
{
zero = Complex.One;
}
else if (Complex.Abs(xk) > Complex.Abs(xi))
{
Complex complex2 = -xi / xk;
complex = 1.0 / Complex.Sqrt(1.0 + complex2 * complex2);
zero = complex * complex2;
}
else
{
Complex complex3 = -xk / xi;
zero = 1.0 / Complex.Sqrt(1.0 + complex3 * complex3);
complex = zero * complex3;
}
return new Complex[2] { zero, complex };
}
public Complex Determinant()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot calc determinant of non-square matrix.");
}
if (columnCount == 1)
{
return this[1, 1];
}
if (IsTrapeze())
{
return DiagProd();
}
Matrix matrix = Clone();
Matrix matrix2 = matrix.LUSafe();
return matrix2.Signum() * matrix.DiagProd();
}
public double ColumnSumNorm()
{
return TaxiNorm();
}
public double RowSumNorm()
{
return MaxNorm();
}
public Complex Permanent()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot compute permanent of non-square matrix.");
}
if (HasZeroRowOrColumn())
{
return Complex.Zero;
}
if (this == Ones(rowCount))
{
return new Complex(Factorial(rowCount));
}
if (IsPermutation())
{
return Complex.One;
}
Complex zero = Complex.Zero;
int minRow = GetMinRow();
int minColumn = GetMinColumn();
if (AbsRowSum(minRow) < AbsColumnSum(minColumn))
{
for (int i = 1; i <= columnCount; i++)
{
if (this[minRow, i] != 0.0)
{
zero += this[minRow, i] * Minor(minRow, i).Permanent();
}
}
}
else
{
for (int j = 1; j <= rowCount; j++)
{
if (this[j, minColumn] != 0.0)
{
zero += this[j, minColumn] * Minor(j, minColumn).Permanent();
}
}
}
return zero;
}
public int GetMinRow()
{
double num = AbsRowSum(1);
int result = 1;
for (int i = 2; i <= rowCount; i++)
{
double num2 = AbsRowSum(i);
if (num2 < num)
{
num = num2;
result = i;
}
}
return result;
}
public int GetMinColumn()
{
double num = AbsColumnSum(1);
int result = 1;
for (int i = 2; i <= columnCount; i++)
{
double num2 = AbsColumnSum(i);
if (num2 < num)
{
num = num2;
result = i;
}
}
return result;
}
private double Factorial(int x)
{
double num = 1.0;
for (int i = 2; i <= x; i++)
{
num *= (double)i;
}
return num;
}
public double Signum()
{
double num = 1.0;
int num2 = rowCount;
int num3 = 0;
for (int i = 1; i < num2; i++)
{
double num4;
for (num4 = 1.0; num4 < (double)num2 && this[i, (int)num4] != Complex.One; num4 += 1.0)
{
num3 = num3++;
}
for (int j = i + 1; j <= num2; j++)
{
double num5;
for (num5 = 1.0; num5 <= (double)num2 && this[j, (int)num5] != Complex.One; num5 += 1.0)
{
}
num *= (num4 - num5) / (double)(i - j);
}
}
return num;
}
public double Condition()
{
return TaxiNorm() * Inverse().TaxiNorm();
}
public double Condition(int p)
{
return PNorm(p) * Inverse().PNorm(p);
}
public double ConditionFro()
{
return FrobeniusNorm() * Inverse().FrobeniusNorm();
}
public double PNorm(double p)
{
if (p <= 0.0)
{
throw new ArgumentException("Argument must be greater than zero.");
}
if (p == 1.0)
{
return TaxiNorm();
}
if (p == double.PositiveInfinity)
{
return MaxNorm();
}
int num = VectorLength();
if (num == 0)
{
throw new InvalidOperationException("Cannot calc p-norm of matrix.");
}
double num2 = 0.0;
for (int i = 1; i <= num; i++)
{
num2 += Math.Pow(Complex.Abs(this[i]), p);
}
return Math.Pow(num2, 1.0 / p);
}
public double Norm()
{
return PNorm(2.0);
}
public double FrobeniusNorm()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot compute frobenius norm of non-square matrix.");
}
int num = columnCount;
double num2 = 0.0;
for (int i = 1; i <= num; i++)
{
for (int j = 1; j <= num; j++)
{
num2 += (this[i, j] * Complex.Conj(this[i, j])).Re;
}
}
return Math.Sqrt(num2);
}
public double TaxiNorm()
{
double num = 0.0;
int num2 = VectorLength();
if (num2 != 0)
{
for (int i = 1; i <= num2; i++)
{
num += Complex.Abs(this[i]);
}
}
else
{
double num3 = 0.0;
for (int j = 1; j <= columnCount; j++)
{
num3 = AbsColumnSum(j);
if (num3 > num)
{
num = num3;
}
}
}
return num;
}
public double MaxNorm()
{
double num = 0.0;
double num2 = 0.0;
int num3 = VectorLength();
if (num3 != 0)
{
for (int i = 1; i <= num3; i++)
{
num2 = Complex.Abs(this[i]);
if (num2 > num)
{
num = num2;
}
}
}
else
{
for (int j = 1; j <= rowCount; j++)
{
num2 = AbsRowSum(j);
if (num2 > num)
{
num = num2;
}
}
}
return num;
}
public Complex ColumnSum(int j)
{
if (j <= 0 || j > columnCount)
{
throw new ArgumentException("Index out of range.");
}
Complex zero = Complex.Zero;
j--;
for (int i = 0; i < rowCount; i++)
{
zero += (Complex)((ArrayList)Values[i])[j];
}
return zero;
}
public double AbsColumnSum(int j)
{
if (j <= 0 || j > columnCount)
{
throw new ArgumentException("Index out of range.");
}
double num = 0.0;
for (int i = 1; i <= rowCount; i++)
{
num += Complex.Abs(this[i, j]);
}
return num;
}
public Complex RowSum(int i)
{
if (i <= 0 || i > rowCount)
{
throw new ArgumentException("Index out of range.");
}
Complex zero = Complex.Zero;
ArrayList arrayList = (ArrayList)Values[i - 1];
for (int j = 0; j < columnCount; j++)
{
zero += (Complex)arrayList[j];
}
return zero;
}
public double AbsRowSum(int i)
{
if (i <= 0 || i > rowCount)
{
throw new ArgumentException("Index out of range.");
}
double num = 0.0;
for (int j = 1; j <= columnCount; j++)
{
num += Complex.Abs(this[i, j]);
}
return num;
}
public Complex DiagProd()
{
Complex one = Complex.One;
int num = Math.Min(rowCount, columnCount);
for (int i = 1; i <= num; i++)
{
one *= this[i, i];
}
return one;
}
public Complex Trace()
{
if (!IsSquare())
{
throw new InvalidOperationException("Cannot calc trace of non-square matrix.");
}
Complex zero = Complex.Zero;
for (int i = 1; i <= rowCount; i++)
{
zero += this[i, i];
}
return zero;
}
private Complex Sqr(Complex x)
{
return x * x;
}
public bool IsNormal()
{
return this * ConjTranspose() == ConjTranspose() * this;
}
public bool IsUnitary()
{
if (!IsSquare())
{
return false;
}
return ConjTranspose() * this == Identity(rowCount);
}
public bool IsHermitian()
{
if (!IsSquare())
{
return false;
}
return ConjTranspose() == this;
}
public bool IsReal()
{
for (int i = 1; i <= rowCount; i++)
{
for (int j = 1; j <= columnCount; j++)
{
if (!this[i, j].IsReal())
{
return false;
}
}
}
return true;
}
public bool IsSymmetricPositiveDefinite()
{
return IsSymmetric() && Definiteness() == DefinitenessType.PositiveDefinite;
}
public bool IsSPD()
{
return IsSymmetricPositiveDefinite();
}
public DefinitenessType Definiteness()
{
if (!IsSquare())
{
throw new InvalidOperationException("Definiteness undefined for non-square matrices.");
}
if (this == Zeros(rowCount, columnCount))
{
return DefinitenessType.Indefinite;
}
if (!IsSymmetric())
{
throw new InvalidOperationException("This test works only for symmetric matrices.");
}
if (!IsReal())
{
throw new InvalidOperationException("This test only works for real matrices.");
}
int num = rowCount;
Matrix[] array = new Matrix[num + 1];
for (int i = 0; i <= num; i++)
{
array[i] = Zeros(num, 1);
}
array[1] = Column(1);
for (int j = 2; j <= num; j++)
{
Matrix matrix = Column(j);
Matrix matrix2 = Zeros(num, 1);
for (int k = 1; k < j; k++)
{
matrix2 += array[k] * Dot(matrix, this * array[k]) / Dot(array[k], this * array[k]);
}
array[j] = matrix - matrix2;
}
bool flag = true;
for (int l = 1; l < num; l++)
{
Complex complex = Dot(array[l], this * array[l]) * Dot(array[l + 1], this * array[l + 1]);
if (complex == 0.0)
{
flag = false;
}
else if (complex.Re < 0.0)
{
return DefinitenessType.Indefinite;
}
}
if (Dot(array[1], this * array[1]).Re >= 0.0)
{
if (flag)
{
return DefinitenessType.PositiveDefinite;
}
return DefinitenessType.PositiveSemidefinite;
}
if (flag)
{
return DefinitenessType.NegativeDefinite;
}
return DefinitenessType.NegativeSemidefinite;
}
public bool HasZeroRowOrColumn()
{
for (int i = 1; i <= rowCount; i++)
{
if (AbsRowSum(i) == 0.0)
{
return true;
}
}
for (int j = 1; j <= columnCount; j++)
{
if (AbsColumnSum(j) == 0.0)
{
return true;
}
}
return false;
}
public bool IsZeroOneMatrix()
{
for (int i = 1; i <= rowCount; i++)
{
for (int j = 1; j <= columnCount; j++)
{
if (this[i, j] != Complex.Zero && this[i, j] != Complex.One)
{
return false;
}
}
}
return true;
}
public bool IsPermutation()
{
return !IsSquare() && IsZeroOneMatrix() && IsInvolutary();
}
public bool IsDiagonal()
{
return Clone() - Diag(DiagVector()) == Zeros(rowCount, columnCount);
}
public int VectorLength()
{
if (columnCount > 1 && rowCount > 1)
{
return 0;
}
return Math.Max(columnCount, rowCount);
}
public bool IsSquare()
{
return columnCount == rowCount;
}
public bool IsInvolutary()
{
return this * this == Identity(rowCount);
}
public bool IsSymmetric()
{
for (int i = 1; i <= rowCount; i++)
{
for (int j = 1; j <= columnCount; j++)
{
if (this[i, j] != this[j, i])
{
return false;
}
}
}
return true;
}
public bool IsOrthogonal()
{
return IsSquare() && this * Transpose() == Identity(rowCount);
}
public bool IsTrapeze()
{
return IsUpperTrapeze() || IsLowerTrapeze();
}
public bool IsTriangular()
{
return IsLowerTriangular() || IsUpperTriangular();
}
public bool IsUpperTriangular()
{
return IsSquare() && IsUpperTrapeze();
}
public bool IsLowerTriangular()
{
return IsSquare() && IsLowerTrapeze();
}
public bool IsUpperTrapeze()
{
for (int i = 1; i <= columnCount; i++)
{
for (int j = i + 1; j <= rowCount; j++)
{
if (this[j, i] != 0.0)
{
return false;
}
}
}
return true;
}
public bool IsLowerTrapeze()
{
for (int i = 1; i <= rowCount; i++)
{
for (int j = i + 1; j <= columnCount; j++)
{
if (this[i, j] != 0.0)
{
return false;
}
}
}
return true;
}
public override string ToString()
{
string text = string.Empty;
for (int i = 1; i <= rowCount; i++)
{
for (int j = 1; j <= columnCount; j++)
{
Complex complex = this[i, j];
text = ((complex.Re != double.PositiveInfinity && complex.Re != double.NegativeInfinity && complex.Im != double.PositiveInfinity && complex.Im != double.NegativeInfinity) ? ((complex.Re != double.NaN && complex.Im != double.NaN) ? (text + complex.ToString()) : (text + "?")) : (text + "oo"));
text += ";\t";
}
text = text + "\\" + Environment.NewLine;
}
return text;
}
public string ToString(string format)
{
string text = string.Empty;
for (int i = 1; i <= rowCount; i++)
{
for (int j = 1; j <= columnCount; j++)
{
Complex complex = this[i, j];
text = ((complex.Re != double.PositiveInfinity && complex.Re != double.NegativeInfinity && complex.Im != double.PositiveInfinity && complex.Im != double.NegativeInfinity) ? ((complex.Re != double.NaN && complex.Im != double.NaN) ? (text + complex.ToString(format)) : (text + "?")) : (text + "oo"));
text += ";\t";
}
text = text + "\\" + Environment.NewLine;
}
return text;
}
public override bool Equals(object obj)
{
return obj.ToString() == ToString();
}
public override int GetHashCode()
{
return -1;
}
public static bool operator ==(Matrix A, Matrix B)
{
if (A.RowCount != B.RowCount || A.ColumnCount != B.ColumnCount)
{
return false;
}
for (int i = 1; i <= A.RowCount; i++)
{
for (int j = 1; j <= A.ColumnCount; j++)
{
if (A[i, j] != B[i, j])
{
return false;
}
}
}
return true;
}
public static bool operator !=(Matrix A, Matrix B)
{
return !(A == B);
}
public static Matrix operator +(Matrix A, Matrix B)
{
if (A.RowCount != B.RowCount || A.ColumnCount != B.ColumnCount)
{
throw new ArgumentException("Matrices must be of the same dimension.");
}
for (int i = 1; i <= A.RowCount; i++)
{
for (int j = 1; j <= A.ColumnCount; j++)
{
A[i, j] += B[i, j];
}
}
return A;
}
public static Matrix operator -(Matrix A, Matrix B)
{
if (A.RowCount != B.RowCount || A.ColumnCount != B.ColumnCount)
{
throw new ArgumentException("Matrices must be of the same dimension.");
}
for (int i = 1; i <= A.RowCount; i++)
{
for (int j = 1; j <= A.ColumnCount; j++)
{
A[i, j] -= B[i, j];
}
}
return A;
}
public static Matrix operator -(Matrix A)
{
for (int i = 1; i <= A.RowCount; i++)
{
for (int j = 1; j <= A.ColumnCount; j++)
{
A[i, j] = -A[i, j];
}
}
return A;
}
public static Matrix operator *(Matrix A, Matrix B)
{
if (A.ColumnCount != B.RowCount)
{
throw new ArgumentException("Inner matrix dimensions must agree.");
}
Matrix matrix = new Matrix(A.RowCount, B.ColumnCount);
for (int i = 1; i <= A.RowCount; i++)
{
for (int j = 1; j <= B.ColumnCount; j++)
{
matrix[i, j] = Dot(A.Row(i), B.Column(j));
}
}
return matrix;
}
public static Matrix operator *(Matrix A, Complex x)
{
Matrix matrix = new Matrix(A.rowCount, A.columnCount);
for (int i = 1; i <= A.RowCount; i++)
{
for (int j = 1; j <= A.ColumnCount; j++)
{
matrix[i, j] = A[i, j] * x;
}
}
return matrix;
}
public static Matrix operator *(Complex x, Matrix A)
{
Matrix matrix = new Matrix(A.RowCount, A.ColumnCount);
for (int i = 1; i <= A.RowCount; i++)
{
for (int j = 1; j <= A.ColumnCount; j++)
{
matrix[i, j] = A[i, j] * x;
}
}
return matrix;
}
public static Matrix operator *(Matrix A, double x)
{
return new Complex(x) * A;
}
public static Matrix operator *(double x, Matrix A)
{
return new Complex(x) * A;
}
public static Matrix operator /(Matrix A, Complex x)
{
return 1.0 / x * A;
}
public static Matrix operator /(Matrix A, double x)
{
return new Complex(1.0 / x) * A;
}
public static Matrix operator ^(Matrix A, int k)
{
if (k < 0)
{
if (A.IsSquare())
{
return A.InverseLeverrier() ^ -k;
}
throw new InvalidOperationException("Cannot take non-square matrix to the power of zero.");
}
switch (k)
{
case 0:
if (A.IsSquare())
{
return Identity(A.RowCount);
}
throw new InvalidOperationException("Cannot take non-square matrix to the power of zero.");
case 1:
if (A.IsSquare())
{
return A;
}
throw new InvalidOperationException("Cannot take non-square matrix to the power of one.");
default:
{
Matrix result = A;
for (int i = 1; i < k; i++)
{
result *= A;
}
return result;
}
}
}
}
}