wolffd@0: function [D,P] = dijk(A,s,t)
wolffd@0: %DIJK Shortest paths from nodes 's' to nodes 't' using Dijkstra algorithm.
wolffd@0: % [D,p] = dijk(A,s,t)
wolffd@0: %     A = n x n node-node weighted adjacency matrix of arc lengths
wolffd@0: %         (Note: A(i,j) = 0   => Arc (i,j) does not exist;
wolffd@0: %                A(i,j) = NaN => Arc (i,j) exists with 0 weight)
wolffd@0: %     s = FROM node indices
wolffd@0: %       = [] (default), paths from all nodes
wolffd@0: %     t = TO node indices
wolffd@0: %       = [] (default), paths to all nodes
wolffd@0: %     D = |s| x |t| matrix of shortest path distances from 's' to 't'
wolffd@0: %       = [D(i,j)], where D(i,j) = distance from node 'i' to node 'j' 
wolffd@0: %     P = |s| x n matrix of predecessor indices, where P(i,j) is the
wolffd@0: %         index of the predecessor to node 'j' on the path from 's(i)' to 'j'
wolffd@0: %         (use PRED2PATH to convert P to paths)
wolffd@0: %       = path from 's' to 't', if |s| = |t| = 1
wolffd@0: %
wolffd@0: %  (If A is a triangular matrix, then computationally intensive node
wolffd@0: %   selection step not needed since graph is acyclic (triangularity is a 
wolffd@0: %   sufficient, but not a necessary, condition for a graph to be acyclic)
wolffd@0: %   and A can have non-negative elements)
wolffd@0: %
wolffd@0: %  (If |s| >> |t|, then DIJK is faster if DIJK(A',t,s) used, where D is now
wolffd@0: %   transposed and P now represents successor indices)
wolffd@0: %
wolffd@0: %  (Based on Fig. 4.6 in Ahuja, Magnanti, and Orlin, Network Flows,
wolffd@0: %   Prentice-Hall, 1993, p. 109.)
wolffd@0: 
wolffd@0: % Copyright (c) 1998-2001 by Michael G. Kay
wolffd@0: % Matlog Version 5 22-Aug-2001
wolffd@0: 
wolffd@0: % Input Error Checking ******************************************************
wolffd@0: error(nargchk(1,3,nargin));
wolffd@0: 
wolffd@0: [n,cA] = size(A);
wolffd@0: 
wolffd@0: if nargin < 2 | isempty(s), s = (1:n)'; else s = s(:); end
wolffd@0: if nargin < 3 | isempty(t), t = (1:n)'; else t = t(:); end
wolffd@0: 
wolffd@0: if ~any(any(tril(A) ~= 0))       % A is upper triangular
wolffd@0:    isAcyclic = 1;
wolffd@0: elseif ~any(any(triu(A) ~= 0))   % A is lower triangular
wolffd@0:    isAcyclic = 2;
wolffd@0: else                             % Graph may not be acyclic
wolffd@0:    isAcyclic = 0;
wolffd@0: end
wolffd@0: 
wolffd@0: if n ~= cA
wolffd@0:    error('A must be a square matrix');
wolffd@0: elseif ~isAcyclic & any(any(A < 0))
wolffd@0:    error('A must be non-negative');
wolffd@0: elseif any(s < 1 | s > n)
wolffd@0:    error(['''s'' must be an integer between 1 and ',num2str(n)]);
wolffd@0: elseif any(t < 1 | t > n)
wolffd@0:    error(['''t'' must be an integer between 1 and ',num2str(n)]);
wolffd@0: end
wolffd@0: % End (Input Error Checking) ************************************************
wolffd@0: 
wolffd@0: A = A';    % Use transpose to speed-up FIND for sparse A
wolffd@0: 
wolffd@0: D = zeros(length(s),length(t));
wolffd@0: if nargout > 1, P = zeros(length(s),n); end
wolffd@0: 
wolffd@0: for i = 1:length(s)
wolffd@0:    j = s(i);
wolffd@0:    
wolffd@0:    Di = Inf*ones(n,1); Di(j) = 0;
wolffd@0:    
wolffd@0:    isLab = logical(zeros(length(t),1));
wolffd@0:    if isAcyclic ==  1
wolffd@0:       nLab = j - 1;
wolffd@0:    elseif isAcyclic == 2
wolffd@0:       nLab = n - j;
wolffd@0:    else
wolffd@0:       nLab = 0;
wolffd@0:       UnLab = 1:n;
wolffd@0:       isUnLab = logical(ones(n,1));
wolffd@0:    end
wolffd@0:    
wolffd@0:    while nLab < n & ~all(isLab)
wolffd@0:       if isAcyclic
wolffd@0:          Dj = Di(j);
wolffd@0:       else	% Node selection
wolffd@0:          [Dj,jj] = min(Di(isUnLab));
wolffd@0:          j = UnLab(jj);
wolffd@0:          UnLab(jj) = [];
wolffd@0:          isUnLab(j) = 0;
wolffd@0:       end
wolffd@0:       
wolffd@0:       nLab = nLab + 1;
wolffd@0:       if length(t) < n, isLab = isLab | (j == t); end
wolffd@0:       
wolffd@0:       [jA,kA,Aj] = find(A(:,j));
wolffd@0:       Aj(isnan(Aj)) = 0;
wolffd@0:             
wolffd@0:       if isempty(Aj), Dk = Inf; else Dk = Dj + Aj; end
wolffd@0:       
wolffd@0:       if nargout > 1, P(i,jA(Dk < Di(jA))) = j; end
wolffd@0:       Di(jA) = min(Di(jA),Dk);
wolffd@0:       
wolffd@0:       if isAcyclic == 1       % Increment node index for upper triangular A
wolffd@0:          j = j + 1;
wolffd@0:       elseif isAcyclic == 2   % Decrement node index for lower triangular A
wolffd@0:          j = j - 1;
wolffd@0:       end
wolffd@0:    end
wolffd@0:    D(i,:) = Di(t)';
wolffd@0: end
wolffd@0: 
wolffd@0: if nargout > 1 & length(s) == 1 & length(t) == 1
wolffd@0:    P = pred2path(P,s,t);
wolffd@0: end