2628 lines
54 KiB
C#
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;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|