using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
public class Matrix {
// Fields
private double[,] data;
private int nRows, nColumns;
// Constructors
public Matrix() { }
public Matrix(int nRows, int nColumns) {
this.data = new double[nRows, nColumns];
this.nRows = nRows;
this.nColumns = nColumns;
}
public Matrix(double[,] matrix) {
nRows = matrix.GetLength(0);
nColumns = matrix.GetLength(1);
this.data = new double[nRows, nColumns];
this.data = matrix;
}
public Matrix(Matrix other) {
nRows = other.GetLength(0);
nColumns = other.GetLength(1);
this.data = new double[nRows, nColumns];
for (int i = 0; i < nRows; ++i) {
for (int j = 0; j < nColumns; ++j) {
this.data[i, j] = other[i, j];
}
}
}
// Indexer
public double this[int i, int j] {
get { return this.data[i, j]; }
set { this.data[i, j] = value; }
}
// Overriding operators
public static bool operator ==(Matrix first, Matrix second) {
for (int i = 0; i < first.GetLength(0); ++i)
for (int j = 0; j < first.GetLength(1); ++j)
if (first[i, j] != second[i, j])
return false;
return true;
}
public static bool operator !=(Matrix first, Matrix second) {
return (first == second) ? false : true;
}
public override bool Equals(object obj) {
obj = (Matrix)obj;
return base.Equals(obj);
}
public override int GetHashCode() {
return base.GetHashCode();
}
public override string ToString() {
string stream = "";
for (int i = 0; i < this.nRows; ++i) {
for (int j = 0; j < this.nColumns; ++j)
stream += String.Format("{0:f4}\t\t", this[i, j]);
stream += "\n";
}
return stream;
}
public static Matrix operator +(Matrix first, Matrix second) {
Matrix retValue = new Matrix(first.GetLength(0), first.GetLength(1));
try {
if (first.GetLength(0) == second.GetLength(0) && first.GetLength(1) == second.GetLength(1)) {
for (int i = 0; i < first.GetLength(0); ++i)
for (int j = 0; j < first.GetLength(1); ++j)
retValue[i, j] = first[i, j] + second[i, j];
}
}
catch (IndexOutOfRangeException err) {
Console.WriteLine(err.Message);
}
return retValue;
}
public static Matrix operator -(Matrix first, Matrix second) {
Matrix retValue = new Matrix(first.GetLength(0), first.GetLength(1));
try {
for (int i = 0; i < first.GetLength(0); ++i)
for (int j = 0; j < first.GetLength(1); ++j)
retValue[i, j] = first[i, j] - second[i, j];
}
catch (IndexOutOfRangeException err) {
Console.WriteLine(err.Message);
}
return retValue;
}
public static Matrix operator *(Matrix matrix, double coeff) {
for (int i = 0; i < matrix.GetLength(0); ++i)
for (int j = 0; j < matrix.GetLength(1); ++j)
matrix[i, j] *= coeff;
return matrix;
}
public static Matrix operator *(double coeff, Matrix matrix) {
for (int i = 0; i < matrix.GetLength(0); ++i)
for (int j = 0; j < matrix.GetLength(1); ++j)
matrix[i, j] *= coeff;
return matrix;
}
public static Matrix operator *(Matrix first, Matrix second) {
Matrix retValue = new Matrix(first.GetLength(0), second.GetLength(1));
double curVal = 0.0;
try {
for (int i = 0; i < first.GetLength(0); ++i) {
for (int j = 0; j < second.GetLength(1); ++j) {
curVal = 0.0;
for (int k = 0; k < Math.Max(first.GetLength(1), second.GetLength(0)); ++k)
curVal += first[i, k] * second[k, j];
retValue[i, j] = curVal;
}
}
}
catch (IndexOutOfRangeException err) {
Console.WriteLine(err.Message);
}
return retValue;
}
// User's functions
public int GetLength(int dimension) {
return (dimension == 0) ? nRows : nColumns;
}
public Matrix GetMinor2(int ind) {
Matrix retValue = new Matrix(this.nRows - ind, this.nColumns - ind);
for (int i = ind; i < this.GetLength(0); ++i) {
for (int j = ind; j < this.GetLength(0); ++j) {
retValue[i - ind, j - ind] = this[i, j];
}
}
return retValue;
}
public Matrix SetMinor2(Matrix minor, int ind) {
for (int i = ind; i < this.GetLength(0); ++i) {
for (int j = ind; j < this.GetLength(0); ++j) {
this[i, j] = minor[i - ind, j - ind];
}
}
return this;
}
public bool isNull() {
for (int i = 0; i < this.GetLength(0); ++i) {
if (this[i, 0] != 0)
return false;
}
return true;
}
public Matrix GetColumn(int indColumn, int begin) {
if (begin == this.GetLength(0))
return new Matrix(0, 0);
Matrix retValue = new Matrix(this.GetLength(0) - begin, 1);
for (int i = begin; i < this.GetLength(0); ++i)
retValue[i - begin, 0] = this[i, indColumn];
return retValue;
}
public Matrix MakeIdentity() {
Matrix retValue = new Matrix(this.nRows, this.nColumns);
for (int i = 0; i < this.nRows; ++i)
retValue[i, i] = 1;
return retValue;
}
public Matrix GetTransposed() {
Matrix retValue = new Matrix(this.nColumns, this.nRows);
for (int i = 0; i < this.nRows; ++i)
for (int j = 0; j < this.nColumns; ++j)
retValue[j, i] = this[i, j];
return retValue;
}
public double GetSqrLength(int indColumn) {
double retValue = 0.0;
for (int i = 0; i < this.nRows; ++i)
retValue += this[i, indColumn] * this[i, indColumn];
return retValue;
}
public double GetSqrtLength(int indColumn) {
return Math.Sqrt(this.GetSqrLength(indColumn));
}
public double GetSummColumn(int indColumn, int start) {
double retValue = 0.0;
for (int i = start; i < this.nRows; ++i)
retValue += this[i, indColumn] * this[i, indColumn];
return retValue;
}
public void SetRow(List<double> values, int ind) {
for (int i = 0; i < values.Count; ++i) {
this[ind, i] = values[i];
}
}
public static Matrix GetNewMinor(Matrix matrix, int r, int c) {
int n = matrix.GetLength(0),
m = matrix.GetLength(1);
int k = 0;
Matrix retValue = new Matrix(n - 1, m - 1);
for (int i = 0; i < n; ++i) {
if (i == r)
continue;
List<double> set = new List<double>();
for (int j = 0; j < m; ++j) {
if (j != c)
set.Add(matrix[i, j]);
}
retValue.SetRow(set, k);
++k;
}
return retValue;
}
public Matrix GetInverse() {
Matrix retValue = new Matrix(this.nRows, this.nColumns);
bool[] rows = new bool[this.nRows],
columns = new bool[this.nColumns];
for (int i = 0; i < this.nRows; ++i) {
rows[i] = true;
columns[i] = true;
}
for (int i = 0; i < this.nRows; ++i) {
for (int j = 0; j < this.nColumns; ++j) {
retValue[i, j] = GetNewMinor(this, i, j).GetDet();
if ((i + j) % 2 != 0)
retValue[i, j] *= -1;
rows[i] = true;
columns[j] = true;
}
}
retValue = retValue.GetTransposed();
retValue *= 1.0 / GetDet();
return retValue;
}
public double GetDet() {
bool[] rows = new bool[this.nRows],
columns = new bool[this.nColumns];
for (int i = 0; i < this.nRows; ++i) {
rows[i] = true;
columns[i] = true;
}
return GetMinor(rows, columns, this.nRows);
}
public double GetMinor(bool[] rows, bool[] columns, int n) {
double minorValue = 0.0;
int i = 0,
j = 0,
order = 1;
i = (Array.IndexOf(rows, true) != -1) ? Array.IndexOf(rows, true) : 0;
j = (Array.IndexOf(columns, true) != -1) ? Array.IndexOf(columns, true) : 0;
if (n == 1) {
return this[i, j];
}
else {
rows[i] = false;
n--;
for (j = 0; j < this.GetLength(0); ++j) {
if (columns[j]) {
columns[j] = false;
if (this[i, j] != 0) {
minorValue += order * this[i, j] * GetMinor(rows, columns, n);
}
columns[j] = true;
order = -order;
}
}
rows[i] = true;
return minorValue;
}
}
public static Matrix QR(Matrix A) {
int n = A.GetLength(0),
m = A.GetLength(1);
Matrix y,
Q = new Matrix(n, m).MakeIdentity();
for (int i = 0; i < n - 1; ++i) {
y = A.GetColumn(i, i);
Matrix curV = new Matrix(n - i, m - i).MakeIdentity(),
curQ = new Matrix(n, m).MakeIdentity(),
omega = new Matrix(n - i, 1),
e = new Matrix(n - i, 1); e[0, 0] = 1;
double yNorm = Math.Sqrt((y.GetTransposed() * y)[0, 0]),
omegaNorm = 0.0;
int sign = Math.Sign(y[0, 0]);
omega = y + sign * (yNorm * e);
omegaNorm = Math.Sqrt((omega.GetTransposed() * omega)[0, 0]);
omega *= (1.0 / omegaNorm);
curV = curV - (omega * omega.GetTransposed()) * 2.0;
curQ.SetMinor2(curV, i);
A = curQ * A;
Q = curQ * Q;
}
return Q.GetTransposed();
}
public static Matrix GetSolve(Matrix A, Matrix b, double eps) {
int n = A.GetLength(0);
Matrix Q = QR(A),
R = Q.GetTransposed() * A,
x = new Matrix(n, 1),
nb = new Matrix(n, 1);
for (int i = 0; i < n; i++) {
nb[i, 0] = 0;
for (int j = 0; j < n; j++)
nb[i, 0] += b[j, 0] * Q[j, i];
}
for (int i = n - 1; i >= 0; --i) {
double s = 0.0;
for (int j = i + 1; j < n; ++j)
s += R[i, j] * x[j, 0];
x[i, 0] = (nb[i, 0] - s) / R[i, i];
}
return x;
}
public static Matrix passMethodSolver(Matrix A, Matrix b) {
int n = A.GetLength(0);
Matrix ans = new Matrix(n, 1);
double[] betas = new double[n];
double[] gammas = new double[n];
betas[0] = -A[0, 1] / A[0, 0];
gammas[0] = b[0, 0] / A[0, 0];
for (int i = 1; i < n - 1; ++i) {
betas[i] = -A[i, i + 1] / (A[i, i] + betas[i - 1] * A[i, i - 1]);
}
betas[n - 1] = 0;
for (int i = 1; i < n; ++i) {
gammas[i] = (b[i, 0] - A[i, i - 1] * gammas[i - 1]) / (A[i, i] + A[i, i - 1] * betas[i - 1]);
}
ans[n - 1, 0] = gammas[n - 1];
for (int i = n - 2; i >= 0; --i) {
ans[i, 0] = betas[i] * ans[i + 1, 0] + gammas[i];
}
return ans;
}
public static void LuDecomposition(Matrix A, ref Matrix L, ref Matrix U) {
if (A.GetLength(0) != A.GetLength(1))
return;
int n = A.GetLength(0);
L = L.MakeIdentity();
for (int i = 0; i < n; ++i) {
for (int j = i; j < n; ++j) {
double temp = A[i,j];
for (int k = 0; k < i; ++k) {
temp = temp - U[k, j] * L[i, k];
}
U[i ,j] = temp;
}
for (int j = i + 1; j < n; ++j) {
double temp = A[j, i];
for (int k = 0; k < i; ++k) {
temp = temp - U[k, i] * L[j, k];
}
temp = temp / U[i, i];
L[j, i] = temp;
}
}
}
public static Matrix inversedLuSolver(Matrix A, Matrix b) {
int n = A.GetLength(1);
Matrix L = new Matrix(n, n);
Matrix U = new Matrix(n, n);
LuDecomposition(A, ref L, ref U);
Matrix inversed = new Matrix(n, n);
for (int j = n - 1; j >= 0; --j) {
double tmp = 1;
for (int k = j + 1; k < n; ++k)
tmp -= U[j, k] * inversed[k, j];
tmp /= U[j, j];
inversed[j, j] = tmp;
for (int i = j - 1; i >= 0; --i) {
tmp = 0;
for (int k = i + 1; k < n; ++k)
tmp += U[i, k] * inversed[k, j];
tmp = -tmp / U[i, i];
inversed[i, j] = tmp;
}
for (int i = j - 1; i >= 0; --i) {
tmp = 0;
for (int k = i + 1; k < n; ++k)
tmp += inversed[j, k] * L[k, i];
inversed[j, i] = - tmp;
}
}
Matrix answer = inversed * b;
return answer;
}
}