#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;
}
Numerical Explorations 