"# H = sp.sparse.lil_matrix((m,m)) # dimension of H is m x m \n",
"# V = sp.sparse.lil_matrix((A.shape[0],m))\n",
"# V[:,0] = v # each column of V is a vector v\n",
"\n",
"# for j in range(m):\n",
"# print(\"j = \", j)\n",
"# w = A @ v \n",
"# for i in range(j):\n",
"# tmp = v.T @ w \n",
"# H[i,j] = tmp[0,0]\n",
"# w = w - H[i,j]*v \n",
" \n",
"# H[j,j-1] = norm(w) \n",
"\n",
"# if H[j,j-1] == 0: \n",
"# print(\"Arnoldi breakdown\")\n",
"# m = j\n",
"# v = 0\n",
"# break\n",
"# else:\n",
"# if j < m-1:\n",
"# v = w/H[j,j-1]\n",
"# # for i in range(A.shape[0]):\n",
"# # V[i,j+1] = v[i,0]\n",
"# V[:,j+1] = v\n",
"\n",
"# print(j, \" iterations completed\")\n",
"# print(\"V = \", V.shape)\n",
"# print(\"H = \", H.shape)\n",
"# print(\"v = \", v.shape)\n",
"# print(\"beta = \", beta)\n",
"\n",
"# return V, H, v, beta, j"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def Arnoldi(A, v, m): \n",
"# review the restarted GMRES method (hereafter referred to as GMRES in short), which is a nonsymmetric Krylov subspace solver based on the Arnoldi decomposition procedure\n",
"\n",
"# return the matrix V, the matrix Hp, the vector v, the scalar beta and the number of iterations j. Where V is a matrix of dimension n x m, Hp is an upper Hessenberg matrix of dimension (m+1) x m. Return also the m+1 vector of the basis V. Return h as the m+1,m element of Hp and the vector v as the m+1 element of the basis V.\n",
"\n",
"\n",
"def ArnoldiGMRES(A,v0,m):\n",
" v = v0\n",
" beta = norm(v)\n",
" v = v/beta \n",
" h = sp.sparse.lil_matrix((m,m)) \n",
" v = v/beta\n",
" H = sp.sparse.lil_matrix((m+1,m)) \n",
" V = sp.sparse.lil_matrix((A.shape[0],m+1))\n",
" V[:,0] = v # each column of V is a vector v\n",
"\n",
" for j in range(m):\n",
" # print(\"j = \", j)\n",
" w = A @ v \n",
" for i in range(j):\n",
" tmp = v.T @ w \n",
" h[i,j] = tmp[0,0]\n",
" w = w - h[i,j]*v \n",
" h[j,j-1] = norm(w) \n",
" tmp = v.T @ w # tmp is a 1x1 matrix, so it's O(1) in memory\n",
" H[i,j] = tmp[0,0] \n",
" w = w - H[i,j]*v \n",
" \n",
" H[j+1,j] = norm(w)\n",
"\n",
" if h[j,j-1] == 0: \n",
" if H[j+1,j] == 0:\n",
" print(\"Arnoldi breakdown\")\n",
" m = j\n",
" v = 0\n",
" break\n",
" else:\n",
" v = w/h[j,j-1] \n",
" \n",
" print(\"Arnoldi finished\\n\")\n",
" print(\"h is \", h.shape)\n",
" print(\"v is \", v.shape)\n",
" print(\"m is \", m)\n",
" print(\"beta is \", beta)\n",
" return v, h, m, beta, j # this is not the output required in the paper, I should return the matrix V and the matrix H\n"
" if j < m-1:\n",
" v = w/H[j+1,j]\n",
" V[:,j+1] = v\n",
"\n",
" print(j, \" iterations completed\")\n",
" print(\"V = \", V.shape)\n",
" print(\"H = \", H.shape)\n",
" print(\"v = \", v.shape)\n",
" print(\"beta = \", beta)\n",
"\n",
" return V, H, v, beta, j \n",
" "
]
},
{
@ -71,35 +129,52 @@
},
{
"cell_type": "code",
"execution_count": 82,
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"m = 100\n",
"A = sp.sparse.rand(m,m, density=0.5, format='lil')\n",
"\u001b[0;32m/tmp/ipykernel_264197/102804730.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 16\u001b[0m \u001b[0mPt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mA\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 17\u001b[0m \u001b[0;31m# run the algorithm\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 18\u001b[0;31m \u001b[0mAlgo4\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mPt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtau\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m/tmp/ipykernel_264197/2668036735.py\u001b[0m in \u001b[0;36mAlgo4\u001b[0;34m(Pt, v, m, a, tau, maxit, x)\u001b[0m\n\u001b[1;32m 83\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 84\u001b[0m \u001b[0;31m# split the result into y and gamma\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 85\u001b[0;31m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 86\u001b[0m \u001b[0mgamma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 87\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mIndexError\u001b[0m: invalid index to scalar variable."
# Refers to Algorithm 2 in the paper, it's needed to implement the algorithm 4. It doesn't work yet. Refer to the file testing.ipynb for more details. This function down here is just a place holder for now