My reference is Numerical Analysis: An Algorithmic Approach 3rd Edition by S. D. Conte and Carl de Boor Chapter 9.1.
Finite Difference Method for Solving Second Order Ordinary Differential Equations by James Pate Williams, Jr.
Category: Matrix Algebra
Determinant Calculations: Bareiss Algorithm and Gaussian Elimination by James Pate Williams, Jr.
Below are three screenshots of two methods of calculating the determinant of a matrix, namely the Bareiss Algorithm and Gaussian Elimination:
using System;
using System.Diagnostics;
using System.Windows.Forms;
namespace MatrixInverseComparison
{
public partial class MainForm : Form
{
public MainForm()
{
InitializeComponent();
}
static private string FormatNumber(double x)
{
string result = string.Empty;
if (x > 0)
result += x.ToString("F5").PadLeft(10);
else
result += x.ToString("F5").PadLeft(10);
return result;
}
private void MultiplyPrintMatricies(
double[,] A, double[,] B, int n)
{
double[,] I = new double[n, n];
textBox1.Text += "Matrix Product:\r\n";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
double sum = 0.0;
for (int k = 0; k < n; k++)
sum += A[i, k] * B[k, j];
I[i, j] = sum;
}
}
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
textBox1.Text += FormatNumber(I[i, j]) + " ";
}
textBox1.Text += "\r\n";
}
}
private void PrintMatrix(string title, double[,] A, int n)
{
textBox1.Text += title + "\r\n";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
textBox1.Text += FormatNumber(A[i, j]) + " ";
}
textBox1.Text += "\r\n";
}
}
private void Compute(double[,] MI, int n)
{
double determinantGE = 1;
double[] b = new double[n];
double[] x = new double[n];
double[,] MB1 = new double[n, n];
double[,] MB2 = new double[n, n];
double[,] MG = new double[n, n];
double[,] MS = new double[n, n];
double[,] IG = new double[n, n];
double[,] IS = new double[n, n];
int[] pivot = new int[n];
Stopwatch sw = new();
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
MB1[i, j] = MB2[i, j] = MG[i, j] = MS[i, j] = MI[i, j];
}
PrintMatrix("Initial Matrix: ", MI, n);
sw.Start();
int flag = DirectMethods.Factor(MG, n, pivot);
if (flag != 0)
{
for (int i = 0; i < n; i++)
determinantGE *= MG[i, i];
determinantGE *= flag;
}
sw.Stop();
PrintMatrix("Gaussian Elimination Final:", MG, n);
for (int i = 0; i < n; i++)
b[i] = 0;
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
b[j] = 1.0;
DirectMethods.Substitute(MG, b, x, n, pivot);
for (int k = 0; k < n; k++)
IG[k, j] = x[k];
b[j] = 0.0;
}
}
PrintMatrix("Gaussian Elimination Inverse:", IG, n);
textBox1.Text += "Determinant: " +
FormatNumber(determinantGE) + "\r\n";
MultiplyPrintMatricies(IG, MI, n);
textBox1.Text += "Runtime (MS) = " +
sw.ElapsedMilliseconds + "\r\n";
sw.Start();
double determinant1 = BareissAlgorithm.Determinant1(MB1, n);
sw.Stop();
PrintMatrix("Bareiss Algorithm Final 1:", MB1, n);
textBox1.Text += "Determinant: " +
FormatNumber(determinant1) + "\r\n";
textBox1.Text += "Runtime (MS) = " +
sw.ElapsedMilliseconds + "\r\n";
sw.Start();
double determinant2 = BareissAlgorithm.Determinant2(MB2, n);
sw.Stop();
PrintMatrix("Bareiss Algorithm Final 2:", MB2, n);
textBox1.Text += "Determinant: " +
FormatNumber(determinant2) + "\r\n";
textBox1.Text += "Runtime (MS) = " +
sw.ElapsedMilliseconds + "\r\n";
}
private void button1_Click(object sender, EventArgs e)
{
int n = (int)numericUpDown1.Value;
int seed = (int)numericUpDown2.Value;
double[,] A = new double[n, n];
Random random = new Random(seed);
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
double q = n * random.NextDouble();
while (q == 0.0)
q = n * random.NextDouble();
A[i, j] = q;
}
}
Compute(A, n);
}
private void button2_Click(object sender, EventArgs e)
{
textBox1.Text = string.Empty;
}
}
}
using System;
namespace MatrixInverseComparison
{
public class DirectMethods
{
// Substitute and Factor translated from FORTRAN 77
// source code found in "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte and Carl de
// Boor. Translator: James Pate Williams, Jr. (c)
// August 14 - 17, 2023
static public void Substitute(
double[,] w,
double[] b,
double[] x,
int n,
int[] pivot)
{
double sum;
int i, j, n1 = n - 1;
if (n == 1)
{
x[0] = b[0] / w[0, 0];
return;
}
// forward substitution
x[0] = b[pivot[0]];
for (i = 0; i < n; i++)
{
for (j = 0, sum = 0.0; j < i; j++)
sum += w[i, j] * x[j];
x[i] = b[pivot[i]] - sum;
}
// backward substitution
x[n1] /= w[n1, n1];
for (i = n - 2; i >= 0; i--)
{
for (j = i + 1, sum = 0.0; j < n; j++)
sum += w[i, j] * x[j];
x[i] = (x[i] - sum) / w[i, i];
}
}
// Factor returns +1 if an even number of exchanges
// Factor returns -1 if an odd number of exchanges
// Factor retrurn 0 if matrix is singular
static public int Factor(
double[,] w, int n, int[] pivot)
// returns 0 if matrix is singular
{
double awikod, col_max, ratio, row_max, temp;
double[] d = new double[n];
int flag = 1, i, i_star, j, k;
for (i = 0; i < n; i++)
{
pivot[i] = i;
row_max = 0;
for (j = 0; j < n; j++)
row_max = Math.Max(row_max, Math.Abs(w[i, j]));
if (row_max == 0)
{
flag = 0;
row_max = 1;
}
d[i] = row_max;
}
if (n <= 1)
return flag;
// factorization
for (k = 0; k < n - 1; k++)
{
// determine pivot row the row i_star
col_max = Math.Abs(w[k, k]) / d[k];
i_star = k;
for (i = k + 1; i < n; i++)
{
awikod = Math.Abs(w[i, k]) / d[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;
i = pivot[i_star];
pivot[i_star] = pivot[k];
pivot[k] = i;
temp = d[i_star];
d[i_star] = d[k];
d[k] = temp;
for (j = 0; j < n; j++)
{
temp = w[i_star, j];
w[i_star, j] = w[k, j];
w[k, j] = temp;
}
}
// eliminate x[k]
for (i = k + 1; i < n; i++)
{
w[i, k] /= w[k, k];
ratio = w[i, k];
for (j = k + 1; j < n; j++)
w[i, j] -= ratio * w[k, j];
}
}
}
if (w[n - 1, n - 1] == 0)
flag = 0;
return flag;
}
}
}
namespace MatrixInverseComparison
{
// One implementation is based on https://en.wikipedia.org/wiki/Bareiss_algorithm
// Another perhaps better implementation is found on the following webpage
// https://cs.stackexchange.com/questions/124759/determinant-calculation-bareiss-vs-gauss-algorithm
class BareissAlgorithm
{
static public double Determinant1(double[,] M, int n)
{
double M00;
for (int k = 0; k < n; k++)
{
if (k - 1 == -1)
M00 = 1;
else
M00 = M[k - 1, k - 1];
for (int i = k + 1; i < n; i++)
{
for (int j = k + 1; j < n; j++)
{
M[i, j] = (
M[i, j] * M[k, k] -
M[i, k] * M[k, j]) / M00;
}
}
}
return M[n - 1, n - 1];
}
static public double Determinant2(double[,] A, int dim)
{
if (dim <= 0)
{
return 0;
}
double sign = 1;
for (int k = 0; k < dim - 1; k++)
{
//Pivot - row swap needed
if (A[k, k] == 0)
{
int m;
for (m = k + 1; m < dim; m++)
{
if (A[m, k] != 0)
{
double[] tempRow = new double[dim];
for (int i = 0; i < dim; i++)
tempRow[i] = A[m, i];
for (int i = 0; i < dim; i++)
A[m, i] = A[k, i];
for (int i = 0; i < dim; i++)
A[k, i] = tempRow[i];
sign = -sign;
break;
}
}
//No entries != 0 found in column k -> det = 0
if (m == dim)
{
return 0;
}
}
//Apply formula
for (int i = k + 1; i < dim; i++)
{
for (int j = k + 1; j < dim; j++)
{
A[i, j] = A[k, k] * A[i, j] - A[i, k] * A[k, j];
if (k >= 1)
{
A[i, j] /= A[k - 1, k - 1];
}
}
}
}
return sign * A[dim - 1, dim - 1];
}
}
}
Blast from the Past 1997 Image of a Matrix by James Pate Williams, Jr.
/*
Author: Pate Williams (c) 1997
"Algorithm 2.3.2 (Image of a Matrix). Given an
m by n matrix M = (m[i][i]) with 1 <= i <= m and
1 <= j <= n having coefficients in the field K,
this algorithm outputs a basis of the image of
M, i. e. vector space spanned by the columns of
M." -Henri Cohen- See "A Course in Computational
Algebraic Number Theory" by Henri Cohen pages
58-59. We specialize the code to the field Zp.
*/
#include <stdio.h>
#include <stdlib.h>
long** create_matrix(long m, long n)
{
long i, ** matrix = (long**)calloc(m, sizeof(long*));
if (!matrix) {
fprintf(stderr, "fatal error\ninsufficient memory\n");
fprintf(stderr, "from create_matrix\n");
exit(1);
}
for (i = 0; i < m; i++) {
matrix[i] = (long*)calloc(n, sizeof(long));
if (!matrix[i]) {
fprintf(stderr, "fatal error\ninsufficient memory\n");
fprintf(stderr, "from create_matrix\n");
exit(1);
}
}
return matrix;
}
void delete_matrix(long m, long** matrix)
{
long i;
for (i = 0; i < m; i++) free(matrix[i]);
free(matrix);
}
void Euclid_extended(long a, long b, long* u,
long* v, long* d)
{
long q, t1, t3, v1, v3;
*u = 1, * d = a;
if (b == 0) {
*v = 0;
return;
}
v1 = 0, v3 = b;
#ifdef DEBUG
printf("----------------------------------\n");
printf(" q t3 *u *d t1 v1 v3\n");
printf("----------------------------------\n");
#endif
while (v3 != 0) {
q = *d / v3;
t3 = *d - q * v3;
t1 = *u - q * v1;
*u = v1, * d = v3;
#ifdef DEBUG
printf("%4ld %4ld %4ld ", q, t3, *u);
printf("%4ld %4ld %4ld %4ld\n", *d, t1, v1, v3);
#endif
v1 = t1, v3 = t3;
}
*v = (*d - a * *u) / b;
#ifdef DEBUG
printf("----------------------------------\n");
#endif
}
long inv(long number, long modulus)
{
long d, u, v;
Euclid_extended(number, modulus, &u, &v, &d);
if (d == 1) return u;
return 0;
}
void image(long m, long n, long p,
long** M, long** X, long* r)
{
int found;
long D, i, j, k, s;
long* c = (long*)calloc(m, sizeof(long));
long* d = (long*)calloc(n, sizeof(long));
long** N = create_matrix(m, n);
if (!c || !d) {
fprintf(stderr, "fatal error\ninsufficient memory\n");
fprintf(stderr, "from kernel\n");
exit(1);
}
for (i = 0; i < m; i++) {
c[i] = -1;
for (j = 0; j < n; j++) N[i][j] = M[i][j];
}
*r = 0;
for (k = 0; k < n; k++) {
found = 0, j = 0;
while (!found && j < m) {
found = M[j][k] != 0 && c[j] == -1;
if (!found) j++;
}
if (found) {
D = p - inv(M[j][k], p);
M[j][k] = p - 1;
for (s = k + 1; s < n; s++)
M[j][s] = (D * M[j][s]) % p;
for (i = 0; i < m; i++) {
if (i != j) {
D = M[i][k];
M[i][k] = 0;
for (s = k + 1; s < n; s++) {
M[i][s] = (M[i][s] + D * M[j][s]) % p;
if (M[i][s] < 0) M[i][s] += p;
}
}
}
c[j] = k;
d[k] = j;
}
else {
*r = *r + 1;
d[k] = -1;
}
}
for (j = 0; j < m; j++) {
if (c[j] != -1) {
for (i = 0; i < n; i++) {
if (i < m) X[i][j] = N[i][c[j]];
else X[i][j] = 0;
}
}
}
delete_matrix(m, N);
free(c);
free(d);
}
void print_matrix(long m, long n, long** a)
{
long i, j;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++)
printf("%2ld ", a[i][j]);
printf("\n");
}
}
int main(void)
{
long i, j, m = 8, n = 8, p = 13, r;
long a[8][8] = { {0, 0, 0, 0, 0, 0, 0, 0},
{2, 0, 7, 11, 10, 12, 5, 11},
{3, 6, 3, 3, 0, 4, 7, 2},
{4, 3, 6, 4, 1, 6, 2, 3},
{2, 11, 8, 8, 2, 1, 3, 11},
{6, 11, 8, 6, 2, 6, 10, 9},
{5, 11, 7, 10, 0, 11, 6, 12},
{3, 3, 12, 5, 0, 11, 9, 11} };
long** M = create_matrix(m, n);
long** X = create_matrix(n, n);
for (i = 0; i < m; i++)
for (j = 0; j < n; j++)
M[i][j] = a[j][i];
printf("the original matrix is as follows:\n");
print_matrix(m, n, M);
image(m, n, p, M, X, &r);
printf("the image of the matrix is as follows:\n");
print_matrix(n, n - r, X);
printf("the rank of the matrix is: %ld\n", n - r);
delete_matrix(m, M);
delete_matrix(n, X);
return 0;
}
First Blast from the Past (1997) Computing the Inverse Image of a Matrix by James Pate Williams, Jr.
/*
Author: Pate Williams (c) 1997
"Algorithm 2.3.5 (Inverse Image Matrix). Let M be
an m by n matrix and V be an m by r matrix, where
n <= m. This algorithm either outputs a message
saying that some column vector of V is not in the
image of M, or outputs an n by r matrix X such
that V = MX." -Henri Cohen- See "A Course in Com-
putational Algebraic Number Theory" by Henri
Cohen pages 60-61. We specialize to the field Q.
*/
#include <stdio.h>
#include <stdlib.h>
double** create_matrix(long m, long n)
{
double** matrix = (double**)calloc(m, sizeof(double*));
long i;
if (!matrix) {
fprintf(stderr, "fatal error\ninsufficient memory\n");
fprintf(stderr, "from create_matrix\n");
exit(1);
}
for (i = 0; i < m; i++) {
matrix[i] = (double*)calloc(n, sizeof(double));
if (!matrix[i]) {
fprintf(stderr, "fatal error\ninsufficient memory\n");
fprintf(stderr, "from create_matrix\n");
exit(1);
}
}
return matrix;
}
void delete_matrix(long m, double** matrix)
{
long i;
for (i = 0; i < m; i++) free(matrix[i]);
free(matrix);
}
void inverse_image_matrix(long m, long n, long r,
double** M, double** V,
double** X)
{
double ck, d, sum, t;
double** B = create_matrix(m, r);
int found;
long i, j, k, l;
for (i = 0; i < m; i++)
for (j = 0; j < r; j++)
B[i][j] = V[i][j];
for (j = 0; j < n; j++) {
found = 0, i = j;
while (!found && i < m) {
found = M[i][j] != 0;
if (!found) i++;
}
if (!found) {
fprintf(stderr, "fatal error\nnot linearly independent\n");
fprintf(stderr, "from inverse_image_matrix\n");
exit(1);
}
if (i > j) {
for (l = 0; l < n; l++)
t = M[i][l], M[i][l] = M[j][l], M[j][l] = t;
for (l = 0; l < r; l++)
t = B[i][l], B[i][l] = B[j][l], B[j][l] = t;
}
d = 1.0 / M[j][j];
for (k = j + 1; k < m; k++) {
ck = d * M[k][j];
for (l = j + 1; l < n; l++)
M[k][l] -= ck * M[j][l];
for (l = 0; l < r; l++)
B[k][l] -= ck * B[j][l];
}
}
for (i = n - 1; i >= 0; i--) {
for (k = 0; k < r; k++) {
sum = 0;
for (j = i + 1; j < n; j++)
sum += M[i][j] * X[j][k];
X[i][k] = (B[i][k] - sum) / M[i][i];
}
}
for (k = n + 1; k < m; k++) {
for (j = 0; j < r; j++) {
sum = 0;
for (i = 0; i < n; i++)
sum += M[k][i] * X[i][j];
if (sum != B[k][j]) {
fprintf(stderr, "fatal error\ncolumn not in image\n");
fprintf(stderr, "from inverse_image_matrix\n");
exit(1);
}
}
}
delete_matrix(m, B);
}
void matrix_multiply(long m, long n, long r,
double** a, double** b,
double** c)
/* c = a * b */
{
double sum;
long i, j, k;
for (i = 0; i < m; i++) {
for (j = 0; j < r; j++) {
sum = 0.0;
for (k = 0; k < n; k++)
sum += a[i][k] * b[k][j];
c[i][j] = sum;
}
}
}
void print_matrix(long m, long n, double** a)
{
long i, j;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++)
printf("%+10.6lf ", a[i][j]);
printf("\n");
}
}
int main(void)
{
long i, j, m = 4, n = 4, r = 4;
double** c = create_matrix(m, n);
double** M = create_matrix(m, n);
double** V = create_matrix(m, r);
double** X = create_matrix(n, r);
for (i = 0; i < m; i++) {
c[i][i] = M[i][i] = 2.0;
if (i > 0)
c[i][i - 1] = M[i][i - 1] = -1.0;
if (i < m - 1)
c[i][i + 1] = M[i][i + 1] = -1.0;
}
for (i = 0; i < m; i++)
for (j = 0; j < r; j++)
V[i][j] = i + j + 1;
printf("M is as follows:\n");
print_matrix(m, n, M);
printf("V is as follows:\n");
print_matrix(m, r, V);
inverse_image_matrix(m, n, r, M, V, X);
printf("X is as follows:\n");
print_matrix(n, r, X);
matrix_multiply(m, n, r, c, X, M);
printf("MX is as follows:\n");
print_matrix(m, r, M);
delete_matrix(m, c);
delete_matrix(m, M);
delete_matrix(m, V);
delete_matrix(n, X);
return 0;
}
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();
}
Part of a New Linear Algebra Package in C++ Implemented by James Pate Williams, Jr.
// 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];
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(size_t n, bool failure);
Vector<double> DblGaussianElimination(
bool& failure);
// The following three 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(
size_t n,
Vector<int>& pivot);
bool DblGaussianSolution(
int n,
Vector<double>& x,
Vector<int>& pivot);
bool DblSubstitution(
size_t n,
Vector<double>& x,
Vector<int>& pivot);
bool DblInverse(
size_t n,
Matrix<double>& A,
Vector<int>& pivot);
};
#include "pch.h"
#include "LinearAlgebra.h"
double LinearAlgebra::DblDeterminant(
size_t n, bool failure)
{
double deter = 1;
Vector<int> pivot(n);
if (!initialized || m != n)
{
failure = true;
return 0.0;
}
if (!DblGaussianFactor(n, 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(
size_t n,
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(
int n,
Vector<double>& x,
Vector<int>& pivot)
{
if (!DblGaussianFactor(n, pivot))
return false;
return DblSubstitution(n, x, pivot);
}
bool LinearAlgebra::DblSubstitution(
size_t n, 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(
size_t n,
Matrix<double>& A,
Vector<int>& pivot)
{
Vector<double> x(n);
if (!DblGaussianFactor(n, 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(n, 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;
}
/*
** 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())
{
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 rows must be >= 1" << endl;
return -101;
}
LinearAlgebra la(m, n);
Matrix<double> copyM(m, n);
Vector<double> copyB(n);
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;
}
}
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)
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(n, X, pivot))
exit(-103);
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(n, failure);
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];
}
}
if (!la.DblInverse(n, A, pivot))
{
cout << "Conte Gaussian inverse matrix failure!" << endl;
exit(-104);
}
else
A.OutputMatrix(outs, ' ', 5, 8);
}
inps.close();
outs.close();
}
2 2 1 1 1 2 7 11 2 2 1 1 1 3 7 11 2 2 6 3 4 8 5 6 2 2 5 3 10 4 8 6 3 3 2 1 -1 -3 -1 2 -2 1 2 8 -11 -3
3 4
3 4
1
2 -1
-1 1
5 2
5 2
2
1.5 -0.5
-0.5 0.5
0.61111 0.44444
0.61111 0.44444
36
0.22222 -0.11111
-0.083333 0.16667
-1.4 5
-1.4 5
-10
-0.4 1
0.3 -0.5
2 3 -1
2 3 -1
1
4 5 -2
3 4 -2
-1 -1 1
Win32 C Double Precision Matrix Operations Including Finding the Determinant and Inverse
Win32 C Real Double Precision and Integer Matrix Add and Subtract by James Pate Williams, Jr.
void DblAdd(int l, double* a, double* b, double* c)
{
_asm
{
mov eax, a; load a matrix address
mov ebx, b; load a matrix address
mov edi, c; load result address
mov edx, l; load matrix length
finit;
ll: fld qword ptr[eax];
fld qword ptr[ebx];
fadd;
fstp qword ptr[edi];
fwait;
add edi, 8;
add eax, 8;
add ebx, 8;
dec edx;
jnz ll;
}
}
void DblSub(int l, double* a, double* b, double* c)
{
_asm
{
mov eax, a; load a matrix address
mov ebx, b; load a matrix address
mov edi, c; load result address
mov edx, l; load matrix length
finit;
ll: fld qword ptr[eax];
fld qword ptr[ebx];
fsub;
fstp qword ptr[edi];
fwait;
add edi, 8;
add eax, 8;
add ebx, 8;
dec edx;
jnz ll;
}
}
C Matrix Operations: Add, Subtract, Multiply, and Determinant Calculation by James Pate Williams, Jr.
Latest Software Development Project Started on Saturday, October 15, 2022, 12:07:19 AM by James Pate Williams, Jr.
I am in the progress of translating (porting) my J. M. Pollard’s algorithm “Factoring with Cubic Integers” C# application to a Free LIP based vanilla C Windows 32-bit console application. The first phase of the method is to generate two factor bases namely a rational prime factor base and an algebraic integer prime factor base. I have included some preliminary results from this fast-moving computer programming task. I generated 2012 algebraic integer primes in about a minute and thirty seconds.
Numerical Explorations 