function [ d ] = diam( A )
%DIAM  calculates the network diameter and average path length
%
%    DIAM(A) returns the diameter  d  of the graph given by
%       its adjacency matrix A
%
%       Note: If A is reducible, a warning is issued.
%
%       A can be a full matrix, but treatment of sparse 
%       matricies is faster
%
%  Last change:  10/2/2005 - Florian Knorn

if nargin ~= 1
    error('Please input adjacency matrix !');
else
	A = logical(sparse(A));
	n = size(A,1);
	if n ~= size(A,2)
        error('Please input SQUARE matrix !');
    end
end

% some initialisation
nn = n^2; d = 0; warning off;

A = logical(A); % make sure A is a logical
Asum2 = logical(speye(n)); Ad = Asum2; Asum1 = sparse(n,n);
while (nnz(Asum2) ~= nn) && (nnz(Asum2)-nnz(Asum1) ~= 0) && (d <= n)
	d = d + 1; Ad = logical(double(Ad)*A);
	Asum1 = Asum2; Asum2 = Asum1 | Ad;
end

% output warning if graph is not connected
if (nnz(Asum2) - nnz(Asum1) == 0) | (d == n)
	disp(' '); warning off backtrace; 
	warning('Graph is not connected !'); d = d - 1;
end
	
if n <= 2000
	[i,j,s] = find(Asum2-Asum1);
	for k=1:length(i)/2
		disp(sprintf('Longest geodesic from %i to %i',j(k),i(k)));
	end
end