{"cells": [{"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["#! /usr/bin/env python3\n", "#\n", "def linpack_bench ( n, verbose = True ):"]}, {"cell_type": "markdown", "metadata": {}, "source": ["****************************************************************************80<br>\n", "<br>\n", " linpack_bench() drives the LINPACK benchmark program.<br>\n", "<br>\n", " Licensing:<br>\n", "<br>\n", "   This code is distributed under the MIT license.<br>\n", "<br>\n", " Modified:<br>\n", "<br>\n", "   06 February 2025<br>\n", "<br>\n", " Author:<br>\n", "<br>\n", "   Original F77 version by Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart.<br>\n", "   This version by John Burkardt.<br>\n", "<br>\n", " Input:<br>\n", "<br>\n", "   integer N, the order of the matrix.<br>\n", ""]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  from time import perf_counter\n", "  import numpy as np\n", "  import platform"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  if ( verbose ):\n", "    print ( '' )\n", "    print ( 'linpack_bench():' )\n", "    print ( '  python version: ' + platform.python_version ( ) )\n", "    print ( '  numpy version:  ' + np.version.version )\n", "    print ( '  The LINPACK benchmark.' )\n", "    print ( '  Datatype: float64' )\n", "    print ( '  Matrix order N = ', n )\n", "#\n", "#  Set the matrix A, a solution x, and a right hand side.\n", "#\n", "  A = 2.0 * np.random.random ( [ n, n ] ) - 1.0\n", "  x_exact = np.ones ( n )\n", "  b = np.matmul ( A, x_exact )\n", "#\n", "#  Factor the matrix.\n", "#\n", "  t = perf_counter ( )\n", "  ALU, ipvt, info = dgefa ( A, n )\n", "  cpu = perf_counter ( ) - t"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  if ( info != 0 ):\n", "    print ( '' )\n", "    print ( 'linpack_bench(): Fatal error!' )\n", "    print ( '  The matrix A is apparently singular.' )\n", "    print ( '  Abnormal end of execution.' )\n", "    raise Exception ( 'linpack_bench(): Fatal error!' )\n", "#\n", "#  Solve the linear system.\n", "#\n", "  t = perf_counter ( )\n", "  x = dgesl ( ALU, n, ipvt, b )\n", "  cpu = cpu + perf_counter ( ) - t\n", "#\n", "#  Compute the residual.\n", "#\n", "  r = np.matmul ( A, x ) - b\n", "#\n", "#  Compute infinity norms.\n", "#\n", "  A_norm = np.linalg.norm ( A, np.inf )\n", "  r_norm = np.linalg.norm ( r, np.inf )\n", "  x_norm = np.linalg.norm ( x, np.inf )\n", "#\n", "#  Report.\n", "#\n", "  eps = np.finfo(float).eps\n", "  ratio = r_norm / A_norm / x_norm / eps"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  ops = ( 2 * n * n * n ) / 3.0 + 2.0 * ( n * n )\n", "  mflops = ops / ( 1000000 * cpu )\n", "  unit = 2.0 / mflops"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  cray = 0.056\n", "  cray_ratio = cpu / cray"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  if ( verbose ):\n", "    print ( \"\" )\n", "    print ( \"  Normalized residual = \", ratio )\n", "    print ( \"  Residual norm       = \", r_norm )\n", "    print ( \"  Machine epsilon     = \", eps )\n", "    print ( \"  First X[]           = \", x[0] )\n", "    print ( \"  Last X[]            = \", x[-1] )\n", "    print ( \"  CPU seconds         = \", cpu )\n", "    print ( \"  MegaFLOPS           = \", mflops )\n", "    print ( \"  Unit                = \", unit )\n", "    print ( \"  Cray ratio          = \", cray_ratio )\n", "#\n", "#  Terminate.\n", "#\n", "  if ( verbose ):\n", "    print ( '' )\n", "    print ( 'linpack_bench():' )\n", "    print ( '  Normal end of execution.' )"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  return mflops"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["def dgefa ( A, n ):"]}, {"cell_type": "markdown", "metadata": {}, "source": ["****************************************************************************80<br>\n", "<br>\n", " dgefa() factors a real matrix in general storage format.<br>\n", "<br>\n", " Discussion:<br>\n", "<br>\n", "   Gaussian elimination with partial pivoting.<br>\n", "<br>\n", " Licensing:<br>\n", "<br>\n", "   This code is distributed under the MIT license.<br>\n", "<br>\n", " Modified:<br>\n", "<br>\n", "   03 June 2024<br>\n", "<br>\n", " Author:<br>\n", "<br>\n", "   Original F77 version by Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart.<br>\n", "   This version by John Burkardt<br>\n", "<br>\n", " Reference:<br>\n", "<br>\n", "   Dongarra, Moler, Bunch and Stewart,<br>\n", "   LINPACK User's Guide,<br>\n", "   SIAM, (Society for Industrial and Applied Mathematics),<br>\n", "   3600 University City Science Center,<br>\n", "   Philadelphia, PA, 19104-2688.<br>\n", "   ISBN 0-89871-172-X<br>\n", "<br>\n", " Input:<br>\n", "<br>\n", "   real A(N,N), the matrix to be factored.<br>\n", "<br>\n", "   integer N, the order of the matrix A.<br>\n", "<br>\n", " Output:<br>\n", "<br>\n", "   real ALU(N,N), an upper triangular matrix and the multipliers <br>\n", "   used to obtain it.  The factorization can be written A=L*U, where L is <br>\n", "   a product of permutation and unit lower triangular matrices, and U is <br>\n", "   upper triangular.<br>\n", "<br>\n", "   integer IPVT(N), the pivot indices.<br>\n", "<br>\n", "   integer INFO, singularity indicator.<br>\n", "   0, normal value.<br>\n", "   K, if U(K-1,K-1) == 0.  This is not an error condition for this subroutine,<br>\n", "   but it does indicate that DGESL or DGEDI will divide by zero if called.<br>\n", "   Use RCOND in DGECO for a reliable indication of singularity.<br>\n", ""]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  import numpy as np"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  ALU = A.copy ( )"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  ipvt = np.zeros ( n, dtype = int )\n", "  info = 0"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  for k in range ( 0, n - 1 ):\n", "#\n", "#  Find L = pivot index,\n", "#  K <= L < N,\n", "#  index of largest entry in column K.\n", "#\n", "    l = -1\n", "    amax = - np.inf\n", "   \n", "    for i in range ( k, n ):\n", "      if ( amax < np.abs ( ALU[i,k] ) ):\n", "        l = i\n", "        amax = np.abs ( ALU[i,k] )\n", "    ipvt[k] = l\n", "#\n", "#  Zero pivot implies this column already triangularized.\n", "#\n", "    if ( ALU[l,k] == 0.0 ):\n", "      info = k + 1\n", "      continue\n", "#\n", "#  Interchange row K and pivot row L if necessary.\n", "#\n", "    if ( l != k ):\n", "      t        = ALU[l,k]\n", "      ALU[l,k] = ALU[k,k]\n", "      ALU[k,k] = t\n", "#\n", "#  Compute multipliers.\n", "#\n", "    ALU[k+1:n,k] = - ALU[k+1:n,k] / ALU[k,k]\n", "#\n", "#  Row elimination with column indexing.\n", "#\n", "    for j in range ( k + 1, n ):\n", "      t = ALU[l,j]\n", "      if ( l != k ):\n", "        ALU[l,j] = ALU[k,j]\n", "        ALU[k,j] = t\n", "      ALU[k+1:n,j] = ALU[k+1:n,j] + t * ALU[k+1:n,k]"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  ipvt[n-1] = n - 1"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  if ( ALU[n-1,n-1] == 0.0 ):\n", "    info = n"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  return ALU, ipvt, info"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["def dgesl ( ALU, n, ipvt, b ):"]}, {"cell_type": "markdown", "metadata": {}, "source": ["****************************************************************************80<br>\n", "<br>\n", " dgesl() solves a real general linear system A * X = B.<br>\n", "<br>\n", " Discussion:<br>\n", "<br>\n", "   The system matrix must have been factored by DGECO or DGEFA.<br>\n", "<br>\n", " Licensing:<br>\n", "<br>\n", "   This code is distributed under the MIT license.<br>\n", "<br>\n", " Modified:<br>\n", "<br>\n", "   03 June 2024<br>\n", "<br>\n", " Author:<br>\n", "<br>\n", "   Original F77 version by Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart.<br>\n", "   This version by John Burkardt.<br>\n", "<br>\n", " Reference:<br>\n", "<br>\n", "   Dongarra, Moler, Bunch and Stewart,<br>\n", "   LINPACK User's Guide,<br>\n", "   SIAM, (Society for Industrial and Applied Mathematics),<br>\n", "   3600 University City Science Center,<br>\n", "   Philadelphia, PA, 19104-2688.<br>\n", "   ISBN 0-89871-172-X<br>\n", "<br>\n", " Input:<br>\n", "<br>\n", "   real ALU(N,N), the LU factorization from DGECO or DGEFA.<br>\n", "<br>\n", "   integer N, the order of the matrix A.<br>\n", "<br>\n", "   integer IPVT(N), the pivot vector from DGECO or DGEFA.<br>\n", "<br>\n", "   real B(N), the right hand side vector.<br>\n", "<br>\n", " Output:<br>\n", "<br>\n", "   real X(N), the solution vector.<br>\n", ""]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  import numpy as np"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  x = b.copy()"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  for k in range ( 0, n - 1 ):\n", "    l = ipvt[k]\n", "    t = x[l]\n", "    if ( l != k ):\n", "      x[l] = x[k]\n", "      x[k] = t\n", "    x[k+1:n] = x[k+1:n] + t * ALU[k+1:n,k]"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  for k in range ( n - 1, -1, -1 ):\n", "    x[k] = x[k] / ALU[k,k]\n", "    t = - x[k]\n", "    x[0:k] = x[0:k] + t * ALU[0:k,k]"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  return x"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["def timestamp ( ):"]}, {"cell_type": "markdown", "metadata": {}, "source": ["****************************************************************************80<br>\n", "<br>\n", " timestamp() prints the date as a timestamp.<br>\n", "<br>\n", " Licensing:<br>\n", "<br>\n", "   This code is distributed under the MIT license. <br>\n", "<br>\n", " Modified:<br>\n", "<br>\n", "   06 April 2013<br>\n", "<br>\n", " Author:<br>\n", "<br>\n", "   John Burkardt<br>\n", ""]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  import time"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  t = time.time ( )\n", "  print ( time.ctime ( t ) )"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["  return"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["if ( __name__ == '__main__' ):\n", "  timestamp ( )\n", "  n = 1000\n", "  linpack_bench ( n )\n", "  timestamp ( )"]}], "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}, "language_info": {"codemirror_mode": {"name": "ipython", "version": 3}, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4"}}, "nbformat": 4, "nbformat_minor": 2}