{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "599913a3",
   "metadata": {},
   "source": [
    "# Solving Linear Equations\n",
    "\n",
    "## Serial Algorithm\n",
    "First, we construct a linear equations system $Ax=b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe3b5b02",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 5\n",
    "A = rand(-10.0:10.0, (n, n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57c912fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = rand(-10.0:10.0, n)\n",
    "b = A * x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4baa0681",
   "metadata": {},
   "source": [
    "The code in the following cell converts the problem $Ax=b$ to the upper triangular equation system $Ux=y$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11d255e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "function convert_to_upper_triangular(A,b)\n",
    "    # Upper Triangularization: convert Ax=b to Ux=y\n",
    "    for k in 1:n\n",
    "        for j in k+1:n\n",
    "            # Divide by pivot\n",
    "            A[k,j] = A[k,j] / A[k,k]\n",
    "        end\n",
    "        b[k] = b[k] / A[k,k]\n",
    "        A[k,k] = 1\n",
    "        # Substract lower rows\n",
    "        for i in k+1:n \n",
    "            for j in k+1:n\n",
    "                A[i,j]=A[i,j] - A[i,k] * A[k,j]\n",
    "            end\n",
    "            b[i] = b[i] - A[i,k] * b[k]\n",
    "            A[i,k] = 0\n",
    "        end\n",
    "    end\n",
    "    return A, b #U,y\n",
    "end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02c47593",
   "metadata": {},
   "source": [
    "The function in the following cell solves the upper triangular equation system using backwards substitution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c19497c",
   "metadata": {},
   "outputs": [],
   "source": [
    "function solve_upper_triangular(U,y)\n",
    "    n = size(U,1)\n",
    "    for step in reverse(1:n)\n",
    "        if U[step,step] == 0\n",
    "            if y[step] != 0\n",
    "                return \"No solution\"\n",
    "            else\n",
    "                return \"Infinity solutions\"\n",
    "            end\n",
    "        else\n",
    "        # Backwards substitution\n",
    "            y[step] = y[step] / U[step,step]\n",
    "        end\n",
    "        for row in reverse(1:step-1)\n",
    "            y[row] -= U[row,step] * y[step]\n",
    "        end\n",
    "    end\n",
    "    return y \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c24d85e",
   "metadata": {},
   "outputs": [],
   "source": [
    "U,y = convert_to_upper_triangular(A,b)\n",
    "sol = solve_upper_triangular(U,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c356b5a",
   "metadata": {},
   "source": [
    "We can test if the obtained solution is correct using `@test`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c71828d",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Test\n",
    "@test sol ≈ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe3c5374",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.9.1",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}