GEPP.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "xoQ5fei_gMa4"
      },
      "outputs": [],
      "source": [
        "import numpy as np"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "##### Utility"
      ],
      "metadata": {
        "id": "Wl0RuoYogU3T"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def swap(B,swaps):\n",
        "  A=np.copy(B)\n",
        "  n=len(swaps)\n",
        "  for i in range(n):\n",
        "    A[[i,swaps[i]],:]=A[[swaps[i],i],:]\n",
        "  return A"
      ],
      "metadata": {
        "id": "VCkP7remgRGw"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def split(A):\n",
        "  U=np.zeros((len(A),len(A)))\n",
        "  L=np.eye(len(A))\n",
        "  for i in range(len(A)):\n",
        "    U[i,i:]=A[i,i:]\n",
        "    L[i+1:,i]=A[i+1:,i]\n",
        "  return L, U"
      ],
      "metadata": {
        "id": "vas3C82dgRET"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def myLU(B):\n",
        "  A = np.copy(B)\n",
        "  p=np.zeros(len(A)-1).astype(np.intc) #indices are int\n",
        "  for i in range(len(A)):\n",
        "    row = np.argmax(np.absolute(A[i:,i])) + i #find the max row\n",
        "    if row != i:\n",
        "      A[[i,row],:]=A[[row,i],:] # swap the two rows\n",
        "    if i < len(p):\n",
        "      p[i] = row #populate p vector as in Problem 1\n",
        "    A[(i+1):,i]=A[(i+1):,i]/A[i,i]\n",
        "    A[(i+1):,(i+1):]-=np.outer(A[(i+1):,i], A[i,(i+1):])\n",
        "  L, U = split(A) #use the split method from Lab1 to split A into L, U\n",
        "  return L, U, p"
      ],
      "metadata": {
        "id": "BqCBKmzcgRCN"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Perform Row-Oriented Forward-Substitution for Lower-Tri System Ax=b\n",
        "# Single for-loop\n",
        "def lowerfSub(A,b):\n",
        "  n=len(A) # get the dimension of A in 0-indexing\n",
        "  x=np.zeros(n)\n",
        "  x[0] = b[0]/A[0][0] # the first entry\n",
        "  for i in range(1,n):\n",
        "    x[i]=(b[i]-np.dot(A[i,:], x))/A[i][i]\n",
        "  return x"
      ],
      "metadata": {
        "id": "egsWknzJgQ_7"
      },
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Example of Row-Oriented Back-Substituion to solve Upper-Tri System Ax=b\n",
        "def upperbSub(A,b):\n",
        "  n=len(A)\n",
        "  x=np.zeros(n)\n",
        "  x[n-1]=b[n-1]/A[n-1,n-1]\n",
        "  for i in range(1,n):\n",
        "    x[n-1-i]=1/A[n-1-i,n-1-i]*(b[n-1-i]-np.dot(A[n-1-i,n-i:],x[n-i:]))\n",
        "  return x"
      ],
      "metadata": {
        "id": "cDNcJ6sQgQ7w"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Implemented swap() to permute a vector\n",
        "def swapfV(b,swaps):\n",
        "  a=np.copy(b)\n",
        "  n=len(swaps)\n",
        "  for i in range(n):\n",
        "    if swaps[i]!=i:\n",
        "      temp = a[i]\n",
        "      a[i]=a[swaps[i]]\n",
        "      a[swaps[i]] = temp\n",
        "  return a"
      ],
      "metadata": {
        "id": "YJT8rT9MgtW-"
      },
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### LU-PA Decomposition"
      ],
      "metadata": {
        "id": "__3UGXorLhn9"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def myLU(B):\n",
        "  A = np.copy(B)\n",
        "  p=np.zeros(len(A)-1).astype(np.intc) #indices are int\n",
        "  for i in range(len(A)):\n",
        "    row = np.argmax(np.absolute(A[i:,i])) + i #find the max row\n",
        "    if row != i:\n",
        "      A[[i,row],:]=A[[row,i],:] # swap the two rows\n",
        "    if i < len(p):\n",
        "      p[i] = row #populate p vector as in Problem 1\n",
        "    A[(i+1):,i]=A[(i+1):,i]/A[i,i]\n",
        "    A[(i+1):,(i+1):]-=np.outer(A[(i+1):,i], A[i,(i+1):])\n",
        "  L, U = split(A) #use the split method from Lab1 to split A into L, U\n",
        "  return L, U, p"
      ],
      "metadata": {
        "id": "aiOc-xF6gtUs"
      },
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Gaussian Elimination with Partial Pivoting"
      ],
      "metadata": {
        "id": "sK3U_CybLkcq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def mySolve(A, b):\n",
        "  x = x=np.zeros(len(b))\n",
        "  L, U, p = myLU(A)\n",
        "  newB = swapfV(b, p)\n",
        "  y = lowerfSub(L, newB)\n",
        "  x = upperbSub(U, y)\n",
        "  return x"
      ],
      "metadata": {
        "id": "FTji9EMsgtSb"
      },
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Testing"
      ],
      "metadata": {
        "id": "ygV-EQw5Lpm8"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# LU-PA Decomposition\n",
        "for i in range(10):\n",
        "  A=np.random.rand(20,20)\n",
        "  L, U, p= myLU(A)\n",
        "  print(np.allclose(L@U, swap(A, p)))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "dQ6YXfBygtQZ",
        "outputId": "57099172-2f74-4f06-afed-f52bad101487"
      },
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "True\n",
            "True\n",
            "True\n",
            "True\n",
            "True\n",
            "True\n",
            "True\n",
            "True\n",
            "True\n",
            "True\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# GEPP\n",
        "for i in range(10):\n",
        "  A=np.random.rand(20,20)\n",
        "  b=np.random.rand(20)\n",
        "  x=mySolve(A,b)\n",
        "  print(np.allclose(A@x,b)) # test if the result is correct\n",
        "  print(np.max(np.abs(np.linalg.solve(A,b)-mySolve(A,b)))) # test accuracy"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "r0XYybV0Lzi8",
        "outputId": "ba954213-da0c-41d5-be92-c71f61fd02ab"
      },
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "True\n",
            "3.9968028886505635e-15\n",
            "True\n",
            "7.105427357601002e-15\n",
            "True\n",
            "2.4424906541753444e-15\n",
            "True\n",
            "1.5543122344752192e-15\n",
            "True\n",
            "2.886579864025407e-15\n",
            "True\n",
            "1.1102230246251565e-15\n",
            "True\n",
            "5.10702591327572e-15\n",
            "True\n",
            "5.773159728050814e-15\n",
            "True\n",
            "5.10702591327572e-15\n",
            "True\n",
            "8.881784197001252e-16\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Note\n",
        "GEPP(myLU()) is more stable compared to outerGEWP() on large-scale linear systems. "
      ],
      "metadata": {
        "id": "Tqd0RexIMDfg"
      }
    }
  ]
}