C++ Linear Algebra Package Extension Implemented by James Pate Williams, Jr.

Input file:

The dimensions of the linear system of equations (m and n, m = n):
2
2
The matrix of the linear system of equations (n by n):
1	1
1	2
The right-hand side of the linear system of equations (n by 1):
7	11
The dimensions of the linear system of equations (m and n, m = 2):
2
2
The matrix of the linear system of equations (n by n):
1	1
1	3
The right-hand side of the linear system of equations (n by 1):
7	11
The dimensions of the linear system of equations (m and n, m = n):
2
2
The matrix of the linear system of equations (n by n):
6	3
4	8
The right-hand side of the linear system of equations (n by 1):
5	6
The dimensions of the linear system of equations (m and n, m = n):
2
2
The matrix of the linear system of equations (n by n):
5	3
10	4
The right-hand side of the linear system of equations (n by 1):
8	6
The dimensions of the linear system of equations (m and n, m = n):
3
3
The matrix of the linear system of equations (n by n):
2	1	-1
-3	-1	2
-2	1	2
The right-hand side of the linear system of equations (n by 1):
8	-11	-3

Output file:

The 1st solution of the linear system of equations:
       3	       4	
The 2nd solution of the linear system of equations:
       3	       4	
The determinant of the linear system of equations:
1
The inverse of the linear system of equations:
       2	      -1	
      -1	       1	
The adjoint of the linear system of equations:
       2	      -0	
      -0	       2	
The characteristic polynomial of the linear system of equations:
       1	       2	
The image of the matrix: 
       1	      -1	
       1	       2	
Rank = 2
The 1st solution of the linear system of equations:
       5	       2	
The 2nd solution of the linear system of equations:
       5	       2	
The determinant of the linear system of equations:
2
The inverse of the linear system of equations:
     1.5	    -0.5	
    -0.5	     0.5	
The adjoint of the linear system of equations:
       3	      -0	
      -0	       3	
The characteristic polynomial of the linear system of equations:
       1	       3	
The image of the matrix: 
       1	      -1	
       1	       3	
Rank = 2
The 1st solution of the linear system of equations:
 0.61111	 0.44444	
The 2nd solution of the linear system of equations:
 0.61111	 0.44444	
The determinant of the linear system of equations:
36
The inverse of the linear system of equations:
 0.22222	-0.11111	
-0.083333	 0.16667	
The adjoint of the linear system of equations:
      48	      -0	
      -0	      48	
The characteristic polynomial of the linear system of equations:
       1	      48	
The image of the matrix: 
       6	      -1	
       4	       8	
Rank = 2
The 1st solution of the linear system of equations:
    -1.4	       5	
The 2nd solution of the linear system of equations:
    -1.4	       5	
The determinant of the linear system of equations:
-10
The inverse of the linear system of equations:
    -0.4	       1	
     0.3	    -0.5	
The adjoint of the linear system of equations:
      20	      -0	
      -0	      20	
The characteristic polynomial of the linear system of equations:
       1	      20	
The image of the matrix: 
       5	      -1	
      10	       4	
Rank = 2
The 1st solution of the linear system of equations:
       2	       3	      -1	
The 2nd solution of the linear system of equations:
       2	       3	      -1	
The determinant of the linear system of equations:
1
The inverse of the linear system of equations:
       4	       5	      -2	
       3	       4	      -2	
      -1	      -1	       1	
The adjoint of the linear system of equations:
      -4	       0	       0	
       0	      -4	       0	
       0	       0	      -4	
The characteristic polynomial of the linear system of equations:
       1	      -7	       4	
The image of the matrix: 
       2	      -1	       2	
      -3	       0	      -1	
      -2	       0	       4	
Rank = 3

// Algorithms from "A Course in Computational
// Algebraic Number Theory" by Henri Cohen
// Implemented by James Pate Williams, Jr.
// Copyright (c) 2023 All Rights Reserved

#pragma once
#include "pch.h"

template<class T> class Matrix
{
public:
	size_t m, n;
	T** data;

	Matrix() { m = 0; n = 0; data = NULL; };
	Matrix(size_t m, size_t n)
	{
		this->m = m;
		this->n = n;
		data = new T*[m];

		if (data == NULL)
			exit(-300);

		for (size_t i = 0; i < m; i++)
		{
			data[i] = new T[n];

			if (data[i] == NULL)
				exit(-301);
		}
	};
	void OutputMatrix(
		fstream& outs, char fill, int precision, int width)
	{
		for (size_t i = 0; i < m; i++)
		{
			for (size_t j = 0; j < n; j++)
			{
				outs << setfill(fill) << setprecision(precision);
				outs << setw(width) << data[i][j] << '\t';
			}

			outs << endl;
		}
	};
};

template<class T> class Vector
{
public:
	size_t n;
	T* data;

	Vector() { n = 0; data = NULL; };
	Vector(size_t n)
	{
		this->n = n;
		data = new T[n];
	};
	void OutputVector(
		fstream& outs, char fill, int precision, int width)
	{
		for (size_t i = 0; i < n; i++)
		{
			outs << setfill(fill) << setprecision(precision);
			outs << setw(width) << data[i] << '\t';
		}

		outs << endl;
	};
};

class LinearAlgebra
{
public:
	bool initialized;
	size_t m, n;
	Matrix<double> M;
	Vector<double> B;

	LinearAlgebra()
	{ 
		initialized = false;
		m = 0; n = 0;
		M.data = NULL;
		B.data = NULL;
	};
	LinearAlgebra(size_t m, size_t n)
	{
		initialized = false;
		this->m = m;
		this->n = n;
		this->M.m = m;
		this->M.n = n;
		this->B.n = n;
		this->M.data = new double*[m];
		this->B.data = new double[n];

		if (M.data != NULL)
		{
			for (size_t i = 0; i < m; i++)
			{
				this->M.data[i] = new double[n];

				for (size_t j = 0; j < n; j++)
					this->M.data[i][j] = 0;
			}
		}

		if (B.data != NULL)
		{
			this->B.data = new double[n];

			for (size_t i = 0; i < n; i++)
				this->B.data[i] = 0;
		}

		initialized = this->B.data != NULL && this->M.data != NULL;
	};
	LinearAlgebra(
		size_t m, size_t n,
		double** M,
		double* B)
	{
		this->m = m;
		this->n = n;
		this->M.m = m;
		this->M.n = n;
		this->M.data = new double*[m];
		this->B.data = new double[n];

		if (M != NULL)
		{
			for (size_t i = 0; i < m; i++)
			{
				this->M.data[i] = new double[n];

				for (size_t j = 0; j < n; j++)
					this->M.data[i][j] = M[i][j];
			}
		}

		if (B != NULL)
		{
			this->B.data = new double[n];

			for (size_t i = 0; i < m; i++)
				this->B.data[i] = B[i];
		}

		initialized = this->B.data != NULL && this->M.data != NULL;
	}
	~LinearAlgebra()
	{
		M.m = m;
		M.n = n;
		B.n = n;

		if (B.data != NULL)
			delete[] B.data;

		for (size_t i = 0; i < m; i++)
			if (M.data[i] != NULL)
				delete[] M.data[i];

		if (M.data != NULL)
			delete[] M.data;
	}
	double DblDeterminant(bool failure);
	Vector<double> DblGaussianElimination(
		bool& failure);
	// The following four methods are from the
	// textbook "Elementary Numerical Analysis
	// An Algorithmic Approach" by S. D. Conte
	// and Carl de Boor Translated from the
	// original FORTRAN by James Pate Williams, Jr.
	// Copyright (c) 2023 All Rights Reserved
	bool DblGaussianFactor(
		Vector<int>& pivot);
	bool DblGaussianSolution(
		Vector<double>& x,
		Vector<int>& pivot);
	bool DblSubstitution(
		Vector<double>& x,
		Vector<int>& pivot);
	bool DblInverse(
		Matrix<double>& A,
		Vector<int>& pivot);
	// Henri Cohen Algorithm 2.2.7
	void DblCharPolyAndAdjoint(
		Matrix<double>& adjoint,
		Vector<double>& a);
	// Henri Cohen Algorithm 2.3.1
	void DblMatrixKernel(
		Matrix<double>& X, size_t& r);
	// Henri Cohen Algorithm 2.3.1
	void DblMatrixImage(
		Matrix<double>& N, size_t& rank);
};

// pch.h: This is a precompiled header file.
// Files listed below are compiled only once, improving build performance for future builds.
// This also affects IntelliSense performance, including code completion and many code browsing features.
// However, files listed here are ALL re-compiled if any one of them is updated between builds.
// Do not add files here that you will be updating frequently as this negates the performance advantage.

#ifndef PCH_H
#define PCH_H
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
using namespace std;
#endif //PCH_H

#include "pch.h"
#include "LinearAlgebra.h"

double LinearAlgebra::DblDeterminant(
    bool failure)
{
    double deter = 1;
    Vector<int> pivot(n);

    if (!initialized || m != n)
    {
        failure = true;
        return 0.0;
    }

    if (!DblGaussianFactor(pivot))
    {
        failure = true;
        return 0.0;
    }

    for (size_t i = 0; i < n; i++)
        deter *= M.data[i][i];

    return deter;
}

Vector<double> LinearAlgebra::DblGaussianElimination(
    bool& failure)
{
    double* C = new double[m];
    Vector<double> X(n);

    X.data = new double[n];

    if (X.data == NULL)
        exit(-200);

    if (!initialized)
    {
        failure = true;
        delete[] C;
        return X;
    }
    
    for (size_t i = 0; i < m; i++)
        C[i] = -1;

    failure = false;

    for (size_t j = 0; j < n; j++)
    {
        bool found = false;
        size_t i = j;

        while (i < n && !found)
        {
            if (M.data[i][j] != 0)
                found = true;
            else
                i++;
        }

        if (!found)
        {
            failure = true;
            break;
        }

        if (i > j)
        {
            for (size_t l = j; l < n; l++)
            {
                double t = M.data[i][l];
                M.data[i][l] = M.data[j][l];
                M.data[j][l] = t;
                t = B.data[i];
                B.data[i] = B.data[j];
                B.data[j] = t;
            }
        }

        double d = 1.0 / M.data[j][j];

        for (size_t k = j + 1; k < n; k++)
            C[k] = d * M.data[k][j];

        for (size_t k = j + 1; k < n; k++)
        {
            for (size_t l = j + 1; l < n; l++)
                M.data[k][l] = M.data[k][l] - C[k] * M.data[j][l];
            
            B.data[k] = B.data[k] - C[k] * B.data[j];
        }
    }

    for (int i = (int)n - 1; i >= 0; i--)
    {
        double sum = 0;

        for (size_t j = i + 1; j < n; j++)
            sum += M.data[i][j] * X.data[j];

        X.data[i] = (B.data[i] - sum) / M.data[i][i];
    }

    delete[] C;
    return X;
}

bool LinearAlgebra::DblGaussianFactor(
    Vector<int>& pivot)
    // returns false if matrix is singular
{
    Vector<double> d(n);
    double awikod, col_max, ratio, row_max, temp;
    int flag = 1;
    size_t i_star, itemp;

    for (size_t i = 0; i < n; i++)
    {
        pivot.data[i] = i;
        row_max = 0;
        for (size_t j = 0; j < n; j++)
            row_max = max(row_max, abs(M.data[i][j]));
        if (row_max == 0)
        {
            flag = 0;
            row_max = 1;
        }
        d.data[i] = row_max;
    }
    if (n <= 1) return flag != 0;
    // factorization
    for (size_t k = 0; k < n - 1; k++)
    {
        // determine pivot row the row i_star
        col_max = abs(M.data[k][k]) / d.data[k];
        i_star = k;
        for (size_t i = k + 1; i < n; i++)
        {
            awikod = abs(M.data[i][k]) / d.data[i];
            if (awikod > col_max)
            {
                col_max = awikod;
                i_star = i;
            }
        }
        if (col_max == 0)
            flag = 0;
        else
        {
            if (i_star > k)
            {
                // make k the pivot row by
                // interchanging with i_star
                flag *= -1;
                itemp = pivot.data[i_star];
                pivot.data[i_star] = pivot.data[k];
                pivot.data[k] = itemp;
                temp = d.data[i_star];
                d.data[i_star] = d.data[k];
                d.data[k] = temp;
                for (size_t j = 0; j < n; j++)
                {
                    temp = M.data[i_star][j];
                    M.data[i_star][j] = M.data[k][j];
                    M.data[k][j] = temp;
                }
            }
            // eliminate x[k]
            for (size_t i = k + 1; i < n; i++)
            {
                M.data[i][k] /= M.data[k][k];
                ratio = M.data[i][k];
                for (size_t j = k + 1; j < n; j++)
                    M.data[i][j] -= ratio * M.data[k][j];
            }
        }

        if (M.data[n - 1][n - 1] == 0) flag = 0;
    }

    if (flag == 0)
        return false;

    return true;
}

bool LinearAlgebra::DblGaussianSolution(
    Vector<double>& x,
    Vector<int>& pivot)
{
    if (!DblGaussianFactor(pivot))
        return false;

    return DblSubstitution(x, pivot);
}

bool LinearAlgebra::DblSubstitution(
    Vector<double>& x,
    Vector<int>& pivot)
{
    double sum;
    size_t j, n1 = n - 1;

    if (n == 1)
    {
        x.data[0] = B.data[0] / M.data[0][0];
        return true;
    }

    // forward substitution

    x.data[0] = B.data[pivot.data[0]];

    for (int i = 1; i < (int)n; i++)
    {
        for (j = 0, sum = 0; j < (size_t)i; j++)
            sum += M.data[i][j] * x.data[j];

        x.data[i] = B.data[pivot.data[i]] - sum;
    }

    // backward substitution

    x.data[n1] /= M.data[n1][n1];

    for (int i = n - 2; i >= 0; i--)
    {
        double sum = 0.0;

        for (j = i + 1; j < n; j++)
            sum += M.data[i][j] * x.data[j];

        x.data[i] = (x.data[i] - sum) / M.data[i][i];
    }

    return true;
}

bool LinearAlgebra::DblInverse(
    Matrix<double>& A,
    Vector<int>& pivot)
{
    Vector<double> x(n);

    if (!DblGaussianFactor(pivot))
        return false;

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            B.data[j] = 0;
    }
    
    for (size_t i = 0; i < n; i++)
    {
        B.data[i] = 1;

        if (!DblSubstitution(x, pivot))
            return false;

        B.data[i] = 0;

        for (size_t j = 0; j < n; j++)
           A.data[i][j] = x.data[pivot.data[j]];
    }

    return true;
}

void LinearAlgebra::DblCharPolyAndAdjoint(
    Matrix<double>& adjoint,
    Vector<double>& a)
{
    Matrix<double> C(n, n);
    Matrix<double> I(n, n);

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            C.data[i][j] = I.data[i][j] = 0;
    }

    for (size_t i = 0; i < n; i++)
        C.data[i][i] = I.data[i][i] = 1;

    a.data[0] = 1;

    for (size_t i = 1; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
            {
                double sum = 0.0;

                for (size_t l = 0; l < n; l++)
                    sum += M.data[j][l] * C.data[l][k];

                C.data[j][k] = sum;
            }
        }

        double tr = 0.0;

        for (size_t j = 0; j < n; j++)
            tr += C.data[j][j];

        a.data[i] = -tr / i;

        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
                C.data[j][k] += a.data[i] * I.data[j][k];
        }
    }

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            double sum = 0.0;

            for (size_t k = 0; k < n; k++)
                sum += M.data[i][k] * C.data[k][j];

            C.data[i][j] = sum;
        }
    }

    double trace = 0.0;

    for (size_t i = 0; i < n; i++)
        trace += C.data[i][i];

    trace /= n;
    a.data[n - 1] = -trace;

    double factor = 1.0;

    if ((n - 1) % 2 != 0)
        factor = -1.0;

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            adjoint.data[i][j] = factor * C.data[i][j];
    }
}

void LinearAlgebra::DblMatrixKernel(Matrix<double>& X, size_t& r)
{
    double D = 0.0;
    Vector<int> c(m);
    Vector<int> d(n);

    r = 0;

    for (size_t i = 0; i < m; i++)
        c.data[i] = -1;

    size_t j, k = 1;
Step2:
    for (j = 0; j < m; j++)
    {
        if (M.data[j][k] != 0 && c.data[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d.data[k] = 0;
        goto Step4;
    }

    D = -1.0 / M.data[j][k];

    M.data[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        M.data[j][s] = D * M.data[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = M.data[i][k];
                M.data[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            M.data[i][s] += D * M.data[j][s];
        }
    }

    c.data[j] = k;
    d.data[k] = j;
Step4:
    if (k < n - 1)
    {
        k++;
        goto Step2;
    }

    X.n = n;

    if (r != 0)
    {
        for (k = 0; k < n; k++)
        {
            if (d.data[k] == 0)
            {
                for (size_t i = 0; i < n; i++)
                {
                    if (d.data[i] > 0)
                        X.data[k][i] = M.data[d.data[i]][k];
                    else if (i == k)
                        X.data[k][i] = 1;
                    else
                        X.data[k][i] = 0;
                }
            }
        }
    }
}

void LinearAlgebra::DblMatrixImage(
    Matrix<double>& N, size_t& rank)
{
    double D = 0.0;
    size_t r = 0;
    Matrix<double> copyM(m, n);
    Vector<int> c(m);
    Vector<int> d(n);

    for (size_t i = 0; i < m; i++)
        c.data[i] = -1;

    size_t j, k = 1;
    N = copyM = M;
Step2:
    for (j = 0; j < m; j++)
    {
        if (copyM.data[j][k] != 0 && c.data[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d.data[k] = 0;
        goto Step4;
    }

    D = -1.0 / copyM.data[j][k];

    copyM.data[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        copyM.data[j][s] = D * copyM.data[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = copyM.data[i][k];
                copyM.data[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            copyM.data[i][s] += D * copyM.data[j][s];
        }
    }

    c.data[j] = k;
    d.data[k] = j;
Step4:
    if (k < n - 1)
    {
        k++;
        goto Step2;
    }

    rank = n - r;

    for (j = 0; j < m; j++)
    {
        if (c.data[j] != 0)
        {
            for (size_t i = 0; i < m; i++)
            {
                N.data[i][c.data[j]] = M.data[i][c.data[j]];
            }
        }
    }
}

/*
** Cohen's linear algebra test program
** Implemented by James Pate Williams, Jr.
** Copyright (c) 2023 All Rights Reserved
*/

#include "pch.h"
#include "LinearAlgebra.h"

double GetDblNumber(fstream& inps)
{
    char ch = inps.get();
    string numberStr;

    while (ch == ' ' || ch == '\t' || ch == '\r' || ch == '\n')
        ch = inps.get();

    while (ch == '+' || ch == '-' || ch == '.' ||
        ch >= '0' && ch <= '9')
    {
        numberStr += ch;
        ch = inps.get();
    }

    double x = atof(numberStr.c_str());
    return x;
}

int GetIntNumber(fstream& inps)
{
    char ch = inps.get();
    string numberStr;

    while (ch == ' ' || ch == '\t' || ch == '\r' || ch == '\n')
        ch = inps.get();

    while (ch == '+' || ch == '-' || ch >= '0' && ch <= '9')
    {
        numberStr += ch;
        ch = inps.get();
    }

    int x = atoi(numberStr.c_str());
    return x;
}

int main()
{
    fstream inps;

    inps.open("CLATestFile.txt", fstream::in);
    
    if (inps.fail())
    {
        cout << "Input file opening error!" << endl;
        return -1;
    }

    fstream outs;

    outs.open("CLAResuFile.txt", fstream::out | fstream::ate);

    if (outs.fail())
    {
        cout << "Output file opening error!" << endl;
        return -2;
    }

    size_t m, n;
    
    while (!inps.eof())
    {
        char buffer[256] = { '\0' };

        inps.getline(buffer, 256);
        m = GetIntNumber(inps);

        if (inps.eof())
            return 0;

        if (m < 1)
        {
            cout << "The number of rows must be >= 1" << endl;
            return -100;
        }

        n = GetIntNumber(inps);

        if (n < 1)
        {
            cout << "The number of colums must be >= 1" << endl;
            return -101;
        }

        LinearAlgebra la(m, n);
        Matrix<double> copyM(m, n);
        Vector<double> copyB(n);

        inps.getline(buffer, 256);

        for (size_t i = 0; i < m; i++)
        {
            for (size_t j = 0; j < n; j++)
            {
                double x = GetDblNumber(inps);

                la.M.data[i][j] = x;
                copyM.data[i][j] = x;
            }
        }

        inps.getline(buffer, 256);

        for (size_t i = 0; i < n; i++)
        {
            la.B.data[i] = GetDblNumber(inps);
            copyB.data[i] = la.B.data[i];
        }

        bool failure = false;
        Vector<double> X = la.DblGaussianElimination(failure);

        if (!failure)
        {
            outs << "The 1st solution of the linear system of equations:" << endl;
            X.OutputVector(outs, ' ', 5, 8);
        }
        else
        {
            cout << "Cohen Gaussian elimination failure!" << endl;
            exit(-102);
        }

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        Matrix<double> A(n, n);
        Vector<int> pivot(n);

        if (!la.DblGaussianSolution(X, pivot))
            exit(-103);

        outs << "The 2nd solution of the linear system of equations:" << endl;

        X.OutputVector(outs, ' ', 5, 8);

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        double deter = la.DblDeterminant(failure);

        outs << "The determinant of the linear system of equations:" << endl;
        outs << deter << endl;

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        outs << "The inverse of the linear system of equations:" << endl;

        if (!la.DblInverse(A, pivot))
        {
            cout << "Conte Gaussian inverse matrix failure!" << endl;
            exit(-104);
        }

        else
            A.OutputMatrix(outs, ' ', 5, 8);

        Matrix<double> adjoint(n, n);
        Vector<double> a(n);

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        la.DblCharPolyAndAdjoint(adjoint, a);
        outs << "The adjoint of the linear system of equations:" << endl;
        adjoint.OutputMatrix(outs, ' ', 5, 8);
        outs << "The characteristic polynomial of the linear system of equations:" << endl;
        a.OutputVector(outs, ' ', 5, 8);

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        Matrix<double> kernel(m, n);
        size_t r = 0;

        la.DblMatrixKernel(kernel, r);

        if (r > 0)
        {
            outs << "The kernel of the matrix: " << endl;
            kernel.OutputMatrix(outs, ' ', 5, 8);
            outs << "Dimension of the kernel: " << r << endl;
        }

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        Matrix<double> N(m, n);
        size_t rank;

        la.DblMatrixImage(N, rank);

        if (rank > 0)
        {
            outs << "The image of the matrix: " << endl;
            N.OutputMatrix(outs, ' ', 5, 8);
            outs << "Rank = " << rank << endl;
        }
    }

    inps.close();
    outs.close();
}
Unknown's avatar

Author: jamespatewilliamsjr

My whole legal name is James Pate Williams, Jr. I was born in LaGrange, Georgia approximately 70 years ago. I barely graduated from LaGrange High School with low marks in June 1971. Later in June 1979, I graduated from LaGrange College with a Bachelor of Arts in Chemistry with a little over a 3 out 4 Grade Point Average (GPA). In the Spring Quarter of 1978, I taught myself how to program a Texas Instruments desktop programmable calculator and in the Summer Quarter of 1978 I taught myself Dayton BASIC (Beginner's All-purpose Symbolic Instruction Code) on LaGrange College's Data General Eclipse minicomputer. I took courses in BASIC in the Fall Quarter of 1978 and FORTRAN IV (Formula Translator IV) in the Winter Quarter of 1979. Professor Kenneth Cooper, a genius poly-scientist taught me a course in the Intel 8085 microprocessor architecture and assembly and machine language. We would hand assemble our programs and insert the resulting machine code into our crude wooden box computer which was designed and built by Professor Cooper. From 1990 to 1994 I earned a Bachelor of Science in Computer Science from LaGrange College. I had a 4 out of 4 GPA in the period 1990 to 1994. I took courses in C, COBOL, and Pascal during my BS work. After graduating from LaGrange College a second time in May 1994, I taught myself C++. In December 1995, I started using the Internet and taught myself client-server programming. I created a website in 1997 which had C and C# implementations of algorithms from the "Handbook of Applied Cryptography" by Alfred J. Menezes, et. al., and some other cryptography and number theory textbooks and treatises.

Leave a comment