KRYLOV .ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ax=b ? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### given a random matrix A and a vector b, we will try to find x st: Ax=b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "n=30\n",
    "A = np.random.rand(n,n)\n",
    "AS=A.T.dot(A)\n",
    "#AS=A.T*A\n",
    "b=np.ones(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exact solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol_ex=np.linalg.lstsq(A,b,rcond=None)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0=np.zeros(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ARNOLDI PROCESS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def arnoldi_iteration(A,r0,M):\n",
    "    n = A.shape[0]\n",
    "    h = np.zeros((M + 1, M))\n",
    "    V= np.zeros((n,M+1))\n",
    "    V[:,0] = r0 / np.linalg.norm(r0)\n",
    "\n",
    "    for j in range(M):\n",
    "        w = A.dot(V[:,j])\n",
    "        for i in range(j + 1):\n",
    "            h[i, j] = np.dot(w,V[:,i])\n",
    "            w = w - h[i, j] * V[:,i]\n",
    "\n",
    "        h[j + 1, j] = np.linalg.norm(w)\n",
    "        if h[j + 1, j] == 0:\n",
    "            return V,h\n",
    "        else:\n",
    "            V[:,j+1] = w / h[j + 1, j]\n",
    "    return V,h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "r0=b-A.dot(x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "V,h=arnoldi_iteration(A,r0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[15.28308974],\n",
       "       [ 1.56989857]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GMRES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GMRES A BASE DE LA FONCTION ARNOLDI DEJA DEFINI\n",
    "def GMRES_R(A, b, x0, nmax_iter,tol):\n",
    "    m=min(nmax_iter, A.shape[0])\n",
    "    n = A.shape[0]\n",
    "    rm=b-A.dot(x0)\n",
    "    for j in range(min(nmax_iter, A.shape[0])):\n",
    "        beta=np.linalg.norm(rm)\n",
    "        V,h=arnoldi_iteration(A,rm,j+1)\n",
    "        s=np.zeros(j+2)\n",
    "        s[0]=beta\n",
    "        y=np.linalg.lstsq(h,s,rcond=None)[0]\n",
    "        xm=x0+V[:,:-1].dot(y)\n",
    "        rm=b-A.dot(xm)\n",
    "        if np.linalg.norm(rm)<=tol:\n",
    "            return xm\n",
    "        else:\n",
    "            x0=xm\n",
    "    return xm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RK:<br>\n",
    "j'ai defini m comme etant le minimum entre nombre d'iteration et la dimension de la matrie A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.06575817,  0.33251067, -0.5729281 , -0.8131979 , -0.54688738,\n",
       "       -0.30520559,  1.13005732,  0.16601282,  1.19785581, -1.1956725 ,\n",
       "        0.07130517,  0.52358126,  0.78286441, -0.12647668,  0.07331804,\n",
       "       -0.24769921, -0.83677894,  1.6340695 ,  0.75501444,  0.49200097,\n",
       "       -0.14564908, -1.70791013,  1.48004496,  0.19575725,  0.07903655,\n",
       "       -0.82935719, -1.21511526, -0.17379711, -0.00354716,  1.63978163])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "GMRES_R(A, b, x0,  n,10e-16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.06575817,  0.33251067, -0.5729281 , -0.8131979 , -0.54688738,\n",
       "       -0.30520559,  1.13005732,  0.16601282,  1.19785581, -1.1956725 ,\n",
       "        0.07130517,  0.52358126,  0.78286441, -0.12647668,  0.07331804,\n",
       "       -0.24769921, -0.83677894,  1.6340695 ,  0.75501444,  0.49200097,\n",
       "       -0.14564908, -1.70791013,  1.48004496,  0.19575725,  0.07903655,\n",
       "       -0.82935719, -1.21511526, -0.17379711, -0.00354716,  1.63978163])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_ex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### THE SYMMETRIC LANCZOS PROCESS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SLP(A,v,m):\n",
    "    n = A.shape[0]\n",
    "    V= np.zeros((n,m+2))\n",
    "    V[:,1]=v\n",
    "    beta=np.zeros(m+2)\n",
    "    alpha=np.zeros(m+2)\n",
    "    for j in range(1,m+1):\n",
    "        w = A.dot(V[:,j])-beta[j]*V[:,j-1]\n",
    "        alpha[j]=np.dot(w,V[:,j])\n",
    "        w=w-alpha[j]*V[:,j]\n",
    "        beta[j+1]=np.linalg.norm(w)\n",
    "        if beta[j+1]==0:\n",
    "            return V,alpha,beta\n",
    "        else:\n",
    "            V[:,j+1]=w/beta[j+1]\n",
    "    Tm=np.zeros([m,m])\n",
    "    np.fill_diagonal(Tm,alpha[1:-1])\n",
    "    np.fill_diagonal(Tm[1:,:-1],beta[2:-1])\n",
    "    np.fill_diagonal(Tm[:-1,1:],beta[2:-1])\n",
    "    return Tm,V[:,1:-1],alpha,beta\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Verification : Tm = Vm.TAVm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = np.random.rand(n)\n",
    "normalized_v = v/np.linalg.norm(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tm,V,_,_=SLP(AS,normalized_v ,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[172.35430877,  93.27587224],\n",
       "       [ 93.27587224,  53.09570218]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Tm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.20262915,  0.08316072],\n",
       "       [ 0.22143518, -0.02320922],\n",
       "       [ 0.23084847,  0.05463803],\n",
       "       [ 0.05292901,  0.25809707],\n",
       "       [ 0.16494616,  0.00771284],\n",
       "       [ 0.27767869, -0.07343001],\n",
       "       [ 0.23406361, -0.05915094],\n",
       "       [ 0.2804858 , -0.11373393],\n",
       "       [ 0.03319008,  0.30429653],\n",
       "       [ 0.1525191 ,  0.05486626],\n",
       "       [ 0.03397658,  0.37057834],\n",
       "       [ 0.1642521 ,  0.14814441],\n",
       "       [ 0.19376546, -0.01276017],\n",
       "       [ 0.08104869,  0.22439775],\n",
       "       [ 0.25713403, -0.12192357],\n",
       "       [ 0.24688544, -0.11235901],\n",
       "       [ 0.11692914,  0.21686359],\n",
       "       [ 0.14437791,  0.11877191],\n",
       "       [ 0.04361118,  0.30455532],\n",
       "       [ 0.17700046,  0.06904763],\n",
       "       [ 0.15265439,  0.0471817 ],\n",
       "       [ 0.09471126,  0.20384811],\n",
       "       [ 0.18791243,  0.0539101 ],\n",
       "       [ 0.04047128,  0.26298254],\n",
       "       [ 0.18779279, -0.03558001],\n",
       "       [ 0.29728957, -0.15696417],\n",
       "       [ 0.27940205, -0.18724161],\n",
       "       [ 0.0313564 ,  0.31830241],\n",
       "       [ 0.22915426, -0.04208545],\n",
       "       [ 0.01359536,  0.36002304]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[172.35430877,  93.27587224],\n",
       "       [ 93.27587224,  53.09570218]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S=np.dot(V.T,np.dot(AS,V))\n",
    "S"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "on constate que Tm~= S, puisque la matrice de depart comporte des valeurs entre (0,1) et du a la trancature, les valeurs trop petites 10-6~=0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### THE SYMMETRIC LANCZOS ALG FOR LINEAR SYSTEMS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SLLS(A, b, x0, nmax_iter,tol):\n",
    "    n=A.shape\n",
    "    rm=b-A.dot(x0)\n",
    "    counter=0\n",
    "    for m in range(min(nmax_iter, A.shape[0])):\n",
    "        beta=np.linalg.norm(rm)\n",
    "        v=rm / beta\n",
    "        Tm,Vm,_,_=SLP(A,v,m+1)\n",
    "        e1=np.zeros(m+1)\n",
    "        e1[0]=1\n",
    "        y=np.dot(beta*np.linalg.inv(Tm),e1)\n",
    "        xm=x0+np.dot(Vm,y)\n",
    "        rm=b-A.dot(xm)\n",
    "        counter+=1\n",
    "        if np.linalg.norm(rm)<=tol:\n",
    "            return xm,counter\n",
    "        else:\n",
    "            x0=xm\n",
    "    return xm,counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "AS=A.T.dot(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ -0.08267342,  -1.75615805,   5.24741237,  10.12015866,\n",
       "          5.2904957 ,   5.88330149,  -6.78756383,   0.56783244,\n",
       "        -12.28698972,  12.76901309,  -1.01982066,  -7.94097067,\n",
       "         -3.88783706,   2.85798572,  -4.38634967,   3.97380469,\n",
       "          8.37136648, -10.28334742,  -7.29027218,  -2.890096  ,\n",
       "          0.6895806 ,  13.27831808, -14.55403273,  -1.55537226,\n",
       "         -1.13266276,   6.80176358,  14.71077431,   3.84063388,\n",
       "         -0.10675035, -15.0046158 ]),\n",
       " 26)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SLLS(AS,b,x0,n,10e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -0.08112645,  -1.75943472,   5.26886144,  10.15649825,\n",
       "         5.31870338,   5.90673874,  -6.83539728,   0.56868664,\n",
       "       -12.32833984,  12.81054416,  -1.01784308,  -7.95701485,\n",
       "        -3.92171129,   2.85982466,  -4.38422581,   3.99231871,\n",
       "         8.40492751, -10.33919565,  -7.3185995 ,  -2.90339488,\n",
       "         0.68540822,  13.33973616, -14.60756679,  -1.56129672,\n",
       "        -1.13837327,   6.83784288,  14.7674138 ,   3.85536289,\n",
       "        -0.10508054, -15.0703437 ])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_ex_sy=np.linalg.lstsq(AS,b,rcond=None)[0]\n",
    "sol_ex_sy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CONSTAT:\n",
    "pour le cas symmetrique, en definissant la matrice symmetrique par A.T.dot(A), l'algorithme converge en 2*n iterations avec n la taille de A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### THE NONSYMMETRIC LANCZOS BIORTHOGONALIZATION PROCESS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def NLBP(A,w,v,m):\n",
    "    n=A.shape[0]\n",
    "    #W=np.column_stack((np.zeros(n),w))\n",
    "    #V=np.column_stack((np.zeros(n),v))\n",
    "    V= np.zeros((n,m+2))\n",
    "    V[:,1]=v\n",
    "    W= np.zeros((n,m+2))\n",
    "    W[:,1]=w\n",
    "    beta=np.zeros(m+2)\n",
    "    gama=np.zeros(m+2)\n",
    "    alpha=np.zeros(m+2)\n",
    "    add_last=np.zeros(m)\n",
    "    for j in range(1,m+1):\n",
    "        alpha[j]=np.dot(A.dot(V[:,j]),W[:,j])\n",
    "        v=A.dot(V[:,j])-alpha[j]*V[:,j]-beta[j]*V[:,j-1]\n",
    "        w=A.T.dot(W[:,j])-alpha[j]*W[:,j]-gama[j]*W[:,j-1]\n",
    "        gama[j+1]=np.sqrt(np.abs(v.dot(w)))\n",
    "        if gama[j+1]!=0:\n",
    "            beta[j+1]=v.dot(w)/gama[j+1]\n",
    "            W[:,j+1]=w/beta[j+1]\n",
    "            V[:,j+1]=v/gama[j+1]\n",
    "    Tm=np.zeros([m,m])\n",
    "    add_last[-1]=gama[-1]\n",
    "    Tm=np.row_stack((Tm,add_last))\n",
    "    np.fill_diagonal(Tm,alpha[1:-1])\n",
    "    np.fill_diagonal(Tm[1:,:-1],gama[2:-1])\n",
    "    np.fill_diagonal(Tm[:-1,1:],beta[2:-1])\n",
    "    return Tm,V[:,1:-1],W[:,1:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "v=np.random.rand(n)\n",
    "w=np.zeros(n)\n",
    "w[0]=1/v[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tm,V,W=NLBP(A,w,v,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Verification Tm=WmTAVm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[7.37475158, 7.28422295],\n",
       "       [7.28422295, 7.8109757 ],\n",
       "       [0.        , 1.66010703]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Tm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[7.37475158, 7.28422295],\n",
       "       [7.28422295, 7.8109757 ]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W.T.dot(A.dot(V))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# THE BICONJUGATE GRADIENT (BCG) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BCG(A,b,x0,max_iter,tol):\n",
    "    n=A.shape\n",
    "    r0=b-A.dot(x0)\n",
    "    p0=r0.copy()\n",
    "    rt0=np.zeros(n)\n",
    "    rt0[0]=1/r0[0]\n",
    "    pt0=rt0.copy()\n",
    "    counter=0\n",
    "    while counter!=max_iter and np.linalg.norm(r0)>tol:\n",
    "        alpha=r0.dot(rt0)/A.dot(p0).dot(pt0)\n",
    "        x0=x0+alpha*p0\n",
    "        r1=r0-alpha*A.dot(p0)\n",
    "        rt1=rt0-alpha*A.T.dot(pt0)\n",
    "        beta=r1.dot(rt1)/r0.dot(rt0)\n",
    "        p1=r1+beta*p0\n",
    "        pt1=rt1+beta*pt0\n",
    "        r0=r1\n",
    "        rt0=rt1\n",
    "        p0=p1\n",
    "        pt0=pt1\n",
    "        counter+=1\n",
    "    return x0,counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-0.71232487, -0.15232019, -0.2480723 , -0.40939757,  0.02078813,\n",
       "         1.01447317, -1.15060546, -1.08452049, -1.4795648 , -0.41534727,\n",
       "         0.25618865,  1.29359381,  0.2837119 ,  0.95949902,  0.75405316,\n",
       "         0.85554373,  0.36217283, -1.04461944,  1.50501175, -0.13241054,\n",
       "         0.85753483, -1.02272784,  0.34989967,  0.31340192,  0.04951765,\n",
       "         0.0582766 ,  0.03200622, -0.75200104,  0.9714533 , -0.50278576]),\n",
       " 100)"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BCG(A,b,x0,100,10e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.06575817,  0.33251067, -0.5729281 , -0.8131979 , -0.54688738,\n",
       "       -0.30520559,  1.13005732,  0.16601282,  1.19785581, -1.1956725 ,\n",
       "        0.07130517,  0.52358126,  0.78286441, -0.12647668,  0.07331804,\n",
       "       -0.24769921, -0.83677894,  1.6340695 ,  0.75501444,  0.49200097,\n",
       "       -0.14564908, -1.70791013,  1.48004496,  0.19575725,  0.07903655,\n",
       "       -0.82935719, -1.21511526, -0.17379711, -0.00354716,  1.63978163])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_ex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QUASI MINIMAL RESIDUAL (QMR) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def QMR(A,b,x0,max_iter):\n",
    "    r0=b-A.dot(x0)\n",
    "    beta=np.linalg.norm(r0)\n",
    "    w=r0/beta\n",
    "    v=r0/beta\n",
    "    for i in range(1,max_iter+1):\n",
    "        Tm,V,_=NLBP(A,w,v,i+1)\n",
    "        s=np.zeros(i+2)\n",
    "        s[0]=beta\n",
    "        y=np.linalg.lstsq(Tm,s,rcond=None)[0]\n",
    "        xm=x0+V.dot(y)\n",
    "    return xm\n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.06575817,  0.33251067, -0.5729281 , -0.8131979 , -0.54688738,\n",
       "       -0.30520559,  1.13005732,  0.16601282,  1.19785581, -1.1956725 ,\n",
       "        0.07130517,  0.52358126,  0.7828644 , -0.12647668,  0.07331804,\n",
       "       -0.24769921, -0.83677894,  1.63406951,  0.75501444,  0.49200097,\n",
       "       -0.14564908, -1.70791012,  1.48004496,  0.19575725,  0.07903655,\n",
       "       -0.82935719, -1.21511526, -0.17379711, -0.00354716,  1.63978163])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "QMR(A,b,x0,2*n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.06575817,  0.33251067, -0.5729281 , -0.8131979 , -0.54688738,\n",
       "       -0.30520559,  1.13005732,  0.16601282,  1.19785581, -1.1956725 ,\n",
       "        0.07130517,  0.52358126,  0.78286441, -0.12647668,  0.07331804,\n",
       "       -0.24769921, -0.83677894,  1.6340695 ,  0.75501444,  0.49200097,\n",
       "       -0.14564908, -1.70791013,  1.48004496,  0.19575725,  0.07903655,\n",
       "       -0.82935719, -1.21511526, -0.17379711, -0.00354716,  1.63978163])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_ex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BICGSTAB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BICGSTAB(A,b,x0,rt0,max_iter,tol):\n",
    "    n=A.shape[0]\n",
    "    r0=b-A.dot(x0)\n",
    "    p0=b-A.dot(x0)\n",
    "    counter=0\n",
    "    while counter!=max_iter and np.linalg.norm(r0)>tol:\n",
    "        alpha=r0.dot(rt0)/A.dot(p0).dot(rt0)\n",
    "        s=r0-alpha*A.dot(p0)\n",
    "        w=A.dot(s).dot(s)/np.dot(A.dot(s),A.dot(s))\n",
    "        x0=x0+alpha*p0+w*s\n",
    "        r1=s-w*(A.dot(s))\n",
    "        beta=(r1.dot(rt0)/r0.dot(rt0))*alpha/w\n",
    "        p1=r1+beta*(p0-w*A.dot(p0))\n",
    "        r0=r1\n",
    "        p0=p1\n",
    "        counter+=1\n",
    "    return x0,counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "rt0=np.random.rand(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-0.06575817,  0.33251067, -0.5729281 , -0.8131979 , -0.54688738,\n",
       "        -0.30520559,  1.13005732,  0.16601282,  1.19785581, -1.1956725 ,\n",
       "         0.07130517,  0.52358126,  0.78286441, -0.12647668,  0.07331804,\n",
       "        -0.24769921, -0.83677894,  1.6340695 ,  0.75501444,  0.49200097,\n",
       "        -0.14564908, -1.70791013,  1.48004496,  0.19575725,  0.07903655,\n",
       "        -0.82935719, -1.21511526, -0.17379711, -0.00354716,  1.63978163]),\n",
       " 208)"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BICGSTAB(A,b,x0,rt0,400,10e-12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.06575817,  0.33251067, -0.5729281 , -0.8131979 , -0.54688738,\n",
       "       -0.30520559,  1.13005732,  0.16601282,  1.19785581, -1.1956725 ,\n",
       "        0.07130517,  0.52358126,  0.78286441, -0.12647668,  0.07331804,\n",
       "       -0.24769921, -0.83677894,  1.6340695 ,  0.75501444,  0.49200097,\n",
       "       -0.14564908, -1.70791013,  1.48004496,  0.19575725,  0.07903655,\n",
       "       -0.82935719, -1.21511526, -0.17379711, -0.00354716,  1.63978163])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_ex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GOLUB KAHAN BIDIAGOBALISATION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "def GKB(A,u,max_iter):\n",
    "    alpha=np.linalg.norm(A.T.dot(u))\n",
    "    v=A.T.dot(u)/alpha\n",
    "    U=np.zeros([n,max_iter])\n",
    "    U[:,0]=u\n",
    "    V=np.zeros([n,max_iter])\n",
    "    V[:,0]=v\n",
    "    for j in range(max_iter-1):\n",
    "        U[:,j+1]=A.dot(V[:,j])-alpha*U[:,j]\n",
    "        beta=np.linalg.norm(U[:,j+1])\n",
    "        U[:,j+1]=U[:,j+1]/beta\n",
    "        V[:,j+1]=A.T.dot(U[:,j+1])-beta*V[:,j]\n",
    "        alpha=np.linalg.norm(V[:,j+1])\n",
    "        V[:,j+1]=V[:,j+1]/alpha\n",
    "    return U,V\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "e1=np.zeros(n)\n",
    "e1[0]=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "U,V=GKB(A,e1,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.17685116,  0.1287993 , -0.00777384],\n",
       "       [ 0.32878172, -0.18121267,  0.21577009],\n",
       "       [ 0.33163416, -0.09810655, -0.28945425],\n",
       "       [ 0.03527031,  0.26608841,  0.01445434],\n",
       "       [ 0.12693343,  0.08042898, -0.26614434],\n",
       "       [ 0.16459293,  0.12723031,  0.33163904],\n",
       "       [ 0.06480777,  0.22750394,  0.24592992],\n",
       "       [ 0.18788437,  0.05995414, -0.06319952],\n",
       "       [ 0.0009993 ,  0.32679516,  0.24646207],\n",
       "       [ 0.27046388, -0.13311787,  0.0827691 ],\n",
       "       [ 0.03382685,  0.34083927, -0.08083464],\n",
       "       [ 0.18568108,  0.10831678, -0.009207  ],\n",
       "       [ 0.25227782, -0.09731564,  0.1653856 ],\n",
       "       [ 0.0494878 ,  0.26118912, -0.16781919],\n",
       "       [ 0.14976813,  0.07114698,  0.24846052],\n",
       "       [ 0.23668948, -0.07134698,  0.01034884],\n",
       "       [ 0.16717796,  0.12193468, -0.06716603],\n",
       "       [ 0.23388424, -0.02916612, -0.11939931],\n",
       "       [ 0.10276196,  0.18513321, -0.04162669],\n",
       "       [ 0.0381822 ,  0.29939531, -0.39786427],\n",
       "       [ 0.00461565,  0.2895454 ,  0.13622166],\n",
       "       [ 0.1053343 ,  0.17260637, -0.08887579],\n",
       "       [ 0.28034919, -0.08956899, -0.15938482],\n",
       "       [ 0.11685762,  0.11518259,  0.30748464],\n",
       "       [ 0.12352601,  0.08034848, -0.01584022],\n",
       "       [ 0.3245597 , -0.17079447,  0.14222365],\n",
       "       [ 0.15317145,  0.04924054,  0.01308501],\n",
       "       [ 0.012511  ,  0.32095517,  0.14342057],\n",
       "       [ 0.09956531,  0.18453464, -0.18800947],\n",
       "       [ 0.24741345, -0.04698808, -0.1753852 ]])"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}