using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace NumMethods
{
class Mtr
{
double[,] data;
int n, k;
public Mtr(int n, int k)
{
this.n = n;
this.k = k;
data = new double[n, k];
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
data[i, j] = 0;
}
public Mtr(Mtr c)
{
this.n = c.n;
this.k = c.k;
data = new double[n, k];
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
data[i, j] = c.data[i,j];
}
public int getColumns()
{
return k;
}
public int getRows()
{
return n;
}
public void makeIdentity()
{
if (n != k)
{
Console.WriteLine("MAKE IDENTITY ERROR: N != K");
return;
}
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
{
data[i, j] = ((i == j) ? 1 : 0);
}
}
public Mtr(double[,] d)
{
n = d.GetLength(0);
k = d.GetLength(1);
data = new double[n, k];
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
data[i, j] = d[i, j];
}
public void setValue(double v, int i, int j)
{
data[i, j] = v;
}
public double getValue(int i, int j)
{
return data[i, j];
}
static public Mtr operator *(Mtr inp1, Mtr inp2)
{
if (inp1.k != inp2.n)
{
Console.WriteLine("MATRIX MULT ERROR");
return new Mtr(1, 1);
}
Mtr rs = new Mtr(inp1.n, inp2.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp2.k; ++j)
{
double tmp = 0;
for (int c = 0; c < inp1.k; ++c)
tmp += inp1.data[i, c] * inp2.data[c, j];
rs.setValue(tmp, i, j);
}
}
return rs;
}
static public Mtr operator *(Mtr inp1, double alpha)
{
Mtr rs = new Mtr(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] * alpha, i, j);
}
}
return rs;
}
static public Mtr operator -(Mtr inp1, double alpha)
{
Mtr rs = new Mtr(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] - alpha, i, j);
}
}
return rs;
}
static public Mtr operator -(Mtr inp1, Mtr inp2)
{
if (inp1.n != inp2.n || inp1.k != inp2.k)
{
Console.WriteLine("Operator- ERROR: Matr size not equal");
return new Mtr(1, 1);
}
Mtr rs = new Mtr(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] - inp2.data[i,j], i, j);
}
}
return rs;
}
static public Mtr operator +(Mtr inp1, Mtr inp2)
{
if (inp1.n != inp2.n || inp1.k != inp2.k)
{
Console.WriteLine("Operator+ ERROR: Matr size not equal");
return new Mtr(1, 1);
}
Mtr rs = new Mtr(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] + inp2.data[i, j], i, j);
}
}
return rs;
}
/// Tensor multiplication of vectors
static public Mtr operator ^(Mtr inp1, Mtr inp2)
{
if (inp1.k != inp2.n || inp1.n != inp2.k)
{
Console.WriteLine("Operator^ ERROR: Matr size not equal");
return new Mtr(1, 1);
}
Mtr rs = new Mtr(inp1.n, inp2.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp2.k; ++j)
{
rs.setValue(inp1.data[i, 0] * inp2.data[0, j], i, j);
}
}
return rs;
}
public Mtr getTranspozed()
{
Mtr ret = new Mtr(k, n);
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < k; ++j)
{
ret.setValue(data[i, j], j, i);
}
}
return ret;
}
public void transpoze()
{
double[,] ret = new double[k, n];
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < k; ++j)
{
ret[j, i] = data[i, j];
}
}
data = ret;
}
public Mtr getSubMatrix(int si, int sj, int ei, int ej)
{
if (si >= n || sj >= k || ei >= n || ej >= k ||
si < 0 || sj < 0 || ei < 0 || ej < 0)
{
Console.WriteLine("getSubMatrix RANGE ERROR");
return new Mtr(1, 1);
}
Mtr retMtr = new Mtr(ei - si + 1, ej - sj + 1);
int fi = 0;
int fj = 0;
for (int i = si; i <= ei; ++i)
{
fj = 0;
for (int j = sj; j <= ej; ++j)
{
retMtr.setValue(data[i, j], fi, fj);
++fj;
}
++fi;
}
return retMtr;
}
public bool isZeroMatrix()
{
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
if (data[i, j] != 0)
return false;
return true;
}
public void printMatrix()
{
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < k; ++j)
{
Console.Write("{0, 3:f} ", data[i, j]);
}
Console.WriteLine();
}
Console.WriteLine(); Console.WriteLine(); Console.WriteLine();
}
public double getMaximumAboveDiag(ref int ii, ref int jj)
{
double val = -1000 * 1000 * 1000;
for (int i = 0; i < n; ++i)
{
for (int j = i + 1; j < k; ++j)
{
if (val < Math.Abs(data[i, j]))
{
val = Math.Abs(data[i, j]);
ii = i;
jj = j;
}
}
}
return val;
}
}
class Program
{
static int mSign(double vl)
{
if (vl >= 0)
return 1;
else
return -1;
}
static bool isNull(Mtr col)
{
for (int i = 1; i < col.getRows(); ++i)
if (col.getValue(i, 0) != 0)
return false;
return true;
}
static Mtr uniteMtr(Mtr A, Mtr B)
{
int bi = A.getRows() - 1;
for (int i = B.getRows() - 1; i >= 0; --i)
{
int bj = A.getColumns() - 1;
for (int j = B.getColumns() - 1; j >= 0; --j)
{
A.setValue(B.getValue(i, j), bi, bj);
--bj;
}
--bi;
}
return A;
}
static void planarRotationMethod(Mtr A)
{
const double EPS = 0.001;
A.printMatrix();
int ii = 0, jj = 0;
Mtr B = A;
int n = A.getRows();
int k = A.getColumns();
Mtr globalT = new Mtr(n, k);
globalT.makeIdentity();
for (int iteration = 0; iteration < 300; ++iteration)
{
Mtr locT = new Mtr(n, k);
locT.makeIdentity();
double aii = B.getMaximumAboveDiag(ref ii, ref jj);
if (Math.Abs(aii) < EPS)
break;
double p = 2 * aii;
double q = B.getValue(ii, ii) - B.getValue(jj, jj);
double d = Math.Sqrt(p * p + q * q);
double r = 0, c = 0, s = 0;
if (Math.Abs(q) > EPS)
{
r = Math.Abs(q) / (2 * d);
c = Math.Sqrt(0.5f + r);
s = Math.Sqrt(0.5f - r) * Math.Sign(p * q);
}
else
{
c = s = Math.Sqrt(2.0f) / 2.0;
}
locT.setValue(c, ii, ii);
locT.setValue(c, jj, jj);
locT.setValue(-s, ii, jj);
locT.setValue(s, jj, ii);
globalT = globalT * locT;
B = locT.getTranspozed() * B * locT;
}
Console.WriteLine("Eigen numbers");
B.printMatrix();
Console.WriteLine("Eigen vectors");
globalT.printMatrix();
}
static Mtr QR_Decomposition(ref Mtr A)
{
int si, sj;
si = A.getRows();
sj = A.getColumns();
Mtr y;
Mtr U = new Mtr(si, sj);
U.makeIdentity();
for (int i = 0; i < A.getColumns() - 1; ++i)
{
y = A.getSubMatrix(i, i, A.getRows() - 1, i);
if (isNull(y))
continue;
Mtr locV = new Mtr(si - i, sj - i);
Mtr locU = new Mtr(si, sj);
locU.makeIdentity();
locV.makeIdentity();
Mtr omega = new Mtr(si - i, 1);
Mtr vecx = new Mtr(si - i, 1);
vecx.setValue(1, 0, 0);
double yNorm = Math.Sqrt((y.getTranspozed() * y).getValue(0, 0));
int sgn = mSign(y.getValue(0, 0));
omega = y + vecx * (sgn * yNorm);
double omegaNorm = Math.Sqrt((omega.getTranspozed() * omega).getValue(0, 0));
omega = omega * (1 / omegaNorm);
locV = locV - (omega ^ omega.getTranspozed()) * 2;
locU = uniteMtr(locU, locV);
A = locU * A;
U = locU * U;
}
return U;
}
static void LuDecomposition(Mtr A, ref Mtr L, ref Mtr U)
{
if (A.getRows() != A.getColumns())
return;
int n = A.getRows();
L.makeIdentity();
for (int i = 0; i < n; ++i)
{
for (int j = i; j < n; ++j)
{
double temp = A.getValue(i, j);
for (int k = 0; k < i; ++k)
{
temp = temp - U.getValue(k, j) * L.getValue(i, k);
}
U.setValue(temp, i, j);
}
for (int j = i + 1; j < n; ++j)
{
double temp = A.getValue(j, i);
for (int k = 0; k < i; ++k)
{
temp = temp - U.getValue(k, i) * L.getValue(j, k);
}
temp = temp / U.getValue(i, i);
L.setValue(temp, j, i);
}
}
}
static void inversedLuSolver(Mtr A, Mtr b)
{
int n = A.getColumns();
Mtr L = new Mtr(n, n);
Mtr U = new Mtr(n, n);
LuDecomposition(A,ref L, ref U);
Mtr inversed = new Mtr(n, n);
for (int j = n - 1; j >= 0; --j)
{
double tmp = 1;
for (int k = j + 1; k < n; ++k)
tmp -= U.getValue(j, k) * inversed.getValue(k, j);
tmp /= U.getValue(j, j);
inversed.setValue(tmp, j, j);
for (int i = j - 1; i >= 0; --i)
{
tmp = 0;
for (int k = i + 1; k < n; ++k)
tmp += U.getValue(i, k) * inversed.getValue(k, j);
tmp = -tmp / U.getValue(i, i);
inversed.setValue(tmp, i, j);
}
for (int i = j - 1; i >= 0; --i)
{
tmp = 0;
for (int k = i + 1; k < n; ++k)
tmp += inversed.getValue(j, k) * L.getValue(k, i);
inversed.setValue(-tmp, j, i);
}
}
Mtr answer = inversed * b;
Console.WriteLine("LU DECOMPOSITION SOLVER: ");
answer.printMatrix();
Console.WriteLine("-------------------------------------");
}
static void QRDecomposeSolver(Mtr A, Mtr b)
{
int SIZE = b.getRows();
Mtr U;
U = QR_Decomposition(ref A);
b = U * b;
double[] results = new double[SIZE];
for (int x = SIZE - 1; x >= 0; --x)
{
double vl = b.getValue(x, 0);
for (int j = x + 1; j < SIZE; ++j)
{
vl -= A.getValue(x, j) * results[j];
}
double dividor = A.getValue(x, x);
results[x] = vl / dividor;
}
Console.WriteLine("QR decomposition solver:");
for (int i = 0; i < SIZE; ++i)
{
Console.WriteLine("X: {0, 3:f} ", results[i]);
}
Console.WriteLine("-------------------------------------");
}
static int CubeDecompositionSolver(Mtr A, Mtr b)
{
int n = A.getRows();
Mtr U = new Mtr(n, n);
double uii;
for (int i = 0; i < n; ++i)
{
double tmp = A.getValue(i, i);
for (int k = 0; k < i; ++k)
tmp -= U.getValue(k, i) * U.getValue(k, i);
if (tmp < 0)
{
Console.WriteLine("Negative tmp in square root!");
return -1;
}
tmp = Math.Sqrt(tmp);
uii = tmp;
if (uii == 0)
{
Console.WriteLine("Zero division in uii!");
return -1;
}
U.setValue(uii, i, i);
for (int j = i + 1; j < n; ++j)
{
tmp = A.getValue(i, j);
for (int k = 0; k < i; ++k)
tmp -= U.getValue(k, i) * U.getValue(k, j);
tmp /= uii;
U.setValue(tmp, i, j);
}
}
double[] y = new double[n];
for (int i = 0; i < n; ++i)
{
double tmp = b.getValue(i, 0);
for (int j = 0; j < i; ++j)
{
tmp -= U.getValue(j, i) * y[j];
}
tmp /= U.getValue(i, i);
y[i] = tmp;
}
Mtr ans = new Mtr(n, 1);
for (int i = n - 1; i >= 0; --i)
{
double tmp = y[i];
for (int k = i + 1; k < n; ++k)
{
tmp -= U.getValue(i, k) * ans.getValue(k, 0);
}
tmp /= U.getValue(i, i);
ans.setValue(tmp, i, 0);
}
Console.WriteLine("Cube decomposition solver:");
ans.printMatrix();
Console.WriteLine("-------------------------------------");
return 1;
}
static void SimpleIteration(Mtr A, Mtr b)
{
int n = A.getRows();
Mtr B = new Mtr(n, n);
Mtr C = new Mtr(n, 1);
Mtr X = new Mtr(n, 1);
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < n; ++j)
{
if (i == j)
continue;
B.setValue(-A.getValue(i, j) / A.getValue(i, i), i, j);
}
C.setValue(b.getValue(i, 0) / A.getValue(i, i), i, 0);
X.setValue(1, i, 0);
}
for (int i = 0; i < 100; ++i)
{
X = B * X + C;
}
Console.WriteLine("Simple iteration solver:");
X.printMatrix();
Console.WriteLine("-------------------------------------");
}
static double getNorm(Mtr a)
{
return Math.Sqrt((a.getTranspozed() * a).getValue(0,0));
}
static int pmAlgorithm(Mtr A, ref double bEigen, ref Mtr vecEigen)
{
int n = A.getRows();
const double pmEPS = 0.0000001;
Mtr y = new Mtr(n, 1);
double[] diffs = new double[n];
for (int i = 0; i < n; ++i)
{ y.setValue(1, i, 0);
diffs[i] = 100;
}
double lastSum = 0;
for (int i = 0; i < 1000; ++i)
{
lastSum = 0;
Mtr x = y * (1 / getNorm(y));
Mtr ny = A * x;
bool stop = true;
for (int j = 0; j < n; ++j)
{
if (Math.Abs(x.getValue(j, 0)) > pmEPS)
{
if (Math.Abs(ny.getValue(j, 0) / x.getValue(j, 0) - diffs[j]) > pmEPS)
stop = false;
diffs[j] = ny.getValue(j, 0) / x.getValue(j, 0);
lastSum += diffs[j];
}
}
y = ny;
if (stop)
break;
}
vecEigen = y * (1 / getNorm(y)); ;
bEigen = lastSum / n;
return 1;
}
static void findMinimalByPm(Mtr A)
{
double maxEigen = 0;
Mtr res = new Mtr(1, 1);
pmAlgorithm(A, ref maxEigen, ref res);
Mtr substr = new Mtr(A.getRows(), A.getColumns());
substr.makeIdentity();
A = A - substr * maxEigen;
double nEigen = 0;
pmAlgorithm(A, ref nEigen, ref res);
Console.WriteLine(nEigen + maxEigen);
res.printMatrix();
}
static void Main(string[] args)
{
int SIZE = 2;
Mtr A = new Mtr(SIZE, SIZE);
Mtr b = new Mtr(SIZE, 1);
Mtr L = new Mtr(SIZE, SIZE);
Mtr P = new Mtr(SIZE, SIZE);
Mtr symmetrical = new Mtr(SIZE, SIZE);
Random rndd = new Random(13335);
for (int i = 0; i < A.getRows(); ++i)
{
for (int j = 0; j < A.getColumns(); ++j)
{
int vl = rndd.Next(100);
A.setValue(vl, i, j);
if (i == j)
A.setValue(vl * 30, i, j);
symmetrical.setValue(vl, i, j);
symmetrical.setValue(vl, j, i);
}
b.setValue(rndd.Next(100), i, 0);
}
inversedLuSolver(A, b);
QRDecomposeSolver(A, b);
SimpleIteration(A, b);
symmetrical.setValue(5, 0, 0);
symmetrical.setValue(1, 0, 1);
// symmetrical.setValue(2, 0, 2);
symmetrical.setValue(1, 1, 0);
symmetrical.setValue(4, 1, 1);
// symmetrical.setValue(1, 1, 2);
// symmetrical.setValue(2, 2, 0);
// symmetrical.setValue(1, 2, 1);
// symmetrical.setValue(3, 2, 2);
planarRotationMethod(symmetrical);
findMinimalByPm(symmetrical);
}
}
}
/*
* Matrices for cube decompose solver
* 5x5
694,00 690,00 140,00 611,00 466,00
690,00 962,00 529,00 693,00 815,00
140,00 529,00 911,00 243,00 638,00
611,00 693,00 243,00 674,00 527,00
466,00 815,00 638,00 527,00 983,00
* 4x4
377,00 207,00 262,00 250,00
207,00 883,00 178,00 1,00
262,00 178,00 823,00 358,00
250,00 1,00 358,00 412,00
*/