The Algorithm Summarized
You are going to manage a matrix of shortest paths. Call it SP, and SP[i][j], at the end of the algorithm, will contain the shortest path from node i to node j. You initialize SP to be the adjacency matrix for a graph. If there is no edge from i to j, then initialize SP[i][j] to be infinity or an appropriately high sentinel value. You should also initialize SP[i][i] to be zero.I'll now summarize the algorithm.
 For each node i in the graph:
 For each pair of nodes x, y in the graph:
 Check to see if SP[x][i] + SP[i][y] is less than SP[x][y].
 If so, what we've discovered is a shorter path from x to y than the one we currently know, and that path goes through node i.
 In that case, update SP[x][y] to equal SP[x][i] + SP[i][y].
 For each pair of nodes x, y in the graph:
I'm not proving here that the FloydWarshall algorithm works, by the way. I'm just telling you that it works, and how to do it.
An Example
This is the same example as in the Wikipedia page (at least as of March, 2016. If Wikipedia changes, go ahead and use the Wayback Machine to make it match up). Here's the graphYou'll note first that it has negative edges. That's ok, as long as there are no negative cycles in the graph (which there aren't). Now, we're going to work through the algorithm, and what I'll do is at each step, show you SP both in graphical form and as a matrix. I'm going to sentinelize the lack of an edge with the value 100. For the weights in this problem, 100 may as well be infinity.
1 2 3 4  1  0 100 2 100 2  4 0 3 100 3  100 100 0 2 4  100 1 100 0 
Step 1: i = 1. We start with i = 1. We next have to iterate through every pair of nodes (x,y), and test to see if the path from x to y through node i is shorter than SP[x][y]. You should see that we can safely any pair where x=i or y=i.
Let's go ahead an enumerate all of the (x,y) pairs, and compare SP[x][y] to SP[x][1] + SP[1][y]. I've colored the "winners" red in each case:
x  y  SP[x][y]  SP[x][1] + SP[1][y] 
2 2 2 3 3 3 4 4 4 
2 3 4 2 3 4 2 3 4 
0 3 100 100 0 2 1 100 0 
104 2 104 200 98 200 200 98 200 
As you can see, there are two places where SP is improved. I'll update the drawing and the matrix below. I've colored the changed edges/values in red, and I colored node 1 in green to show that it was the intermediate node in these paths:
1 2 3 4  1  0 100 2 100 2  4 0 2 100 3  100 100 0 2 4  100 1 98 0 
Step 2: i = 2. Let's make the same table as before, only now with i = 2:
x  y  SP[x][y]  SP[x][2] + SP[2][y] 
1 1 1 3 3 3 4 4 4 
1 3 4 1 3 4 1 3 4 
0 2 100 100 0 2 100 98 0 
104 102 200 104 102 200 3 1 99 
There are three places where SP is improved. I'll update the drawing and the matrix below:
1 2 3 4  1  0 100 2 100 2  4 0 2 100 3  100 100 0 2 4  3 1 1 0 
Step 3: i = 3. Here's our table:
x  y  SP[x][y]  SP[x][3] + SP[3][y] 
1 1 1 2 2 2 4 4 4 
1 2 4 1 2 4 1 2 4 
0 100 100 4 0 100 3 1 0 
98 98 0 102 102 4 101 101 3 
Once again there are two places where SP is improved. I'll update the drawing and the matrix below:
1 2 3 4  1  0 98 2 0 2  4 0 2 4 3  100 100 0 2 4  3 1 1 0 
The Last Step: i = 4. Here's our table:
x  y  SP[x][y]  SP[x][4] + SP[4][y] 
1 1 1 2 2 2 3 3 3 
1 2 3 1 2 3 1 2 3 
0 98 2 4 0 2 100 100 0 
3 1 1 7 3 5 5 1 3 
We're done  the final drawing and matrix are below. As you can see, three values were changed, and there are no more big values on the graph at all.
1 2 3 4  1  0 1 2 0 2  4 0 2 4 3  5 1 0 2 4  3 1 1 0 
The matrix now has your allpairs shortest paths. If any of the diagonal entries are negative, then your graph has negative cycles, and the matrix is invalid.
Exploration #1: TheTips problem from Topcoder, FloydWarshall vs DFS.
Statement
This is the 500point problem from Qualification Round 1B of the Topcoder Open competition from 2015. Here is the problem statement. Here's a summary, in case Topcoder's surveys are down.
 Fred's parents have hidden N easter eggs in their house. The eggs are numbered 0 to N1, because Fred's parents understand that the world is zeroindexed.
 You are given a vector of integers called probability. Its size is N.
 Without any clues, Fred's probability of finding egg i is probability[i] percent.
 Inside of some eggs, Fred's parents have written on paper the location of some other eggs.
 The information about these eggs is stored in a vector of strings called clues.
 For each i, j, if clues[i][j] is 'Y', then egg i contains the location of egg j. If clues[i][j] is 'N', then egg i does not contain the location of egg j.
 N is a number between 1 and 50.
Examples
#  clues  probability  Answer  Comments 
0  { "Y" } 
{ 50 } 
0.5  Pretty straightforward. One egg. 50% chance of finding it. 
1 
{ "YN", "NY" } 
{ 100, 50 } 
1.5  The clues are worthless. 100% chance of egg 0, and 50% chance of egg 1. 
2 
{ "YYY", "NYY", "NNY" } 
{ 100, 0, 0 } 
3.0  Fred is finding egg 0, and with the clues, will also find eggs 1 and 2 
3 
{ "NNN", "NNN", "NNN" } 
{ 0, 0, 0 } 
0.0  Fred's parents are mean. 
4 
{ "NYYNYYNNNN", "NNNNYNNNYN", "YNNYYYYYNN", "YYNYNNNNYN", "NYNNNNNNNY", "YNYYNNYNNY", "NYNNYYYYYY", "NYYYYNNNNN", "YYNYNNYYYN", "NNYYNYNYYY"} 
{ 11, 66, 99, 37, 64, 45, 21, 67, 71, 62 } 
9.999891558057332  I'm not working through this one by hand. 
3 
{ "NNY", "NNY", "NNN" } 
{ 50, 50, 1 } 
1.7525  I'll go over this one below. 
Solution
This problem has a solution based on connectivity, that you can hack up with DFS or with FloydWarshall. The solution works as follows: Each clue i will have an expected probability of being found. Call that p_{i}.
 Start by setting all p_{i} to zero.
 For each node j, and each node i reachable from j set p_{i} to p_{i} + (1  p_{i})*probability[j].
 At the end, sum up all p_{i}, and that's your answer.
 We start with p = { 0, 0, 0 }.
 We note that nodes 0 and 2 are reachable from node 0, whose probability is 0.50. This makes p become { 0.5, 0, 0.5 }.
 Nodes 1 and 2 are reachable from node 1, whose probability is also 0.50. This makes p become { 0.5, 0.5, 0.75 }. The 0.75 is calculated as .5 + (1  .5)*.5.
 Finally, node 2 is reachable from node 2, whose probability is 0.01. That means that p_{2} is set to 0.75 + (10.75)*0.01 = 0.7525, and the final value of p is { 0.5, 0.5, 0.7525 }.
 So the final expected value is the sum of the values in p, which is 1.7525.
I have also programmed a FloydWarshall solution for this in TheTipsFloydIf.cpp. It turns the input into an adjacency matrix named C of zero's and ones, and then uses FloydWarshall to turn the matrix into a connectivity matrix  if C[i][j] is equal to one, then node j is reachable from node i. The connectivity matrix is then used to calculate the values of p above. Here's the FloydWarshall code, the calculation of p, and the calculation of the final return value:
/* The FloydWarshall Calculation */ for (v = 0; v < C.size(); v++) { for (i = 0; i < C.size(); i++) { for (j = 0; j < C.size(); j++) { if (!C[i][j] && C[i][v] && C[v][j]) C[i][j] = 1; } } } /* Calculate the values of p from the probabilities and reachability: */ p.resize(C.size(), 0); for (i = 0; i < C.size(); i++) { x = probability[i]; x /= 100.0; for (j = 0; j < C.size(); j++) { if (C[i][j]) p[j] += ((1  p[j]) * x); } } /* Calculate the final return value */ x = 0; for (i = 0; i < C.size(); i++) x += p[i]; return x; } 
I have also programmed up a solution that uses bit arithmetic instead of the if statement. This solution doesn't do anything fancy with the bits; it simply replaces the if statement with bit arithmetic. This is in TheTipsFloydBits.cpp. Here's the FloydWarshall part:
for (v = 0; v < C.size(); v++) { for (i = 0; i < C.size(); i++) { for (j = 0; j < C.size(); j++) { C[i][j] = (C[i][v] & C[v][j]); } } } 
Think about the tradeoffs with the two pieces of code. With the if statement, you don't evaluate C[i][v] or C[j][v] when C[i][j] is one. That saves some work. On the flip side, if statements involve comparisons and jumps, which require more instructions than the simple bit arithmetic.
Here's the relevant assembly code for each (Intel core i5, compiled with O2). I find it hard to read, but it does back up the comments above.
if (!C[i][j] && C[i][v] && C[v][j]) C[i][j] = 1;  C[i][j] = (C[i][v] & C[v][j]); 
One last tweak  suppose I want to move the c[i][v] lookup out of the inner loop. Why? Because if it's false, I won't have to do the inner loop at all. This is in TheTipsFloydBits2.cpp:
for (v = 0; v < C.size(); v++) { for (i = 0; i < C.size(); i++) { if (C[i][v]) { for (j = 0; j < C.size(); j++) { C[i][j] = C[v][j]; } } } } 
Remember this  it will come in handy for your lab. All of these work on the topcoder examples:
UNIX> sh sh3.2$ for i in 0 1 2 3 4 5 ; do > tipfwif $i  tail n 1 > tipfwbits $i  tail n 1 > tipfwbits2 $i  tail n 1 > done 0.5000000000 0.5000000000 0.5000000000 1.5000000000 1.5000000000 1.5000000000 3.0000000000 3.0000000000 3.0000000000 0.0000000000 0.0000000000 0.0000000000 9.9998915581 9.9998915581 9.9998915581 1.7525000000 1.7525000000 1.7525000000 sh3.2$ exit UNIX>I've timed these on my MacBook (2.4 Ghz Core i5), and the results make sense. Because these matrices are dense, C[i][j] is true a lot of the time, and that makes the if statement exit early a lot. On the flip side, C[i][v] is also true a lot, so moving it out of the inner loop doesn't eliminate the inner loop very much.
Using DFS Instead of FloydWarshall
Instead of FloydWarshall, you can simply do a DFS for every starting node. To do that, I have programmed four variants. They all make an adjacency list representation of the graph in Adj, and then do a DFS from each node v to calculate which nodes are reachable from v. Here are the variants: TheTipsDFSR.cpp: This is a standard recursive
DFS, which you call on every starting node (clearing out the visited vector each time):
void DFS::DoDFS(int n, int from) { int j; if (visited[n]) return; p[n] += ((1  p[n]) * pr[from]); visited[n] = 1; for (j = 0; j < Adj[n].size(); j++) { DoDFS(Adj[n][j], from); } }

TheTipsDFSR2.cpp: This does an extra level of
pruning. It checks visited *before* you make the DFS call:
for (j = 0; j < Adj[n].size(); j++) { if (!visited[Adj[n][j]]) DoDFS(Adj[n][j], from); }
 TheTipsDFSStack1.cpp: This one manages a stack explicitly instead of doing recursion. It calls push_back() to push a value onto the stack and pop_back() to pop it.
 TheTipsDFSStack2.cpp: This one also manages a stack, but it keeps the current index in sp, rather than calling push_back() and pop_back().
The timing results may surprise you:
They surprise me, as I thought that the stack version with the index would be the best, but instead, it's the pruned recursion. I don't have time to explore it.
Packing Bits Together for InstructionLevel Parallelism
Now, recall this version of the FloydWarshall algorithm, which used bit arithmetic, and moved the checking of C[i][v] out of the inner loop (in TheTipsFloydBits2.cpp):for (v = 0; v < C.size(); v++) { for (i = 0; i < C.size(); i++) { if (C[i][v]) { for (j = 0; j < C.size(); j++) { C[i][j] = C[v][j]; } } } } 
Instead of storing one bit per char, what if we packed each char with 8 bits? In other words, we can store an entire row C with ceil(C.size()/8) char's. Let's use an example:
01234567  0  10001001 1  01100000 2  10101000 3  00010000 4  10001000 5  00100100 6  00000010 7  00100001In our previous FloydWarshall implementations, each entry of the adjacency matrix consumed a char. What if instead we had each row of the adjacency matrix be a single byte, and each entry consumes a bit? Look and see what that lets us do in that inner loop. Suppose v equals 0 and i equals 2. You'll note that C[i][v] equals one. That means that for each value of j, we will OR C[i][j] with C[v][j] and set the result in C[i][j]. Here's row i (2) of C:
2  10101000And here's row v (0):
0  10001001Since both rows are bytes, we can do the OR in one instruction, which will set row 2 to be:
2  10101001Do you see how that can improve performance? Instead of doing 8 OR operations in the inner loop, you only have to do one!
The graph below shows how we use this trick for 1, 2, 4, 8 and 16byte data types. In that latter case, I'm using the Microsoft Intrinsic instruction _mm_or_si128(), which compiles to code that uses the 128bit OR instruction from the Intel SSE2 instruction set. Make sure you pay attention to the units on the Y axis  the savings are HUGE.
This last graph shows the original FloydWarshall implementation, the fastest DFS implementation, and the 128bit SIMD implementation  that's an amazing difference, is it not?
An aside for 2017  a student, Gregory Rouleau, got a bit enthused pursuing the understanding of performance in The Tips. I'm including his writeups. I love it when students lock into a research question!
All Pairs Max Flow Paths
This is too much fun  read on.The problem that we're going to solve is, given an adjacency matrix of a weighted, directed graph, compute the flow of the maximum flow path between every pair of nodes. We're going to constrain flow values so that they are between 0 and 255. I've drawn an example graph on the left, and its allpairs maxflow path graph on the right. I've drawn the edges that have been updated with better flow values in red:
As you can imagine, we'll solve the problem with the FloydWarshall algorithm. To do so, I've written a driver program in APFlowMain.cpp. I have three programs that solve the problem:
 apflowfw, implemented in APFlowFW.cpp, solves it with the FloydWarshall algorithm.
 apflowd, implemented in APFlowDijkstra.cpp, solves it by applying Dijkstra's algorithm to every starting node (this is similar to my Network Flow lecture notes in CS302, if you remember).
 apflowsimd, solves it with FloydWarshall, but it uses SIMD instructions to speed up the inner loop. You'll get to write this in your lab.
usage: APFlow nodes seed( to read from stdin) print(YNH) 
If you give it a seed, it will simply create a random graph. Otherwise, it reads the graph from standard input. If print is "Y", then it will print the graph's adjacency matrix before the flow calculation, and the allpairs maxflow path matrix after the flow calculation. For example, the graph above was generated with a seed of 12.
UNIX> apflowfw 3 12 Y Adjacency Matrix: 255 37 64 93 255 52 98 62 255 Flow Matrix: 255 62 64 93 255 64 98 62 255 UNIX>There is a header file in APFlow.h:
class APFlow { public: int N; uint8_t *Adj; uint8_t *Flow; void CalcFlow(); }; 
N is the number of nodes. Adj is the adjacency matrix, where Adj[i*N+j] stores the weight of the edge from i to j. Flow will be the flow matrix. You will create it by calling CalcFlow(). When CalcFlow() finishes then Flow[i*N+j] will contain the flow of the maximum flow path from i to j. APFlowMain.cpp creates an instance of the class, and initializes N and Adj. It then calls CalcFlow(). If print is 'Y', then it prints out the two matrices. If print is 'H', then it prints out the MD5 hash of Flow. That allows us to check correctness without having to look at giant matrices:
UNIX> apflowfw 5 1 Y Adjacency Matrix: 255 8 41 52 19 44 255 1 11 5 27 44 255 49 60 29 12 108 255 115 53 29 11 29 255 Flow Matrix: 255 44 52 52 52 44 255 44 44 44 53 44 255 52 60 53 44 108 255 115 53 44 52 52 255 UNIX> apflowfw 5 1 H 3B65A77BC185B3A15603EF2268873233 UNIX> apflowd 5 1 H 3B65A77BC185B3A15603EF2268873233 UNIX>BTW, the diagonal elements are always given a weight of 255.
The FloydWarshall solution in APFlowFW.cpp is straightforward, and is very much like our other FloydWarshall solutions. What we're doing now, is saying that if the path from i to j through v has a higher flow value than what we currently know (which is in Flow[i*N+j]), then update it:
#include "APFlow.h" void APFlow::CalcFlow() { int i, j, v, f; for (i = 0; i < N*N; i++) Flow[i] = Adj[i]; for (v = 0; v < N; v++) { for (i = 0; i < N; i++) { for (j = 0; j < N; j++) { /* f is the flow from i to j through v */ f = (Flow[i*N+v] < Flow[v*N+j]) ? Flow[i*N+v] : Flow[v*N+j]; if (f > Flow[i*N+j]) Flow[i*N+j] = f; } } } } 
I won't go through the Dijkstra solution, but suffice it to say that it is slower than FloydWarshall:
UNIX> sh sh3.2$ for i in 160 320 480 960 ; do csh c "time apflowfw $i 1 N"; done 0.016u 0.000s 0:00.01 100.0% 0+0k 0+0io 0pf+0w 0.118u 0.000s 0:00.12 91.6% 0+0k 0+0io 0pf+0w 0.381u 0.001s 0:00.38 100.0% 0+0k 0+0io 0pf+0w 2.905u 0.004s 0:02.91 99.6% 0+0k 0+0io 0pf+0w sh3.2$ for i in 160 320 480 960 ; do csh c "time apflowd $i 1 N"; done 0.039u 0.000s 0:00.04 75.0% 0+0k 0+0io 0pf+0w 0.193u 0.000s 0:00.19 100.0% 0+0k 0+0io 0pf+0w 0.516u 0.001s 0:00.51 100.0% 0+0k 0+0io 0pf+0w 3.036u 0.004s 0:03.04 99.6% 0+0k 0+0io 0pf+0w sh3.2$ exit UNIX>Now, as with the connectivity problem above, we can use SIMD to make this faster. I'm going to illustrate by annotating the FloydWarshall solution:
#include "APFlow.h" void APFlow::CalcFlow() { int i, j, v, f; for (i = 0; i < N*N; i++) Flow[i] = Adj[i]; for (v = 0; v < N; v++) { for (i = 0; i < N; i++) { /* Create a vector alli, which is 16 instances of Flow[i*N+v] */ for (j = 0; j < N; j += 16) { /* Load Flow[v*N+j] through Flow[v*N+j+15] to vector vv */ /* Create fv, which is the flow from i to each of j through j+15 through v. This is simply the min of alli and vv. */ /* Load Flow[i*N+j] through Flow[i*N+j+15] to vector iv */ /* Create rv, which is the max of each value of fv and iv. */ /* Store rv into Flow[i*N+j] through Flow[i*N+j+15] */ } } } } 
Let me illustrate. Suppose that N is 16, and suppose that row i is:
30 95 101 255 104 106 69 106 11 109 73 75 108 7 15 37Let's suppose that v is 2, and let's also suppose that row v is:
119 66 255 62 80 4 47 123 48 99 22 35 100 31 13 99Now, the flow from i to v is 101. The current flow from i to 0 is 30. However, I can get from i to v with a flow of 101, and from v to 0 with a flow of 119. That means that I can get from i to 0 through v with a flow of 101, and I should update entry zero in row i to be 101.
On the flip side, the flow from i to 1 is 95, and the flow from v to 1 is 66, so I can't improve my flow by going through v.
This should give you a flavor of how row i gets updated using the flow from i to v and row v. Now, let's do it with SIMD:
alli 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101 vv 119 66 255 62 80 4 47 123 48 99 22 35 100 31 13 99 fv 101 66 101 62 80 4 47 101 48 99 22 35 100 31 13 99 iv 30 95 101 255 104 106 69 106 11 109 73 75 108 7 15 37 rv 101 95 101 255 104 106 69 106 48 109 73 75 108 31 15 99You'll store rv as the new value of row i.
That is your job  to write APFlowSIMD.cpp, which solves the allpairs maxflow paths problem using FloydWarshall and SIMD. You may assume that N is always a multiple of 16. My program simply exits if it is not:
UNIX> apflowsimd 17 1 H For SIMD, N must be a multiple of 16 UNIX>The only SIMD routines that you need to solve this problem are _mm_set1_epi8(), _mm_min_epu8() and _mm_max_epu8().
How much faster is this? Oh my.....