Category: Matrix Algebra
Vector Analysis by James Pate Williams, Jr. Exercises and Supplementary Problems from Introduction to Vector Analysis Fourth Edition© 1979 by Harry F. Davis and Arthur David Snider Selected Exercises from Chapter 1 Pages 48 and 50
// VectorAnalysis.cpp © Tuesday, April 14, 2026
// by James Pate Williams, Jr.
// Reference: Introduction to Vector Analysis Fourth Edition
// © 1979 by Harry F. Davis and Arthur David Snider
#include <iostream>
#include <vector>
static double InnerProduct(
std::vector<double> A,
std::vector<double> B,
int n)
{
double sum = 0.0;
for (int i = 0; i < n; i++)
sum += A[i] * B[i];
return sum;
}
static void VectorProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double>&C)
{
C.resize(3);
C[0] = A[1] * B[2] - A[2] * B[1];
C[1] = A[2] * B[0] - A[0] * B[2];
C[2] = A[0] * B[1] - A[1] * B[0];
}
static double TripleProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double> C)
{
double sum0 = 0.0, sum1 = 0.0;
sum0 += A[0] * B[1] * C[2] + A[1] * B[2] * C[0] + A[2] * B[0] * C[1];
sum1 += A[1] * B[0] * C[2] + A[0] * B[2] * C[1] + A[2] * B[1] * C[0];
return sum0 - sum1;
}
static void Exercises_Section_1_13_1_Triple_Products()
{
// triple products (a)
std::vector<double> A1 = { 2, 0, 0 };
std::vector<double> B1 = { 0, 3, 0 };
std::vector<double> C1 = { 0, 0, 5 };
double tp_1 = TripleProduct(A1, B1, C1);
std::cout << "Exercise 1 (a)" << '\t';
std::cout << tp_1 << std::endl;
// triple products (b)
std::vector<double> A2 = { 1, 1, 1 };
std::vector<double> B2 = { 3, 1, 0 };
std::vector<double> C2 = { 0, -1, 5 };
std::cout << "Exercise 1 (b)" << '\t';
double tp_2 = TripleProduct(A2, B2, C2);
std::cout << tp_2 << std::endl;
// triple products (c)
std::vector<double> A3 = { 2, -1, 1 };
std::vector<double> B3 = { 1, 1, 1 };
std::vector<double> C3 = { 2, 0, 3 };
double tp_3 = TripleProduct(A3, B3, C3);
std::cout << "Exercise 1 (c)" << '\t';
std::cout << tp_3 << std::endl;
// triple products (d)
std::vector<double> A4 = { 0, 0, 1 };
std::vector<double> B4 = { 1, 0, 0 };
std::vector<double> C4 = { 0, 1, 0 };
double tp_4 = TripleProduct(A4, B4, C4);
std::cout << "Exercise 1 (d)" << '\t';
std::cout << tp_4 << std::endl;
// volume of a parallelpipped
std::vector<double> A5 = { 3, 4, 1 };
std::vector<double> B5 = { 2, 3, 4 };
std::vector<double> C5 = { 0, 0, 5 };
double tp_5 = TripleProduct(A5, B5, C5);
std::cout << "Exercise 2" << '\t';
std::cout << tp_5 << std::endl;
// volume of a parallelpipped
std::vector<double> A6 = { 3, 2, 1 };
std::vector<double> B6 = { 4, 2, 1 };
std::vector<double> C6 = { 0, 1, 4 };
std::vector<double> D6 = { 0, 0, 7 };
std::vector<double> AB = { -1, 0, 0 };
std::vector<double> AC = { 3, 1, -3 };
std::vector<double> AD = { 3, 2, -6 };
std::cout << "Exercise 3" << '\t';
double tp_6 = TripleProduct(AB, AC, AD);
std::cout << tp_6 << std::endl;
// volume of a tetrahedron
std::vector<double> AB1 = { 1, 1, 0 };
std::vector<double> AC1 = { 1, -1, 0 };
std::vector<double> AD1 = { 0, 0, 2 };
std::cout << "Exercise 4" << '\t';
double tp_7 = TripleProduct(AB1, AC1, AD1) / 6.0;
std::cout << fabs(tp_7) << std::endl;
std::vector<double> P1 = { 0, 0, 0 };
std::vector<double> P2 = { 1, 1, 0 };
std::vector<double> P3 = { 3, 4, 0 };
std::vector<double> P4 = { 4, 5, 0 };
std::vector<double> P5 = { 0, 0, 1 };
std::vector<double> Q1 = { 1, 1, 0 };
std::vector<double> Q2 = { 2, 3, 0 };
std::vector<double> Q3 = { 1, 1, 1 };
std::cout << "Exercise 5" << '\t';
double tp_8 = TripleProduct(Q1, Q2, Q3);
std::cout << fabs(tp_8) << std::endl;
std::vector<double> A10 = { 1, 1, 1 };
std::vector<double> B10 = { 2, 4, -1 };
std::vector<double> C10 = { 1, 1, 3 };
double tp_10 = TripleProduct(A10, B10, C10);
std::vector<double> D10 = { 0 };
VectorProduct(A10, B10, D10);
double magnitude = sqrt(InnerProduct(D10, D10, 3));
std::cout << "Exercise 10" << '\t';
std::cout << tp_10 / magnitude << '\t';
std::cout << 2.0 * sqrt(38.0) / 19.0 << std::endl;
std::vector<double> A11 = { 1, 1, 1 };
std::vector<double> B11 = { 2, 4, -1 };
std::vector<double> C11 = { 1, 1, 3 };
std::vector<double> D11 = { 3, 2, 1 };
std::vector<double> AB11(3), BC11(3), CA11(3), BCCA11(3);
VectorProduct(A11, B11, AB11);
VectorProduct(B11, C11, BC11);
VectorProduct(C11, A11, CA11);
VectorProduct(BC11, CA11, BCCA11);
double Q11 = InnerProduct(AB11, BCCA11, 3);
double A_x = A11[0], A_y = A11[1], A_z = A11[2];
double B_x = B11[0], B_y = B11[1], B_z = B11[2];
double C_x = C11[0], C_y = C11[1], C_z = C11[2];
double term1 = +(A_y * B_z - A_z * B_y) * (B_z * C_x - B_x * C_z) * (C_x * A_y - C_y * A_x);
double term2 = -(A_y * B_z - A_z * B_y) * (B_x * C_y - B_y * C_x) * (C_z * A_x - C_x * A_z);
double term3 = +(A_z * B_x - A_x * B_z) * (B_x * C_y - B_y * C_x) * (C_y * A_z - C_z * A_y);
double term4 = -(A_z * B_x - A_x * B_z) * (B_y * C_z - B_z * C_y) * (C_x * A_y - C_y * A_x);
double term5 = +(A_x * B_y - A_y * B_x) * (B_y * C_z - B_z * C_y) * (C_z * A_x - C_x * A_z);
double term6 = -(A_x * B_y - A_y * B_x) * (B_z * C_x - B_x * C_z) * (C_y * A_z - C_z * A_y);
double P11 = term1 + term2 + term3 + term4 + term5 + term6;
std::cout << "Q = (A x B) . (B x C) x (C x A) = " << Q11 << std::endl;
std::cout << "P = (A x B) . (B x C) x (C x A) = " << P11 << std::endl;
}
static void Exercises_Section_1_14_Vector_Identities()
{
std::vector<double> A11 = { 1, 1, 1 };
std::vector<double> B11 = { 2, 4, -1 };
std::vector<double> C11 = { 1, 1, 3 };
std::vector<double> D11 = { 3, 2, 1 };
std::cout << "Section 1.14 page 50 Exercises Exercise 1" << std::endl;
std::cout << "A = " << A11[0] << '\t' << A11[1] << '\t' << A11[2] << std::endl;
std::cout << "B = " << B11[0] << '\t' << B11[1] << '\t' << B11[2] << std::endl;
std::cout << "C = " << C11[0] << '\t' << C11[1] << '\t' << C11[2] << std::endl;
std::cout << "D = " << D11[0] << '\t' << D11[1] << '\t' << D11[2] << std::endl;
std::cout << "TPI1 = (A x B) x (C x D) = [A, B, D]C - [A, B, C]D = " << std::endl;
double TP1411a = TripleProduct(A11, B11, D11);
double TP1411b = TripleProduct(A11, B11, C11);
std::cout << "[A, B, D] = " << TP1411a << std::endl;
std::cout << "[A, B, C] = " << TP1411b << std::endl;
std::cout << "TPI1_x = " << TP1411a * C11[0] << std::endl;
std::cout << "TPI1_y = " << TP1411a * C11[1] << std::endl;
std::cout << "TPI1_z = " << TP1411a * C11[2] << std::endl;
std::cout << "TPI2_x = " << TP1411b * D11[0] << std::endl;
std::cout << "TPI2_y = " << TP1411b * D11[1] << std::endl;
std::cout << "TPI2_z = " << TP1411b * D11[2] << std::endl;
std::cout << "RHS1 = [A, B, D]C - [A, B, C]D = " << std::endl;
std::vector<double> RHS1(3);
RHS1[0] = TP1411a * C11[0] - TP1411b * D11[0];
RHS1[1] = TP1411a * C11[1] - TP1411b * D11[1];
RHS1[2] = TP1411a * C11[2] - TP1411b * D11[2];
std::cout << "RHS1_x = " << RHS1[0] << std::endl;
std::cout << "RHS1_y = " << RHS1[1] << std::endl;
std::cout << "RHS1_z = " << RHS1[2] << std::endl;
std::vector<double> CD11(3), TD11(3);
std::vector<double> AB11(3), BC11(3), CA11(3), BCCA11(3);
VectorProduct(A11, B11, AB11);
VectorProduct(B11, C11, BC11);
VectorProduct(C11, A11, CA11);
VectorProduct(BC11, CA11, BCCA11);
VectorProduct(A11, B11, AB11);
VectorProduct(C11, D11, CD11);
VectorProduct(AB11, CD11, TD11);
std::cout << "A = " << A11[0] << '\t' << A11[1] << '\t' << A11[2] << std::endl;
std::cout << "B = " << B11[0] << '\t' << B11[1] << '\t' << B11[2] << std::endl;
std::cout << "C = " << C11[0] << '\t' << C11[1] << '\t' << C11[2] << std::endl;
std::cout << "D = " << D11[0] << '\t' << D11[1] << '\t' << D11[2] << std::endl;
std::cout << "A x B = " << AB11[0] << '\t' << AB11[1] << '\t' << AB11[2] << std::endl;
std::cout << "C x D = " << CD11[0] << '\t' << CD11[1] << '\t' << CD11[2] << std::endl;
std::cout << "TD11 = (A x B) x (C x D) = " << std::endl;
std::cout << "TD11_x = " << TD11[0] << std::endl;
std::cout << "TD11_y = " << TD11[1] << std::endl;
std::cout << "TD11_z = " << TD11[2] << std::endl;
VectorProduct(B11, C11, BC11);
VectorProduct(C11, A11, CA11);
VectorProduct(BC11, CA11, D11);
VectorProduct(A11, B11, AB11);
double ip12 = InnerProduct(AB11, D11, 3);
std::cout << "2. Inner Product = " << ip12 << std::endl;
double tp12 = TripleProduct(A11, B11, C11);
std::cout << "2. Triple Product ^ 2 = " << tp12 * tp12 << std::endl;
std::vector<double> ABC11(3), BAC11(3), CAB11(3);
VectorProduct(A11, BC11, ABC11);
VectorProduct(B11, CA11, BAC11);
VectorProduct(C11, AB11, CAB11);
double zx = ABC11[0] + BAC11[0] + CAB11[0];
double zy = ABC11[1] + BAC11[1] + CAB11[1];
double zz = ABC11[2] + BAC11[2] + CAB11[2];
std::cout << "3. Zero Vector = " << zx << ' ' << zy << ' ' << zz;
std::cout << std::endl;
}
int main()
{
Exercises_Section_1_13_1_Triple_Products();
Exercises_Section_1_14_Vector_Identities();
return 0;
}
Blog Entry © Tuesday, April 14, 2026, by James Pate Williams, Jr. Exercises and Supplementary Problems from Introduction to Vector Analysis Fourth Edition © 1979 by Harry F. Davis and Arthur David Snider Selected Exercises from Chapter 1 Page 48
// VectorAnalysis.cpp © Tuesday, April 14, 2026
// by James Pate Williams, Jr.
// Reference: Introduction to Vector Analysis Fourth Edition
// © 1979 by Harry F. Davis and Arthur David Snider
#include <iostream>
#include <vector>
static double InnerProduct(
std::vector<double> A,
std::vector<double> B,
int n)
{
double sum = 0.0;
for (int i = 0; i < n; i++)
sum += A[i] * B[i];
return sum;
}
static void VectorProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double>&C)
{
C.resize(3);
C[0] = A[1] * B[2] - A[2] * B[1];
C[1] = A[0] * B[2] - A[2] * B[0];
C[2] = A[0] * B[1] - A[1] * B[0];
}
static double TripleProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double> C)
{
double sum0 = 0.0, sum1 = 0.0;
sum0 += A[0] * B[1] * C[2] + A[1] * B[2] * C[0] + A[2] * B[0] * C[1];
sum1 += A[1] * B[0] * C[2] + A[0] * B[2] * C[1] + A[2] * B[1] * C[0];
return sum0 - sum1;
}
static void Exercises_Section_1_13_1_Triple_Products()
{
// triple products (a)
std::vector<double> A1 = { 2, 0, 0 };
std::vector<double> B1 = { 0, 3, 0 };
std::vector<double> C1 = { 0, 0, 5 };
double tp_1 = TripleProduct(A1, B1, C1);
std::cout << "Exercise 1 (a)" << '\t';
std::cout << tp_1 << std::endl;
// triple products (b)
std::vector<double> A2 = { 1, 1, 1 };
std::vector<double> B2 = { 3, 1, 0 };
std::vector<double> C2 = { 0, -1, 5 };
std::cout << "Exercise 1 (b)" << '\t';
double tp_2 = TripleProduct(A2, B2, C2);
std::cout << tp_2 << std::endl;
// triple products (c)
std::vector<double> A3 = { 2, -1, 1 };
std::vector<double> B3 = { 1, 1, 1 };
std::vector<double> C3 = { 2, 0, 3 };
double tp_3 = TripleProduct(A3, B3, C3);
std::cout << "Exercise 1 (c)" << '\t';
std::cout << tp_3 << std::endl;
// triple products (d)
std::vector<double> A4 = { 0, 0, 1 };
std::vector<double> B4 = { 1, 0, 0 };
std::vector<double> C4 = { 0, 1, 0 };
double tp_4 = TripleProduct(A4, B4, C4);
std::cout << "Exercise 1 (d)" << '\t';
std::cout << tp_4 << std::endl;
// volume of a parallelpipped
std::vector<double> A5 = { 3, 4, 1 };
std::vector<double> B5 = { 2, 3, 4 };
std::vector<double> C5 = { 0, 0, 5 };
double tp_5 = TripleProduct(A5, B5, C5);
std::cout << "Exercise 2" << '\t';
std::cout << tp_5 << std::endl;
// volume of a parallelpipped
std::vector<double> A6 = { 3, 2, 1 };
std::vector<double> B6 = { 4, 2, 1 };
std::vector<double> C6 = { 0, 1, 4 };
std::vector<double> D6 = { 0, 0, 7 };
std::vector<double> AB = { -1, 0, 0 };
std::vector<double> AC = { 3, 1, -3 };
std::vector<double> AD = { 3, 2, -6 };
std::cout << "Exercise 3" << '\t';
double tp_6 = TripleProduct(AB, AC, AD);
std::cout << tp_6 << std::endl;
// volume of a tetrahedron
std::vector<double> AB1 = { 1, 1, 0 };
std::vector<double> AC1 = { 1, -1, 0 };
std::vector<double> AD1 = { 0, 0, 2 };
std::cout << "Exercise 4" << '\t';
double tp_7 = TripleProduct(AB1, AC1, AD1) / 6.0;
std::cout << fabs(tp_7) << std::endl;
std::vector<double> P1 = { 0, 0, 0 };
std::vector<double> P2 = { 1, 1, 0 };
std::vector<double> P3 = { 3, 4, 0 };
std::vector<double> P4 = { 4, 5, 0 };
std::vector<double> P5 = { 0, 0, 1 };
std::vector<double> Q1 = { 1, 1, 0 };
std::vector<double> Q2 = { 2, 3, 0 };
std::vector<double> Q3 = { 1, 1, 1 };
std::cout << "Exercise 5" << '\t';
double tp_8 = TripleProduct(Q1, Q2, Q3);
std::cout << fabs(tp_8) << std::endl;
std::vector<double> A10 = { 1, 1, 1 };
std::vector<double> B10 = { 2, 4, -1 };
std::vector<double> C10 = { 1, 1, 3 };
double tp_10 = TripleProduct(A10, B10, C10);
std::vector<double> D10 = { 0 };
VectorProduct(A10, B10, D10);
double magnitude = sqrt(InnerProduct(D10, D10, 3));
std::cout << "Exercise 10" << '\t';
std::cout << tp_10 / magnitude << '\t';
std::cout << 2.0 * sqrt(38.0) / 19.0 << std::endl;
}
int main()
{
Exercises_Section_1_13_1_Triple_Products();
return 0;
}
Blog Entry © Thursday, April 2, 2026, by James Pate Williams, Jr., Matrix Inverses that Exist and Their Characteristic Polynomials
#pragma once
#include <vector>
class Inverse
{
public:
static bool Compute(
int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& X);
static void CharacteristicPolynomial(
int n, std::vector<double>& a,
std::vector<std::vector<double>>& C,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& Madj);
};
// Reference: 5.2: The Characteristic Polynomial - Mathematics LibreTexts
#include "Inverse.h"
#include <vector>
bool Inverse::Compute(
int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& X) {
double d = 0.0;
int i = 0, j = 0, k = 0, l = 0;
std::vector<double> C(n);
std::vector<std::vector<double>> B;
B.resize(n);
X.resize(n);
for (i = 0; i < n; i++)
{
B[i].resize(n);
X[i].resize(n);
for (j = 0; j < n; j++)
{
B[i][j] = 0.0;
X[i][j] = 0.0;
}
B[i][i] = 1.0;
}
j = -1;
Step2:
j++;
if (j == n)
goto Step6;
for (i = j; i < n; i++) {
if (M[i][j] != 0)
break;
}
if (i == n)
return false;
if (i > j) {
for (l = j; l < n; l++) {
double temp = M[i][l];
M[i][l] = M[j][l];
M[j][l] = temp;
}
for (l = 0; l < n; l++) {
double temp = B[i][l];
B[i][l] = B[j][l];
B[j][l] = temp;
}
}
d = 1.0 / M[j][j];
for (k = j + 1; k < n; k++) {
C[k] = d * M[k][j];
}
for (k = j + 1; k < n; k++) {
for (l = j + 1; l < n; l++) {
M[k][l] -= C[k] * M[j][l];
}
}
for (k = j + 1; k < n; k++) {
for (l = 0; l < n; l++) {
B[k][l] -= C[k] * B[j][l];
}
}
goto Step2;
Step6:
for (i = n - 1; i >= 0; i--) {
std::vector<double> sum(n, 0);
for (j = i + 1; j < n; j++) {
for (k = 0; k < n; k++) {
sum[k] += M[i][j] * X[j][k];
}
}
for (j = 0; j < n; j++)
X[i][j] = (B[i][j] - sum[j]) / M[i][i];
}
return true;
}
double delta(int i, int j) {
return i == j ? 1.0 : 0.0;
}
void Inverse::CharacteristicPolynomial(
int n, std::vector<double>& a,
std::vector<std::vector<double>>& C,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& Madj) {
int i = 0, j = 0, k = 0, l = 0;
std::vector<double> X(n);
std::vector<std::vector<double>> MC;
a.resize(n + 1LL);
C.resize(n);
M.resize(n);
MC.resize(n);
Madj.resize(n);
for (i = 0; i < n; i++) {
C[i].resize(n);
M[i].resize(n);
MC[i].resize(n);
Madj[i].resize(n);
for (j = 0; j < n; j++) {
C[i][j] = delta(i, j);
}
}
a[0] = 1.0;
i = 0;
Step2:
i++;
if (i == n) {
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
double sum = 0.0;
for (l = 0; l < n; l++) {
sum += M[j][l] * C[l][k];
}
MC[j][k] = sum;
}
}
double trace = 0.0;
for (j = 0; j < n; j++) {
trace += MC[j][j];
}
a[n] = -trace / n;
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
Madj[j][k] = pow(-1, n - 1) * C[j][k];
}
}
return;
}
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
double sum = 0.0;
for (l = 0; l < n; l++) {
sum += M[j][l] * C[l][k];
}
MC[j][k] = sum;
}
}
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
C[j][k] = MC[j][k];
}
}
double tr = 0.0;
for (j = 0; j < n; j++) {
tr += C[j][j];
}
a[i] = -tr / i;
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
C[j][k] += a[i] * delta(j, k);
}
}
goto Step2;
}
Blog Entry © Sunday, March 29, 2026, James Pate Williams, Jr. Properties of Determinants
// Reference: "A Course in Computational
// Algebraic Number Theory by Henri Cohen
// Chapter 2 Algorithm 2.2.1 and Algorithm
// 2.2.3
#pragma once
#include <vector>
class Determinants
{
public:
static bool GaussianElimnation(
int n,
std::vector<std::vector<double>>& M,
std::vector<double> B,
std::vector<double> X);
static bool Determinant(
double& det, int n,
std::vector<std::vector<double>>& M);
};
#include "Determinants.h"
bool Determinants::GaussianElimnation(
int n,
std::vector<std::vector<double>>& M,
std::vector<double> B,
std::vector<double> X) {
double d = 0.0;
int i = 0, j = -1, k = 0, l = 0;
std::vector<double> C(n);
Step2:
j++;
if (j == n)
goto Step6;
for (i = j; i < n; i++) {
if (M[i][j] > 0)
break;
}
if (i == n)
return false;
if (i > j) {
for (l = j; j < n; l++) {
double temp = M[i][l];
M[i][l] = M[j][l];
M[j][l] = temp;
temp = B[i];
B[i] = B[j];
B[j] = temp;
}
}
d = 1.0 / M[j][j];
for (k = j + 1; k < n; k++) {
C[k] = d * M[k][j];
}
for (k = j + 1; k < n; k++) {
for (l = j + 1; l < n; l++) {
M[k][l] -= C[k] * M[j][l];
}
}
for (k = j + 1; k < n; k++) {
B[k] -= C[k] * B[j];
}
goto Step2;
Step6:
for (i = n - 1; i >= 0; i--) {
double sum = 0.0;
for (j = i + 1; j < n; j++) {
sum += M[i][j] * X[j];
}
X[i] = (B[i] - sum) / M[i][i];
}
return true;
}
bool Determinants::Determinant(
double& det, int n,
std::vector<std::vector<double>>& M) {
double d = 0.0, x = 1.0;
int i = 0, j = -1, k = 0, l = 0;
std::vector<double> C(n, 0);
Step2:
j++;
if (j == n) {
det = x;
return true;
}
for (i = j; i < n; i++) {
if (M[i][j] != 0.0)
break;
}
if (i == n) {
det = 0.0;
return false;
}
if (i > j) {
for (l = j; j < n; l++) {
double temp = M[i][l];
M[i][l] = M[j][l];
M[j][l] = temp;
}
x = -x;
}
d = 1.0 / M[j][j];
for (k = j + 1; k < n; k++) {
C[k] = d * M[k][j];
}
for (k = j + 1; k < n; k++) {
for (l = j + 1; l < n; l++) {
M[k][l] -= C[k] * M[j][l];
}
}
x *= M[j][j];
goto Step2;
}
// ComputeRealDeterminant.cpp : Defines the entry point for the application.
//
#include "framework.h"
#include "ComputeRealDeterminant.h"
#include "Determinants.h"
#include <string>
#include <vector>
#define MAX_LOADSTRING 100
// Global Variables:
HINSTANCE hInst; // current instance
WCHAR szTitle[MAX_LOADSTRING]; // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING]; // the main window class name
WCHAR line[128], numWstr[128];
WCHAR inpStr[16384], outStr[16384];
// Forward declarations of functions included in this code module:
ATOM MyRegisterClass(HINSTANCE hInstance);
BOOL InitInstance(HINSTANCE, int);
LRESULT CALLBACK WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK About(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK DataInputDialog(HWND, UINT, WPARAM, LPARAM);
int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
_In_opt_ HINSTANCE hPrevInstance,
_In_ LPWSTR lpCmdLine,
_In_ int nCmdShow)
{
UNREFERENCED_PARAMETER(hPrevInstance);
UNREFERENCED_PARAMETER(lpCmdLine);
// TODO: Place code here.
// Initialize global strings
LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
LoadStringW(hInstance, IDC_COMPUTEREALDETERMINANT, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_COMPUTEREALDETERMINANT));
MSG msg;
// Main message loop:
while (GetMessage(&msg, nullptr, 0, 0))
{
if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
{
TranslateMessage(&msg);
DispatchMessage(&msg);
}
}
return (int) msg.wParam;
}
//
// FUNCTION: MyRegisterClass()
//
// PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
WNDCLASSEXW wcex = { 0 };
wcex.cbSize = sizeof(WNDCLASSEX);
wcex.style = CS_HREDRAW | CS_VREDRAW;
wcex.lpfnWndProc = WndProc;
wcex.cbClsExtra = 0;
wcex.cbWndExtra = 0;
wcex.hInstance = hInstance;
wcex.hIcon = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_COMPUTEREALDETERMINANT));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_COMPUTEREALDETERMINANT);
wcex.lpszClassName = szWindowClass;
wcex.hIconSm = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));
return RegisterClassExW(&wcex);
}
//
// FUNCTION: InitInstance(HINSTANCE, int)
//
// PURPOSE: Saves instance handle and creates main window
//
// COMMENTS:
//
// In this function, we save the instance handle in a global variable and
// create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
hInst = hInstance; // Store instance handle in our global variable
HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);
if (!hWnd)
{
return FALSE;
}
ShowWindow(hWnd, nCmdShow);
UpdateWindow(hWnd);
return TRUE;
}
//
// FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
// PURPOSE: Processes messages for the main window.
//
// WM_COMMAND - process the application menu
// WM_PAINT - Paint the main window
// WM_DESTROY - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
switch (message)
{
case WM_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDM_START:
DialogBox(hInst, MAKEINTRESOURCE(IDD_DATA_INPUT_DIALOG), hWnd,
DataInputDialog);
break;
case IDM_EXIT:
DestroyWindow(hWnd);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
}
break;
case WM_PAINT:
{
PAINTSTRUCT ps;
HDC hdc = BeginPaint(hWnd, &ps);
// TODO: Add any drawing code that uses hdc here...
EndPaint(hWnd, &ps);
}
break;
case WM_DESTROY:
PostQuitMessage(0);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
return 0;
}
// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
switch (message)
{
case WM_INITDIALOG:
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
INT_PTR CALLBACK DataInputDialog(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam) {
static int n = 0;
static std::wstring wstr;
static std::vector < std::vector<double>> A;
static HFONT hFont = nullptr;
static HWND hEditMultiline1 = nullptr;
static HWND hEditMultiline2 = nullptr;
UNREFERENCED_PARAMETER(lParam);
switch (message) {
case WM_INITDIALOG:
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"2");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"3");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"4");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"5");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"6");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"7");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"8");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_SETCURSEL, 0, 0);
hFont = CreateFont(
-MulDiv(10, GetDeviceCaps(GetDC(hDlg), LOGPIXELSY), 72),
0, 0, 0, FW_BOLD, FALSE, FALSE, FALSE,
DEFAULT_CHARSET, OUT_DEFAULT_PRECIS,
CLIP_DEFAULT_PRECIS, DEFAULT_QUALITY,
FIXED_PITCH | FF_MODERN,
TEXT("Courier New")
);
hEditMultiline1 = CreateWindowEx(
WS_EX_CLIENTEDGE, // Extended style for sunken border
TEXT("EDIT"), // Class name
TEXT(""), // Initial text (can be blank)
WS_CHILD | WS_VISIBLE | WS_VSCROLL |
ES_LEFT | ES_MULTILINE | ES_AUTOVSCROLL,
10, 100, 250, 250, // Position and size
hDlg, // Parent window handle
(HMENU)IDC_EDIT_MULTILINE1, // Unique control ID
hInst, // Application instance
NULL // Extra parameter
);
hEditMultiline2 = CreateWindowEx(
WS_EX_CLIENTEDGE, // Extended style for sunken border
TEXT("EDIT"), // Class name
TEXT(""), // Initial text (can be blank)
WS_CHILD | WS_VISIBLE | ES_LEFT | ES_READONLY,
120, 360, 120, 30, // Position and size
hDlg, // Parent window handle
(HMENU)IDC_EDIT_MULTILINE2, // Unique control ID
hInst, // Application instance
NULL // Extra parameter
);
break;
case WM_COMMAND:
if (LOWORD(wParam) == IDC_BUTTON_COMPUTE) {
GetDlgItemText(hDlg, IDC_COMBO_N, inpStr, 16384);
WCHAR* endptr = nullptr;
n = (int)wcstol(inpStr, &endptr, 10);
A.resize(n);
GetDlgItemText(hDlg, IDC_EDIT_MULTILINE1, inpStr, 16384);
int k = 0;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
A[i].resize(n);
int l = 0;
while (inpStr[k] != L' ') {
numWstr[l++] = inpStr[k++];
}
k++;
numWstr[l++] = L'\0';
A[i][j] = (double)wcstof(numWstr, &endptr);
}
}
double det = 0.0;
bool result = Determinants::Determinant(det, n, A);
swprintf_s(line, 128, L"%lf", det);
SetDlgItemText(hDlg, IDC_EDIT_MULTILINE2, line);
break;
}
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL) {
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
//{{NO_DEPENDENCIES}}
// Microsoft Visual C++ generated include file.
// Used by ComputeRealDeterminant.rc
#define IDS_APP_TITLE 103
#define IDR_MAINFRAME 128
#define IDD_COMPUTEREALDETERMINANT_DIALOG 102
#define IDD_ABOUTBOX 103
#define IDM_ABOUT 104
#define IDM_EXIT 105
#define IDI_COMPUTEREALDETERMINANT 107
#define IDI_SMALL 108
#define IDC_COMPUTEREALDETERMINANT 109
#define IDC_MYICON 2
#ifndef IDC_STATIC
#define IDC_STATIC -1
#endif
// Next default values for new objects
//
#ifdef APSTUDIO_INVOKED
#ifndef APSTUDIO_READONLY_SYMBOLS
#define _APS_NO_MFC 130
#define _APS_NEXT_RESOURCE_VALUE 129
#define _APS_NEXT_COMMAND_VALUE 32771
#define _APS_NEXT_CONTROL_VALUE 1000
#define _APS_NEXT_SYMED_VALUE 110
#endif
#endif
#define IDD_DATA_INPUT_DIALOG 1000
#define IDC_COMBO_N 1010
#define IDC_EDIT_MULTILINE1 1020
#define IDC_EDIT_MULTILINE2 1030
#define IDC_BUTTON_COMPUTE 1040
#define IDM_START 32771
//Microsoft Visual C++ generated resource script.
//
#include "resource.h"
#define APSTUDIO_READONLY_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//
#ifndef APSTUDIO_INVOKED
#include "targetver.h"
#endif
#define APSTUDIO_HIDDEN_SYMBOLS
#include "windows.h"
#undef APSTUDIO_HIDDEN_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
#undef APSTUDIO_READONLY_SYMBOLS
#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_ENU)
LANGUAGE 9, 1
/////////////////////////////////////////////////////////////////////////////
//
// Icon
//
// Icon with lowest ID value placed first to ensure application icon
// remains consistent on all systems.
IDI_COMPUTEREALDETERMINANT ICON "ComputeRealDeterminant.ico"
IDI_SMALL ICON "small.ico"
/////////////////////////////////////////////////////////////////////////////
//
// Menu
//
IDC_COMPUTEREALDETERMINANT MENU
BEGIN
POPUP "&Begin"
BEGIN
MENUITEM "&Start", IDM_START
MENUITEM SEPARATOR
MENUITEM "E&xit", IDM_EXIT
END
POPUP "&Help"
BEGIN
MENUITEM "&About ...", IDM_ABOUT
END
END
/////////////////////////////////////////////////////////////////////////////
//
// Accelerator
//
IDC_COMPUTEREALDETERMINANT ACCELERATORS
BEGIN
"?", IDM_ABOUT, ASCII, ALT
"/", IDM_ABOUT, ASCII, ALT
END
/////////////////////////////////////////////////////////////////////////////
//
// Dialog
//
IDD_ABOUTBOX DIALOGEX 0, 0, 170, 62
STYLE DS_SETFONT | DS_MODALFRAME | DS_FIXEDSYS | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "About ComputeRealDeterminant"
FONT 8, "MS Shell Dlg"
BEGIN
ICON IDI_COMPUTEREALDETERMINANT,IDC_STATIC,14,14,21,20
LTEXT "ComputeRealDeterminant, Version 1.0",IDC_STATIC,42,14,114,8,SS_NOPREFIX
LTEXT "Copyright (c) 2026",IDC_STATIC,42,26,114,8
DEFPUSHBUTTON "OK",IDOK,113,41,50,14,WS_GROUP
END
IDD_DATA_INPUT_DIALOG DIALOGEX 0, 0, 280, 260
STYLE DS_SETFONT | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "Data Input Dialog"
FONT 10, "Courier New", 700
BEGIN
LTEXT "n:", IDC_STATIC, 10, 10, 30, 10
COMBOBOX IDC_COMBO_N, 70, 10, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT "input matrix:", IDC_STATIC, 10, 30, 200, 10
LTEXT "determinant:", IDC_STATIC, 10, 185, 200, 10
PUSHBUTTON "Compute", IDC_BUTTON_COMPUTE, 10, 220, 50, 16
PUSHBUTTON "Cancel", IDCANCEL, 220, 220, 50, 16
END
/////////////////////////////////////////////////////////////////////////////
//
// DESIGNINFO
//
#ifdef APSTUDIO_INVOKED
GUIDELINES DESIGNINFO
BEGIN
IDD_ABOUTBOX, DIALOG
BEGIN
LEFTMARGIN, 7
RIGHTMARGIN, 163
TOPMARGIN, 7
BOTTOMMARGIN, 55
END
END
#endif // APSTUDIO_INVOKED
#ifdef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// TEXTINCLUDE
//
1 TEXTINCLUDE
BEGIN
"resource.h\0"
END
2 TEXTINCLUDE
BEGIN
"#ifndef APSTUDIO_INVOKED\r\n"
"#include ""targetver.h""\r\n"
"#endif\r\n"
"#define APSTUDIO_HIDDEN_SYMBOLS\r\n"
"#include ""windows.h""\r\n"
"#undef APSTUDIO_HIDDEN_SYMBOLS\r\n"
"\0"
END
3 TEXTINCLUDE
BEGIN
"\r\n"
"\0"
END
#endif // APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// String Table
//
STRINGTABLE
BEGIN
IDC_COMPUTEREALDETERMINANT "COMPUTEREALDETERMINANT"
IDS_APP_TITLE "ComputeRealDeterminant"
END
#endif
/////////////////////////////////////////////////////////////////////////////
#ifndef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//
/////////////////////////////////////////////////////////////////////////////
#endif // not APSTUDIO_INVOKED
Blog Entry © Wednesday, September 3, 2025, by James Pate Williams, Jr. Solution of Exercise 1.11 in Modern Quantum Chemistry an Introduction to Advanced Electronic Structure Theory © 1996 by Attila Szabo and Neil S. Ostlund
// EigenVV2x2.cpp (c) Wednesday, September 3, 2025
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Solution to Exercise 1.11 in "Quantum Chemistry
// an Introduction to Advanced Electronic Structure
// Theory (c) 1996 by Attila Szabo and Neil S. Ostlund
#include "pch.h"
#include "framework.h"
#include "EigenVV2x2.h"
#define MAX_LOADSTRING 100
// Global Variables:
HINSTANCE hInst; // current instance
WCHAR szTitle[MAX_LOADSTRING]; // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING]; // the main window class name
std::wstring text; // output wide string
// Forward declarations of functions included in this code module:
ATOM MyRegisterClass(HINSTANCE hInstance);
BOOL InitInstance(HINSTANCE, int);
LRESULT CALLBACK WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK About(HWND, UINT, WPARAM, LPARAM);
int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
_In_opt_ HINSTANCE hPrevInstance,
_In_ LPWSTR lpCmdLine,
_In_ int nCmdShow)
{
UNREFERENCED_PARAMETER(hPrevInstance);
UNREFERENCED_PARAMETER(lpCmdLine);
// TODO: Place code here.
// Initialize global strings
LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
LoadStringW(hInstance, IDC_EIGENVV2X2, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_EIGENVV2X2));
MSG msg;
// Main message loop:
while (GetMessage(&msg, nullptr, 0, 0))
{
if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
{
TranslateMessage(&msg);
DispatchMessage(&msg);
}
}
return (int) msg.wParam;
}
//
// FUNCTION: MyRegisterClass()
//
// PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
WNDCLASSEXW wcex;
wcex.cbSize = sizeof(WNDCLASSEX);
wcex.style = CS_HREDRAW | CS_VREDRAW;
wcex.lpfnWndProc = WndProc;
wcex.cbClsExtra = 0;
wcex.cbWndExtra = 0;
wcex.hInstance = hInstance;
wcex.hIcon = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_EIGENVV2X2));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_EIGENVV2X2);
wcex.lpszClassName = szWindowClass;
wcex.hIconSm = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));
return RegisterClassExW(&wcex);
}
//
// FUNCTION: InitInstance(HINSTANCE, int)
//
// PURPOSE: Saves instance handle and creates main window
//
// COMMENTS:
//
// In this function, we save the instance handle in a global variable and
// create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
hInst = hInstance; // Store instance handle in our global variable
HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);
if (!hWnd)
{
return FALSE;
}
ShowWindow(hWnd, nCmdShow);
UpdateWindow(hWnd);
return TRUE;
}
#define IDC_STATIC1 1000
#define IDC_STATIC2 1010
#define IDC_STATIC3 1020
#define IDC_STATIC4 1030
#define IDC_EDIT_A11 2000
#define IDC_EDIT_A12 2010
#define IDC_EDIT_A21 2020
#define IDC_EDIT_A22 2030
#define IDC_EDIT_MULTILINE 3000
#define IDC_BUTTON_COMPUTE 4000
#define IDC_BUTTON_CANCEL 4010
//
// FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
// PURPOSE: Processes messages for the main window.
//
// WM_COMMAND - process the application menu
// WM_PAINT - Paint the main window
// WM_DESTROY - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
static HFONT hFont = { };
static HWND hEditMultiline = { };
static HWND hEditA11 = { };
static HWND hEditA12 = { };
static HWND hEditA21 = { };
static HWND hEditA22 = { };
switch (message)
{
case WM_CREATE:
{
hFont = CreateFont(16, 0, 0, 0, FW_NORMAL, FALSE, FALSE, FALSE,
ANSI_CHARSET, OUT_DEFAULT_PRECIS, CLIP_DEFAULT_PRECIS,
DEFAULT_QUALITY, DEFAULT_PITCH | FF_SWISS, L"Courier New Bold");
CreateWindowEx(0, L"STATIC", L"A[1][1]:", WS_CHILD | WS_VISIBLE,
10, 10, 80, 20, hWnd, (HMENU)IDC_STATIC1, hInst, NULL);
CreateWindowEx(0, L"STATIC", L"A[1][2]:", WS_CHILD | WS_VISIBLE,
10, 40, 80, 20, hWnd, (HMENU)IDC_STATIC2, hInst, NULL);
CreateWindowEx(0, L"STATIC", L"A[2][1]:", WS_CHILD | WS_VISIBLE,
10, 70, 80, 20, hWnd, (HMENU)IDC_STATIC3, hInst, NULL);
CreateWindowEx(0, L"STATIC", L"A[2][2]:", WS_CHILD | WS_VISIBLE,
10, 100, 80, 20, hWnd, (HMENU)IDC_STATIC4, hInst, NULL);
hEditA11 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
100, 10, 200, 20, hWnd, (HMENU)IDC_EDIT_A11, hInst, NULL);
hEditA12 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
100, 40, 200, 20, hWnd, (HMENU)IDC_EDIT_A12, hInst, NULL);
hEditA21 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
100, 70, 200, 20, hWnd, (HMENU)IDC_EDIT_A21, hInst, NULL);
hEditA22 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
100, 100, 200, 20, hWnd, (HMENU)IDC_EDIT_A22, hInst, NULL);
hEditMultiline = CreateWindowEx(
WS_EX_CLIENTEDGE, // Extended style for sunken border
TEXT("EDIT"), // Class name
TEXT(""), // Initial text (can be blank)
WS_CHILD | WS_VISIBLE | WS_VSCROLL | ES_AUTOHSCROLL |
ES_LEFT | ES_MULTILINE | ES_AUTOVSCROLL | WS_HSCROLL | WS_VSCROLL,
310, 10, 300, 300, // Position and size
hWnd, // Parent window handle
(HMENU)IDC_EDIT_MULTILINE, // Unique control ID
hInst, // Application instance
NULL // Extra parameter
);
CreateWindowEx(0, L"BUTTON", L"Compute", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
10, 130, 80, 30, hWnd, (HMENU)IDC_BUTTON_COMPUTE, hInst, NULL);
CreateWindowEx(0, L"BUTTON", L"Cancel", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
220, 130, 80, 30, hWnd, (HMENU)IDC_BUTTON_CANCEL, hInst, NULL);
SendMessage(hEditMultiline, WM_SETFONT, (WPARAM)hFont, TRUE);
SetDlgItemText(hWnd, IDC_EDIT_A11, L"3");
SetDlgItemText(hWnd, IDC_EDIT_A12, L"1");
SetDlgItemText(hWnd, IDC_EDIT_A21, L"1");
SetDlgItemText(hWnd, IDC_EDIT_A22, L"3");
}
break;
case WM_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDC_BUTTON_COMPUTE:
{
WCHAR line[128] = L"";
text = L"";
// eigenvalues
std::vector<double> omega(3);
// matrix
std::vector<std::vector<double>> A(3,
std::vector<double>(3));
std::vector<std::vector<double>> B(3,
std::vector<double>(3));
// eigenvectors
std::vector<std::vector<double>> c(3,
std::vector<double>(3));
GetWindowText(hEditA11, line, 128);
std::wstring A11Str(line);
A[1][1] = std::stod(A11Str);
GetWindowText(hEditA12, line, 128);
std::wstring A12Str(line);
A[1][2] = std::stod(A12Str);
GetWindowText(hEditA21, line, 128);
std::wstring A21Str(line);
A[2][1] = std::stod(A21Str);
GetWindowText(hEditA22, line, 128);
std::wstring A22Str(line);
A[2][2] = std::stod(A22Str);
double term1 = A[1][1] + A[2][2];
double term2 = pow(A[2][2] - A[1][1], 2.0);
double term3 = 4.0 * A[1][2] * A[2][1];
// compute eigenvalues
omega[1] = 0.5 * (term1 - sqrt(term2 + term3));
omega[2] = 0.5 * (term1 + sqrt(term2 + term3));
swprintf_s(line, L"Eigenvalues:\r\n\r\n");
text += std::wstring(line);
swprintf_s(line, L"omega 1 = %13.10lf\r\n", omega[1]);
text += std::wstring(line);
swprintf_s(line, L"omega 2 = %13.10lf\r\n\r\n", omega[2]);
text += std::wstring(line);
// compute eigenvalues using a unitary transformation
// matrix A must be symmetric
double theta0 = 0.0;
if (A[1][1] == A[2][2])
{
theta0 = 0.5 * acos(0.0);
}
else
{
theta0 = 0.5 * atan(2.0 * A[1][2] /
(A[1][1] - A[2][2]));
}
// compute the eigenvalues
omega[1] = A[1][1] * pow(cos(theta0), 2.0) +
A[2][2] * pow(sin(theta0), 2.0) +
A[1][2] * sin(2.0 * theta0);
omega[2] = A[1][1] * pow(cos(theta0), 2.0) +
A[2][2] * pow(sin(theta0), 2.0) -
A[1][2] * sin(2.0 * theta0);
swprintf_s(line, L"Eigenvalues:\r\n\r\n");
text += std::wstring(line);
swprintf_s(line, L"omega 1 = %13.10lf\r\n", omega[1]);
text += std::wstring(line);
swprintf_s(line, L"omega 2 = %13.10lf\r\n\r\n", omega[2]);
text += std::wstring(line);
// compute eigenvectors
c[1][1] = cos(theta0);
c[1][2] = sin(theta0);
c[2][1] = sin(theta0);
c[2][2] = -cos(theta0);
swprintf_s(line, L"Eigenvectors:\r\n\r\n");
text += std::wstring(line);
swprintf_s(line, L"c 11 = %13.10lf\r\n", c[1][1]);
text += std::wstring(line);
swprintf_s(line, L"c 12 = %13.10lf\r\n", c[1][2]);
text += std::wstring(line);
swprintf_s(line, L"c 21 = %13.10lf\r\n", c[2][1]);
text += std::wstring(line);
swprintf_s(line, L"c 22 = %13.10lf\r\n", c[2][2]);
text += std::wstring(line);
SetWindowText(hEditMultiline, text.c_str());
break;
}
case IDC_BUTTON_CANCEL:
{
PostQuitMessage(0);
break;
}
break;
case IDM_EXIT:
DestroyWindow(hWnd);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
}
break;
case WM_PAINT:
{
PAINTSTRUCT ps;
HDC hdc = BeginPaint(hWnd, &ps);
// TODO: Add any drawing code that uses hdc here...
EndPaint(hWnd, &ps);
}
break;
case WM_DESTROY:
PostQuitMessage(0);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
return 0;
}
// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
switch (message)
{
case WM_INITDIALOG:
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
Blog Entry © Tuesday, September 2, 2025, by James Pate Williams, Jr. Testing of a Backpropagation Neural Network Function Approximator
Blog Entry © Friday, August 29, 2025, by James Pate Williams, Jr., Eigenvalue and Eigenvector Calculators and Linear System of Equations Solver
Blog Entry © Friday, August 22, 2025, by James Pate Williams, Jr. New Quantum Chemical Total Molecular Ground-State Energies for the Helium Hydride Cation (a Hetero Nuclear molecule) and the Hydrogen Molecule (a Homo Nuclear Molecule)
Blog Entry © Tuesday, August 19, 2025, by James Pate Williams, Jr., Continuation of Answers to the Exercises in Chapter 1 of Modern Quantum Chemistry an Introduction to Advanced Electronic Structure Theory by Attila Szabo and Neil S. Ostlund
Note: Later on, Tuesday, August 19, 2025, I added five C++ source code files.
#include <vector>
#include <random>
class DblLinearAlgebra
{
public:
static void DblPrintMatrix(
int m, int n, std::vector<std::vector<double>>& A);
static void DblAddition(
size_t m, size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblSubtraction(
size_t m, size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblMultiply(
size_t m, size_t n, size_t p,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblAnticommutator(
size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblCommutator(
size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static double DblDeterminant(
int n, int row, int col,
std::vector<std::vector<double>>& A);
static bool DblGaussianElimination(
int m, int n, std::vector<std::vector<double>>& A,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot);
static bool DblGaussianFactor(
int n, std::vector<std::vector<double>>& M,
std::vector<size_t>& pivot);
static bool DblGaussianSolution(
int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot);
static bool DblSubstitution(
int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot);
static bool DblInverse(
int n, std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& A);
static void DblCharPolyAndAdjoint(
int n,
std::vector<std::vector<double>>& C,
std::vector<std::vector<double>>& I,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& adjoint,
std::vector<double>& a);
static void DblMatrixKernel(
int m, int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& X,
size_t& r);
static void DblMatrixImage(
int m, int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& N,
std::vector<std::vector<double>>& X,
int rank);
static void DblGenerateNonSingular(
double scale, double& determinant,
int n, unsigned int seed,
std::vector<std::vector<double>>& Mr);
};
#include "DblLinearAlgebra.h"
#include <iomanip>
#include <iostream>
void DblLinearAlgebra::DblPrintMatrix(
int m, int n, std::vector<std::vector<double>>& A)
{
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
std::cout << std::setprecision(6) << std::setw(9);
if (fabs(A[i][j]) > 1.0e-12)
{
std::cout << A[i][j] << ' ';
}
else
{
std::cout << 0 << ' ';
}
}
std::cout << std::endl;
}
}
void DblLinearAlgebra::DblAddition(
size_t m, size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
C[i][j] = A[i][j] + B[i][j];
}
}
}
void DblLinearAlgebra::DblSubtraction(
size_t m, size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
C[i][j] = A[i][j] - B[i][j];
}
}
}
void DblLinearAlgebra::DblMultiply(
size_t m, size_t n, size_t p,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
double sum = { 0 };
for (size_t k = 0; k < p; k++)
{
sum += A[i][k] * B[k][j];
}
C[i][j] = sum;
}
}
}
void DblLinearAlgebra::DblAnticommutator(
size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C)
{
std::vector<std::vector<double>> D(n,
std::vector<double>(n));
std::vector<std::vector<double>> E(n,
std::vector<double>(n));
DblMultiply(n, n, n, A, B, D);
DblMultiply(n, n, n, B, A, E);
DblAddition(n, n, D, E, C);
}
void DblLinearAlgebra::DblCommutator(
size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C)
{
std::vector<std::vector<double>> D(n,
std::vector<double>(n));
std::vector<std::vector<double>> E(n,
std::vector<double>(n));
DblMultiply(n, n, n, A, B, D);
DblMultiply(n, n, n, B, A, E);
DblSubtraction(n, n, D, E, C);
}
// https://www.geeksforgeeks.org/dsa/determinant-of-a-matrix/
double getDet(std::vector<std::vector<double>>& mat, int n) {
// Base case: if the matrix is 1x1
if (n == 1) {
return mat[0][0];
}
// Base case for 2x2 matrix
if (n == 2) {
return mat[0][0] * mat[1][1] -
mat[0][1] * mat[1][0];
}
// Recursive case for larger matrices
double res = 0;
for (int col = 0; col < n; ++col) {
// Create a submatrix by removing the first
// row and the current column
std::vector<std::vector<double>> sub(n - 1,
std::vector<double>(n - 1));
for (int i = 1; i < n; ++i) {
int subcol = 0;
for (int j = 0; j < n; ++j) {
// Skip the current column
if (j == col) continue;
// Fill the submatrix
sub[i - 1LL][subcol++] = mat[i][j];
}
}
// Cofactor expansion
int sign = (col % 2 == 0) ? 1 : -1;
res += sign * mat[0][col] * getDet(sub, n - 1);
}
return res;
}
double DblLinearAlgebra::DblDeterminant(
int n, int row, int col,
std::vector<std::vector<double>>& A)
{
return getDet(A, A.size());
}
bool DblLinearAlgebra::DblGaussianElimination(
int m, int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot)
{
bool failure = false;
std::vector<double> c(m);
b.resize(n);
x.resize(n);
for (size_t i = 0; i < m; i++)
c[i] = -1;
for (size_t j = 0; j < n; j++)
{
bool found = false;
size_t i = j;
while (i < n && !found)
{
if (M[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[i][l];
M[i][l] = M[j][l];
M[j][l] = t;
t = b[i];
b[i] = b[j];
b[j] = t;
}
}
double d = 1.0 / M[j][j];
for (size_t k = j + 1; k < n; k++)
c[k] = d * M[k][j];
for (size_t k = j + 1; k < n; k++)
{
for (size_t l = j + 1; l < n; l++)
M[k][l] = M[k][l] - c[k] * M[j][l];
b[k] = b[k] - c[k] * b[j];
}
}
for (long long i = (long long)n - 1; i >= 0; i--)
{
double sum = 0;
for (size_t j = i + 1; j < n; j++)
sum += M[i][j] * x[j];
x[i] = (b[i] - sum) / M[i][i];
}
return failure;
}
bool DblLinearAlgebra::DblGaussianFactor(
int n, std::vector<std::vector<double>>& M,
std::vector<size_t>& pivot)
{
// returns false if matrix is singular
std::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[i] = i;
row_max = 0;
for (size_t j = 0; j < n; j++)
row_max = fmax(row_max, fabs(M[i][j]));
if (row_max == 0)
{
flag = 0;
row_max = 1.0;
}
d[i] = row_max;
}
if (n <= 1) return flag != 0;
// factorization
for (size_t k = 0; k < (size_t)(n - 1LL); k++)
{
// determine pivot row the row i_star
col_max = fabs(M[k][k]) / d[k];
i_star = k;
for (size_t i = k + 1; i < n; i++)
{
awikod = fabs(M[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;
itemp = pivot[i_star];
pivot[i_star] = pivot[k];
pivot[k] = itemp;
temp = d[i_star];
d[i_star] = d[k];
d[k] = temp;
for (size_t j = 0; j < n; j++)
{
temp = M[i_star][j];
M[i_star][j] = M[k][j];
M[k][j] = temp;
}
}
// eliminate x[k]
for (size_t i = k + 1; i < n; i++)
{
M[i][k] /= M[k][k];
ratio = M[i][k];
for (size_t j = k + 1; j < n; j++)
M[i][j] -= ratio * M[k][j];
}
}
if (M[n - 1LL][n - 1LL] == 0) flag = 0;
}
if (flag == 0)
return false;
return true;
}
bool DblLinearAlgebra::DblGaussianSolution(
int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot)
{
if (!DblGaussianFactor(n, M, pivot))
return false;
return DblSubstitution(n, M, b, x, pivot);
}
bool DblLinearAlgebra::DblSubstitution(
int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot)
{
double sum = 0.0;
size_t n1 = n - 1LL;
if (n == 1)
{
x[0] = b[0] / M[0][0];
return true;
}
// forward substitution
x[0] = b[pivot[0]];
for (size_t i = 1; i < n; i++)
{
double sum = 0.0;
for (size_t j = 0; j < i; j++)
sum += M[i][j] * x[j];
x[i] = b[pivot[i]] - sum;
}
// backward substitution
x[n1] /= M[n1][n1];
for (long long i = n - 2LL; i >= 0; i--)
{
double sum = 0.0;
for (size_t j = i + 1; j < n; j++)
sum += M[i][j] * x[j];
x[i] = (x[i] - sum) / M[i][i];
}
return true;
}
bool DblLinearAlgebra::DblInverse(
int n, std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& Mi)
{
std::vector<double> b(n);
std::vector<double> x(n);
std::vector<size_t> pivot(n);
std::vector<std::vector<double>> Mc(n,
std::vector<double>(n));
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
{
Mc[i][j] = M[i][j];
}
}
if (!DblGaussianFactor(n, Mc, pivot))
return false;
for (size_t i = 0; i < n; i++)
{
b[i] = 0;
}
for (size_t i = 0; i < n; i++)
{
b[i] = 1;
if (!DblSubstitution(n, Mc, b, x, pivot))
return false;
b[i] = 0;
for (size_t j = 0; j < n; j++)
Mi[j][i] = x[j];
}
return true;
}
void DblLinearAlgebra::DblCharPolyAndAdjoint(
int n,
std::vector<std::vector<double>>& C,
std::vector<std::vector<double>>& I,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& adjoint,
std::vector<double>& a)
{
C.resize(n, std::vector<double>(n));
I.resize(n, std::vector<double>(n));
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
C[i][j] = I[i][j] = 0;
}
for (size_t i = 0; i < n; i++)
C[i][i] = I[i][i] = 1;
a[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[j][l] * C[l][k];
C[j][k] = sum;
}
}
double tr = 0.0;
for (size_t j = 0; j < n; j++)
tr += C[j][j];
a[i] = -tr / i;
for (size_t j = 0; j < n; j++)
{
for (size_t k = 0; k < n; k++)
C[j][k] += a[i] * I[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[i][k] * C[k][j];
C[i][j] = sum;
}
}
double trace = 0.0;
for (size_t i = 0; i < n; i++)
trace += C[i][i];
trace /= n;
a[n - 1LL] = -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[i][j] = factor * C[i][j];
}
}
void DblLinearAlgebra::DblMatrixKernel(
int m, int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& X,
size_t& r)
{
double D = 0.0;
std::vector <int> c(m);
std::vector <int> d(n);
r = 0;
for (size_t i = 0; i < m; i++)
c[i] = -1;
size_t j, k = 1;
Step2:
for (j = 0; j < m; j++)
{
if (M[j][k] != 0 && c[j] == -1)
break;
}
if (j == m)
{
r++;
d[k] = 0;
goto Step4;
}
D = -1.0 / M[j][k];
M[j][k] = -1;
for (size_t s = k + 1; s < n; s++)
{
M[j][s] = D * M[j][s];
for (size_t i = 0; i < m; i++)
{
if (i != j)
{
D = M[i][k];
M[i][k] = 0;
}
}
}
for (size_t s = k + 1; s < n; s++)
{
for (size_t i = 0; i < m; i++)
{
M[i][s] += D * M[j][s];
}
}
c[j] = (int)k;
d[k] = (int)j;
Step4:
if (k < n - 1)
{
k++;
goto Step2;
}
X.resize(n, std::vector<double>(n));
if (r != 0)
{
for (k = 0; k < n; k++)
{
if (d[k] == 0)
{
for (size_t i = 0; i < n; i++)
{
if (d[i] > 0)
X[k][i] = M[d[i]][k];
else if (i == k)
X[k][i] = 1;
else
X[k][i] = 0;
}
}
}
}
}
void DblLinearAlgebra::DblMatrixImage(
int m, int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& N,
std::vector<std::vector<double>>& X,
int rank)
{
double D = 0.0;
size_t r = 0;
std::vector<std::vector<double>> copyM(
m, std::vector<double>(n));
std::vector <int> c(m);
std::vector <int> d(n);
for (size_t i = 0; i < m; i++)
c[i] = -1;
size_t j = 0, k = 1;
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
N[i][j] = copyM[i][j] = M[i][j];
}
}
Step2:
for (size_t j = 0; j < m; j++)
{
if (copyM[j][k] != 0 && c[j] == -1)
break;
}
if (j == m)
{
r++;
d[k] = 0;
goto Step4;
}
D = -1.0 / copyM[j][k];
copyM[j][k] = -1;
for (size_t s = k + 1; s < n; s++)
{
copyM[j][s] = D * copyM[j][s];
for (size_t i = 0; i < m; i++)
{
if (i != j)
{
D = copyM[i][k];
copyM[i][k] = 0;
}
}
}
for (size_t s = k + 1; s < n; s++)
{
for (size_t i = 0; i < m; i++)
{
copyM[i][s] += D * copyM[j][s];
}
}
c[j] = (int)k;
d[k] = (int)j;
Step4:
if (k < (size_t)(n - 1LL))
{
k++;
goto Step2;
}
rank = (int)(n - r) ;
for (j = 0; j < m; j++)
{
if (c[j] != 0)
{
for (size_t i = 0; i < m; i++)
{
N[i][c[j]] = M[i][c[j]];
}
}
}
}
void DblLinearAlgebra::DblGenerateNonSingular(
double scale, double& determinant,
int n, unsigned int seed,
std::vector<std::vector<double>>& Mr)
{
bool failure = false;
std::mt19937 rng(seed);
std::uniform_real_distribution<double> dist(0.0, 1.0);
while (true)
{
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
Mr[i][j] = scale * dist(rng);
}
}
determinant = DblDeterminant(n, 0, 0, Mr);
failure = determinant == 0;
if (!failure)
return;
}
}
#include <complex>
#include <vector>
class CmpLinearAlgebra
{
public:
static void CmpPrintMatrix(
int m, int n,
std::vector<std::vector<std::complex<double>>>& Ac);
static void CmpAddition(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpSubtraction(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpMultiply(
size_t m, size_t n, size_t p,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpAnticommutator(
size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpCommutator(
size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static std::complex<double> CmpDeterminant(
int n,
std::vector<std::vector<std::complex<double>>>& Ac);
static void CmpAdjoint(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& Ac,
std::vector<std::vector<std::complex<double>>>& Ad);
static bool CmpGaussianElimination(
int m, int n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot);
static bool CmpGaussianFactor(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<size_t>& pivot);
static bool CmpGaussianSolution(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot);
static bool CmpSubstitution(
int m, int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot);
static bool CmpInverse(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& Mi);
static void CmpCharPolyAndAdjoint(
int n,
std::vector<std::vector<std::complex<double>>>& C,
std::vector<std::vector<std::complex<double>>>& I,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& adjoint,
std::vector<std::complex<double>>& a);
static void CmpMatrixKernel(
int m, int n,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& X,
size_t& r);
static void CmpMatrixImage(
int m, int n,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& N,
std::vector<std::vector<std::complex<double>>>& X,
int rank);
static void CmpGenerateNonSingular(
double scale, std::complex<double>& determinant,
int n, unsigned int seed,
std::vector<std::vector<std::complex<double>>>& Mc);
};
#include "CmpLinearAlgebra.h"
#include <iomanip>
#include <iostream>
#include <random>
void CmpLinearAlgebra::CmpPrintMatrix(
int m, int n,
std::vector<std::vector<std::complex<double>>>& Ac)
{
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
if (Ac[i][j]._Val[0] >= 0)
std::cout << '+';
else
std::cout << '-';
std::cout << std::setprecision(6) << std::setw(9);
if (fabs(Ac[i][j]._Val[0]) > 1.0e-12)
{
std::cout << fabs(Ac[i][j]._Val[0]) << ' ';
}
else
{
std::cout << 0 << ' ';
}
if (Ac[i][j]._Val[1] >= 0)
std::cout << '+';
else
std::cout << '-';
std::cout << std::setprecision(6) << std::setw(9);
if (fabs(Ac[i][j]._Val[1]) > 1.0e-12)
{
std::cout << fabs(Ac[i][j]._Val[1]) << "i\t";
}
else
{
std::cout << 0 << "i\t";
}
}
std::cout << std::endl;
}
}
void CmpLinearAlgebra::CmpAddition(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
C[i][j] = A[i][j] + B[i][j];
}
}
}
void CmpLinearAlgebra::CmpSubtraction(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
C[i][j] = A[i][j] - B[i][j];
}
}
}
void CmpLinearAlgebra::CmpMultiply(
size_t m, size_t n, size_t p,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
std::complex<double> sum = 0;
for (size_t k = 0; k < p; k++)
{
sum += A[i][k] * B[k][j];
}
C[i][j] = sum;
}
}
}
void CmpLinearAlgebra::CmpAnticommutator(
size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C)
{
std::vector<std::vector<std::complex<double>>> D(n,
std::vector<std::complex<double>>(n));
std::vector<std::vector<std::complex<double>>> E(n,
std::vector<std::complex<double>>(n));
CmpMultiply(n, n, n, A, B, D);
CmpMultiply(n, n, n, B, A, E);
CmpAddition(n, n, D, E, C);
}
void CmpLinearAlgebra::CmpCommutator(
size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C)
{
std::vector<std::vector<std::complex<double>>> D(n,
std::vector<std::complex<double>>(n));
std::vector<std::vector<std::complex<double>>> E(n,
std::vector<std::complex<double>>(n));
CmpMultiply(n, n, n, A, B, D);
CmpMultiply(n, n, n, B, A, E);
CmpSubtraction(n, n, D, E, C);
}
void CmpLinearAlgebra::CmpAdjoint(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& Ac,
std::vector<std::vector<std::complex<double>>>& Ad)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
Ad[j][i] = std::conj(Ac[i][j]);
}
}
}
// https://www.geeksforgeeks.org/dsa/determinant-of-a-matrix/
std::complex<double> getDet(
std::vector<std::vector<std::complex<double>>>& mat, int n) {
// Base case: if the matrix is 1x1
if (n == 1) {
return mat[0][0];
}
// Base case for 2x2 matrix
if (n == 2) {
return mat[0][0] * mat[1][1] -
mat[0][1] * mat[1][0];
}
// Recursive case for larger matrices
std::complex<double> res = 0;
for (int col = 0; col < n; ++col) {
// Create a submatrix by removing the first
// row and the current column
std::vector<std::vector<std::complex<double>>> sub(n - 1,
std::vector<std::complex<double>>(n - 1));
for (int i = 1; i < n; ++i) {
int subcol = 0;
for (int j = 0; j < n; ++j) {
// Skip the current column
if (j == col) continue;
// Fill the submatrix
sub[i - 1LL][subcol++] = mat[i][j];
}
}
// Cofactor expansion
int sign = (col % 2 == 0) ? 1 : -1;
std::complex<double> csign(sign, 0.0);
res = res + csign * mat[0][col] * getDet(sub, n - 1);
}
return res;
}
std::complex<double> CmpLinearAlgebra::CmpDeterminant(
int n, std::vector<std::vector<std::complex<double>>>& A)
{
return getDet(A, A.size());
}
bool CmpLinearAlgebra::CmpGaussianElimination(
int m, int n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot)
{
bool failure = false;
std::vector<std::complex<double>> c(m);
b.resize(n);
x.resize(n);
for (size_t i = 0; i < m; i++)
c[i] = -1;
for (size_t j = 0; j < n; j++)
{
bool found = false;
size_t i = j;
while (i < n && !found)
{
if (abs(A[i][j]) != 0)
found = true;
else
i++;
}
if (!found)
{
failure = true;
break;
}
if (i > j)
{
for (size_t l = j; l < n; l++)
{
std::complex<double> t = A[i][l];
A[i][l] = A[j][l];
A[j][l] = t;
t = b[i];
b[i] = b[j];
b[j] = t;
}
}
std::complex<double> d = 1.0 / A[j][j];
for (size_t k = j + 1; k < n; k++)
c[k] = d * A[k][j];
for (size_t k = j + 1; k < n; k++)
{
for (size_t l = j + 1; l < n; l++)
A[k][l] = A[k][l] - c[k] * A[j][l];
b[k] = b[k] - c[k] * b[j];
}
}
for (long long i = (long long)n - 1; i >= 0; i--)
{
std::complex<double> sum = 0;
for (size_t j = i + 1; j < n; j++)
sum += A[i][j] * x[j];
x[i] = (b[i] - sum) / A[i][i];
}
return failure;
}
bool CmpLinearAlgebra::CmpSubstitution(
int m, int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot)
{
std::complex<double> sum = 0;
size_t n1 = n - 1LL;
if (n == 1)
{
x[0] = b[0] / M[0][0];
return true;
}
// forward substitution
x[0] = b[pivot[0]];
for (size_t i = 1; i < n; i++)
{
std::complex<double> sum = 0;
for (size_t j = 0; j < i; j++)
sum += M[i][j] * x[j];
x[i] = b[pivot[i]] - sum;
}
// backward substitution
x[n1] /= M[n1][n1];
for (long long i = n - 2LL; i >= 0; i--)
{
std::complex<double> sum = 0;
for (size_t j = i + 1; j < n; j++)
sum += M[i][j] * x[j];
x[i] = (x[i] - sum) / M[i][i];
}
return true;
}
static std::complex<double> complex_max(
std::complex<double> a, std::complex<double> b) {
return (std::abs(a) > std::abs(b)) ? a : b;
}
bool CmpLinearAlgebra::CmpGaussianFactor(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<size_t>& pivot)
{
// returns false if matrix is singular
std::vector<std::complex<double>> d(n);
std::complex<double> awikod = 0, col_max = 0, ratio = 0, row_max = 0, temp = 0;
int flag = 1;
size_t i_star, itemp;
for (size_t i = 0; i < n; i++)
{
pivot[i] = i;
row_max = 0;
for (size_t j = 0; j < n; j++)
row_max = complex_max(row_max, abs(M[i][j]));
if (abs(row_max) == 0)
{
flag = 0;
row_max = 1;
}
d[i] = row_max;
}
if (n <= 1) return flag != 0;
// factorization
for (size_t k = 0; k < (size_t)n - 1LL; k++)
{
// determine pivot row the row i_star
col_max = abs(M[k][k]) / d[k];
i_star = k;
for (size_t i = k + 1; i < n; i++)
{
awikod = abs(M[i][k]) / d[i];
if (abs(awikod) > abs(col_max))
{
col_max = awikod;
i_star = i;
}
}
if (abs(col_max) == 0)
flag = 0;
else
{
if (i_star > k)
{
// make k the pivot row by
// interchanging with i_star
flag *= -1;
itemp = pivot[i_star];
pivot[i_star] = pivot[k];
pivot[k] = itemp;
temp = d[i_star];
d[i_star] = d[k];
d[k] = temp;
for (size_t j = 0; j < n; j++)
{
temp = M[i_star][j];
M[i_star][j] = M[k][j];
M[k][j] = temp;
}
}
// eliminate x[k]
for (size_t i = k + 1; i < n; i++)
{
M[i][k] /= M[k][k];
ratio = M[i][k];
for (size_t j = k + 1; j < n; j++)
M[i][j] -= ratio * M[k][j];
}
}
if (abs(M[n - 1LL][n - 1LL]) == 0) flag = 0;
}
if (flag == 0)
return false;
return true;
}
bool CmpLinearAlgebra::CmpGaussianSolution(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot)
{
if (!CmpGaussianFactor(n, M, pivot))
return false;
return CmpSubstitution(n, n, M, b, x, pivot);
}
bool CmpLinearAlgebra::CmpInverse(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& Mi)
{
std::vector<std::complex<double>> b(n);
std::vector<std::complex<double>> x(n);
std::vector<size_t> pivot(n);
std::vector<std::vector<std::complex<double>>> Mc(n,
std::vector<std::complex<double>>(n));
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
{
Mc[i][j] = M[i][j];
}
}
if (!CmpGaussianFactor(n, Mc, pivot))
return false;
for (size_t i = 0; i < n; i++)
{
b[i] = 0;
}
for (size_t i = 0; i < n; i++)
{
b[i] = 1;
if (!CmpSubstitution(n, n, Mc, b, x, pivot))
return false;
b[i] = 0;
for (size_t j = 0; j < n; j++)
Mi[j][i] = x[j];
}
return true;
}
void CmpLinearAlgebra::CmpCharPolyAndAdjoint(
int n,
std::vector<std::vector<std::complex<double>>>& C,
std::vector<std::vector<std::complex<double>>>& I,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& adjoint,
std::vector<std::complex<double>>& a)
{
C.resize(n, std::vector<std::complex<double>>(n));
I.resize(n, std::vector<std::complex<double>>(n));
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
C[i][j] = I[i][j] = 0;
}
for (size_t i = 0; i < n; i++)
C[i][i] = I[i][i] = 1;
a[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++)
{
std::complex<double> sum = 0.0;
for (size_t l = 0; l < n; l++)
sum += M[j][l] * C[l][k];
C[j][k] = sum;
}
}
std::complex<double> tr = 0.0;
for (size_t j = 0; j < n; j++)
tr += C[j][j];
std::complex<double> ci = 0;
ci._Val[0] = (double)i;
a[i] = -tr / ci;
for (size_t j = 0; j < n; j++)
{
for (size_t k = 0; k < n; k++)
C[j][k] += a[i] * I[j][k];
}
}
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
{
std::complex<double> sum = 0.0;
for (size_t k = 0; k < n; k++)
sum += M[i][k] * C[k][j];
C[i][j] = sum;
}
}
std::complex<double> trace = 0.0;
for (size_t i = 0; i < n; i++)
trace += C[i][i];
trace /= n;
a[n - 1LL] = -trace;
std::complex<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[i][j] = factor * C[i][j];
}
}
void CmpLinearAlgebra::CmpMatrixKernel(
int m, int n,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& X,
size_t& r)
{
std::complex<double> D = 0;
std::vector<int> c(m);
std::vector<int> d(n);
r = 0;
for (size_t i = 0; i < m; i++)
c[i] = -1;
size_t j = 0, k = 1;
Step2:
for (j = 0; j < m; j++)
{
if (abs(M[j][k]) != 0 && c[j] == -1)
break;
}
if (j == m)
{
r++;
d[k] = 0;
goto Step4;
}
D = -1.0 / M[j][k];
M[j][k] = -1;
for (size_t s = k + 1; s < n; s++)
{
M[j][s] = D * M[j][s];
for (size_t i = 0; i < m; i++)
{
if (i != j)
{
D = M[i][k];
M[i][k] = 0;
}
}
}
for (size_t s = k + 1; s < n; s++)
{
for (size_t i = 0; i < m; i++)
{
M[i][s] += D * M[j][s];
}
}
c[j] = (int)k;
d[k] = (int)j;
Step4:
if (k < n - 1)
{
k++;
goto Step2;
}
X.resize(n, std::vector<std::complex<double>>(n));
if (r != 0)
{
for (k = 0; k < n; k++)
{
if (d[k] == 0)
{
for (size_t i = 0; i < n; i++)
{
if (d[i] > 0)
X[k][i] = M[d[i]][k];
else if (i == k)
X[k][i] = 1;
else
X[k][i] = 0;
}
}
}
}
}
void CmpLinearAlgebra::CmpMatrixImage(
int m, int n,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& N,
std::vector<std::vector<std::complex<double>>>& X,
int rank)
{
std::complex<double> D = 0.0;
size_t r = 0;
std::vector<std::vector<std::complex<double>>> copyM(
m, std::vector<std::complex<double>>(n));
std::vector<int> c(m);
std::vector<int> d(n);
for (size_t i = 0; i < m; i++)
c[i] = -1;
size_t j = 0, k = 1;
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
N[i][j] = copyM[i][j] = M[i][j];
}
}
Step2:
for (size_t j = 0; j < m; j++)
{
if (abs(copyM[j][k]) != 0 && c[j] == -1)
break;
}
if (j == m)
{
r++;
d[k] = 0;
goto Step4;
}
D = -1.0 / copyM[j][k];
copyM[j][k] = -1;
for (size_t s = k + 1; s < n; s++)
{
copyM[j][s] = D * copyM[j][s];
for (size_t i = 0; i < m; i++)
{
if (i != j)
{
D = copyM[i][k];
copyM[i][k] = 0;
}
}
}
for (size_t s = k + 1; s < n; s++)
{
for (size_t i = 0; i < m; i++)
{
copyM[i][s] += D * copyM[j][s];
}
}
c[j] = (int)k;
d[k] = (int)j;
Step4:
if (k < (size_t)(n - 1LL))
{
k++;
goto Step2;
}
rank = (int)(n - r);
for (j = 0; j < m; j++)
{
if (c[j] != 0)
{
for (size_t i = 0; i < m; i++)
{
N[i][c[j]] = M[i][c[j]];
}
}
}
}
void CmpLinearAlgebra::CmpGenerateNonSingular(
double scale, std::complex<double>& cDeterminant,
int n, unsigned int seed,
std::vector<std::vector<std::complex<double>>>& Mc)
{
bool failure = false;
std::mt19937 rng(seed);
std::uniform_real_distribution<double> dist(0.0, 1.0);
while (true)
{
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
Mc[i][j]._Val[0] = scale * dist(rng);
Mc[i][j]._Val[1] = scale * dist(rng);
}
}
cDeterminant = CmpDeterminant(n, Mc);
if (cDeterminant._Val[0] != 0 || cDeterminant._Val[1] != 0)
break;
}
}
// Exercises from "Modern Quantum Chemistry an Introduction to Advanced
// Electronic Structure Theory" by Attila Szabo and Neil S. Ostlund
// https://chemistlibrary.wordpress.com/wp-content/uploads/2015/02/modern-quantum-chemistry.pdf
// Program (c) Tuesday, August 12, 2025 by James Pate Williams, Jr.
// Program (c) Tuesday, August 19, 2025 by James Pate Williams, Jr.
#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>
#include "DblLinearAlgebra.h"
#include "CmpLinearAlgebra.h"
int main()
{
// static data matrices
double AArcb[3][3] = { { 2, 3, -1 }, { 4, 4, -3 }, { -2, 3, -1 } };
double AArso[3][3] = { { 1, 1, 0 }, { 1, 2, 2 }, { 0, 2, -1 } };
double BBrso[3][3] = { { 1, -1, 1 }, { -1, 0, 0 }, { 1, 0, 1} };
double BBr[3][3] = { { 1, -1, 1 }, { -1 , 0, 0 }, { 1, 0, 1 } };
double AAcr[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
double AAci[3][3] = { { 1, 1, 2 }, { 3, 0, 1 }, { 0, 2, 4 } };
double BBcr[3][3] = { { 1, 0, 1 }, { 1 , 1, 0 }, { 0, 1, 1 } };
double BBci[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
// some array dimensions
int m = 3, n = 3, p = 3;
// a couple of 3x1 vectors
std::vector<double> br(3);
std::vector<size_t> pivot(3);
// 3x3 real matrices
std::vector<std::vector<double>> Arcb(3, std::vector<double>(3));
std::vector<std::vector<double>> Arso(3, std::vector<double>(3));
std::vector<std::vector<double>> Brso(3, std::vector<double>(3));
std::vector<std::vector<double>> Br(3, std::vector<double>(3));
std::vector<std::vector<double>> Cr(3, std::vector<double>(3));
std::vector<std::vector<double>> Ai(3, std::vector<double>(3));
std::vector<std::vector<double>> Ari(3, std::vector<double>(3));
// a 4x4 real matrix
std::vector<std::vector<double>> Mr(4, std::vector<double>(4));
// 3x3 complex matrices
std::vector<std::vector<std::complex<double>>> Ac(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Bc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Cc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Dc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Ec(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Fc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Gc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Mc(4,
std::vector<std::complex<double>>(4));
// copy static real matrices to dynamic matrices
for (int i = 0; i < m; i++)
{
for (int j = 0; j < p; j++)
{
Arcb[i][j] = AArcb[i][j];
Arso[i][j] = AArso[i][j];
Brso[i][j] = BBrso[i][j];
Ac[i][j]._Val[0] = AAcr[i][j];
Ac[i][j]._Val[1] = AAci[i][j];
}
}
// copy static complex matrices to dynamic matrices
for (int i = 0; i < p; i++)
{
for (int j = 0; j < n; j++)
{
Br[i][j] = BBr[i][j];
Bc[i][j]._Val[0] = BBcr[i][j];
Bc[i][j]._Val[1] = BBci[i][j];
}
}
// See "Elementary Numerical Analysis an
// Algorithmic Approach" (c) 1980 by S. D. Conte
// and Carl de Boor
DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Br, Cr);
std::cout << "Ar * Br = Cr Conte & de Boor" << std::endl;
DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
std::cout << std::endl;
// complex matrix multiplication
CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
std::cout << "Ac * Bc = Cc" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
std::cout << std::endl;
// Exercise 1.2 from Szabo and Ostlund
std::cout << "Exercise 1.2 page 5 Commutator" << std::endl;
DblLinearAlgebra::DblCommutator(3, Arso, Brso, Cr);
DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
std::cout << std::endl;
std::cout << "Exercise 1.2 page 5 Anticommutator" << std::endl;
DblLinearAlgebra::DblAnticommutator(3, Arso, Brso, Cr);
DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
std::cout << std::endl;
CmpLinearAlgebra::CmpAdjoint(3, 3, Cc, Dc);
std::cout << "Exercise 1.3 page 6 Cc adjoint" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Dc);
std::cout << std::endl;
CmpLinearAlgebra::CmpAdjoint(3, 3, Ac, Ec);
CmpLinearAlgebra::CmpAdjoint(3, 3, Bc, Fc);
CmpLinearAlgebra::CmpMultiply(3, 3, 3, Fc, Ec, Gc);
std::cout << "Exercise 1.3 page 6 Bc adjoint * Ac adjoint"
<< std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Gc);
std::cout << std::endl;
std::cout << "Ar matrix" << std::endl;
DblLinearAlgebra::DblPrintMatrix(3, 3, Arcb);
bool inv = DblLinearAlgebra::DblInverse(n, Arcb, Ai);
std::cout << std::endl;
std::cout << "Ar Conte & de Boor inverse flag = "
<< inv << std::endl;
DblLinearAlgebra::DblPrintMatrix(3, 3, Ai);
std::cout << std::endl;
std::cout << "Ar * Ar inverse" << std::endl;
DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Ai, Ari);
DblLinearAlgebra::DblPrintMatrix(3, 3, Ari);
std::cout << std::endl;
double rDeterminant = 0;
DblLinearAlgebra::DblGenerateNonSingular(
2.0, rDeterminant, 4, 1, Mr);
std::cout << "rDeterminant = ";
std::cout << rDeterminant << std::endl;
std::cout << std::endl;
std::cout << "Ac" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Ac);
std::cout << std::endl;
inv = CmpLinearAlgebra::CmpInverse(3, Ac, Bc);
std::cout << "Ac inverse flag = " << inv << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Bc);
CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
std::cout << std::endl;
std::cout << "Ac * Ac inverse" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
std::complex<double> cDeterminant = 0;
CmpLinearAlgebra::CmpGenerateNonSingular(
2.0, cDeterminant, 4, 1, Mc);
std::cout << std::endl;
std::cout << "complex determinant = ";
std::cout << cDeterminant << std::endl;
double rDeterminantA = 0;
std::vector<std::vector<double>> A44r(4,
std::vector<double>(4));
DblLinearAlgebra::DblGenerateNonSingular(
2.0, rDeterminantA, 4, 2, A44r);
double rDeterminantB = 0;
std::vector<std::vector<double>> B44r(4,
std::vector<double>(4));
DblLinearAlgebra::DblGenerateNonSingular(
2.0, rDeterminantB, 4, 3, B44r);
std::cout << std::endl;
std::vector<std::vector<double>> C44r(4,
std::vector<double>(4));
DblLinearAlgebra::DblMultiply(4, 4, 4, A44r, B44r, C44r);
std::cout << "|A| = " << rDeterminantA << std::endl;
std::cout << "|B| = " << rDeterminantB << std::endl;
bool failure = false;
double rDeterminantC =
DblLinearAlgebra::DblDeterminant(4, 0, 0, C44r);
std::cout << "|AB| = " << rDeterminantC << std::endl;
std::cout << "|A||B| = " << rDeterminantA *
rDeterminantB << std::endl;
// Exercise 1.6 with 4x4 complex determinants
std::vector<std::vector<std::complex<double>>> A44c(4,
std::vector<std::complex<double>>(4));
std::cout << std::endl;
std::complex<double> cDeterminantA = 0;
CmpLinearAlgebra::CmpGenerateNonSingular(
2.0, cDeterminantA, 4, 2, A44c);
std::vector<std::vector<std::complex<double>>> B44c(4,
std::vector<std::complex<double>>(4));
std::complex<double> cDeterminantB = 0;
CmpLinearAlgebra::CmpGenerateNonSingular(
2.0, cDeterminantB, 4, 3, B44c);
std::vector<std::vector<std::complex<double>>> C44c(4,
std::vector<std::complex<double>>(4));
CmpLinearAlgebra::CmpMultiply(4, 4, 4, A44c, B44c, C44c);
std::cout << "|A| = " << cDeterminantA << std::endl;
std::cout << "|B| = " << cDeterminantB << std::endl;
failure = false;
std::complex<double> cDeterminantC =
CmpLinearAlgebra::CmpDeterminant(4, C44c);
std::cout << "|AB| = " << cDeterminantC << std::endl;
std::cout << "|A||B| = " << cDeterminantA *
cDeterminantB << std::endl;
std::cout << "\nEnter any key to halt: ";
char line[128] = "";
std::cin.getline(line, 128);
}
Numerical Explorations 