Blog Entry Tuesday, June 16, 2026, by James Pate Williams, Jr., Linear Systems of Equations Solutions

Blog Entry Wednesday June 17 2026Download
#pragma once
#include <vector>

class DirectMethods
{
public:
    static void Substitute(
        int n,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x,
        std::vector<int>& pivot);
    static bool GaussianElimimation(
        int n,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x);
    static void LUDecomposition(
        int n,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x);
};


#include "DirectMethods.h"
#include <vector>

void DirectMethods::Substitute(
    int n,
    std::vector<std::vector<double>>& w,
    std::vector<double>& b,
    std::vector<double>& x,
    std::vector<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 = 1; i < n; i++)
    {
        for (j = 0, sum = 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; j < n; j++)
            sum += w[i][j] * x[j];
        x[i] = (x[i] - sum) / w[i][i];
    }
}

bool DirectMethods::GaussianElimimation(
    int n,
    std::vector<std::vector<double>>& w,
    std::vector<double>& b,
    std::vector<double>& x)
    // returns false if matrix is singular
{
    double awikod, col_max, ratio, row_max, temp;
    std::vector<double> d(n);
    int flag = 1, i, i_star, j, k;
    std::vector<int> pivot(n);

    for (i = 0; i < n; i++)
    {
        pivot[i] = i;
        row_max = 0;
        for (j = 0; j < n; j++)
            row_max = fmax(row_max, fabs(w[i][j]));

        if (row_max == 0)
        {
            flag = 0;
            row_max = 1;
        }

        d[i] = row_max;
    }

    if (n <= 1) return flag != 0;

    // factorization
    
    for (k = 0; k < n - 1; k++)
    {
        // determine pivot row the row i_star
        col_max = fabs(w[k][k]) / d[k];
        i_star = k;
        
        for (i = k + 1; i < n; i++)
        {
            awikod = fabs(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;

    if (flag == 0)
        return false;

    Substitute(n, w, b, x, pivot);
    return true;
}

void DirectMethods::LUDecomposition(
    int n,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x)
{
    std::vector<double> y(n);
    std::vector<std::vector<double>> l(n);
    std::vector<std::vector<double>> u(n);

    for (int k = 0; k < n; k++)
    {
        l[k].resize(n);
        u[k].resize(n);
    }

    for (int k = 0; k < n; k++)
    {
        for (int j = k; j < n; j++)
        {
            double sum = 0.0;

            for (int p = 0; p <= k - 1; p++)
                sum += l[k][p] * u[p][j];

            u[k][j] = a[k][j] - sum;
        }

        for (int i = k + 1; i < n; i++)
        {
            double sum = 0.0;

            for (int p = 0; p < k - 1; p++)
                sum += l[i][p] * u[p][k];

            l[i][k] = (a[i][k] - sum) / u[k][k];
        }
    }

    // forward substitution

    for (int k = 0; k < n; k++)
    {
        double sum = 0.0;

        for (int j = 0; j < k; j++)
            sum += l[k][j] * y[j];

        y[k] = b[k] - sum;
    }

    // backward substitution

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

        for (int j = k + 1; j < n; j++)
            sum += u[k][j] * x[j];

        x[k] = (y[k] - sum) / u[k][k];
    }
}


#pragma once
#include <vector>

class ClassicalIterativeMethods
{
public:
    static double Jacobi(
        int n,
        int maxIterations,
        double tolerance,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x0);
    static double GaussSeidel(
        int n,
        int maxIterations,
        double tolerance,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x0);
    static double SOR(
        int n,
        int maxIterations,
        double omega,
        double tolerance,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x0);
    static double InnerProduct(
        int n,
        std::vector<double>& a,
        std::vector<double>& b);
    static double GradientMethod(
        int n,
        int maxIterations,
        double tolerance,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x0);
};


#include "ClassicalIterativeMethods.h"
#include <vector>

double ClassicalIterativeMethods::Jacobi(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

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

            for (int j = 0; j != i && j < n; j++)
                sum += a[i][j] * x0[j];

            x[i] = (b[i] - sum) / a[i][i];
        }

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::GaussSeidel(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

        for (int i = 0; i < n; i++)
        {
            double sum0 = 0;
            double sum1 = 0;

            for (int j = 0; j <= i - 1; j++)
                sum0 += a[i][j] * x[j];

            for (int j = i + 1; j < n; j++)
                if (i != n)
                    sum1 += a[i][j] * x0[j];

            x[i] = (b[i] - sum0 - sum1) / a[i][i];
        }

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::SOR(
    int n,
    int maxIterations,
    double omega,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

        for (int i = 0; i < n; i++)
        {
            double sum0 = 0;
            double sum1 = 0;

            for (int j = 0; j <= i - 1; j++)
                sum0 += a[i][j] * x[j];

            for (int j = i + 1; j < n; j++)
                if (i != n)
                    sum1 += a[i][j] * x0[j];

            x[i] = (1.0 - omega) * x0[i] +
                omega * (b[i] - sum0 - sum1) / a[i][i];
        }

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::InnerProduct(
    int n, 
    std::vector<double>& a,
    std::vector<double>& b)
{
    double sum = 0.0;

    for (int i = 0; i < n; i++)
        sum += a[i] * b[i];

    return sum;
}

double ClassicalIterativeMethods::GradientMethod(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> r(n);
    std::vector<double> x(n);
    std::vector<double> y(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

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

            for (int j = 0; j < n; j++)
                sum += a[i][j] * x0[j];

            r[i] = b[i] - sum;
        }

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

            for (int j = 0; j < n; j++)
                sum += a[i][j] * r[j];

            y[i] = sum;
        }

        double alpha = InnerProduct(n, r, r);

        alpha /= InnerProduct(n, r, y);

        for (int i = 0; i < n; i++)
            x[i] = x0[i] + alpha * r[i];

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}


#include "ClassicalIterativeMethods.h"
#include <vector>

double ClassicalIterativeMethods::Jacobi(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

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

            for (int j = 0; j != i && j < n; j++)
                sum += a[i][j] * x0[j];

            x[i] = (b[i] - sum) / a[i][i];
        }

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::GaussSeidel(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

        for (int i = 0; i < n; i++)
        {
            double sum0 = 0;
            double sum1 = 0;

            for (int j = 0; j <= i - 1; j++)
                sum0 += a[i][j] * x[j];

            for (int j = i + 1; j < n; j++)
                if (i != n)
                    sum1 += a[i][j] * x0[j];

            x[i] = (b[i] - sum0 - sum1) / a[i][i];
        }

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::SOR(
    int n,
    int maxIterations,
    double omega,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

        for (int i = 0; i < n; i++)
        {
            double sum0 = 0;
            double sum1 = 0;

            for (int j = 0; j <= i - 1; j++)
                sum0 += a[i][j] * x[j];

            for (int j = i + 1; j < n; j++)
                if (i != n)
                    sum1 += a[i][j] * x0[j];

            x[i] = (1.0 - omega) * x0[i] +
                omega * (b[i] - sum0 - sum1) / a[i][i];
        }

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::InnerProduct(
    int n, 
    std::vector<double>& a,
    std::vector<double>& b)
{
    double sum = 0.0;

    for (int i = 0; i < n; i++)
        sum += a[i] * b[i];

    return sum;
}

double ClassicalIterativeMethods::GradientMethod(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> r(n);
    std::vector<double> x(n);
    std::vector<double> y(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

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

            for (int j = 0; j < n; j++)
                sum += a[i][j] * x0[j];

            r[i] = b[i] - sum;
        }

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

            for (int j = 0; j < n; j++)
                sum += a[i][j] * r[j];

            y[i] = sum;
        }

        double alpha = InnerProduct(n, r, r);

        alpha /= InnerProduct(n, r, y);

        for (int i = 0; i < n; i++)
            x[i] = x0[i] + alpha * r[i];

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}


#include "Utility.h"
#include <ctime>
#include <cstring>
#include <windows.h>

wchar_t* Utility::ConvertToWide(
    char* input,
    UINT codePage = CP_UTF8)
{
    if (!input) 
    {
        return nullptr;
    }

    // Determine required wide-character size
    int requiredSize = MultiByteToWideChar(
        codePage,
        MB_ERR_INVALID_CHARS,
        input,
        -1,       // Input is null-terminated
        nullptr,
        0
    );

    if (requiredSize == 0) {
        return nullptr;
    }

    // Allocate buffer (requiredSize includes null terminator)
    wchar_t* output = new (std::nothrow) wchar_t[requiredSize];

    if (!output) {
        return nullptr;
    }

    // Perform conversion
    int result = MultiByteToWideChar(
        codePage,
        MB_ERR_INVALID_CHARS,
        input,
        -1,
        output,
        requiredSize
    );

    if (result == 0) {
        delete[] output;
        output = nullptr;
        return output;
    }

    return output;
}

double Utility::Factorial(int n)
{
    if (n <= 1)
        return 1.0;
    else
        return n * Factorial(n - 1);
}

double Utility::Binomial(int m, int n)
{
    return Factorial(m) / Factorial(m - n) / Factorial(n);
}

std::wstring Utility::Format(double x)
{
    wchar_t line[128] = L"";
    std::wstring result = L"";

    if (x > 0.0)
        result += L"+";
    else if (x == 0.0)
        result += L" ";
    else
        result += L"-";

    x = fabs(x);
    swprintf_s(line, 32, L"%14.8e", x);
    return result + line;
}

bool Utility::FormatTime(
    char buffer[], size_t bufferSize, std::time_t t)
{
    bool result = false;
    const char format[] = "%Y-%m-%d %H:%M:%S";
    
    if (buffer == nullptr || format == nullptr || bufferSize == 0)
    {
        return false;
    }

    std::tm local0 = {};
    
    localtime_s(&local0, &t);

    size_t written = std::strftime(buffer, bufferSize, format, &local0);
    result = written > 0;

    if (!result)
        return false;

    return true;
}

void Utility::Matrix(
    int matrix, int n,
    std::vector<std::vector<double>>& a,
    std::vector<std::vector<double>>& c,
    std::vector<double>& b)
{
    if (matrix == 1)
    {
        // Pascal matrix

        for (int i = 0; i < n; i++)
            b[i] = pow(2, i + 1);

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                a[i][j] = Binomial(i + j, j);
    }

    else if (matrix == 2)
    {
        b[0] = -0.5;

        for (int i = 1; i < n - 1; i++)
            b[i] = -1.5;

        b[n - 1] = +0.5;

        for (int i = 0; i < n; i++)
        {
            a[i][i] = -2.0;

            if (i < n - 1)
                a[i][i + 1LL] = 1.0;

            if (i > 0)
                a[i][i - 1LL] = 1.0;
        }
    }

    else if (matrix == 3)
    {
        // Cauchy matrix

        for (int i = 0; i < n; i++)
        {
            double xi = 2 * i;

            for (int j = 0; j < n; j++)
            {
                double yj = 2 * j + 1;

                a[i][j] = 1.0 / (xi + yj);
            }

            b[i] = i + 1;
        }
    }

    else if (matrix == 4)
    {
        // Lehmer matrix

        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
                a[i][j] = fmin(i + 1, j + 1) + fmax(i + 1, j + 1);

            a[i][i] = i + 1;
            b[i] = i;
        }
    }

    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            c[i][j] = a[i][j];
}


// LinearSystems.cpp : Defines the entry point for the application.
//

#include "framework.h"
#include "LinearSystems.h"
#include "DirectMethods.h"
#include "ClassicalIterativeMethods.h"
#include "Utility.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
WCHAR line[256], text[16384];                   // wide character buffers

// 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_LINEARSYSTEMS, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

    // Perform application initialization:
    if (!InitInstance (hInstance, nCmdShow))
    {
        return FALSE;
    }

    HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_LINEARSYSTEMS));

    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_LINEARSYSTEMS));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_LINEARSYSTEMS);
    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)
{
    static char ascLine[256] = "";
    static wchar_t* unicodePtr = nullptr;
    static int matrix = 0, method = 0;
    static HWND hEditN = NULL;
    static HWND hEditMultiline = NULL;
    static HFONT hFont = NULL;
    static WCHAR matrixStr[32] = L"";
    static std::time_t t0 = {}, t1 = {};

    switch (message)
    {
    case WM_CREATE:
        CreateWindowEx(0, L"STATIC", L"n", WS_CHILD | WS_VISIBLE,
            10, 10, 60, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
        hEditN = CreateWindowEx(
            WS_EX_CLIENTEDGE,                       // Extended style for sunken border
            TEXT("EDIT"),                           // Class name
            TEXT("2"),                              // Initial text (can be blank)
            WS_CHILD | WS_VISIBLE | ES_LEFT | ES_AUTOHSCROLL,
            100, 10, 60, 20,                         // Position and size
            hWnd,                                   // Parent window handle
            (HMENU)IDC_EDIT_N,                      // Unique control ID
            hInst,                                  // Application instance
            NULL                                    // Extra parameter
        );

        hFont = CreateFont(
            16,                // Height of font
            0,                 // Width of font (0 = default)
            0,                 // Escapement angle
            0,                 // Orientation angle
            FW_BOLD,           // Weight (FW_NORMAL, FW_BOLD, etc.)
            FALSE,             // Italic
            FALSE,             // Underline
            FALSE,             // Strikeout
            UNICODE,           // Character set
            OUT_DEFAULT_PRECIS,// Output precision
            CLIP_DEFAULT_PRECIS,// Clipping precision
            DEFAULT_QUALITY,   // Output quality
            FIXED_PITCH | FF_MODERN, // Pitch and family
            TEXT("Courier New")// Typeface name
        );

        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,
            10, 50, 990, 400,                      // 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,
            200, 10, 80, 30, hWnd, (HMENU)IDC_COMPUTE_BUTTON, hInst, NULL);
        CreateWindowEx(0, L"BUTTON", L"Cancel", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
            300, 10, 80, 30, hWnd, (HMENU)IDC_CANCEL_BUTTON, hInst, NULL);
        CreateWindowEx(0, L"BUTTON", L"Clear", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
            400, 10, 80, 30, hWnd, (HMENU)IDC_CLEAR_BUTTON, hInst, NULL);

        SendMessage(hEditMultiline, WM_SETFONT, (WPARAM)hFont, 0);
        return 0;
    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_CAUCHY:
                matrix = 1;
                break;
            case IDM_LEHMER:
                matrix = 2;
                break;
            case IDM_OTHER:
                matrix = 3;
                break;
            case IDM_PASCAL:
                matrix = 4;
                break;
            case IDM_GAUSSIAN:
                method = 1;
                break;
            case IDM_GAUSS_SEIDEL:
                method = 2;
                break;
            case IDM_GRADIENT:
                method = 3;
                break;
            case IDM_JACOBI:
                method = 4;
                break;
            case IDM_LUD:
                method = 5;
                break;
            case IDM_SOR:
                method = 6;
                break;
            case IDC_COMPUTE_BUTTON:
            {
                t0 = time(NULL);
                int n = GetDlgItemInt(hWnd, IDC_EDIT_N, FALSE, FALSE);
                int maxIterations = 20000 * n;
                double omega = 1.5;
                double tolerance = n <= 5 ? 1.0e-12 : 1.0e-8;
                std::vector<double> b(n);
                std::vector<double> x(n);
                std::vector<std::vector<double>> a(n);
                std::vector<std::vector<double>> c(n);

                for (int i = 0; i < n; i++)
                {
                    a[i].resize(n);
                    c[i].resize(n);
                }
                
                if (matrix == 4)
                {
                    // Pascal matrix

                    for (int i = 0; i < n; i++)
                        b[i] = pow(2, i + 1);

                    for (int i = 0; i < n; i++)
                        for (int j = 0; j < n; j++)
                            a[i][j] = Utility::Binomial(i + j, j);

                    wcscpy_s(matrixStr, L"Matrix Type Pascal");
                }

                else if (matrix == 3)
                {
                    wcscpy_s(matrixStr, L"Matrix Type Other");

                    b[0] = -0.5;

                    for (int i = 1; i < n - 1; i++)
                        b[i] = -1.5;

                    b[n - 1LL] = +0.5;

                    for (int i = 0; i < n; i++)
                    {
                        a[i][i] = -2.0;

                        if (i < n - 1)
                            a[i][i + 1LL] = 1.0;

                        if (i > 0)
                            a[i][i - 1LL] = 1.0;
                    }
                }

                else if (matrix == 1)
                {
                    // Cauchy matrix

                    for (int i = 0; i < n; i++)
                    {
                        double xi = 2 * i;

                        for (int j = 0; j < n; j++)
                        {
                            double yj = 2 * j + 1;

                            a[i][j] = 1.0 / (xi + yj);
                        }

                        b[i] = i + 1;
                    }

                    wcscpy_s(matrixStr, L"Matrix Type Cauchy");
                }

                else if (matrix == 2)
                {
                    // Lehmer matrix

                    for (int i = 0; i < n; i++)
                    {
                        for (int j = 0; j < n; j++)
                            a[i][j] = fmin(i + 1, j + 1) + fmax(i + 1, j + 1);

                        a[i][i] = i + 1;
                        b[i] = i;
                    }

                    wcscpy_s(matrixStr, L"Matrix Type Lehmer");
                }

                for (int i = 0; i < n; i++)
                    for (int j = 0; j < n; j++)
                        c[i][j] = a[i][j];

                text[0] = L'\0';
                wcscpy_s(text, matrixStr);
                wcscat_s(text, L"\r\n");

                switch (method)
                {
                case 1:
                    wcscat_s(text, L"Gaussian Elimination:\r\n");
                    DirectMethods::GaussianElimimation(n, a, b, x);
                    break;
                case 2:
                    wcscat_s(text, L"Gauss-Seidel Method:\r\n");
                    ClassicalIterativeMethods::GaussSeidel(
                        n, maxIterations, tolerance, a, b, x);
                    break;
                case 3:
                    wcscat_s(text, L"Gradient Method:\r\n");
                    ClassicalIterativeMethods::GradientMethod(
                        n, maxIterations, tolerance, a, b, x);
                    break;
                case 4:
                    wcscat_s(text, L"Jacobi Method:\r\n");
                    ClassicalIterativeMethods::Jacobi(
                        n, maxIterations, tolerance, a, b, x);
                    break;
                case 5:
                    wcscat_s(text, L"LU Decomposition:\r\n");
                    DirectMethods::LUDecomposition(n, a, b, x);
                    break;
                case 6:
                    wcscat_s(text, L"Successive Over Relaxation:\r\n");
                    ClassicalIterativeMethods::SOR(
                        n, maxIterations, omega, tolerance, a, b, x);
                    break;
                }

                for (int i = 0; i < n; i++)
                {
                    std::wstring format = Utility::Format(x[i]) + L"\r\n";
                    wcscat_s(text, format.c_str());
                }

                double error = 0.0;

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

                    for (int j = 0; j < n; j++)
                        sum += c[i][j] * x[j];

                    error += pow(b[i] - sum, 2.0);
                }

                error = sqrt(error / n);
                swprintf_s(line, 256, L"RMS Error = %14.8e\r\n", error);
                wcscat_s(text, line);

                t1 = time(NULL);
                Utility::FormatTime(ascLine, 256, t0);
                unicodePtr = Utility::ConvertToWide(ascLine, CP_UTF8);
                wcscat_s(text, unicodePtr);
                wcscat_s(text, L"\r\n");
                Utility::FormatTime(ascLine, 256, t1);
                unicodePtr = Utility::ConvertToWide(ascLine, CP_UTF8);
                wcscat_s(text, unicodePtr);
                wcscat_s(text, L"\r\n");
                SetWindowText(hEditMultiline, text);
                return 1;
            }
            case IDC_CANCEL_BUTTON:
            case IDM_EXIT:
                DestroyWindow(hWnd);
                break;
            case IDC_CLEAR_BUTTON:
                text[0] = L'\0';
                SetWindowText(hEditMultiline, text);
                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;
}
Unknown's avatar

Author: jamespatewilliamsjr

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

Leave a comment