--- a/main/queueing/inst/qnopenmulti.m	Wed Oct 17 20:38:24 2012 +0000
+++ b/main/queueing/inst/qnopenmulti.m	Thu Oct 18 10:03:41 2012 +0000
@@ -39,7 +39,7 @@
 ## @item S
 ## @code{@var{S}(c,k)} is the mean service time of class @math{c}
 ## customers on the service center @math{k} (@code{@var{S}(c,k)>0}).
-## For FCFS nodes, average service times must be class-independent.
+## For FCFS nodes, mean service times must be class-independent.
 ##
 ## @item V
 ## @code{@var{V}(c,k)} is the average number of visits of class @math{c}
@@ -47,10 +47,11 @@
 ##
 ## @item m
 ## @code{@var{m}(k)} is the number of servers at service center
-## @math{k}. Valid values are @code{@var{m}(k) < 1} to denote a delay
-## center (@math{-/G/\infty}), and @code{@var{m}(k)==1} to denote
-## a single server queueing center (@math{M/M/1}--FCFS,
-## @math{-/G/1}--LCFS-PR or @math{-/G/1}--PS).
+## @math{k}. If @code{@var{m}(k) < 1}, center @math{k} is a delay center
+## (@math{-/G/\infty}); if @code{@var{m}(k)==1} then center @math{k} is
+## either (@math{M/M/1}--FCFS,
+## @math{-/G/1}--LCFS-PR or @math{-/G/1}--PS). Finally, if @code{@var{m}(k)>1},
+## center math{k} is @math{M/M/m} with @code{@var{m}(k)} identical servers.
 ##
 ## @end table
 ##
@@ -89,7 +90,7 @@
 ## Author: Moreno Marzolla <marzolla(at)cs.unibo.it>
 ## Web: http://www.moreno.marzolla.name/
 function [U R Q X] = qnopenmulti( lambda, S, V, m )
-  if ( nargin < 3 || nargin > 4 )
+  if ( nargin < 2 || nargin > 4 )
     print_usage();
   endif
   isvector(lambda) && all(lambda > 0) || \
@@ -97,32 +98,64 @@
   lambda = lambda(:)'; # make lambda a row vector
   C = length(lambda);
   K = columns(S);
-  [C,K] == size(S) || \
-      error( "rows(S) must be equal to length(lambda)" );
+  (ismatrix(S) && [C,K] == size(S) ) || \
+      error( "S must be a %d x %d matrix", C, K);
   all(all( S > 0 )) || \
       error( "S(c,k) must be > 0" );
-  [C,K] == size(V) || \
-      error( "V must be a matrix of the same size as S" );
-  all( all(V>= 0) ) || \
-      error( "V must be >= 0 " );
+  if ( nargin < 3 || isempty(V) )
+    V = ones(size(S));
+  else
+    (ismatrix(V) && [C,K] == size(V)) || \
+	error( "V must be a %d x %d matrix", C, K);
+    all( all(V>=0) ) || \
+	error( "V must be >= 0 " );
+  endif
 
-  if ( nargin < 4 )
+  D = S .* V;  # Service demands: D(c,k) = S(c,k) * V(c,k)
+
+  if ( nargin < 4 || isempty(m) )
     m = ones(1,K);
   else
-    ( isvector( m ) && length(m) == K && all( m <= 1 ) ) || \
-        error( "m must be less than or equal to ones(1,K)" );
+    ( isvector( m ) && length(m) == K ) || \
+        error( "m must be a vector wiht %d elements", K);
     m = m(:)'; # make m a row vector
   endif
 
-  D = S .* V; # Service demands: D(c,k) = S(c,k) * V(c,k)
-  max(lambda * D) < 1 || \
-      error( "processing capacity exceeded" );
+  ## If there are M/M/k servers with k>=1, compute the maximum
+  ## processing capacity
+  m(m<1) = -1; # avoids division by zero in next line
+  [Umax kmax] = max(lambda * D ./ m);
+  (Umax < 1) || \
+      error( "Processing capacity exceeded at center %d", kmax );
 
   U = R = Q = X = zeros(C,K);
+  X = diag(lambda)*V; # X(c,k) = lambda(c)*V(c,k);
+
+  ## Compute utilizations (together with resp. time and queue lenghts,
+  ## but only for IS nodes)
+  for k=1:K
+    for c=1:C
+      if ( m(k) > 1 ) # M/M/m-FCFS
+	[U(c,k)] = qnmmm( X(c,k), 1/S(c,k), m(k) );
+      elseif ( m(k) == 1 ) # M/M/1-FCFS
+	[U(c,k)] = qnmm1( X(c,k), 1/S(c,k) );
+      else # -/G/1-PS
+  	[U(c,k) R(c,k) Q(c,k)] = qnmminf( X(c,k), 1/S(c,k) );
+      endif
+    endfor
+  endfor
+  ## Adjust response times and queue lengths for FCFS queues
+  k_fcfs = find(m>=1);
+  for c=1:C
+    Q(c,k_fcfs) = U(c,k_fcfs) ./ ( 1 - sum(U(:,k_fcfs),1) );
+    R(c,k_fcfs) = Q(c,k_fcfs) ./ X(c,k_fcfs); # Use Little's law
+  endfor
+
+#{
   i_delay  = find(m<1);
   i_single = find(m==1);
   U = diag(lambda)*D; # U(c,:) = lambda(c)*D(c,:);
-  X = diag(lambda)*V; # X(c,:) = lambda(c)*V(c,:);
+
   
   ## delay centers
   R(:,i_delay) = S(:,i_delay);
@@ -133,8 +166,17 @@
     R(c,i_single) = S(c,i_single) ./ ( 1 - sum(U(:,i_single),1) );
     Q(c,i_single) = U(c,i_single) ./ ( 1 - sum(U(:,i_single),1) );
   endfor
+#}
 endfunction
 %!test
+%! fail( "qnopenmulti([1 1], [.9; 1.0])", "exceeded at center 1");
+%! fail( "qnopenmulti([1 1], [0.9 .9; 0.9 1.0])", "exceeded at center 2");
+%! qnopenmulti([1 1], [.9; 1.0],[],2); # should not fail, M/M/2-FCFS
+%! qnopenmulti([1 1], [.9; 1.0],[],-1); # should not fail, -/G/1-PS
+%! fail( "qnopenmulti(1./[2 3], [1.9 1.9 0.9; 2.9 3.0 2.9])", "exceeded at center 2");
+%! qnopenmulti(1./[2 3], [1 1.9 0.9; 0.3 3.0 1.5],[],[1 2 1]); # should not fail
+
+%!test
 %! V = [1 1; 1 1];
 %! S = [1 3; 2 4];
 %! lambda = [3/19 2/19];

--- a/main/queueing/inst/qnopensingle.m	Wed Oct 17 20:38:24 2012 +0000
+++ b/main/queueing/inst/qnopensingle.m	Thu Oct 18 10:03:41 2012 +0000
@@ -101,31 +101,33 @@
     print_usage();
   endif
   ( isscalar(lambda) && lambda>0 ) || \
-      error( "lambda must be a positive number" );
-  lambda = lambda(:)';
+      error( "lambda must be a positive scalar" );
   ( isvector( S ) && all(S>0) ) || \
       error( "S must be a vector >0" );
   S = S(:)';
   K = length(S);
-  ( isvector( V ) && length(V)==K && all(V>=0) ) || \
-      error( "V must be a vector >=0 and of the same length as S" );
-  V = V(:)';
+  if ( nargin < 3 || isempty(V) )
+    V = ones(size(S));
+  else
+    ( isvector( V ) && length(V)==K && all(V>=0) ) || \
+	error( "V must be a vector with %d elements >=0", K );
+    V = V(:)';
+  endif
 
-  if ( nargin < 4 )
-    m = ones(1,K);
+  if ( nargin < 4 || isempty(m) )
+    m = ones(size(S));
   else
     (isvector(m) && (length(m) == K)) || \
 	error( "m must be a vector of %d elements", K);
     m = m(:)';
-    [err m] = common_size(m,S);
-    ( err == 0 ) || \
-        error( "m and S are not of common size" );
   endif
 
-  ## Compute maximum processing capacity
-  lambda_sat = m ./ max( S .* V );
-  (lambda <= lambda_sat) || \
-      error( "Processing capacity exceeded (lambda must be less than %f)", lambda_sat );
+  ## If there are M/M/k servers with k>=1, compute the maximum
+  ## processing capacity
+  m(m<1) = -1; # avoids division by zero in next line
+  [Umax kmax] = max( lambda * S .* V ./ m );
+  (Umax < 1) || \
+      error( "Processing capacity exceeded at center %d", kmax );
 
   l = lambda*V; # arrival rates
 
@@ -156,7 +158,11 @@
 %! m = [1 1];
 %! fail( "qnopensingle(lambda,S,V,m)","m must be a vector");
 %! V = [1 1 1 1];
-%! fail( "qnopensingle(lambda,S,V)","same length as S");
+%! fail( "qnopensingle(lambda,S,V)","3 elements");
+%! fail( "qnopensingle(1.0, [0.9 1.2], [1 1])", "exceeded at center 2");
+%! fail( "qnopensingle(1.0, [0.9 2.0], [1 1], [1 2])", "exceeded at center 2");
+%! qnopensingle(1.0, [0.9 1.9], [1 1], [1 2]); # should not fail
+%! qnopensingle(1.0, [0.9 1.9], [1 1], [1 0]); # should not fail
  
 %!test
 %! # Example 34.1 p. 572 Bolch et al.
@@ -208,11 +214,14 @@
 
 ## Check if processing capacity is properly accounted for
 %!test
-%! lambda = [1.1];
+%! lambda = 1.1;
 %! V = 1;
 %! m = [2];
 %! S = [1];
-%! [U1 R1 Q1 X1] = qnopensingle(sum(lambda),S,V,m); 
+%! [U1 R1 Q1 X1] = qnopensingle(lambda,S,V,m); 
+%! m = [-1];
+%! lambda = 90.0;
+%! [U1 R1 Q1 X1] = qnopensingle(lambda,S,V,m); 
 
 %!demo
 %! lambda = 3;