lukefleed
/
ShfitedPowGMRES
mirror of https://github.com/lukefleed/ShfitedPowGMRES
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

algo 2 still not working. idem for algo4

Browse Source
pull/1/head
Luca Lombardo 2 years ago
parent 2d101f402f
commit 679b17b3ce
2 changed files with 106 additions and 27 deletions
Show Stats Download Patch File Download Diff File

101
src/algo1_testing.ipynb vendored
Unescape Escape View File

@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 35,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@ -15,7 +15,9 @@
"import plotly.express as px\n", "import plotly.express as px\n",
"import plotly.graph_objs as go\n", "import plotly.graph_objs as go\n",
"from scipy.sparse import *\n", "from scipy.sparse import *\n",
"from scipy.sparse.linalg import norm" "from scipy import linalg\n",
"from scipy.sparse.linalg import norm\n",
"from scipy.optimize import least_squares"
] ]
}, },
{ {
@ -204,7 +206,7 @@
" start_time = time.time()\n", " start_time = time.time()\n",
"\n", "\n",
" u = Pt.dot(v) - v \n", " u = Pt.dot(v) - v \n",
" mv = 1 # number of iteration\n", " mv = 1 # number of matrix vector products\n",
" r = sp.sparse.lil_matrix((n,1)) \n", " r = sp.sparse.lil_matrix((n,1)) \n",
" Res = sp.sparse.lil_matrix((len(a),1))\n", " Res = sp.sparse.lil_matrix((len(a),1))\n",
" x = sp.sparse.lil_matrix((n,1)) \n", " x = sp.sparse.lil_matrix((n,1)) \n",
@ -326,7 +328,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 30,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@ -338,32 +340,49 @@
" h = sp.sparse.lil_matrix((m,m))\n", " h = sp.sparse.lil_matrix((m,m))\n",
" print(\"C\")\n", " print(\"C\")\n",
"\n", "\n",
" for j in range(1,m):\n", " for j in range(m):\n",
" w = A.dot(v)\n", " w = A.dot(v)\n",
" for i in range(1,j):\n", " print(\"D\")\n",
" for i in range(j):\n",
" h[i,j] = v.T.dot(w)\n", " h[i,j] = v.T.dot(w)\n",
" w = w - h[i,j]*v\n", " print(\"E\")\n",
" w = w - h[i,j]*v[i]\n",
" print(\"F\")\n",
"\n", "\n",
" h[j+1,j] = norm(w)\n", " h[j+1,j] = norm(w)\n",
" print(\"G\")\n",
"\n", "\n",
" if h[j+1,j] == 0:\n", " if h[j+1,j] == 0:\n",
" print(\"The algorithm didn't converge\")\n",
" m = j\n", " m = j\n",
" v[m+1] = 0\n", " v[m+1] = 0\n",
" break\n", " break\n",
" else:\n", " else:\n",
" v = w**h[j+1,j]\n", " print(\"H\")\n",
" v[j+1] = w**h[j+1,j]\n",
" print(\"I\")\n",
"\n", "\n",
" return v, h, m, beta, j" " return v, h, m, beta, j"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 29,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"A = sp.sparse.rand(100,100, density=0.5, format='lil')\n", "A = sp.sparse.rand(100,100, density=0.5, format='lil')\n",
"v = sp.sparse.rand(100,1, density=1, format='lil')" "v = sp.sparse.rand(100,1, density=1, format='lil')\n",
"m = 100"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"v, h, m, beta, j = Arnoldi(A, v, m)"
] ]
}, },
{ {
@ -372,7 +391,67 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"v, h, m, beta, j = Arnoldi(A, v, 100)" "def Algo4(Pt, v, m, a: list, tau, maxit: int, x):\n",
" \n",
" iter = 1\n",
" mv = 0\n",
" e1 = sp.sparse.lil_matrix((1,n))\n",
" e1[0,0] = 1\n",
" x = sp.sparse.lil_matrix((len(a),1))\n",
" I = sp.sparse.eye(n, n, format='lil')\n",
" res = sp.sparse.lil_matrix((len(a),1))\n",
" r = sp.sparse.lil_matrix((n,1))\n",
" y = sp.sparse.lil_matrix((n,1))\n",
"\n",
" for i in range(len(a)):\n",
" r = ((1-a[i])**a[i])*v - ((1**a[i])*I - Pt).dot(x)\n",
" res[i] = a[i]*norm(r)\n",
"\n",
" def Find_k(res, maxit):\n",
" k = 0\n",
" for i in range(len(a)):\n",
" if res[i] == max(res):\n",
" k = i\n",
" break\n",
" return k\n",
"\n",
" def Find_gamma(res, a, k):\n",
" gamma = sp.sparse.lil_matrix((len(a),1))\n",
" for i in range(len(a)):\n",
" if i != k:\n",
" gamma[i] = (res[i]*a[k])/(res[k]*a[i])\n",
" else:\n",
" gamma[i] = 0\n",
" return gamma\n",
"\n",
"\n",
" while max(res) > tau and iter < maxit:\n",
" k = Find_k(res, maxit)\n",
" gamma = Find_gamma(res, a, k)\n",
" v, h, m, beta, j = Arnoldi((1**a[k])*I - Pt, r, m)\n",
" Hbar = sp.sparse.lil_matrix((m+1,m))\n",
" Hbar[0:m,0:m] = h\n",
" Hbar[m+1,0:m] = e1\n",
"\n",
" mv += j\n",
"\n",
" # solve the least squares problem for Hbar*x = beta*e1\n",
" y = sp.sparse.linalg.least_squares(Hbar, beta*e1)\n",
" res[k] = a[k]*norm(beta*e1 - Hbar*y)\n",
" x[k] = x[k] + v*y[k]\n",
"\n",
" for i in range(len(a)):\n",
" if i != k:\n",
" if res[i] >= tau:\n",
" Hbar[i] = Hbar[k] + ((1-a[i])/a[i] - (1-a[k])/a[k])*I\n",
" z = beta*e1 - Hbar*y\n",
" y = sp.sparse.linalg.solve(Hbar, gamma*beta*e1)\n",
" x = x + v*y\n",
" res[i] = a[i]**a[k]*gamma[i]*res[k]\n",
" \n",
" iter += 1\n",
" \n",
" return x, res, mv\n"
] ]
} }
], ],

32
src/main.py
Unescape Escape View File

@ -10,15 +10,15 @@ import scipy as sp
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import networkx as nx import networkx as nx
import plotly.graph_objs as go from os.path import exists
from scipy.sparse import * from scipy.sparse import *
from scipy.sparse.linalg import norm from scipy.sparse.linalg import norm
from os.path import exists import plotly.graph_objs as go
warnings.simplefilter(action='ignore', category=FutureWarning) warnings.simplefilter(action='ignore', category=FutureWarning)
# some stupid pandas function that doesn't work # some stupid pandas function that doesn't work
class utilities: class Utilities:
# Importing the dataset # Importing the dataset
def load_data(): def load_data():
# Loading the dataset # Loading the dataset
@ -124,8 +124,7 @@ class utilities:
return a return a
class Plotting: class Plotting:
def tau_over_iterations(dataframe): def tau_over_iterations(df):
dataframe = df
x = df['tau'][::-1].tolist() x = df['tau'][::-1].tolist()
y = df['iterations'].tolist() y = df['iterations'].tolist()
@ -154,7 +153,7 @@ class Algorithms:
print("STARTING ALGORITHM 1...") print("STARTING ALGORITHM 1...")
u = Pt.dot(v) - v u = Pt.dot(v) - v
mv = 1 # number of iteration mv = 1 # number of matrix-vector multiplications
r = sp.sparse.lil_matrix((n,1)) r = sp.sparse.lil_matrix((n,1))
Res = sp.sparse.lil_matrix((len(a),1)) Res = sp.sparse.lil_matrix((len(a),1))
x = sp.sparse.lil_matrix((n,1)) x = sp.sparse.lil_matrix((n,1))
@ -168,7 +167,7 @@ class Algorithms:
x = r + v x = r + v
while max(Res) > tau and mv < max_mv: while max(Res) > tau and mv < max_mv:
u = Pt*u # should it be the same u of the beginning? u = Pt*u
mv += 1 mv += 1
for i in range(len(a)): for i in range(len(a)):
@ -182,7 +181,7 @@ class Algorithms:
if mv == max_mv: if mv == max_mv:
print("The algorithm didn't converge in ", max_mv, " iterations") print("The algorithm didn't converge in ", max_mv, " iterations")
else: else:
print("The algorithm converged in ", mv, " iterations") print("The algorithm converged with ", mv, " matrix-vector multiplications executed")
total_time = time.time() - start_time total_time = time.time() - start_time
total_time = round(total_time, 2) total_time = round(total_time, 2)
@ -217,16 +216,16 @@ df = pd.DataFrame(columns=['alpha', 'iterations', 'tau', 'time'])
# Main # Main
if __name__ == "__main__": if __name__ == "__main__":
dataset = utilities.load_data() dataset = Utilities.load_data()
# maximum number of iterations, asked to the user # maximum number of iterations, asked to the user
max_mv = int(input("Insert the maximum number of iterations: ")) max_mv = int(input("\nInsert the maximum number of iterations: "))
G, n = utilities.create_graph(dataset) G, n = Utilities.create_graph(dataset)
P = utilities.create_matrix(G) P = Utilities.create_matrix(G)
d = utilities.dangling_nodes(P,n) d = Utilities.dangling_nodes(P,n)
v = utilities.probability_vector(n) v = Utilities.probability_vector(n)
Pt = utilities.transition_matrix(P, v, d) Pt = Utilities.transition_matrix(P, v, d)
a = utilities.alpha() a = Utilities.alpha()
# run the algorithm for different values of tau from 10^-5 to 10^-9 with step 10^-1 # run the algorithm for different values of tau from 10^-5 to 10^-9 with step 10^-1
for i in range(5,10): for i in range(5,10):
@ -245,4 +244,5 @@ if __name__ == "__main__":
Plotting.tau_over_time(df) Plotting.tau_over_time(df)
# print in the terminal the columns of the dataframe iterations, tau and time # print in the terminal the columns of the dataframe iterations, tau and time
print("Computations done. Here are the results:")
print("\n", df[['iterations', 'tau', 'time']]) print("\n", df[['iterations', 'tau', 'time']])

Write Preview
Loading…
Cancel
Save