SqRoot.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Linq.Expressions;
using System.Diagnostics;

namespace SquareRoot
{
    class SystemCount
    {

        public static void GetSMatrix(double[,] MatrixA, double[,] MatrixS,int Dimension)
        {
            double summ;    

            for(int i = 0; i < Dimension; i++){     //columns
                for(int j = i; j < Dimension; j++){     //rows

                    summ = 0.0;

                    if (i == j){
                        
                        for(int k = 0; k < i; k++)
                            summ += MatrixS[k,i]*MatrixS[k,i];

                        if ((MatrixA[i,i] - summ) < 0.0){
                                Console.WriteLine("Matrix minor is less the zero");
                                return;}

                        MatrixS[i,i] = Math.Sqrt(MatrixA[i,i] - summ);}
                    else{
                        for(int k = 0; k < i; k++)
                            summ += MatrixS[k,i]*MatrixS[k,j];

                        if (MatrixS[i,i] == 0.0){
                                Console.WriteLine("MatrixS contains zero diagonal values");
                                return;}
                                
                        MatrixS[i,j] = (MatrixA[i,j]-summ)/MatrixS[i,i];}
                }    
            }
        }

        public static void GetVectorY(double[,] MatrixS, double[] VectorB, double[] VectorY, int Dimension)
        {
            double summ = 0.0;

            for(int i = 0; i < Dimension; i++)
            {
                for(int k = 0; k < i; k++)
                    summ += MatrixS[k,i]*VectorY[k];
            
                VectorY[i] = (VectorB[i]-summ)/MatrixS[i,i];

            }

            return;
        }

        public static void GetSolution(double[,] MatrixS, double[] Solution, double[] VectorY, int Dimension)
        {
            double summ = 0.0;
            for(int i = Dimension-1; i >= 0; i--)
            {

                for(int k = i+1; k < Dimension; k++)
                    summ =+ MatrixS[i,k]*Solution[k];

                Solution[i] = (VectorY[i] - summ)/MatrixS[i,i];

            }

            return;
        }

        public static bool ChechMartrixA(double[,] MatrixA, int Dimension)
        {
            for(int i =0;i<Dimension;i++){
                for(int k =0;k<Dimension;k++){
                    if(MatrixA[i,k] != MatrixA[k,i]){
                    Console.WriteLine("Matrix is not symmetrical!");
                    return false;    
                    }
                }
            }
            return true;
        }

        public static void TestMethod()
        {
            const int TestDimension = 3;

            double[ , ] MatrixA = new double[TestDimension,TestDimension] {{1,0,-1},{0,4,4},{-1,4,14}};
            double[ , ] MatrixS = new double[TestDimension,TestDimension];

            double[] Solution = new double[TestDimension];
            double[] VectorB = new double[] {2,12,19};
            double[] VectorY = new double[TestDimension];

            double[,] MatrixSCheck = new double[TestDimension,TestDimension] {{1,0,-1},{0,2,2},{0,0,3}};
            GetSMatrix(MatrixA,MatrixS,TestDimension);
            for(int i =0;i<TestDimension;i++){
                for(int k =0;k<TestDimension;k++){
                    Debug.Assert(Math.Abs(MatrixS[i,k] - MatrixSCheck[i,k]) < 0.001);}}

            double[] VectorYCheck = new double [TestDimension] {2,6,3};
            GetVectorY(MatrixS,VectorB,VectorY,TestDimension);   
            for (int i = 0; i < TestDimension; i++)
                Debug.Assert(Math.Abs(VectorY[i]-VectorYCheck[i]) < 0.001); 
                
            double[] SolutionCheck = new double[TestDimension] {3,2,1};
            GetSolution(MatrixS,Solution,VectorY,TestDimension);
            for (int i = 0; i < TestDimension; i++)
                Debug.Assert(Math.Abs(Solution[i]-SolutionCheck[i]) < 0.001);                          

            return;
        }

        static void Main(string[] args)
        {
            const int Dimension = 3;

            // Ax = y
            double[ , ] MatrixA = new double[Dimension,Dimension] {{1,0,-1},{0,4,4},{-1,4,14}};
            double[ , ] MatrixS = new double[Dimension,Dimension];

            double[] Solution = new double[Dimension];
            double[] VectorB = new double[] {2,12,19};
            double[] VectorY = new double[Dimension];

            TestMethod();

            if(ChechMartrixA(MatrixA,Dimension)){
                GetSMatrix(MatrixA,MatrixS,Dimension);
                GetVectorY(MatrixS,VectorB,VectorY,Dimension);
                Console.WriteLine("Solution vector is: ");
                GetSolution(MatrixS,Solution,VectorY,Dimension);}
                for(int j = 0; j < Dimension; j++)
                    Console.WriteLine("{0} ",Solution[j]);

        }   
    }
}