function retval = gamm_rnd(a,b,m,n)
% Usage: xgamm = gamm_rnd(a,b,m,n)
% Description: This program accepts scalar arguments a
% and b, and index arguments m and n  and returns an array
% of size mxn of random numbers generated from a gamma
% distribution with parameters:
%   a = alpha
%   b = beta
% in the standard description so that the distribution has
% mean alpha*beta and variance alpha*beta^2. The
% last index n is optional, with a default value of 1.
% The builtin uniform random number generator rand is
% used so resetting the seed is accomplished via rand.

% Contact: T. Shores at tshores1@math.unl.edu

if (nargin < 3)
  error('gamm_rnd: need three arguments');
end
if ((a < 0) | (b < 0))
  error('gamm_rnd: parameters a, b must be positive');
end
if (m < 1)
  error('gamm_rnd: row number m must be positive integer');
end
if (nargin == 4)
  if (n < 1)
    error('gamm_rnd: column number must be positive integer');
  end
else
  n = 1;
end
retval = rand(m,n);
retval = gamm_inv(retval,a,b);