simponic/lizfcm
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

hw 5

Browse Source
This commit is contained in:
Elizabeth Hunt 2023-10-30 19:07:43 -06:00
parent 81979f09cf
commit 562ba9a9b6
Signed by: simponic
GPG Key ID: 52B3774857EB24B1
14 changed files with 717 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

189
doc/software_manual.org
View File

@ -5,7 +5,7 @@
#+STARTUP: entitiespretty fold inlineimages
* Design
The LIZFCM static library (at [[https://github.com/Simponic/math-4610][[https://github.com/Simponic/math-4610]] is a successor to my
The LIZFCM static library (at [[https://github.com/Simponic/math-4610]]) is a successor to my
attempt at writing codes for the Fundamentals of Computational Mathematics course in Common
Lisp, but the effort required to meet the requirement of creating a static library became
too difficult to integrate outside of the ~ASDF~ solution that Common Lisp already brings
@ -14,13 +14,14 @@ to the table.
All of the work established in ~deprecated-cl~ has been painstakingly translated into
the C programming language. I have a couple tenets for its design:
+ Implemntations of routines should all be done immutably in respect to arguments.
+ Implementations of routines should all be done immutably in respect to arguments.
+ Functional programming is good (it's... rough in C though).
+ Routines are separated into "module" c files, and not individual files per function.
+ Routines are separated into "modules" that follow a form of separation of concerns
in files, and not individual files per function.
* Compilation
A provided ~Makefile~ is added for convencience. It has been tested on an M1 machine running MacOS as
well as Arch Linux.
A provided ~Makefile~ is added for convencience. It has been tested on an ~arm~-based M1 machine running
MacOS as well as ~x86~ Arch Linux.
1. ~cd~ into the root of the repo
2. ~make~
@ -317,6 +318,39 @@ void free_vector(Array_double *v) {
}
#+END_SRC
*** ~add_element~
+ Author: Elizabeth Hunt
+ Name: ~add_element~
+ Location: ~src/vector.c~
+ Input: a pointer to an ~Array_double~
+ Output: a new ~Array_double~ with element ~x~ appended.
#+BEGIN_SRC c
Array_double *add_element(Array_double *v, double x) {
Array_double *pushed = InitArrayWithSize(double, v->size + 1, 0.0);
for (size_t i = 0; i < v->size; ++i)
pushed->data[i] = v->data[i];
pushed->data[v->size] = x;
return pushed;
}
#+END_SRC
*** ~slice_element~
+ Author: Elizabeth Hunt
+ Name: ~slice_element~
+ Location: ~src/vector.c~
+ Input: a pointer to an ~Array_double~
+ Output: a new ~Array_double~ with element ~x~ sliced.
#+BEGIN_SRC c
Array_double *slice_element(Array_double *v, size_t x) {
Array_double *sliced = InitArrayWithSize(double, v->size - 1, 0.0);
for (size_t i = 0; i < v->size - 1; ++i)
sliced->data[i] = i >= x ? v->data[i + 1] : v->data[i];
return sliced;
}
#+END_SRC
*** ~copy_vector~
+ Author: Elizabeth Hunt
+ Name: ~copy_vector~
@ -370,25 +404,24 @@ void format_vector_into(Array_double *v, char *s) {
matrix $L$, $U$, respectively such that $LU = m$.
+ Output: a pointer to the location in memory in which two ~Matrix_double~'s reside: the first
representing $L$, the second, $U$.
+ Errors: Exits and throws a status code of ~-1~ when encountering a matrix that cannot be
+ Errors: Fails assertions when encountering a matrix that cannot be
decomposed
#+BEGIN_SRC c
Matrix_double **lu_decomp(Matrix_double *m) {
Matrix_double **lu_decomp(Matrix_double *m) {
assert(m->cols == m->rows);
Matrix_double *u = copy_matrix(m);
Matrix_double *l_empt = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
Matrix_double *l = put_identity_diagonal(l_empt);
free(l_empt);
free_matrix(l_empt);
Matrix_double **u_l = malloc(sizeof(Matrix_double *) * 2);
for (size_t y = 0; y < m->rows; y++) {
if (u->data[y]->data[y] == 0) {
printf("ERROR: a pivot is zero in given matrix\n");
exit(-1);
assert(false);
}
}
@ -399,7 +432,7 @@ void format_vector_into(Array_double *v, char *s) {
if (denom == 0) {
printf("ERROR: non-factorable matrix\n");
exit(-1);
assert(false);
}
double factor = -(u->data[y]->data[x] / denom);
@ -418,7 +451,7 @@ void format_vector_into(Array_double *v, char *s) {
u_l[0] = u;
u_l[1] = l;
return u_l;
}
}
#+END_SRC
*** ~bsubst~
+ Author: Elizabeth Hunt
@ -467,7 +500,7 @@ Array_double *fsubst(Matrix_double *l, Array_double *b) {
}
#+END_SRC
*** ~solve_matrix~
*** ~solve_matrix_lu_bsubst~
+ Author: Elizabeth Hunt
+ Location: ~src/matrix.c~
+ Input: a pointer to a ~Matrix_double~ $m$ and a pointer to an ~Array_double~ $b$
@ -480,7 +513,7 @@ Here we make use of forward substitution to first solve $Ly = b$ given $L$ as th
Then, $LUx = b$, thus $x$ is a solution.
#+BEGIN_SRC c
Array_double *solve_matrix(Matrix_double *m, Array_double *b) {
Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b) {
assert(b->size == m->rows);
assert(m->rows == m->cols);
@ -495,11 +528,97 @@ Array_double *solve_matrix(Matrix_double *m, Array_double *b) {
free_matrix(u);
free_matrix(l);
free(u_l);
return x;
}
#+END_SRC
*** ~gaussian_elimination~
+ Author: Elizabeth Hunt
+ Location: ~src/matrix.c~
+ Input: a pointer to a ~Matrix_double~ $m$
+ Output: a pointer to a copy of $m$ in reduced echelon form
This works by finding the row with a maximum value in the column $k$. Then, it uses that as a pivot, and
applying reduction to all other rows. The general idea is available at [[https://en.wikipedia.org/wiki/Gaussian_elimination]].
#+BEGIN_SRC c
Matrix_double *gaussian_elimination(Matrix_double *m) {
uint64_t h = 0;
uint64_t k = 0;
Matrix_double *m_cp = copy_matrix(m);
while (h < m_cp->rows && k < m_cp->cols) {
uint64_t max_row = 0;
double total_max = 0.0;
for (uint64_t row = h; row < m_cp->rows; row++) {
double this_max = c_max(fabs(m_cp->data[row]->data[k]), total_max);
if (c_max(this_max, total_max) == this_max) {
max_row = row;
}
}
if (max_row == 0) {
k++;
continue;
}
Array_double *swp = m_cp->data[max_row];
m_cp->data[max_row] = m_cp->data[h];
m_cp->data[h] = swp;
for (uint64_t row = h + 1; row < m_cp->rows; row++) {
double factor = m_cp->data[row]->data[k] / m_cp->data[h]->data[k];
m_cp->data[row]->data[k] = 0.0;
for (uint64_t col = k + 1; col < m_cp->cols; col++) {
m_cp->data[row]->data[col] -= m_cp->data[h]->data[col] * factor;
}
}
h++;
k++;
}
return m_cp;
}
#+END_SRC
*** ~solve_matrix_gaussian~
+ Author: Elizabeth Hunt
+ Location: ~src/matrix.c~
+ Input: a pointer to a ~Matrix_double~ $m$ and a target ~Array_double~ $b$
+ Output: a pointer to a vector $x$ being the solution to the equation $mx = b$
We first perform ~gaussian_elimination~ after augmenting $m$ and $b$. Then, as $m$ is in reduced echelon form, it's an upper
triangular matrix, so we can perform back substitution to compute $x$.
#+BEGIN_SRC c
Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b) {
assert(b->size == m->rows);
assert(m->rows == m->cols);
Matrix_double *m_augment_b = add_column(m, b);
Matrix_double *eliminated = gaussian_elimination(m_augment_b);
Array_double *b_gauss = col_v(eliminated, m->cols);
Matrix_double *u = slice_column(eliminated, m->rows);
Array_double *solution = bsubst(u, b_gauss);
free_matrix(m_augment_b);
free_matrix(eliminated);
free_matrix(u);
free_vector(b_gauss);
return solution;
}
#+END_SRC
*** ~m_dot_v~
+ Author: Elizabeth Hunt
+ Location: ~src/matrix.c~
@ -535,6 +654,48 @@ Matrix_double *put_identity_diagonal(Matrix_double *m) {
}
#+END_SRC
*** ~slice_column~
+ Author: Elizabeth Hunt
+ Location: ~src/matrix.c~
+ Input: a pointer to a ~Matrix_double~
+ Output: a pointer to a copy of the given ~Matrix_double~ with column at ~x~ sliced
#+BEGIN_SRC c
Matrix_double *slice_column(Matrix_double *m, size_t x) {
Matrix_double *sliced = copy_matrix(m);
for (size_t row = 0; row < m->rows; row++) {
Array_double *old_row = sliced->data[row];
sliced->data[row] = slice_element(old_row, x);
free_vector(old_row);
}
sliced->cols--;
return sliced;
}
#+END_SRC
*** ~add_column~
+ Author: Elizabet Hunt
+ Location: ~src/matrix.c~
+ Input: a pointer to a ~Matrix_double~ and a new vector representing the appended column ~x~
+ Output: a pointer to a copy of the given ~Matrix_double~ with a new column ~x~
#+BEGIN_SRC c
Matrix_double *add_column(Matrix_double *m, Array_double *v) {
Matrix_double *pushed = copy_matrix(m);
for (size_t row = 0; row < m->rows; row++) {
Array_double *old_row = pushed->data[row];
pushed->data[row] = add_element(old_row, v->data[row]);
free_vector(old_row);
}
pushed->cols++;
return pushed;
}
#+END_SRC
*** ~copy_matrix~
+ Author: Elizabeth Hunt
+ Location: ~src/matrix.c~

BIN
doc/software_manual.pdf
View File

Binary file not shown.

96
doc/software_manual.tex
View File

@ -1,4 +1,4 @@
% Created 2023-10-18 Wed 13:06
% Created 2023-10-30 Mon 09:44
% Intended LaTeX compiler: pdflatex
\documentclass[11pt]{article}
\usepackage[utf8]{inputenc}
@ -31,8 +31,8 @@
\setlength\parindent{0pt}
\section{Design}
\label{sec:org7204e7c}
The LIZFCM static library (at \href{https://github.com/Simponic/math-4610}{[https://github.com/Simponic/math-4610} is a successor to my
\label{sec:org7f9c526}
The LIZFCM static library (at \url{https://github.com/Simponic/math-4610}) is a successor to my
attempt at writing codes for the Fundamentals of Computational Mathematics course in Common
Lisp, but the effort required to meet the requirement of creating a static library became
too difficult to integrate outside of the \texttt{ASDF} solution that Common Lisp already brings
@ -42,13 +42,13 @@ All of the work established in \texttt{deprecated-cl} has been painstakingly tra
the C programming language. I have a couple tenets for its design:
\begin{itemize}
\item Implemntations of routines should all be done immutably in respect to arguments.
\item Implementations of routines should all be done immutably in respect to arguments.
\item Functional programming is good (it's\ldots{} rough in C though).
\item Routines are separated into "module" c files, and not individual files per function.
\end{itemize}
\section{Compilation}
\label{sec:org16cc307}
\label{sec:org911a41e}
A provided \texttt{Makefile} is added for convencience. It has been tested on an M1 machine running MacOS as
well as Arch Linux.
@ -74,11 +74,11 @@ Which is then bundled into a static library in \texttt{lib/lizfcm.a} which can b
in the standard method.
\section{The LIZFCM API}
\label{sec:org832532a}
\label{sec:orgd74cd2d}
\subsection{Simple Routines}
\label{sec:org540b602}
\label{sec:org66bba13}
\subsubsection{\texttt{smaceps}}
\label{sec:org4d03b6e}
\label{sec:orgeae9531}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{smaceps}
@ -104,7 +104,7 @@ float smaceps() {
\end{verbatim}
\subsubsection{\texttt{dmaceps}}
\label{sec:org2603bfc}
\label{sec:org237c904}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{dmaceps}
@ -130,9 +130,9 @@ double dmaceps() {
\end{verbatim}
\subsection{Derivative Routines}
\label{sec:org95c28e9}
\label{sec:org9cf9027}
\subsubsection{\texttt{central\_derivative\_at}}
\label{sec:org950de62}
\label{sec:org1fcd333}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{central\_derivative\_at}
@ -163,7 +163,7 @@ double central_derivative_at(double (*f)(double), double a, double h) {
\end{verbatim}
\subsubsection{\texttt{forward\_derivative\_at}}
\label{sec:org832eda6}
\label{sec:org6a768fc}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{forward\_derivative\_at}
@ -194,7 +194,7 @@ double forward_derivative_at(double (*f)(double), double a, double h) {
\end{verbatim}
\subsubsection{\texttt{backward\_derivative\_at}}
\label{sec:org591836d}
\label{sec:org610ce76}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{backward\_derivative\_at}
@ -225,9 +225,9 @@ double backward_derivative_at(double (*f)(double), double a, double h) {
\end{verbatim}
\subsection{Vector Routines}
\label{sec:org5254fe4}
\label{sec:orgfd176e6}
\subsubsection{Vector Arithmetic: \texttt{add\_v, minus\_v}}
\label{sec:orgf802d61}
\label{sec:org2dbc55f}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name(s): \texttt{add\_v}, \texttt{minus\_v}
@ -258,7 +258,7 @@ Array_double *minus_v(Array_double *v1, Array_double *v2) {
\end{verbatim}
\subsubsection{Norms: \texttt{l1\_norm}, \texttt{l2\_norm}, \texttt{linf\_norm}}
\label{sec:orgc56e22d}
\label{sec:org53a2ffc}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name(s): \texttt{l1\_norm}, \texttt{l2\_norm}, \texttt{linf\_norm}
@ -292,7 +292,7 @@ double linf_norm(Array_double *v) {
\end{verbatim}
\subsubsection{\texttt{vector\_distance}}
\label{sec:orgb54922f}
\label{sec:org1fb3f8c}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{vector\_distance}
@ -313,7 +313,7 @@ double vector_distance(Array_double *v1, Array_double *v2,
\end{verbatim}
\subsubsection{Distances: \texttt{l1\_distance}, \texttt{l2\_distance}, \texttt{linf\_distance}}
\label{sec:orgf22f8e0}
\label{sec:org4a25a94}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name(s): \texttt{l1\_distance}, \texttt{l2\_distance}, \texttt{linf\_distance}
@ -339,7 +339,7 @@ double linf_distance(Array_double *v1, Array_double *v2) {
\end{verbatim}
\subsubsection{\texttt{sum\_v}}
\label{sec:org4593341}
\label{sec:org035a547}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{sum\_v}
@ -359,7 +359,7 @@ double sum_v(Array_double *v) {
\subsubsection{\texttt{scale\_v}}
\label{sec:org3123f61}
\label{sec:org12b0853}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{scale\_v}
@ -378,7 +378,7 @@ Array_double *scale_v(Array_double *v, double m) {
\end{verbatim}
\subsubsection{\texttt{free\_vector}}
\label{sec:org983efcf}
\label{sec:org70ba90c}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{free\_vector}
@ -396,7 +396,7 @@ void free_vector(Array_double *v) {
\end{verbatim}
\subsubsection{\texttt{copy\_vector}}
\label{sec:orgde05d32}
\label{sec:org57afc74}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{copy\_vector}
@ -416,7 +416,7 @@ Array_double *copy_vector(Array_double *v) {
\end{verbatim}
\subsubsection{\texttt{format\_vector\_into}}
\label{sec:org2e779f3}
\label{sec:orgc346c3c}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{format\_vector\_into}
@ -446,9 +446,9 @@ void format_vector_into(Array_double *v, char *s) {
\end{verbatim}
\subsection{Matrix Routines}
\label{sec:org2354147}
\label{sec:org3b053ab}
\subsubsection{\texttt{lu\_decomp}}
\label{sec:org3690faa}
\label{sec:org5553968}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{lu\_decomp}
@ -509,7 +509,7 @@ Matrix_double **lu_decomp(Matrix_double *m) {
}
\end{verbatim}
\subsubsection{\texttt{bsubst}}
\label{sec:orgdeba296}
\label{sec:org253efdc}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{bsubst}
@ -534,7 +534,7 @@ Array_double *bsubst(Matrix_double *u, Array_double *b) {
}
\end{verbatim}
\subsubsection{\texttt{fsubst}}
\label{sec:org60d3435}
\label{sec:orge0c7bc6}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{fsubst}
@ -562,7 +562,7 @@ Array_double *fsubst(Matrix_double *l, Array_double *b) {
\end{verbatim}
\subsubsection{\texttt{solve\_matrix}}
\label{sec:org914121f}
\label{sec:orgbcd445a}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -598,7 +598,7 @@ Array_double *solve_matrix(Matrix_double *m, Array_double *b) {
\end{verbatim}
\subsubsection{\texttt{m\_dot\_v}}
\label{sec:orgae0f4c9}
\label{sec:orga9b1f68}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -620,7 +620,7 @@ Array_double *m_dot_v(Matrix_double *m, Array_double *v) {
\end{verbatim}
\subsubsection{\texttt{put\_identity\_diagonal}}
\label{sec:org6d84f6a}
\label{sec:org33ead5e}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -639,7 +639,7 @@ Matrix_double *put_identity_diagonal(Matrix_double *m) {
\end{verbatim}
\subsubsection{\texttt{copy\_matrix}}
\label{sec:orge750c56}
\label{sec:org34b3f5b}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -659,7 +659,7 @@ Matrix_double *copy_matrix(Matrix_double *m) {
\end{verbatim}
\subsubsection{\texttt{free\_matrix}}
\label{sec:org4ebcf85}
\label{sec:org9c91101}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -678,7 +678,7 @@ void free_matrix(Matrix_double *m) {
\end{verbatim}
\subsubsection{\texttt{format\_matrix\_into}}
\label{sec:org308ee0d}
\label{sec:org51f3e27}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{format\_matrix\_into}
@ -705,9 +705,9 @@ void format_matrix_into(Matrix_double *m, char *s) {
}
\end{verbatim}
\subsection{Root Finding Methods}
\label{sec:org8981156}
\label{sec:org0e83d47}
\subsubsection{\texttt{find\_ivt\_range}}
\label{sec:orga5835b0}
\label{sec:org3e4e34e}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{find\_ivt\_range}
@ -737,7 +737,7 @@ double *find_ivt_range(double (*f)(double), double start_x, double delta,
}
\end{verbatim}
\subsubsection{\texttt{bisect\_find\_root}}
\label{sec:orgb118fc7}
\label{sec:org48f0967}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name(s): \texttt{bisect\_find\_root}
@ -765,7 +765,7 @@ double bisect_find_root(double (*f)(double), double a, double b,
}
\end{verbatim}
\subsubsection{\texttt{bisect\_find\_root\_with\_error\_assumption}}
\label{sec:orgf5124e7}
\label{sec:org15e3c2d}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{bisect\_find\_root\_with\_error\_assumption}
@ -789,9 +789,9 @@ double bisect_find_root_with_error_assumption(double (*f)(double), double a,
\end{verbatim}
\subsection{Linear Routines}
\label{sec:org6f4fce5}
\label{sec:org98cb54b}
\subsubsection{\texttt{least\_squares\_lin\_reg}}
\label{sec:orge810f5f}
\label{sec:org0c0c5d7}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{least\_squares\_lin\_reg}
@ -821,12 +821,12 @@ Line *least_squares_lin_reg(Array_double *x, Array_double *y) {
}
\end{verbatim}
\subsection{Appendix / Miscellaneous}
\label{sec:org85d2eae}
\label{sec:orge34af18}
\subsubsection{Data Types}
\label{sec:org198ca2d}
\label{sec:org0f2f877}
\begin{enumerate}
\item \texttt{Line}
\label{sec:org1866885}
\label{sec:org3f27166}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/types.h}
@ -839,7 +839,7 @@ typedef struct Line {
} Line;
\end{verbatim}
\item The \texttt{Array\_<type>} and \texttt{Matrix\_<type>}
\label{sec:org4a1c956}
\label{sec:org83fc1f3}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/types.h}
@ -871,10 +871,10 @@ typedef struct {
\end{enumerate}
\subsubsection{Macros}
\label{sec:org1976330}
\label{sec:org2bf9bf0}
\begin{enumerate}
\item \texttt{c\_max} and \texttt{c\_min}
\label{sec:org208b148}
\label{sec:orgcaa569e}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/macros.h}
@ -888,7 +888,7 @@ typedef struct {
\end{verbatim}
\item \texttt{InitArray}
\label{sec:orgccc4528}
\label{sec:org5805999}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/macros.h}
@ -909,7 +909,7 @@ typedef struct {
\end{verbatim}
\item \texttt{InitArrayWithSize}
\label{sec:org7e87550}
\label{sec:org264d6b7}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/macros.h}
@ -930,7 +930,7 @@ typedef struct {
\end{verbatim}
\item \texttt{InitMatrixWithSize}
\label{sec:orge6ec2b1}
\label{sec:org310d41a}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/macros.h}

59
homeworks/hw-5.org Normal file
View File

@ -0,0 +1,59 @@
#+TITLE: Homework 5
#+AUTHOR: Elizabeth Hunt
#+LATEX_HEADER: \notindent \notag \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
#+LATEX: \setlength\parindent{0pt}
#+OPTIONS: toc:nil
* Question One
See LIZFCM \rightarrow Matrix Routines \rightarrow ~lu decomp~ & ~bsubst~.
The test ~UTEST(matrix, lu_decomp)~ is a unit test for the ~lu_decomp~ routine,
and ~UTEST(matrix, bsubst)~ verifies back substitution on an upper triangular
3 \times 3 matrix with a known solution that can be verified manually.
Both can be found in ~tests/matrix.t.c~.
* Question Two
Unless the following are met, the resulting solution will be garbage.
1. The matrix $U$ must be not be singular.
2. $U$ must be square (or it will fail the ~assert~).
3. The system created by $Ux = b$ must be consistent.
4. $U$ is in (obviously) upper-triangular form.
Thus, the actual calculation performing the $LU$ decomposition
(in ~lu_decomp~) does a sanity
check for 1-3 will fail an assert, should a point along the diagonal (pivot) be
zero, or the matrix be non-factorable.
* Question Three
See LIZFCM \rightarrow Matrix Routines \rightarrow ~fsubst~.
~UTEST(matrix, fsubst)~ verifies forward substitution on a lower triangular 3 \times 3
matrix with a known solution that can be verified manually.
* Question Four
See LIZFCM \rightarrow Matrix Routines \rightarrow ~gaussian_elimination~ and ~solve_gaussian_elimination~.
* Question Five
See LIZFCM \rightarrow Matrix Routines \rightarrow ~m_dot_v~, and the ~UTEST(matrix, m_dot_v)~ in
~tests/matrix.t.c~.
* Question Six
See ~UTEST(matrix, solve_gaussian_elimination)~ in ~tests/matrix.t.c~, which generates a diagonally dominant 10 \times 10 matrix
and shows that the solution is consistent with the initial matrix, according to the steps given. Then,
we do a dot product between each row of the diagonally dominant matrix and the solution vector to ensure
it is near equivalent to the input vector.
* Question Seven
See ~UTEST(matrix, solve_matrix_lu_bsubst)~ which does the same test in Question Six with the solution according to
~solve_matrix_lu_bsubst~ as shown in the Software Manual.
* Question Eight
No, since the time complexity for Gaussian Elimination is always less than that of the LU factorization solution by $O(n^2)$ operations
(in LU factorization we perform both backwards and forwards substitutions proceeding the LU decomp, in Gaussian Elimination we only need
back substitution).
* Question Nine, Ten
See LIZFCM Software manual and shared library in ~dist~ after compiling.

BIN
homeworks/hw-5.pdf Normal file
View File

Binary file not shown.

95
homeworks/hw-5.tex Normal file
View File

@ -0,0 +1,95 @@
% Created 2023-10-30 Mon 19:05
% Intended LaTeX compiler: pdflatex
\documentclass[11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{graphicx}
\usepackage{longtable}
\usepackage{wrapfig}
\usepackage{rotating}
\usepackage[normalem]{ulem}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{capt-of}
\usepackage{hyperref}
\notindent \notag \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
\author{Elizabeth Hunt}
\date{\today}
\title{Homework 5}
\hypersetup{
pdfauthor={Elizabeth Hunt},
pdftitle={Homework 5},
pdfkeywords={},
pdfsubject={},
pdfcreator={Emacs 28.2 (Org mode 9.7-pre)},
pdflang={English}}
\begin{document}
\maketitle
\setlength\parindent{0pt}
\section{Question One}
\label{sec:org88abf18}
See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{lu decomp} \& \texttt{bsubst}.
The test \texttt{UTEST(matrix, lu\_decomp)} is a unit test for the \texttt{lu\_decomp} routine,
and \texttt{UTEST(matrix, bsubst)} verifies back substitution on an upper triangular
3 \texttimes{} 3 matrix with a known solution that can be verified manually.
Both can be found in \texttt{tests/matrix.t.c}.
\section{Question Two}
\label{sec:org098a7f1}
Unless the following are met, the resulting solution will be garbage.
\begin{enumerate}
\item The matrix \(U\) must be not be singular.
\item \(U\) must be square (or it will fail the \texttt{assert}).
\item The system created by \(Ux = b\) must be consistent.
\item \(U\) is in (obviously) upper-triangular form.
\end{enumerate}
Thus, the actual calculation performing the \(LU\) decomposition
(in \texttt{lu\_decomp}) does a sanity
check for 1-3 will fail an assert, should a point along the diagonal (pivot) be
zero, or the matrix be non-factorable.
\section{Question Three}
\label{sec:org40d5983}
See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{fsubst}.
\texttt{UTEST(matrix, fsubst)} verifies forward substitution on a lower triangular 3 \texttimes{} 3
matrix with a known solution that can be verified manually.
\section{Question Four}
\label{sec:orgf7d23bb}
See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{gaussian\_elimination} and \texttt{solve\_gaussian\_elimination}.
\section{Question Five}
\label{sec:org54e966c}
See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{m\_dot\_v}, and the \texttt{UTEST(matrix, m\_dot\_v)} in
\texttt{tests/matrix.t.c}.
\section{Question Six}
\label{sec:org413b527}
See \texttt{UTEST(matrix, solve\_gaussian\_elimination)} in \texttt{tests/matrix.t.c}, which generates a diagonally dominant 10 \texttimes{} 10 matrix
and shows that the solution is consistent with the initial matrix, according to the steps given. Then,
we do a dot product between each row of the diagonally dominant matrix and the solution vector to ensure
it is near equivalent to the input vector.
\section{Question Seven}
\label{sec:orgd3d7443}
See \texttt{UTEST(matrix, solve\_matrix\_lu\_bsubst)} which does the same test in Question Six with the solution according to
\texttt{solve\_matrix\_lu\_bsubst} as shown in the Software Manual.
\section{Question Eight}
\label{sec:orgf8ac9bf}
No, since the time complexity for Gaussian Elimination is always less than that of the LU factorization solution by \(O(n^2)\) operations
(in LU factorization we perform both backwards and forwards substitutions proceeding the LU decomp, in Gaussian Elimination we only need
back substitution).
\section{Question Nine, Ten}
\label{sec:orgb270171}
See LIZFCM Software manual and shared library in \texttt{dist} after compiling.
\end{document}

8
inc/lizfcm.h
View File

@ -25,6 +25,8 @@ extern double l2_distance(Array_double *v1, Array_double *v2);
extern double l1_distance(Array_double *v1, Array_double *v2);
extern double linf_distance(Array_double *v1, Array_double *v2);
extern Array_double *copy_vector(Array_double *v1);
extern Array_double *add_element(Array_double *v, double x);
extern Array_double *slice_element(Array_double *v, size_t x);
extern void free_vector(Array_double *v);
extern void format_vector_into(Array_double *v, char *s);
extern int vector_equal(Array_double *a, Array_double *b);
@ -33,11 +35,15 @@ extern Matrix_double *put_identity_diagonal(Matrix_double *m);
extern Matrix_double **lu_decomp(Matrix_double *m);
extern Array_double *bsubst(Matrix_double *u, Array_double *b);
extern Array_double *fsubst(Matrix_double *l, Array_double *b);
extern Array_double *solve_matrix(Matrix_double *m, Array_double *b);
extern Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b);
extern Matrix_double *gaussian_elimination(Matrix_double *m);
extern Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b);
extern Array_double *m_dot_v(Matrix_double *m, Array_double *v);
extern Matrix_double *m_dot_m(Matrix_double *a, Matrix_double *b);
extern Array_double *col_v(Matrix_double *m, size_t x);
extern Matrix_double *copy_matrix(Matrix_double *m);
extern Matrix_double *add_column(Matrix_double *m, Array_double *col);
extern Matrix_double *slice_column(Matrix_double *m, size_t col);
extern void free_matrix(Matrix_double *m);
extern void format_matrix_into(Matrix_double *m, char *s);
extern int matrix_equal(Matrix_double *a, Matrix_double *b);

4
inc/macros.h
View File

@ -52,7 +52,7 @@
#define c_max(x, y) (((x) >= (y)) ? (x) : (y))
#define c_min(x, y) (((x) <= (y)) ? (x) : (y))
#define true 1;
#define false 0;
#define true 1
#define false 0
#endif // MACROS_H

26
notes/Oct-27.org Normal file
View File

@ -0,0 +1,26 @@
Use a bisection criterion for a start
Hybrid Method: combine Bisection and Higher Order Method:
- Newton's Method
- Secant Method (Newton's method with secant approx.)
#+BEGIN_SRC c
fa = f(a)
fb = f(b)
if (fa * fb >= 0) return
error = 10 * tol
iter = 0
while (error > tol && iter < maxiter) {
x0 = 0.5 * (a + b)
x1 = x0 - f(x0) / f'(x0)
if (abs(x1 - x0) > 0.5 * (b - a)) {
// do bisection
} else{
// do newton's method
}
}
#+END_SRC

34
notes/Oct-30.org Normal file
View File

@ -0,0 +1,34 @@
* Power Method for computing the largest eigenvalue of a square matrix
An eigenvector, v \in R^n is a nonzero vector such that for some number, \lambda \in C, Av = \lambda v
\Rightarrow || v || = 1
Suppose we start with some vector v and assume, v = \alpha_0 v_0 + \alpha_1 v_1 + \cdots + \alpha_n v_n, where {v_1, \cdots, v_n}
are the eigenvectors of A. Assume {v_1, \cdots, v_n} is a basis for R^n
We can order the eigenvalues such that \lambda_1 \ge \lambda_2 \ge \lambda_3 \ge \cdots \ge \lambda_n
Compute u = Av
= A(\alpha_1 v_1 + \cdots + \alpha_n v_n)
= \alpha_1 Av_1 + A(\cdots) + \alpha_n A v_n
= \alpha_1 \lambda_1 v_1 + \alpha_2 \lambda_2 v_2 + \cdots + \alpha_n \lambda_n v_n
w = A (Av)
= \alpha_1 \lambda_1^2 v_1 + \alpha_2 \lambda_2^2 v_2 + \cdots + \alpha_n \lambda_n^2 v_n
Thus,
A^k v = \alpha_1 \lambda_1^k v_1 + \alpha_2 \lambda_2^k v_2 + \cdots + \alpha_n \lambda_n^k v_n
= \lambda_1^k ( \alpha_1 v_1 + \alpha_2 \frac{\lambda_2^k}{\lambda_1^k} v_2 + \cdots + \alpha_n \frac{\lambda_3^k}{\lambda_1^k} v_n)
As k \rightarrow \infty
A^k v = \lambda_1^k (\alpha_1 v_1) + \text{negligble terms}
Algorithm:
v \ne 0 with v \in R^n
y = Av = \alpha_1 v_1 + \cdots + \alpha_n v_n
w = \frac{1}{||y||} \cdot y
Rayleigh Quotient:
If $v$ is an eigenvector of A with eigenvalue \lambda then \frac{v^T A v}{v^T v} = \lambda

96
src/matrix.c
View File

@ -1,5 +1,6 @@
#include "lizfcm.h"
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
@ -71,7 +72,7 @@ Matrix_double **lu_decomp(Matrix_double *m) {
for (size_t y = 0; y < m->rows; y++) {
if (u->data[y]->data[y] == 0) {
printf("ERROR: a pivot is zero in given matrix\n");
exit(-1);
assert(false);
}
}
@ -82,7 +83,7 @@ Matrix_double **lu_decomp(Matrix_double *m) {
if (denom == 0) {
printf("ERROR: non-factorable matrix\n");
exit(-1);
assert(false);
}
double factor = -(u->data[y]->data[x] / denom);
@ -129,7 +130,7 @@ Array_double *fsubst(Matrix_double *l, Array_double *b) {
return x;
}
Array_double *solve_matrix(Matrix_double *m, Array_double *b) {
Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b) {
assert(b->size == m->rows);
assert(m->rows == m->cols);
@ -144,10 +145,99 @@ Array_double *solve_matrix(Matrix_double *m, Array_double *b) {
free_matrix(u);
free_matrix(l);
free(u_l);
return x;
}
Matrix_double *gaussian_elimination(Matrix_double *m) {
uint64_t h = 0;
uint64_t k = 0;
Matrix_double *m_cp = copy_matrix(m);
while (h < m_cp->rows && k < m_cp->cols) {
uint64_t max_row = 0;
double total_max = 0.0;
for (uint64_t row = h; row < m_cp->rows; row++) {
double this_max = c_max(fabs(m_cp->data[row]->data[k]), total_max);
if (c_max(this_max, total_max) == this_max) {
max_row = row;
}
}
if (max_row == 0) {
k++;
continue;
}
Array_double *swp = m_cp->data[max_row];
m_cp->data[max_row] = m_cp->data[h];
m_cp->data[h] = swp;
for (uint64_t row = h + 1; row < m_cp->rows; row++) {
double factor = m_cp->data[row]->data[k] / m_cp->data[h]->data[k];
m_cp->data[row]->data[k] = 0.0;
for (uint64_t col = k + 1; col < m_cp->cols; col++) {
m_cp->data[row]->data[col] -= m_cp->data[h]->data[col] * factor;
}
}
h++;
k++;
}
return m_cp;
}
Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b) {
assert(b->size == m->rows);
assert(m->rows == m->cols);
Matrix_double *m_augment_b = add_column(m, b);
Matrix_double *eliminated = gaussian_elimination(m_augment_b);
Array_double *b_gauss = col_v(eliminated, m->cols);
Matrix_double *u = slice_column(eliminated, m->rows);
Array_double *solution = bsubst(u, b_gauss);
free_matrix(m_augment_b);
free_matrix(eliminated);
free_matrix(u);
free_vector(b_gauss);
return solution;
}
Matrix_double *slice_column(Matrix_double *m, size_t x) {
Matrix_double *sliced = copy_matrix(m);
for (size_t row = 0; row < m->rows; row++) {
Array_double *old_row = sliced->data[row];
sliced->data[row] = slice_element(old_row, x);
free_vector(old_row);
}
sliced->cols--;
return sliced;
}
Matrix_double *add_column(Matrix_double *m, Array_double *v) {
Matrix_double *pushed = copy_matrix(m);
for (size_t row = 0; row < m->rows; row++) {
Array_double *old_row = pushed->data[row];
pushed->data[row] = add_element(old_row, v->data[row]);
free_vector(old_row);
}
pushed->cols++;
return pushed;
}
void free_matrix(Matrix_double *m) {
for (size_t y = 0; y < m->rows; ++y)
free_vector(m->data[y]);

15
src/vector.c
View File

@ -88,6 +88,21 @@ Array_double *copy_vector(Array_double *v) {
return copy;
}
Array_double *add_element(Array_double *v, double x) {
Array_double *pushed = InitArrayWithSize(double, v->size + 1, 0.0);
for (size_t i = 0; i < v->size; ++i)
pushed->data[i] = v->data[i];
pushed->data[v->size] = x;
return pushed;
}
Array_double *slice_element(Array_double *v, size_t x) {
Array_double *sliced = InitArrayWithSize(double, v->size - 1, 0.0);
for (size_t i = 0; i < v->size - 1; ++i)
sliced->data[i] = i >= x ? v->data[i + 1] : v->data[i];
return sliced;
}
void free_vector(Array_double *v) {
free(v->data);
free(v);

114
test/matrix.t.c
View File

@ -7,6 +7,38 @@ UTEST(matrix, free) {
EXPECT_NE(data_addr, (uint64_t)(m->data));
}
UTEST(matrix, add_column) {
Matrix_double *m = InitMatrixWithSize(double, 5, 5, 0.0);
Array_double *col = InitArray(double, {1.0, 2.0, 3.0, 4.0, 5.0});
Matrix_double *new_m = add_column(m, col);
for (size_t row = 0; row < m->rows; row++)
EXPECT_EQ(new_m->data[row]->data[m->cols], col->data[row]);
EXPECT_EQ(new_m->cols, m->cols + 1);
free_matrix(m);
free_matrix(new_m);
free_vector(col);
}
UTEST(matrix, slice_column) {
size_t slice = 1;
Matrix_double *m = InitMatrixWithSize(double, 5, 5, 1.0 * (rand() % 10));
Matrix_double *new_m = slice_column(m, slice);
for (size_t row = 0; row < m->rows; row++) {
Array_double *sliced_row = slice_element(m->data[row], slice);
EXPECT_TRUE(vector_equal(new_m->data[row], sliced_row));
free_vector(sliced_row);
}
EXPECT_EQ(new_m->cols, m->cols - 1);
free_matrix(m);
free_matrix(new_m);
}
UTEST(matrix, put_identity_diagonal) {
Matrix_double *m = InitMatrixWithSize(double, 8, 8, 0.0);
Matrix_double *ident = put_identity_diagonal(m);
@ -47,11 +79,53 @@ UTEST(matrix, m_dot_v) {
free_vector(dotted);
}
UTEST(matrix, bsubst) {
Matrix_double *u = InitMatrixWithSize(double, 3, 3, 0.0);
u->data[0]->data[0] = 1.0;
u->data[0]->data[1] = 2.0;
u->data[0]->data[2] = 3.0;
u->data[1]->data[1] = 4.0;
u->data[1]->data[2] = 5.0;
u->data[2]->data[2] = 6.0;
Array_double *b = InitArray(double, {14.0, 29.0, 30.0});
Array_double *solution = bsubst(u, b);
EXPECT_NEAR(solution->data[0], -3.0, 0.0001);
EXPECT_NEAR(solution->data[1], 1.0, 0.0001);
EXPECT_NEAR(solution->data[2], 5.0, 0.0001);
free_matrix(u);
free_vector(b);
free_vector(solution);
}
UTEST(matrix, fsubst) {
Matrix_double *l = InitMatrixWithSize(double, 3, 3, 0.0);
l->data[0]->data[0] = 1.0;
l->data[1]->data[0] = 2.0;
l->data[1]->data[1] = 3.0;
l->data[2]->data[0] = 4.0;
l->data[2]->data[1] = 5.0;
l->data[2]->data[2] = 6.0;
Array_double *b = InitArray(double, {14.0, 13.0, 32.0});
Array_double *solution = fsubst(l, b);
EXPECT_NEAR(solution->data[0], 14.0, 0.0001);
EXPECT_NEAR(solution->data[1], -5.0, 0.0001);
EXPECT_NEAR(solution->data[2], 0.16667, 0.0001);
free_matrix(l);
free_vector(b);
free_vector(solution);
}
UTEST(matrix, lu_decomp) {
Matrix_double *m = InitMatrixWithSize(double, 8, 8, 0.0);
Matrix_double *m = InitMatrixWithSize(double, 10, 10, 0.0);
for (size_t y = 0; y < m->rows; ++y) {
for (size_t x = 0; x < m->cols; ++x)
m->data[y]->data[x] = x == y ? 5.0 : (5.0 - rand() % 10 + 1);
m->data[y]->data[x] = x == y ? 20.0 : (100.0 - rand() % 100) / 100.0;
}
Matrix_double **ul = lu_decomp(m);
@ -75,15 +149,40 @@ UTEST(matrix, lu_decomp) {
free(ul);
}
UTEST(matrix, solve_matrix) {
Matrix_double *m = InitMatrixWithSize(double, 8, 8, 0.0);
UTEST(matrix, solve_gaussian_elimination) {
Matrix_double *m = InitMatrixWithSize(double, 10, 10, 0.0);
for (size_t y = 0; y < m->rows; ++y) {
for (size_t x = 0; x < m->cols; ++x)
m->data[y]->data[x] = x == y ? 10.0 : (5.0 - rand() % 10 + 1);
m->data[y]->data[x] = x == y ? 20.0 : (100.0 - rand() % 100) / 100.0;
}
Array_double *b = InitArray(double, {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0});
Array_double *solution = solve_matrix(m, b);
Array_double *b_1 = InitArrayWithSize(double, m->rows, 1.0);
Array_double *b = m_dot_v(m, b_1);
Array_double *solution = solve_matrix_gaussian(m, b);
for (size_t y = 0; y < m->rows; y++) {
double dot = v_dot_v(m->data[y], solution);
EXPECT_NEAR(b->data[y], dot, 0.0001);
}
free_vector(b_1);
free_matrix(m);
free_vector(b);
free_vector(solution);
}
UTEST(matrix, solve_matrix_lu_bsubst) {
Matrix_double *m = InitMatrixWithSize(double, 10, 10, 0.0);
for (size_t y = 0; y < m->rows; ++y) {
for (size_t x = 0; x < m->cols; ++x)
m->data[y]->data[x] = x == y ? 20.0 : (100.0 - rand() % 100) / 100.0;
}
Array_double *b_1 = InitArrayWithSize(double, m->rows, 1.0);
Array_double *b = m_dot_v(m, b_1);
Array_double *solution = solve_matrix_lu_bsubst(m, b);
for (size_t y = 0; y < m->rows; y++) {
double dot = v_dot_v(m->data[y], solution);
@ -92,6 +191,7 @@ UTEST(matrix, solve_matrix) {
free_matrix(m);
free_vector(b);
free_vector(b_1);
free_vector(solution);
}

22
test/vector.t.c
View File

@ -10,6 +10,28 @@ UTEST(vector, copy_vector) {
free_vector(w);
}
UTEST(vector, add_element) {
Array_double *v = InitArray(double, {3, 1, -4});
Array_double *w = add_element(v, -2);
Array_double *w_expect = InitArray(double, {3, 1, -4, -2});
EXPECT_TRUE(vector_equal(w, w_expect));
free_vector(v);
free_vector(w);
free_vector(w_expect);
}
UTEST(vector, slice_element) {
Array_double *v = InitArray(double, {3, 1, -4});
Array_double *w = slice_element(v, 1);
Array_double *w_expect = InitArray(double, {3, -4});
EXPECT_TRUE(vector_equal(w, w_expect));
free_vector(v);
free_vector(w);
free_vector(w_expect);
}
UTEST(vector, free_vector) {
Array_double *v = InitArray(double, {3, 1, -4});
uint64_t arr_addr = (uint64_t)v->data;