- Commit
- c3244248e2ea438014ad0173abe3700567a03ed4
- Parent
- 7bc7ca3114c0e9e4e103d21fba2886e33177b47b
- Author
- Pablo <pablo-pie@riseup.net>
- Date
Created a program to benchmark the getrf LAPACK function
numeric-linalg
Educational material on the SciPy implementation of numerical linear algebra algorithms
Diffstat
7 files changed, 316 insertions, 0 deletions
| Status | File Name | N° Changes | Insertions | Deletions |
| Added | getrf/benchmark/.gitignore | 2 | 2 | 0 |
| Added | getrf/benchmark/benchmark.ipynb | 86 | 86 | 0 |
| Added | getrf/benchmark/build.sh | 10 | 10 | 0 |
| Added | getrf/benchmark/config.h | 19 | 19 | 0 |
| Added | getrf/benchmark/histogram.bin | 0 | 0 | 0 |
| Added | getrf/benchmark/main.c | 163 | 163 | 0 |
| Added | getrf/benchmark/progress-bar.h | 36 | 36 | 0 |
diff --git a/getrf/benchmark/benchmark.ipynb b/getrf/benchmark/benchmark.ipynb @@ -0,0 +1,86 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "9decf754-a174-4338-9c54-5eef69cc91bb", + "metadata": {}, + "source": [ + "To test the complexity of the `getrf` algorithm we plot the log of the excution time on progressively larger inputs: `getrf` is $O(n^d)$, where $d$ is the slope of the result line in the plot. We can thus confirm that `getrf` is $O(n^3)$!" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "30318f00-d2d5-444f-9e8a-abd78af4b74f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[<matplotlib.lines.Line2D at 0x7f8412999580>]" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAGdCAYAAADT1TPdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLIklEQVR4nO3de1xVVcLG8R8gHhEBUQLUhNRMSx3vmbfE62iNZU1WY5M2U41OYpo2yE2RkIvZaGlp2aSVvb068+pMt5nyQuDd1LTMNDNNNMW74A0Q2O8fa5QsNS8c9jmH5/v5nE/txQGeOsh5XHvvtbwsy7IQERERcVHedgcQERERuRyVFREREXFpKisiIiLi0lRWRERExKWprIiIiIhLU1kRERERl6ayIiIiIi5NZUVERERcWhW7A1yv0tJS9u3bR0BAAF5eXnbHERERkStgWRYnTpygbt26eHtffu7E7cvKvn37qF+/vt0xRERE5Brs2bOHG2+88bLPcfuyEhAQAJj/2MDAQJvTiIiIyJXIz8+nfv3659/HL8fty8q5Uz+BgYEqKyIiIm7mSi7hcOoFtunp6bRv356AgABCQ0MZMGAA33zzzQXPeeyxx/Dy8rrgcccddzgzloiIiLgRp5aV7Oxshg8fzpo1a1i8eDHFxcX06dOHU6dOXfC8vn37sn///vOPf//7386MJSIiIm7EqaeBPv744wuO58yZQ2hoKBs2bODOO+88P+5wOAgPD3dmFBEREXFTFbrOSl5eHgC1atW6YDwrK4vQ0FBuueUWnnzySQ4ePHjJr1FYWEh+fv4FDxEREfFcXpZlWRXxjSzL4t577+XYsWMsX778/Pj8+fOpUaMGkZGR7Nq1i3HjxlFcXMyGDRtwOBw/+zoTJkwgOTn5Z+N5eXm6wFZERMRN5OfnExQUdEXv3xVWVoYPH85HH33EihUrLns/9f79+4mMjGTevHncf//9P/t4YWEhhYWF54/P3fqksiIiIuI+rqasVMityyNGjOD9999n2bJlv7jwS506dYiMjOTbb7+96McdDsdFZ1xERETEMzm1rFiWxYgRI/jnP/9JVlYWDRo0+MXPOXLkCHv27KFOnTrOjCYiIiJuwqkX2A4fPpx33nmHd999l4CAAHJzc8nNzeXMmTMAnDx5kmeffZbVq1fz/fffk5WVRf/+/QkJCeG+++5zZjQRERFxE069ZuVSq9LNmTOHxx57jDNnzjBgwAA2btzI8ePHqVOnDt27dyclJeWK9/u5mnNeIiIi4hpc5pqVX+pBfn5+fPLJJ86MICIiIm6uQtdZEREREblaKisiIiJupNQq5d3N7zJg3gBKrVK741QIt991WUREpDKwLIuPd3xM3NI4vjjwBQDzvprHoBaDbE7mfCorIiIiLm7N3jXELokle3c2AIGOQGI6xXBPk3tsTlYxVFZERERc1LbD24hfGs8/t/0TgKo+VYluH01c1zhCqofYnK7iqKyIiIi4mL35e5mQNYE5m+ZQapXi7eXN4JaDSY5KJiIowu54FU5lRURExEUcPXOUjBUZTP9sOgXFBQDc0+Qe0nqk0Sy0mc3p7KOyIiIiYrPTZ08zbe00Jq2cxPGC4wB0iehCRs8MOkd0tjecC1BZERERscnZkrPM2TSH5Oxk9p3YB0Dz0Oak90zn7sZ3X3Il+MpGZUVERKSCWZbFgq0LSMhMYPuR7QBEBkWS0j2FQS0G4ePtY3NC16KyIiIiUoGW7lxK7NJY1u9bD0BI9RASuyYyrN0wHFUcNqdzTSorIiIiFeDz/Z8TuySWxTsXA+Dv68+YjmMY02kMgQ5txHs5KisiIiJOtOPoDhIzE5m/ZT4Avt6+DG07lMQ7EwmrEWZzOvegsiIiIuIEuSdzeS77OV7//HWKS4sBGNRiECndU2gY3NDmdO5FZUVERKQc5RXkMXnVZKaumcrps6cB6HtzX9J7ptMqvJW94dyUyoqIiEg5KCguYMa6GaQtT+PImSMA3F7vdib1mkTUTVH2hnNzKisiIiLXoaS0hLe/eJukrCT25O8BoEntJqT1TOO+pvdprZRyoLIiIiJyDSzL4v1v3ic+M56vD30NQL2AeiRHJTOk1RCqeOsttrzo/6SIiMhVWr57ObFLY1m1ZxUAwdWCiesSR/Tt0fj5+tmczvOorIiIiFyhzQc2E7c0jo++/QgAvyp+jOwwkpjOMQT7BducznOprIiIiPyC749/z/hPx/POl+9gYeHj5cPjrR8nKSqJugF17Y7n8VRWRERELuHQqUOkLk9l5vqZFJUUAfDAbQ8wsftEmoQ0sTld5aGyIiIi8hMni04yZfUUXlj1AieKTgDQo0EPMnpm0L5ee5vTVT4qKyIiIv9VVFLErA2zSFmWwsFTBwFoHd6ajF4Z9G7YW7ch20RlRUREKr1Sq5R5X81j3Kfj2HlsJwCNghsxscdEHmz2IN5e3jYnrNxUVkREpNKyLIuPd3xM3NI4vjjwBQBh/mEkdUviiTZP4Ovja3NCAZUVERGppNbsXUPskliyd2cDEOgIJKZTDKPuGIV/VX+b08mPqayIiEilsu3wNuKXxvPPbf8EoKpPVaLbRxPXNY6Q6iE2p5OLUVkREZFKYW/+XiZkTWDOpjmUWqV4e3kzuOVgkqOSiQiKsDueXIbKioiIeLSjZ46SsSKD6Z9Np6C4AIB7mtxDWo80moU2szmdXAmVFRER8Uinz55m2tppZKzIIK8wD4AuEV3I6JlB54jONqeTq6GyIiIiHuVsyVnmbJpDcnYy+07sA6B5aHMyemZwV+O7tFaKG1JZERERj2BZFgu2LiAhM4HtR7YDEBkUSUr3FAa1GISPt4/NCeVaqayIiIjbW7pzKbFLY1m/bz0AIdVDSOyayLB2w3BUcdicTq6XyoqIiLitz/d/TuySWBbvXAyAv68/YzqOYUynMQQ6Am1OJ+VFZUVERNzOjqM7SMxMZP6W+QD4evsytO1QEu9MJKxGmM3ppLw5dbOD9PR02rdvT0BAAKGhoQwYMIBvvvnmgudYlsWECROoW7cufn5+REVFsWXLFmfGEhERN5V7MpenPnqKW1+59XxRGdRiENuitzH9rukqKh7KqWUlOzub4cOHs2bNGhYvXkxxcTF9+vTh1KlT55/z/PPPM2XKFF5++WXWrVtHeHg4vXv35sSJE86MJiIibiSvII/EzEQaTWvEzPUzKS4tpu/Nfdk4dCP/c///0DC4od0RxYm8LMuyKuqbHTp0iNDQULKzs7nzzjuxLIu6desyatQoxo4dC0BhYSFhYWFMmjSJoUOH/uLXzM/PJygoiLy8PAIDdX5SRMSTFBQXMGPdDNKWp3HkzBEAOtTrQEavDKJuirI3nFyXq3n/rtBrVvLyzKI8tWrVAmDXrl3k5ubSp0+f889xOBx069aNVatWXbSsFBYWUlhYeP44Pz/fyalFRKSilZSWMPfLuSRlJZGTlwNA05CmpPVIY0DTAVorpZKpsLJiWRajR4+mS5cuNG/eHIDc3FwAwsIuPMcYFhbG7t27L/p10tPTSU5Odm5YERGxhWVZfLD9A+KXxrPlkLl+sV5APZKjkhnSaghVvHVfSGVUYa96dHQ0X375JStWrPjZx37akC3LumRrjouLY/To0eeP8/PzqV+/fvmGFRGRCrciZwWxS2JZuWclAMHVgonrEkf07dH4+frZnE7sVCFlZcSIEbz//vssW7aMG2+88fx4eHg4YGZY6tSpc3784MGDP5ttOcfhcOBwaIEfERFPsfnAZuIz4/lw+4cA+FXxY2SHkYztMpaa1WraG05cglPvBrIsi+joaBYuXEhmZiYNGjS44OMNGjQgPDycxYsXnx8rKioiOzubTp06OTOaiIjYbPfx3Qz51xBavtqSD7d/iI+XD39q8yd2PL2D9F7pKipynlNnVoYPH867777Le++9R0BAwPlrVIKCgvDz88PLy4tRo0aRlpZG48aNady4MWlpaVSvXp1BgwY5M5qIiNjk0KlDpC5PZeb6mRSVFAHwwG0PMLH7RJqENLE5nbgip5aVmTNnAhAVFXXB+Jw5c3jssccAiImJ4cyZMzz11FMcO3aMDh06sGjRIgICApwZTUREKtjJopNMWT2FF1a9wIkis5ZWjwY9yOiZQft67W1OJ66sQtdZcQatsyIi4tqKSoqYtWEWKctSOHjqIABt6rQho2cGvRr20m3IlZTLrrMiIiKVR6lVyryv5jHu03HsPLYTgEbBjUjtkcrAZgPx9nLqZZPiQVRWRESkXFmWxSfffULc0jg25W4CIMw/jKRuSTzR5gl8fXztDShuR2VFRETKzdq9a4ldGkvW91kABDoCiekUw6g7RuFf1d/ecOK2VFZEROS6bTu8jYTMBBZuXQhAVZ+qRLePJr5rPLWr17Y5nbg7lRUREblmP+T/wISsCczeNJtSqxRvL28GtxxMclQyEUERdscTD6GyIiIiV+3YmWNkrMhg2mfTKCguAODeJveS2iOVZqHNbE4nnkZlRURErtjps6eZvnY6GSszOF5wHICuEV3J6JVBp/paeVycQ2VFRER+UXFpMXM2zmFC9gT2ndgHQIvQFqT3TOeuxndprRRxKpUVERG5JMuyWLh1IfGZ8Ww/sh2AyKBIUrqnMKjFIHy8fWxOKJWByoqIiFxU5q5MYpfEsm7fOgBCqoeQ2DWRYe2G4ajisDmdVCYqKyIicoGN+zcSuzSWRd8tAsDf158xHccwptMYAh3a1kQqnsqKiIgA8N3R70j8NJF5X80DwNfbl6Fth5J4ZyJhNcJsTieVmcqKiEgll3syl5TsFGZ9Povi0mIABrUYREr3FBoGN7Q5nYjKiohIpZVfmM/klZOZsmYKp8+eBqDfzf1I65lGq/BW9oYT+RGVFRGRSqaguICZ62aSujyVI2eOANChXgcm9ZpEt5u62ZxO5OdUVkREKomS0hLmfjmXpKwkcvJyAGga0pS0HmkMaDpAa6WIy1JZERHxcJZl8cH2D4hfGs+WQ1sAqBdQj+SoZIa0GkIVb70ViGvTT6iIiAdbkbOC2CWxrNyzEoDgasHEdYkj+vZo/Hz9bE4ncmVUVkREPNDmA5uJz4znw+0fAuBXxY+RHUYytstYalaraW84kauksiIi4kF2H9/N+KzxzP1iLhYWPl4+PN76cZKikqgbUNfueCLXRGVFRMQDHD59mNRlqcxYP4OikiIABt42kJTuKTQJaWJzOpHro7IiIuLGThadZOrqqUxeNZkTRScA6NmgJ+k902lfr73N6UTKh8qKiIgbKiop4vUNr5OyLIUDpw4A0KZOGzJ6ZtC7UW+b04mUL5UVERE3UmqVMu+reYz7dBw7j+0E4OZaNzOx+0QGNhuIt5e3zQlFyp/KioiIG7Asi0+++4S4pXFsyt0EQHiNcMbfOZ4n2jyBr4+vvQFFnEhlRUTExa3du5bYpbFkfZ8FQKAjkJhOMYy6YxT+Vf3tDSdSAVRWRERc1LbD20jITGDh1oUAOHwcRN8eTVyXOGpXr21zOpGKo7IiIuJifsj/gQlZE5i9aTalVineXt4MaTmECVETiAiKsDueSIVTWRERcRHHzhwjY0UG0z6bRkFxAQD3NrmX1B6pNAttZnM6EfuorIiI2Oz02dNMWzuNSSsncbzgOABdI7qS0SuDTvU72RtOxAWorIiI2KS4tJjZG2eTnJ3MvhP7AGgR2oL0nunc1fguvLy8bE4o4hpUVkREKphlWSzYuoCEzAS2H9kOQGRQJCndUxjUYhA+3j42JxRxLSorIiIVKHNXJrFLYlm3bx0AIdVDSOyayLB2w3BUcdicTsQ1qayIiFSAjfs3Ers0lkXfLQLA39efZzs9y5iOYwhwBNicTsS1qayIiDjRjqM7GPfpOOZ9NQ8AX29fhrUbRkLXBMJqhNmcTsQ9qKyIiDhB7slcUrJTmPX5LIpLi/HCi0EtBvFc9+doGNzQ7ngibsWpO14tW7aM/v37U7duXby8vPjXv/51wccfe+wxvLy8LnjccccdzowkIuJUeQV5jMscR6NpjZixfgbFpcX0u7kfG4du5J3731FREbkGTp1ZOXXqFC1btuQPf/gDv/3tby/6nL59+zJnzpzzx1WrVnVmJBERpygoLmDGuhmkLU/jyJkjAHSo14FJvSbR7aZuNqcTcW9OLSv9+vWjX79+l32Ow+EgPDzcmTFERJympLSEuV/OJSkriZy8HACahjQlrUcaA5oO0FopIuXA9mtWsrKyCA0NpWbNmnTr1o3U1FRCQ0PtjiUiclmWZfHB9g+IXxrPlkNbAKgXUI/kqGSGtBpCFW/bf72KeAxb/zT169ePgQMHEhkZya5duxg3bhw9evRgw4YNOBwXX2+gsLCQwsLC88f5+fkVFVdEBIAVOSuIXRLLyj0rAQiuFkx813iGtx+On6+fzelEPI+tZeWhhx46/+/NmzenXbt2REZG8tFHH3H//fdf9HPS09NJTk6uqIgiIudtPrCZ+Mx4Ptz+IQB+VfwY2WEkY7uMpWa1mvaGE/FgLjVPWadOHSIjI/n2228v+Zy4uDhGjx59/jg/P5/69etXRDwRqaR2H9/N+KzxzP1iLhYWPl4+PNHmCcZ3G0/dgLp2xxPxeC5VVo4cOcKePXuoU6fOJZ/jcDgueYpIRKQ8HT59mNRlqcxYP4OikiIABt42kIk9JnJL7VtsTidSeTi1rJw8eZIdO3acP961axebNm2iVq1a1KpViwkTJvDb3/6WOnXq8P333xMfH09ISAj33XefM2OJiFzWyaKTTF09lcmrJnOi6AQAPRv0JL1nOu3rtbc5nUjl49Sysn79erp3737++NzpmyFDhjBz5kw2b97M22+/zfHjx6lTpw7du3dn/vz5BARonwwRqXhFJUXM2jCLlGUpHDx1EIA2ddqQ0TOD3o1625xOpPLysizLsjvE9cjPzycoKIi8vDwCAwPtjiMibqjUKmXeV/MY9+k4dh7bCcDNtW5mYveJDGw2EG8vpy72LVIpXc37t0tdsyIiUpEsy+KT7z4hbmkcm3I3ARBeI5zxd47niTZP4Ovja29AEQFUVkSkklq7dy2xS2PJ+j4LgEBHIDGdYhh1xyj8q/rbG05ELqCyIiKVyrbD20jITGDh1oUAOHwcRN8eTVyXOGpXr21zOhG5GJUVEakU9ubvJTkrmdmbZlNqleLt5c2QlkOYEDWBiKAIu+OJyGWorIiIRzt65igZKzKY/tl0CooLALi3yb2k9kilWWgzm9OJyJVQWRERj3T67GmmrZ3GpJWTOF5wHICuEV3J6JVBp/qd7A0nIldFZUVEPEpxaTGzN84mOTuZfSf2AdAitAXpPdO5q/FdeHl52ZxQRK6WyoqIeATLsliwdQEJmQlsP7IdgMigSFK6pzCoxSB8vH1sTigi10plRUTcXuauTGKXxLJu3zoAQqqHkNg1kWHthuGoor3ERNydyoqIuK3P939O3NI4Fn23CAB/X3/GdBzDmE5jCHRoRWsRT6GyIiJuZ8fRHSRmJjJ/y3wAfL19Gdp2KIl3JhJWI8zmdCLup7AQfH3B20V3lnDRWCIiP5d7MpfhHw3n1lduPV9UBrUYxLbobUy/a7qKishVOnUKRo+GGjXMY9gwKC62O9XPaWZFRFxeXkEek1dNZuqaqZw+exqAvjf3Jb1nOq3CW9kbTsQNnT0LGRkwdSocO2bGiovhtddg+XJo3RqmTIHQUHtznqOyIiIuq6C4gBnrZpC2PI0jZ44AcHu925nUaxJRN0XZG07ETWVnm9mUzz83xw0bwksvmVNBgwbB11+bxxdfQFYW1HaBXShUVkTE5ZSUljD3y7kkZSWRk5cDQJPaTUjrmcZ9Te/TWiki1+CHH2DECPjnP81xcDBMnw4PPww+/72zf8sWWL0axo6Fr76CAQNg8WKoVs222IDKioi4EMuy+GD7B8QvjWfLoS0A1AuoR3JUMkNaDaGKt35liVyt3FwYPx7eftvMnlSpAn/6EyQmQp06Fz735pvNo21b6NQJVqyAtDR47jl7sp/jZVmWZW+E65Ofn09QUBB5eXkEBupWRRF3tSJnBbFLYlm5ZyUAwdWCiesSR/Tt0fj5+tmcTsT9lJTAW2/BX/4CR4+asU6dYOZM+NWvfvnz//53eOghCAqC3bvNP8vT1bx/668pImKrzQc2E58Zz4fbPwTAr4ofIzuMJKZzDMF+wTanE3FPa9ea2ZMvvzTHrVub61K6dIErPYv6wANw223m+pVXXoH4eOfl/SUqKyJii++Pf8/4T8fzzpfvYGHh4+XD460fJykqiboBde2OJ+KWjh6FiRNNMSkthZo1TckYNcqso3I1vL0hIQF+/3szs2InlRURqVCHTh0idXkqM9fPpKikCIAHbnuAid0n0iSkic3pRNzTmTPmYtn0dDh+3Iw98ogpLddzN8+DD5pZmVtvLZeY10xlRUQqxMmik0xZPYUXVr3AiaITAPRo0IOMnhm0r9fe5nQi7qmgAN54w5SUH34wYy1awPPPQ9++1//1q1Sxv6iAyoqIOFlRSRGzNswiZVkKB08dBKB1eGsyemXQu2Fv3YYsco3WrYNHH4VvvjHH9etDSoo5bePjYZuMq6yIiFOUWqXM+2oe4z4dx85jOwFoFNyIiT0m8mCzB/H20m4fIteiuNjMpCQnmzt+wsPNrcl//CM4PHSTcZUVESlXlmXxyXefELc0jk25mwAI8w9jfLfxPNHmCar6VLU3oIgb+/Zb+MMfYKW5w5+HHoIZM6BWLXtzOZvKioiUm7V71xK7NJas77MACHQEEtMphpF3jKRG1Rr2hhNxY/v3w9/+ZhZoKyiAgACzXsojj9idrGKorIjIddt2eBsJmQks3LoQgKo+VYluH01c1zhCqofYnE7EfR0/bmZS3nsPzi3h2rs3zJoFN91kZ7KKpbIiItdsb/5ekrOSmb1pNqVWKd5e3gxuOZjkqGQigiLsjifi1g4fhl//umzDwY4dYfhws9lgZbsuXWVFRK7a0TNHyViRwfTPplNQXADAPU3uIa1HGs1Cm9mcTsT9bd0Kv/2t+ecNN8B//mP266msVFZE5IqdPnuaaWunMWnlJI4XHAegS0QXMnpm0Dmis73hRDzE22/Dn/8Mp09D3bqwZIlrrHViJ5UVEflFxaXFzN44m+TsZPad2AdA89DmpPdM5+7Gd2utFJFyUFRklsWfOdMc9+gB77zz852RKyOVFRG5JMuyWLB1AQmZCWw/sh2AyKBInuv+HI+0eAQfbw9beUrEJnv3wsMPm1uSvbzMuinjxnne4m7XSmVFRC4qc1cmsUtiWbdvHQAh1UNI6JrAn9v9GUcVD115SsQG//gHDB0Kx45BYCD8z//Ab35jdyrXorIiIhfYuH8jsUtjWfTdIgD8ff0Z03EMYzqNIdARaHM6Ec+xfz888wzMn2+O27c3RaVxY3tzuSKVFREBYMfRHYz7dBzzvpoHgK+3L0PbDiXxzkTCaoTZnE7Ec5SUwGuvQXw85OWBtzfExUFSEvj62p3ONamsiFRyuSdzSclOYdbnsyguLQZgUItBpHRPoWFwQ5vTiXiW7Gx4+mn48ktz3K6dKS5t2tiby9WprIhUUnkFeUxeNZmpa6Zy+uxpAPre3Jf0num0Cm9lbzgRD7Ntm5k9+de/zHFwMDz3nLlFWRfR/jKVFZFKpqC4gBnrZpC2PI0jZ44AcHu925nUaxJRN0XZG07Ew+zaZXZInj3bnP7x9oY//QlSUiBEO1FcMafu0b5s2TL69+9P3bp18fLy4l/nKuV/WZbFhAkTqFu3Ln5+fkRFRbFlyxZnRhKptEpKS3hz05s0ebkJYxaN4ciZIzSp3YQFDy5gzeNrVFREytGOHTB4sLlY9vXXTVHp3x82bzbrqKioXB2nlpVTp07RsmVLXn755Yt+/Pnnn2fKlCm8/PLLrFu3jvDwcHr37s2JEyecGUukUrEsi/e/eZ+Wr7bkD+/9gZy8HOoF1OP1/q/z1VNfcf+t92tRN5FyUlhoTu80bw5z55qS8utfw/Ll8P77cNttdid0T049DdSvXz/69et30Y9ZlsWLL75IQkIC999/PwBvvfUWYWFhvPvuuwwdOtSZ0UQqhRU5K4hdEsvKPSsBCK4WTFyXOKJvj8bP18/mdCKe5ZNPzMWz2836ifTuDWlp5iJauT62XbOya9cucnNz6dOnz/kxh8NBt27dWLVq1SXLSmFhIYWFheeP8/PznZ5VxN1sPrCZ+Mx4Ptz+IQB+VfwY2WEkMZ1jCPYLtjmdiGf5+msYMwY+/tgch4fDSy/BwIGVb3dkZ7GtrOTm5gIQFnbh+g1hYWHs3r37kp+Xnp5OcnKyU7OJuKvdx3czPms8c7+Yi4WFj5cPj7d+nKSoJOoG1LU7nohH+eEHc8rnjTfM6R5fXxgxwiyVHxRkdzrPYvvdQD89V25Z1mXPn8fFxTF69Ojzx/n5+dSvX99p+UTcweHTh0ldlsqM9TMoKikC4IHbHmBi94k0CWliczoRz3LsGEyaZGZPCgrM2H33mTGtPusctpWV8PBwwMyw1PnRlpIHDx782WzLjzkcDhwO7UsiAnCy6CRTV09l8qrJnCgyF6b3aNCDjJ4ZtK/X3uZ0Ip7lzBmYNg0yMuD4cTPWubO5NblrV1ujeTzbykqDBg0IDw9n8eLFtG7dGoCioiKys7OZNGmSXbFE3EJRSRGzNswiZVkKB08dBKB1eGsyemXQu2Fv3d0jUo5KSmDOHLMc/r59Zqx5c3Px7G9+o+tSKoJTy8rJkyfZsWPH+eNdu3axadMmatWqRUREBKNGjSItLY3GjRvTuHFj0tLSqF69OoMGDXJmLBG3VWqVMu+reYz7dBw7j+0EoFFwIyb2mMiDzR7E28upqxGIVDrLlpk7fL74whxHRJgF3R55RCvPViSnlpX169fTvXv388fnrjUZMmQIb775JjExMZw5c4annnqKY8eO0aFDBxYtWkRAQIAzY4m4Hcuy+OS7T4hbGsem3E0AhPmHkdQtiSfaPIGvj3Y/EylPu3dDTAz8/e/muGZNGDcOnnoKqlWzNVql5GVZlmV3iOuRn59PUFAQeXl5BAZq+3rxPGv3riV2aSxZ32cBEOgIJKZTDKPuGIV/VX97w4l4mIMHzTUoM2ZAUZGWx3emq3n/tv1uIBG5uG2Ht5GQmcDCrQsBqOpTlej20cR1jSOkun5ripSnvDz4619h6lQ4edKMde9ujlu2tDebqKyIuJwf8n9gQtYEZm+aTalVireXN4NbDiY5KpmIoAi744l4lPx8c4fPlCnmlmQwK86mpUGvXrp41lWorIi4iGNnjpGxIoNpn02joNgs3nBPk3tI65FGs9BmNqcT8SwnTsDLL8MLL8DRo2bs1lth4kSzZopKimtRWRGx2emzp5m+djoZKzM4XnAcgC4RXcjomUHniM72hhPxMAUF8OqrkJoKhw+bsSZNzG3JDz6oO3xclcqKiE2KS4uZs3EOE7InsO+EWbyheWhzMnpmcFfju7RWikg5KikxuyAnJUFOjhm7+WZz/LvfqaS4OpUVkQpmWRYLti4gITOB7UfM9qyRQZGkdE9hUItB+Hjrt6ZIebEs+OADiIszGw4C1KsHEybAY49BFb0LugW9TCIVKHNXJrFLYlm3bx0AIdVDSOyayLB2w3BU0TYSIuVp3Tp49lmzsBtAcDDEx8Pw4eDnZ282uToqKyIVYOP+jcQujWXRd4sA8Pf1Z0zHMYzpNIZAh9YHEilPu3aZUjJvnjmuVg1GjYKxY83ibuJ+VFZEnGjH0R2M+3Qc874yvzV9vX0Z2nYoiXcmElbj0ht2isjVO3DAbDJ4bkE3Ly8YPNgs6Fa/vt3p5HqorIg4Qe7JXFKyU5j1+SyKS4sBGNRiECndU2gY3NDmdCKe5cgRmDwZpk+H06fNWK9eZqxVK1ujSTlRWREpR3kFebyw6gWmrJnC6bPmt2bfm/uS3jOdVuGt7A0n4mHy8swKs1OmmHVTAG6/3cyk9O6ttVI8icqKSDkoKC5gxroZpC1P48iZIwB0qNeBjF4ZRN0UZW84EQ9TWGhO9UycWLag269+ZUpK//4qKZ5IZUXkOpSUljD3y7kkZSWRk2cWb2ga0pS0HmkMaDpAa6WIlCPLgv/7P4iNhZ07zVjTppCcDA88YDYdFM+ksiJyDSzL4oPtHxC/NJ4th7YAUC+gHslRyQxpNYQq3vqjJVKeVq0ytyGvXm2Ow8PNTIrWSqkc9BKLXKUVOSuIXRLLyj0rAQiuFkxclziib4/Gz1eLN4iUp+3bITER/vEPc1y9OvzlL6a41KhhbzapOCorIldo84HNxGfG8+H2DwHwq+LHyA4jGdtlLDWr1bQ3nIiH2bHDzJy88w6UlppTPH/8oznlU7eu3emkoqmsiPyC3cd3Mz5rPHO/mIuFhY+XD4+3fpykqCTqBui3pkh52rnTXDj79ttmPx8wF82mpkKLFvZmE/uorIhcwuHTh0ldlsqM9TMoKikCYOBtA0npnkKTkCY2pxPxLN9/b0rKW29BsVmaiLvvNnv4tGtnZzJxBSorIj9xsugkU1dPZfKqyZwoMos39GjQg4yeGbSv197mdCKe5bvvzKqzb75ZVlL69jUlpUMHO5OJK1FZEfmvopIiZm2YRcqyFA6eOghAmzptyOiZQa+GvXQbskg52roV0tPh3XfLTvf07m2uSenY0d5s4npUVqTSK7VKmffVPMZ9Oo6dx8ziDY2CG5HaI5WBzQbi7aXFG0TKyxdfmOtP/u//zLopAP36QUICdO5sbzZxXSorUmlZlsUn331C3NI4NuVuAiDMP4ykbkk80eYJfH187Q0o4kE++8yUlPffLxsbMMDclty2rW2xxE2orEiltHbvWmKXxpL1fRYAgY5AYjrFMOqOUfhX9bc3nIiHsCzIzIRJk2DxYjPm5QUPPQTx8bq7R66cyopUKtsObyMhM4GFWxcCUNWnKtHto4nvGk/t6rVtTifiGUpKYOFCU1I2bDBjPj7w+99DXBw00c10cpVUVqRS2Ju/l+SsZGZvmk2pVYq3lzeDWw4mOSqZiKAIu+OJeISCAnPr8QsvmEXdAPz84PHHYfRoaNDA3nzivlRWxKMdPXOUjBUZTP9sOgXFBQDc2+ReUnuk0iy0mc3pRDzD8eMwcya89BIcOGDGatWC6GjzuOEGW+OJB1BZEY90+uxppq2dxqSVkzhecByArhFdyeiVQaf6newNJ+Ih9u2DqVPhtdfghFmSiPr1YcwYM5uivXukvKisiEcpLi1m9sbZJGcns+/EPgBahLYgvWc6dzW+S2uliJSDb7+F5583S+IXmcWdad4cYmLg4YfBVzfSSTlTWRGPYFkWC7YuICEzge1HtgMQGRRJSvcUBrUYhI+3j80JRdzf55+b1WZ/vEZKly4QGwt33WXu9BFxBpUVcXuZuzKJXRLLun3rAAipHkJi10SGtRuGo4rD5nQi7s2yIDvbrDa7aFHZ+N13m5LSpYt92aTyUFkRt/X5/s+JWxrHou/Mb1B/X3/GdBzDmE5jCHQE2pxOxL2VlsIHH5iSsnatGfPxMad5YmLgV7+yN59ULior4nZ2HN1BYmYi87fMB8DX25ehbYeSeGciYTXCbE4n4t7OnjX79UyaZPbvAXA4zAWzY8ZAw4b25pPKSWVF3EbuyVxSslOY9fksikvN9qyDWgwipXsKDYP1G1TkehQUwOzZpqTk5JixwEAYPhxGjoQw/T1AbKSyIi4vryCPyasmM3XNVE6fPQ1Av5v7kdYzjVbhrewNJ+LmTp+GWbPM3T3795uxsDB45hkYNgyCguzNJwIqK+LCCooLmLFuBmnL0zhy5ggAHep1YFKvSXS7qZvN6UTc28mTMGMG/PWvcPCgGbvxRhg71pzy8fOzN5/Ij6msiMspKS1h7pdzScpKIifPzEc3DWlKWo80BjQdoLVSRK5DXh68/DJMmQJHj5qxm24ye/YMGWKuTxFxNd52B5gwYQJeXl4XPMLDw+2OJTawLIv3v3mflq+25A/v/YGcvBzqBdTjb/3/xuY/b+a+W+9TURG5RkePQlKSKSaJiea4cWOYMwe2b4c//UlFRVyXS8ysNGvWjCVLlpw/9vHRAl6VzYqcFcQuiWXlnpUABFcLJq5LHNG3R+Pnq/lokWt16JCZRXnllbIl8W+91RSWBx+EKi7xLiByeS7xY1qlShXNplRSmw9sJj4zng+3fwiAXxU/RnYYydguY6lZraa94UTcWG6u2f145kxzES2YtVHGjYP77wdv2+fVRa6cS5SVb7/9lrp16+JwOOjQoQNpaWk0vMTN/IWFhRQWFp4/zs/Pr6iYUo52H9/N+KzxzP1iLhYWPl4+PN76cZKikqgbUNfueCJuybJg0yZzauf1183tyABt25qS0r+/Soq4J9vLSocOHXj77be55ZZbOHDgABMnTqRTp05s2bKF2rVr/+z56enpJCcn25BUysOhU4dIXZ7KzPUzKSoxO6ANvG0gKd1TaBLSxOZ0Iu7n7FmzHP5778H775etkQLQsaMpKX37at8ecW9elnVuOyrXcOrUKRo1akRMTAyjR4/+2ccvNrNSv3598vLyCAzUEuuu6mTRSaasnsILq17gRJE5cd6jQQ8yembQvl57m9OJuJf8fPj4Y1NQPvrI3OFzTvXq0KcPREdDjx4qKeK68vPzCQoKuqL3b9tnVn7K39+fFi1a8O2331704w6HA4cuWXcbRSVFzNowi5RlKRw8ZRZzaFOnDRk9M+jVsJfu7hG5Qj/8YGZO3nsPMjPNjMo5oaHmFM+990KvXlojRTyPy5WVwsJCtm7dSteuXe2OIteh1Cpl3lfzGPfpOHYe2wlAo+BGpPZIZWCzgXh76cS5yOVYFnz1lSkn770H69df+PEmTUw5ufde6NDBbDIo4qlsLyvPPvss/fv3JyIigoMHDzJx4kTy8/MZMmSI3dHkGliWxSfffULc0jg25W4CIMw/jKRuSTzR5gl8fXztDSjiwoqLYeXKsoKyc2fZx7y84I47ygpK06b25RSpaLaXlb179/K73/2Ow4cPc8MNN3DHHXewZs0aIiMj7Y4mV2nt3rXELo0l6/ssAAIdgcR0imHUHaPwr+pvbzgRF3XqFHzyiSknH35YtqosmEXaevc25aR/f20mKJWX7WVl3rx5dkeQ67Tt8DYSMhNYuHUhAFV9qhLdPpq4rnGEVA+xOZ2Ia8rLg1Gj4H//F350zwC1asFvfgMDBpgLZf3V80XsLyvivvbm7yU5K5nZm2ZTapXi7eXN4JaDSY5KJiIowu54Ii5rwwazeuy50zwNG5ad3uncWavKivyU/kjIVTt65igZKzKY/tl0CorNqlP3NLmHtB5pNAttZnM6EddlWWYTwWefhaIiiIyEuXOhSxfdYixyOSorcsVOnz3NtLXTmLRyEscLjgPQJaILGT0z6BzR2d5wIi7u2DF4/HH45z/N8YABMHs2BAfbGkvELaisyC8qLi1m9sbZJGcns+/EPgCahzYno2cGdzW+S2uliPyCtWvh4Yfh+++halWzZ090tGZTRK6UyopckmVZLNi6gITMBLYf2Q5AZFAkKd1TGNRiED7eWthB5HIsC6ZOhbFjzW3JDRvC3/9u9uoRkSunsiIXlbkrk9glsazbtw6AkOohJHZNZFi7YTiqaAVhkV9y5Ag89pi5HRlg4ECzuWBQkK2xRNySyopcYOP+jcQujWXRd4sA8Pf1Z0zHMYzpNIZAh/ZeErkSq1aZ0z579pi1Ul58EYYO1WkfkWulsiIA7Di6g3GfjmPeV2bdG19vX4a2HUrinYmE1dBKVCJXorQUJk+GhAQoKYHGjc1pn1at7E4m4t5UViq53JO5pGSnMOvzWRSXFgMwqMUgUrqn0DC4oc3pRNzHoUMweLDZDRlg0CB49VUICLA3l4gnUFmppPIK8pi8ajJT10zl9NnTAPS9uS/pPdNpFd7K3nAibmbZMvjd72DfPqhWzayl8sc/6rSPSHlRWalkCooLmLFuBmnL0zhy5ggAt9e7nUm9JhF1U5S94UTcTEkJpKdDUpI5BdS0KfzjH9C8ud3JRDyLykolUVJawtwv55KUlUROXg4ATWo3Ia1nGvc1vU9rpYhcpQMH4Pe/hyVLzPGQIfDKK9rLR8QZVFY8nGVZfLD9A+KXxrPl0BYA6gXUY0LUBB5r9RhVvPUjIHK1MjPNNSkHDkD16jBjhikrIuIceqfyYCtyVhC7JJaVe1YCULNaTeK6xDHi9hH4+frZnE7E/ZSUwHPPQUqKWfCteXOYPx9uu83uZCKeTWXFA20+sJn4zHg+3G5Wo/Kr4sfIDiOJ6RxDsJ82IhG5Fvv2wSOPQFaWOX7iCXjpJTOzIiLOpbLiQXYf3834rPHM/WIuFhY+Xj483vpxkqKSqBtQ1+54Im5r0SJzfcqhQ1CjBrz2mjkNJCIVQ2XFAxw+fZjUZanMWD+DopIiAB647QEmdp9Ik5AmNqcTcV/FxTB+vLnjB6BlS7PI2y232JtLpLJRWXFjJ4tOMnX1VCavmsyJohMA9GjQg4yeGbSv197mdCLube9es3bKihXm+M9/hilTzDoqIlKxVFbcUFFJEbM2zCJlWQoHTx0EoHV4azJ6ZdC7YW/dhixynf79b7Ma7ZEjZgXav/0NHnzQ7lQilZfKihsptUqZ99U8xn06jp3HdgLQKLgRE3tM5MFmD+Lt5W1zQhH3dvas2ddn8mRz3LatudunUSN7c4lUdiorbsCyLD757hPilsaxKXcTAGH+YSR1S+KJNk/g6+Nrb0ARD7B7t9kpec0aczxihCktDoe9uUREZcXlrd27ltilsWR9nwVAoCOQmE4xjLxjJDWq1rA3nIiHeO89+MMf4NgxqFkTZs+G++6zO5WInKOy4qK2Hd5GQmYCC7cuBKCqT1Wi20cT1zWOkOohNqcT8QxFRRATY9ZLAbj9dnPa56abbI0lIj+hsuJi9ubvJTkrmdmbZlNqleLt5c3gloNJjkomIijC7ngiHmPnTnjoIVi/3hyPGQNpaVC1qr25ROTnVFZcxNEzR8lYkcH0z6ZTUFwAwD1N7iGtRxrNQpvZnE7EsyxYAH/8I+TnQ61a8Oab0L+/3alE5FJUVmx2+uxppq2dxqSVkzhecByALhFdyOiZQeeIzvaGE/EwBQXw7LNmd2SATp3gf/8XIjRpKeLSVFZsUlxazOyNs0nOTmbfiX0ANA9tTnrPdO5ufLfWShEpZ599BsOGwcaN5njsWLMhoa9uphNxeSorFcyyLBZsXUBCZgLbj2wHIDIokpTuKQxqMQgfbx+bE4p4ljVrIDkZPv7YHIeEwNy50LevvblE5MqprFSgzF2ZxC6JZd2+dQCEVA8hsWsiw9oNw1FFizmIlKfVq01J+eQTc+zjA48+ChMnQr169mYTkaujslIBNu7fSOzSWBZ9twgAf19/xnQcw5hOYwh0BNqcTsSzXKykDB5sVqbVSrQi7kllxYl2HN3BuE/HMe+reQD4evsytO1QEu9MJKxGmM3pRDzLqlWmpCwyfyfAxweGDIH4eJUUEXensuIEuSdzSclOYdbnsyguLQZgUItBpHRPoWFwQ5vTiXiWlStNSVm82BxXqVJWUhrqj5uIR1BZKUd5BXlMXjWZqWumcvrsaQD63tyX9J7ptApvZW84EQ+zYoUpKUuWmOMqVeCxx0xJadDA1mgiUs5UVspBQXEBM9bNIG15GkfOHAHg9nq3M6nXJKJuirI3nIiHWbECJkyApUvNsUqKiOdTWbkOJaUlzP1yLklZSeTk5QDQpHYT0nqmcV/T+7RWikg5Wr7clJTMTHNcpYrZfDA+Xnv5iHg6lZVrYFkWH2z/gPil8Ww5tAWAegH1SI5KZkirIVTx1v9WkfKybJkpKZ9+ao59fU1JiYtTSRGpLLztDgAwY8YMGjRoQLVq1Wjbti3Lly+3O9IlrchZQdc5Xbl33r1sObSF4GrBPN/reb4d8S2Pt3lcRUWknGRnQ48e0K2bKSq+vjB0KHz7Lbz2moqKSGVi+zvr/PnzGTVqFDNmzKBz58689tpr9OvXj6+//poIF9qwY/OBzcRnxvPh9g8B8Kvix8gOI4npHEOwX7DN6UQ8R1aWuXA2K8sc+/rC44+bmRQX+pUgIhXIy7Isy84AHTp0oE2bNsycOfP82K233sqAAQNIT0//xc/Pz88nKCiIvLw8AgPLf4G13cd3Mz5rPHO/mIuFhY+XD4+3fpykqCTqBtQt9+8nUlllZZnTPdnZ5tjXF554AmJjVVJEPNHVvH/bOrNSVFTEhg0biI2NvWC8T58+rFq1yqZUxqFTh0hdnsrM9TMpKikC4IHbHmBi94k0CWliazYRT2FZZSVl2TIzVrVqWUmpX9/OdCLiKmwtK4cPH6akpISwsAtXcw0LCyM3N/ein1NYWEhhYeH54/z8fKdky1iRwUtrXwKgR4MeZPTMoH299k75XiKVjWWZ61CSk1VSROSX2X7NCvCzW3wty7rkbb/p6ekkJyc7PdPYLmNZv389CV0T6N2wt25DFikHlmVuPU5ONrcigykpTz5pSsqNN9qbT0Rck61lJSQkBB8fn5/Nohw8ePBnsy3nxMXFMXr06PPH+fn51HfCX8NC/UPJfiy73L+uSGVkWWYRt+Rks6gbgMNhSsrYsSopInJ5tt66XLVqVdq2bcvic5t6/NfixYvp1KnTRT/H4XAQGBh4wUNEXJNlmeXwu3aF3r1NUXE4YMQI+O47mD5dRUVEfpntp4FGjx7No48+Srt27ejYsSOzZs0iJyeHYcOG2R1NRK7RuZIyYYLZDRlMSRk61Myk1NWNdCJyFWwvKw899BBHjhzhueeeY//+/TRv3px///vfREZG2h1NRK6SZZndjydMgNWrzVi1aqakxMSopIjItbF9nZXr5ex1VkTkl1kWLFpkSsqaNWasWjUYNsyUlDp1bI0nIi7IbdZZERH3ZlnwySfmwlmVFBFxFpUVEblqlgUff2xKytq1ZszPr6ykhIfbm09EPIvKiohcMcuC//zHlJTPPjNjfn7w5z/DX/6ikiIizqGyIiK/6FxJmTAB1q0zY35+8NRTpqRcYlkkEZFyobIiIpdkWfDvf5uSsn69GfPzg+HD4dlnVVJEpGKorIjIz1gWfPSROd1zrqRUr15WUkJD7c0nIpWLyoqInGdZ8OGHpqRs2GDGqleH6GgYM0YlRUTsobIiIlgWfPCBKSmff27G/P3LZlJuuMHefCJSuamsiFRilgXvv29KysaNZszfv2wmRSVFRFyByopIJXSxklKjRllJCQmxN5+IyI+prIhUIpYF771nSsqmTWasRg2zC/Lo0SopIuKaVFZEKoHS0rKS8sUXZqxGDXj6aVNSate2N5+IyOWorIh4sNJS+Ne/TEn58kszFhBgSsozz6ikiIh7UFkR8UClpfDPf5qSsnmzGQsIgJEjTUmpVcvefCIiV0NlRcSDlJbCwoXw3HMqKSLiOVRWRDxAaSksWGBKyldfmbHAQFNSRo1SSRER96ayIuLGzpWU5GTYssWMBQaagjJqFAQH25lORKR8qKyIuKHSUvi//zMzKedKSlCQKSgjR6qkiIhnUVkRcSMlJWUl5euvzVhQkLkeZeRIqFnT1ngiIk6hsiLiBkpK4B//MCVl61YzVrOmKSlPP62SIiKeTWVFxIWVlMDf/25KyrZtZqxmTbOQ29NPm1kVERFPp7Ii4oJKSmD+fEhJUUkREVFZEXEhJSUwb54pKd98Y8aCg01JGTFCJUVEKieVFREXcLGSUqtWWUkJDLQ3n4iInVRWRGxUXFxWUrZvN2O1asGYMRAdrZIiIgIqKyK2KC6G//1fU1K+/daM1aoFzz5rSkpAgL35RERcicqKSAUqLoZ33zUlZccOM1a7tikpw4erpIiIXIzKikgFKC6G//kfmDjxwpLyl7/AU0+ppIiIXI7KiogTFRfDO++YkvLdd2YsJKRsJqVGDXvziYi4A5UVESc4e9aUlNTUC0vKuZkUlRQRkSunsiJSjs6VlIkTYedOM3bDDaak/PnPKikiItdCZUWkHJw9C3PnmpKya5cZu+EGiIkxJcXf3958IiLuTGVF5DqcPQtvv21O95wrKaGhpqQMG6aSIiJSHlRWRK7B2bPw1lumpHz/vRkLDYWxY01JqV7d1ngiIh5FZUXkKhQVmZKSllZWUsLCTEkZOlQlRUTEGVRWRK5AURG8+aYpKbt3mzGVFBGRiqGyInIZRUUwZ44pKTk5Ziw83JSUP/1JJUVEpCJ42/nNb7rpJry8vC54xMbG2hlJBDAl5bXXoHFjcw1KTg7UqQMvvmhuSR41SkVFRKSi2D6z8txzz/Hkk0+eP66hhSjERoWFZTMpe/aYsTp1IDYWnnwS/PzszSciUhnZXlYCAgIIDw+3O4ZUcoWFMHu2KSl795qxOnUgLg6eeEIlRUTETraeBgKYNGkStWvXplWrVqSmplJUVHTZ5xcWFpKfn3/BQ+RaFRbCjBlw881mGfy9e6FuXZg+3ZzuGTFCRUVExG62zqyMHDmSNm3aEBwczGeffUZcXBy7du3ib3/72yU/Jz09neTk5ApMKZ6ooADeeAMyMspmUurVMzMpjz8O1arZm09ERMp4WZZllecXnDBhwi+WiXXr1tGuXbufjS9YsIAHHniAw4cPU7t27Yt+bmFhIYWFheeP8/PzqV+/Pnl5eQQGBl5fePF4BQXwt7+ZkvLDD2ZMJUVEpOLl5+cTFBR0Re/f5T6zEh0dzcMPP3zZ59x0000XHb/jjjsA2LFjxyXLisPhwOFwXFdGqXzOlZT0dNi3z4zdeGNZSdGPlIiI6yr3shISEkJISMg1fe7GjRsBqFOnTnlGkkqsoABef93MpJwrKfXrm5Lyxz+qpIiIuAPbrllZvXo1a9asoXv37gQFBbFu3TqeeeYZ7rnnHiIiIuyKJR7izJmykrJ/vxmrXx/i4+EPf1BJERFxJ7aVFYfDwfz580lOTqawsJDIyEiefPJJYmJi7IokHuDMGZg1y5SU3FwzVr8+JCTAY4+ppIiIuCPbykqbNm1Ys2aNXd9ePMyZM2bF2UmTykpKRERZSala1dZ4IiJyHWxfFE7kepw+bUrK88+XlZTISFNShgxRSRER8QQqK+KWTp+GV181JeXAATOmkiIi4plUVsStnCspkybBwYNm7KabTEkZPFglRUTEE6msiFs4dapsJuVcSWnQoKyk+Pram09ERJxHZUVc2qlTMHOmKSmHDpmxBg0gMREefVQlRUSkMlBZEZd06pTZYHDy5LKS0rChKSm//71KiohIZaKyIi7l5MmyknL4sBlr1MiUlEceUUkREamMVFbEJZw8Ca+8Ai+8cGFJGTfOlJQq+kkVEam09BYgtjpxoqykHDlixm6+uWwmRSVFRET0ViC2uFhJadzYlJRBg1RSRESkjN4SpEKdOAEvv2xKytGjZqxxY3O653e/U0kREZGf01uDVIj8fFNS/vrXspJyyy2mpDz8sEqKiIhcmt4ixKny82H6dFNSjh0zY7fcAuPHm5Li42NvPhERcX0qK+IU+fkwbRpMmVJWUpo0MSXloYdUUkRE5MqprEi5ysszJWXq1LKS0rSpKSkPPqiSIiIiV09lRcpFXh689JIpKcePm7FbbzXXpKikiIjI9VBZkety/HjZTMqPS8r48TBwoEqKiIhcP5UVuSbHj5fNpOTlmbHbbjMl5YEHVFJERKT8qKzIVTl+HF580Tx+XFKSkkxJ8fa2MZyIiHgklRW5IseOlZWU/Hwz1qyZKSm//a1KioiIOI/KilzWsWPmVM9LL5WVlObNTUm5/36VFBERcT6VFbmoo0dNSZk2rayktGhhSsp996mkiIhIxVFZkQscPWoWcps2zezjA/CrX5kLZ1VSRETEDiorApidj8/NpPy4pCQlwYABKikiImIflZVK7siRspmUkyfNWMuWpqTce69KioiI2E9lpZI6fNiUlOnTy0pKq1ampNxzj0qKiIi4DpWVSubwYbMD8vTpcOqUGWvVCiZMMCXFy8vOdCIiIj+nslJJHDpkSsrLL5eVlNatTUnp318lRUREXJfKioc7dAheeAFeeaWspLRpY0rKb36jkiIiIq5PZcVDHTxYVlJOnzZjbduaa1JUUkRExJ2orHiYi5WUdu1MSbn7bpUUERFxPyorHuLgQZg8GWbMKCsp7dubknLXXSopIiLivlRW3NyBA2Ul5cwZM9a+vbkmpV8/lRQREXF/KituKjfXlJSZM8tKyu23m5LSt69KioiIeA6VFTeTmwvPPw+vvlpWUjp0MCXl179WSREREc+jsuIm9u8vKykFBWasQwdIToY+fVRSRETEczl1UfXU1FQ6depE9erVqVmz5kWfk5OTQ//+/fH39yckJISnn36aoqIiZ8ZyK/v3wzPPQMOG8OKLpqjccQd8/DGsXq3ZFBER8XxOnVkpKipi4MCBdOzYkTfeeONnHy8pKeHuu+/mhhtuYMWKFRw5coQhQ4ZgWRbTp093ZjSXt38/TJoEr71WNpPSsaM53dO7twqKiIhUHk4tK8nJyQC8+eabF/34okWL+Prrr9mzZw9169YF4K9//SuPPfYYqampBAYGOjOeS9q3r6ykFBaasU6dTEnp1UslRUREKh9b99ZdvXo1zZs3P19UAH79619TWFjIhg0bLvo5hYWF5OfnX/DwBD/8AE8/bU73TJtmikrnzrB4MaxYodkUERGpvGwtK7m5uYSFhV0wFhwcTNWqVcnNzb3o56SnpxMUFHT+Ub9+/YqI6jQ//AAjRkCjRmYn5MJC6NIFliyB5cs1myIiInLVZWXChAl4eXld9rF+/for/npeF3kntizrouMAcXFx5OXlnX/s2bPnav8TXMLevRAdbWZSXn7ZlJSuXWHpUli2DHr2VEkRERGBa7hmJTo6mocffviyz7npppuu6GuFh4ezdu3aC8aOHTvG2bNnfzbjco7D4cDhcFzR13dFe/ZARgb87W9w7qanrl3NLchRUSooIiIiP3XVZSUkJISQkJBy+eYdO3YkNTWV/fv3U6dOHcBcdOtwOGjbtm25fA9XcbGScued5sJZlRQREZFLc+rdQDk5ORw9epScnBxKSkrYtGkTADfffDM1atSgT58+3HbbbTz66KNMnjyZo0eP8uyzz/Lkk096zJ1Ae/ZAejq88UZZSenWraykiIiIyOU5tayMHz+et9566/xx69atAfj000+JiorCx8eHjz76iKeeeorOnTvj5+fHoEGDeOGFF5wZq0Lk5JSVlLNnzVhUlNkFWSVFRETkynlZlmXZHeJ65OfnExQURF5enkvMxuzebUrK7NllJaV7d1NSunWzN5uIiIiruJr3b+0NVE5274a0NJgzp6yk9OhhSsqdd9qbTURExJ2prFyn7783JeXNN8tKSs+epqR07WpnMhEREc+gsnKNdu0qKynFxWasVy9TUrp0sTWaiIiIR1FZuUoqKSIiIhVLZeUK7dxpSspbb5WVlN69TUnp3NnebCIiIp5MZeUX7NwJqammpJSUmLE+fUxJ6dTJ3mwiIiKVgcrKJXz3nSkpb79dVlJ+/WtTUjp2tDebiIhIZaKycglvvmluQwbo29eUlDvusDWSiIhIpaSycgnPPANffw0xMdChg91pREREKi+VlUuoVQsWLLA7hYiIiHjbHUBERETkclRWRERExKWprIiIiIhLU1kRERERl6ayIiIiIi5NZUVERERcmsqKiIiIuDSVFREREXFpKisiIiLi0lRWRERExKWprIiIiIhLU1kRERERl6ayIiIiIi7N7XddtiwLgPz8fJuTiIiIyJU697597n38cty+rJw4cQKA+vXr25xERERErtaJEycICgq67HO8rCupNC6stLSUffv2ERAQgJeXV7l+7fz8fOrXr8+ePXsIDAws168t5Uevk3vQ6+Qe9Dq5B094nSzL4sSJE9StWxdv78tfleL2Myve3t7ceOONTv0egYGBbvvDUJnodXIPep3cg14n9+Dur9MvzaicowtsRURExKWprIiIiIhLU1m5DIfDQVJSEg6Hw+4ochl6ndyDXif3oNfJPVS218ntL7AVERERz6aZFREREXFpKisiIiLi0lRWRERExKWprIiIiIhLU1m5hBkzZtCgQQOqVatG27ZtWb58ud2R5EfS09Np3749AQEBhIaGMmDAAL755hu7Y8kvSE9Px8vLi1GjRtkdRS7ihx9+4Pe//z21a9emevXqtGrVig0bNtgdS36kuLiYxMREGjRogJ+fHw0bNuS5556jtLTU7mhOpbJyEfPnz2fUqFEkJCSwceNGunbtSr9+/cjJybE7mvxXdnY2w4cPZ82aNSxevJji4mL69OnDqVOn7I4ml7Bu3TpmzZrFr371K7ujyEUcO3aMzp074+vry3/+8x++/vpr/vrXv1KzZk27o8mPTJo0iVdffZWXX36ZrVu38vzzzzN58mSmT59udzSn0q3LF9GhQwfatGnDzJkzz4/deuutDBgwgPT0dBuTyaUcOnSI0NBQsrOzufPOO+2OIz9x8uRJ2rRpw4wZM5g4cSKtWrXixRdftDuW/EhsbCwrV67ULLKL+81vfkNYWBhvvPHG+bHf/va3VK9enblz59qYzLk0s/ITRUVFbNiwgT59+lww3qdPH1atWmVTKvkleXl5ANSqVcvmJHIxw4cP5+6776ZXr152R5FLeP/992nXrh0DBw4kNDSU1q1b8/rrr9sdS36iS5cuLF26lO3btwPwxRdfsGLFCu666y6bkzmX229kWN4OHz5MSUkJYWFhF4yHhYWRm5trUyq5HMuyGD16NF26dKF58+Z2x5GfmDdvHp9//jnr1q2zO4pcxs6dO5k5cyajR48mPj6ezz77jKeffhqHw8HgwYPtjif/NXbsWPLy8mjatCk+Pj6UlJSQmprK7373O7ujOZXKyiV4eXldcGxZ1s/GxDVER0fz5ZdfsmLFCrujyE/s2bOHkSNHsmjRIqpVq2Z3HLmM0tJS2rVrR1paGgCtW7dmy5YtzJw5U2XFhcyfP5933nmHd999l2bNmrFp0yZGjRpF3bp1GTJkiN3xnEZl5SdCQkLw8fH52SzKwYMHfzbbIvYbMWIE77//PsuWLePGG2+0O478xIYNGzh48CBt27Y9P1ZSUsKyZct4+eWXKSwsxMfHx8aEck6dOnW47bbbLhi79dZbWbBggU2J5GL+8pe/EBsby8MPPwxAixYt2L17N+np6R5dVnTNyk9UrVqVtm3bsnjx4gvGFy9eTKdOnWxKJT9lWRbR0dEsXLiQzMxMGjRoYHckuYiePXuyefNmNm3adP7Rrl07HnnkETZt2qSi4kI6d+78s9v/t2/fTmRkpE2J5GJOnz6Nt/eFb90+Pj4ef+uyZlYuYvTo0Tz66KO0a9eOjh07MmvWLHJychg2bJjd0eS/hg8fzrvvvst7771HQEDA+ZmwoKAg/Pz8bE4n5wQEBPzsOiJ/f39q166t64tczDPPPEOnTp1IS0vjwQcf5LPPPmPWrFnMmjXL7mjyI/379yc1NZWIiAiaNWvGxo0bmTJlCn/84x/tjuZcllzUK6+8YkVGRlpVq1a12rRpY2VnZ9sdSX4EuOhjzpw5dkeTX9CtWzdr5MiRdseQi/jggw+s5s2bWw6Hw2ratKk1a9YsuyPJT+Tn51sjR460IiIirGrVqlkNGza0EhISrMLCQrujOZXWWRERERGXpmtWRERExKWprIiIiIhLU1kRERERl6ayIiIiIi5NZUVERERcmsqKiIiIuDSVFREREXFpKisiIiLi0lRWRERExKWprIiIiIhLU1kRERERl6ayIiIiIi7t/wERDR/SZr5RPQAAAABJRU5ErkJggg==", + "text/plain": [ + "<Figure size 640x480 with 1 Axes>" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "MAX_N = 5000\n", + "STEP = 100\n", + "\n", + "# Used to avoid errors when taking the log of the histogram\n", + "EPS = 0.00001\n", + "\n", + "x = np.arange(1, MAX_N, STEP)\n", + "x_ref = np.arange(np.log(x[0]), np.log(x[-1]))\n", + "\n", + "histogram = np.fromfile(\"histogram.bin\", dtype=np.uint32)\n", + "plt.plot(np.log(x), np.log(histogram + EPS), color=\"blue\")\n", + "plt.plot(x_ref, 3 * x_ref, color=\"green\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63849432-9c7b-4866-a647-eb051161c590", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.12.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}
diff --git a/getrf/benchmark/build.sh b/getrf/benchmark/build.sh @@ -0,0 +1,10 @@ +#!/bin/sh + +set -xe + +EXE_NAME="benchmark" +C_FLAGS="-Wall -Wextra -llapacke" + +ctags -R . +cc ./main.c -o "./target/$EXE_NAME-debug" $C_FLAGS -g -DDEBUG # debug version +cc ./main.c -o "./target/$EXE_NAME" $C_FLAGS -O2 # release version
diff --git a/getrf/benchmark/config.h b/getrf/benchmark/config.h @@ -0,0 +1,19 @@ +/* + * The number of tests for each size of matrix + */ +#define N_TESTS 1 + +/* + * The biggest possible size of the input matrix + */ +#define MAX_N 10000 + +/* + * The number of threads in which the benchmarks will be distributed + */ +#define N_THREADS 8 + +/* + * The interval between successive terms in the histogram + */ +#define STEP 100
diff --git a/getrf/benchmark/histogram.bin b/getrf/benchmark/histogram.bin Binary files differ.
diff --git a/getrf/benchmark/main.c b/getrf/benchmark/main.c @@ -0,0 +1,163 @@ +#include <stdio.h> +#include <lapacke.h> +#include <pthread.h> + +#include <stdint.h> +#include <stdbool.h> +#include <time.h> + +#include <string.h> +#include <errno.h> + +#include "config.h" + +#define HISTOGRAM_SIZE (MAX_N/STEP) +#define CLOCKS_PER_MILLIS (CLOCKS_PER_SEC/1000) +#define LOGISTICS_INITIAL_CONDITION (-800.) + +#define PROGRESS_BAR_TOTAL HISTOGRAM_SIZE +#include "progress-bar.h" + +uint32_t histogram[HISTOGRAM_SIZE]; + +// .data is a pointer because it should be allocated in the heap +// (.data DOES NOT fit in the stack 🤡) +typedef struct { + double *data; + double *ref_data; + int32_t ipiv[MAX_N]; + size_t id; +} Thread; + +typedef struct { + Thread threads[N_THREADS]; +} Benchmarker; + +uint32_t thread_run_benchmark(Thread *thread, size_t n) +{ + int64_t total_time = 0; + clock_t start, end; + + for (size_t i = 0; i < N_TESTS; i++) { + // Reinitializes the values in the .data array: this avoids progressively + // moving larger values to the beginning of the array, which would decrise + // decrease the number of row interchanges required for computations + memcpy(thread->data, thread->ref_data, sizeof(double)*n*n); + + start = clock(); + (void)LAPACKE_dgetrf(LAPACK_ROW_MAJOR, + (lapack_int)n, (lapack_int)n, + thread->data, (lapack_int)n, thread->ipiv); + end = clock(); + total_time += (end - start)/CLOCKS_PER_MILLIS; + } + + return (uint32_t)(total_time/N_TESTS); +} + +void *thread_benchmark(void *arg) +{ + Thread *thread = (Thread*) arg; + + // Computations are distributed evenly across threads + for (size_t n = 1 + thread->id*STEP; n < MAX_N; n += STEP*N_THREADS) { + uint32_t duration = thread_run_benchmark(thread, n); + + histogram[(n - 1)/STEP] = duration; + progress_bar_inc(); + } + + return NULL; +} + +Benchmarker benchmarker_new(double *ref_data) +{ + Benchmarker bench = {0}; + + // These arrays will live for the entire duration of the program, + // so we just leak them 🤡 + double *data = malloc(sizeof(double)*MAX_N*MAX_N*N_THREADS); + if (data == NULL) { + fprintf(stderr, "ERROR: Buy more RAM!\n"); + exit(EXIT_FAILURE); + } + + for (size_t i = 0; i < N_THREADS; i++) { + bench.threads[i].data = data + i*MAX_N*MAX_N; + bench.threads[i].ref_data = ref_data; + bench.threads[i].id = i; + } + + return bench; +} + +void benchmarker_run(Benchmarker *bench) +{ + pthread_t handles[N_THREADS] = {0}; + for (size_t i = 0; i < N_THREADS; i++) { + if (pthread_create(&handles[i], NULL, + thread_benchmark, + &bench->threads[i]) != 0) { + fprintf(stderr, "ERROR: Failed to spawn thread %lu!\n", i); + exit(EXIT_FAILURE); + } + } + + for (size_t i = 0; i < N_THREADS; i++) { + pthread_join(handles[i], NULL); + } +} + +int main(void) +{ + const char *output_path = "histogram.bin"; + FILE *output = fopen(output_path, "w"); + + if (output == NULL) { + fprintf(stderr, + "ERROR: Coundn't open output file \"%s\": %s", + output_path, strerror(errno)); + return EXIT_FAILURE; + } + + // ======================================================================== + printf("INFO: Initializing random input data... "); + + double *ref_data = malloc(sizeof(double)*MAX_N*MAX_N); + if (ref_data == NULL) { + fprintf(stderr, "ERROR: Buy more RAM!\n"); + exit(EXIT_FAILURE); + } + + // Pseudorandom data given by the logistic map + double acc = LOGISTICS_INITIAL_CONDITION; + for (size_t i = 0; i < MAX_N*MAX_N; i++) { + ref_data[i] = acc; + acc = 1000. - acc*acc/500.; + } + + printf("done!\n"); + + // ======================================================================== + printf("INFO: Benchmarking the dgetrf function on %ux%u matrices... (using %u threads)\n", + MAX_N, MAX_N, N_THREADS); + + Benchmarker bench = benchmarker_new(ref_data); + benchmarker_run(&bench); + + // ======================================================================== + printf("INFO: Done with benchmarking! Saving histogram to disk...\n"); + + size_t written = fwrite(histogram, sizeof(uint32_t), HISTOGRAM_SIZE, output); + + if (written < HISTOGRAM_SIZE) { + fprintf(stderr, + "ERROR: Failed to write histogram to output file \"%s\": %s\n", + output_path, strerror(errno)); + return EXIT_FAILURE; + } + + printf("INFO: Done!\n"); + + return EXIT_SUCCESS; +}
diff --git a/getrf/benchmark/progress-bar.h b/getrf/benchmark/progress-bar.h @@ -0,0 +1,36 @@ +#ifndef PROGRESS_BAR_H_ +#define PROGRESS_BAR_H_ + +#include <stdio.h> +#include <stdint.h> +#include <pthread.h> + +#define PROGRESS_BAR_LENGTH 50 + +size_t progress_last_printed_count = 0; +size_t progress_count = 0; +pthread_mutex_t progress_mutex = PTHREAD_MUTEX_INITIALIZER; + +void progress_bar_inc(void) +{ + pthread_mutex_lock(&progress_mutex); + + progress_count++; + size_t filled_length = (PROGRESS_BAR_LENGTH*progress_count)/PROGRESS_BAR_TOTAL; + size_t empty_length = PROGRESS_BAR_LENGTH - filled_length; + + printf("\r["); + for (size_t i = 0; i < filled_length; i++) printf("="); + for (size_t i = 0; i < empty_length; i++) printf(" "); + + size_t percent = (100 * progress_count) / PROGRESS_BAR_TOTAL; + if (percent == 100) printf("] %3zu%%\n", percent); + else printf("] %3zu%%", percent); + + progress_last_printed_count = progress_count; + + fflush(stdout); + pthread_mutex_unlock(&progress_mutex); +} + +#endif // PROGRESS_BAR_H_