using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace QR {
class Program {
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);
Console.WriteLine("Матрица Q:");
Console.WriteLine(Q);
Console.WriteLine("Матрица R:");
Console.WriteLine(R);
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];
}
Console.WriteLine("Решение СЛАУ:");
Console.WriteLine(x);
return x;
}
public static Matrix TryNewSolve(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];
}
x[A.GetLength(0) - 1, 0] = 1.0;
for (int i = n - 2; 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];
}
for (int i = 0; i < n; ++i)
x[i, 0] /= Math.Sqrt(x.GetSummColumn(0, 0));
return x;
}
public static List<double> GetEigenvalues(Matrix A, double eps) {
Matrix reserve = new Matrix(A);
List<double> retValue = new List<double>();
int countIter = 0;
int indColumn = A.GetLength(0) - 2;
while (true) {
bool stop = true;
for (int i = indColumn; i >= 0; --i)
if (A.GetColumn(i, i + 1).GetNorm() >= eps)
stop = false;
if (stop)
break;
if (countIter > 100) {
countIter = 0;
indColumn--;
A = reserve;
}
Matrix Q = QR(A),
R = Q.GetTransposed() * A;
A = R * Q;
countIter++;
}
for (int i = 0; i < A.GetLength(0) - 1; ++i) {
if (A.GetColumn(i, i + 1).GetNorm() < eps) {
retValue.Add(A[i, i]);
}
if (retValue.Count == A.GetLength(0) - 1)
retValue.Add(A[A.GetLength(0) - 1, A.GetLength(1) - 1]);
}
return retValue;
}
public static List<Matrix> GetEigenvectors(Matrix A, List<double> eigenvalues) {
List<Matrix> retValue = new List<Matrix>();
Matrix b = new Matrix(A.GetLength(0), 1);
for (int i = 0; i < A.GetLength(0); ++i)
b[i, 0] = 0;
Matrix reserve = new Matrix(A);
for (int i = 0; i < eigenvalues.Count; ++i) {
A = new Matrix(reserve);
A = A - A.MakeIdentity() * eigenvalues[i];
retValue.Add(TryNewSolve(A, b, 0.001));
}
return retValue;
}
static void Main(string[] args) {
Random rand = new Random();
Matrix A = new Matrix(new double[,] {{21, 18, 6},
{18, 17, 5},
{6, 5, 2}});
Matrix b = new Matrix(new double[,] {{1},
{2},
{3}});
Console.WriteLine("Исходная матрица А:");
Console.WriteLine(A);
Console.WriteLine("Вектор свободных членов:");
Console.WriteLine(b);
GetSolve(A, b, 0.00001);
List<double> eigenValues = GetEigenvalues(A, 0.0000001);
List<Matrix> eigenVectors = GetEigenvectors(A, eigenValues);
for (int i = 0; i < eigenVectors.Count; ++i) {
Console.WriteLine("Собственное значение : {0:f2}", eigenValues[i]);
Console.WriteLine("Собственный вектор :");
Console.WriteLine(eigenVectors[i]);
Console.WriteLine();
}
}
}
}