import numpy as np
a=np.zeros((100,100),dtype=int)
cost=np.zeros((100,100),dtype=int)
w=np.zeros((100,100),dtype=int)
def out(A,n):
for i in range(1,n+1):
for j in range(1,n+1):
print "%3d " %A[i][j],
print ""
def edge(i,j):
k=w[i][j]
if k!=0:
edge(i,k)
print "(%d,%d)" %(i,j)
edge(k,j)
else:
print "(%d,%d)" %(i,j)
if __name__ == "__main__":
n=6
a[1][1]=0; a[1][2]=1; a[1][3]=3; a[1][4]=3; a[1][5]=7; a[1][6]=1;
a[2][1]=1; a[2][2]=0; a[2][3]=1; a[2][4]=4; a[2][5]=6; a[2][6]=1;
a[3][1]=1; a[3][2]=1; a[3][3]=0; a[3][4]=1; a[3][5]=3; a[3][6]=8;
a[4][1]=1; a[4][2]=1; a[4][3]=1; a[4][4]=0; a[4][5]=1; a[4][6]=2;
a[5][1]=1; a[5][2]=1; a[5][3]=1; a[5][4]=1; a[5][5]=0; a[5][6]=1;
a[6][1]=1; a[6][2]=1; a[6][3]=1; a[6][4]=1; a[6][5]=1; a[6][6]=0;
for i in range(1,n+1):
cost[i][i]=0
for i in range(1,n):
cost[i][i+1]=a[i][i+1]
cost[1][n]=a[1][n]
for t in range(2,n):
for i in range(1,n-t+1):
j=i+t
x=99l
for k in range(i+1,j):
if x>cost[i][k]+cost[k][j]+a[i][j]:
x=cost[i][k]+cost[k][j]+a[i][j]
w[i][j]=k
cost[i][j]=x
print "a="
out(a,n)
print "cost="
out(cost,n)
print "witness="
out(w,n)
edge(1,n)
#!/usr/bin/python
import sys
import time
import random
import math
def floyd(D, n):
for k in range(n):
for i in range(n):
for j in range(n):
D[i][j] = min(D[i][j], D[i][k] + D[k][j])
def floyd_sub_cubic(D, n, m):
M = n / m
for K in range(M):
baseK = K * m
# run traditional Floyd over a "cross" of m-by-m matrices
for k in range(baseK, baseK + m):
for i in range(baseK, baseK + m):
for j in range(n):
D[i][j] = min(D[i][j], D[i][k] + D[k][j])
for i in range(n):
for j in range(baseK, baseK + m):
D[i][j] = min(D[i][j], D[i][k] + D[k][j])
# compute (a_ir - a_is) for all m-by-l sub-matrices on row K
# compute (b_sj - b_rj) for all l-by-m sub-matrices on column K
# Use lookup table to compute DMM for all other m-by-m sub-matrices
for i in range(M):
baseI = i * m
for j in range(M):
baseJ = j * m
if (baseI != baseK and baseJ != baseK):
DMM(D, m, baseI, baseJ, baseK)
def DMM(D, m, baseI, baseJ, baseK):
for i in range(baseI, baseI + m):
for j in range(baseJ, baseJ + m):
for k in range(baseK, baseK + m):
D[i][j] = min(D[i][j], D[i][k] + D[k][j])
def main():
# if (len(sys.argv) < 4):
# print("Usage: %s <n> <m> <r>" % sys.argv[0])
# return
# n = int(sys.argv[1])
# m = int(sys.argv[2])
# r = int(sys.argv[3])
n=16
m=4
r=4
if (n % m) != 0:
print("n must be divisibe by m")
return
seed = int(time.time())
#seed = 1463481178
random.seed(seed)
print("n = %d, m = %d, M = n / m = %d, r = %d, seed = %d" % (n, m, n / m, r, seed))
for r in range(r):
D1 = [[0 for i in range(n)] for j in range(n)]
D2 = [[0 for i in range(n)] for j in range(n)]
# Full graph with random edge costs
for i in range(n):
for j in range(n):
if not (i == j):
D1[i][j] = random.randint(0, n - 1)
D2[i][j] = D1[i][j]
# Solve using traditional Floyd
floyd(D1, n)
# Solve using sub-cubic Floyd
floyd_sub_cubic(D2, n, m)
# Check the results
for i in range(n):
for j in range(n):
if (D1[i][j] != D2[i][j]):
print("!!! ERROR !!! D1[%d][%d] = %d, D2[%d][%d] = %d" % (i, j, D1[i][j], i, j, D2[i][j]))
return
if __name__ == "__main__":
main()