[neltmx,factor,ierr,n,nelt,isym,ia,ja,a,f,soln,dsum,itmp,idiag]=srmgen(neltmx,factor,ierr,n,nelt,isym,ia,ja,a,f,soln,dsum,itmp,idiag);

function [neltmx,factor,ierr,n,nelt,isym,ia,ja,a,f,soln,dsum,itmp,idiag]=srmgen(neltmx,factor,ierr,n,nelt,isym,ia,ja,a,f,soln,dsum,itmp,idiag);
persistent dummy i icol inum irow iseed k nl ; 

if isempty(dummy), dummy=0; end;
if isempty(i), i=0; end;
if isempty(icol), icol=0; end;
if isempty(inum), inum=0; end;
if isempty(irow), irow=0; end;
if isempty(iseed), iseed=0; end;
if isempty(k), k=0; end;
if isempty(nl), nl=0; end;
dummy = 16381.0;
[iseed ,dummy]=rand(dummy);
ierr = 0;
for i = 1 : n;
idiag(i) = 0;
dsum(i) = -1.0e0;
end; i = fix(n+1);
dummy = 0.0;
nelt = 0;
for icol = 1 : n;
nl = fix(n + 1 - icol);
inum =fix(fix((fix(rand(dummy).*nl)+1)./3));
[nl,inum,itmp]=ismpl(nl,inum,itmp);
for irow = 1 : inum;
nelt = fix(nelt + 1);
if( nelt>neltmx )
ierr = 1;
return;
end;
ia(nelt) = fix(n + 1 - itmp(irow));
ja(nelt) = fix(icol);
if( ia(nelt)==icol )
idiag(icol) = fix(nelt);
else;
a(nelt) = -rand(dummy);
dsum(icol) = dsum(icol) + a(nelt);
if( isym==0 )
nelt = fix(nelt + 1);
if( nelt>neltmx )
ierr = 1;
return;
end;
ia(nelt) = fix(icol);
ja(nelt) = fix(ia(nelt-1));
a(nelt) = a(nelt-1).*factor;
dsum(ja(nelt)) = dsum(ja(nelt)) + a(nelt);
else;
dsum(ia(nelt)) = dsum(ia(nelt)) + a(nelt);
end;
end;
end; irow = fix(inum+1);
if( idiag(icol)==0 )
nelt = fix(nelt + 1);
if( nelt>neltmx )
ierr = 1;
return;
end;
idiag(icol) = fix(nelt);
a(nelt) = 0.0e0;
ia(nelt) = fix(icol);
ja(nelt) = fix(icol);
end;
end; icol = fix(n+1);
for i = 1 : n;
a(idiag(i)) = -1.0001e0.*dsum(i);
end; i = fix(n+1);
for i = 1 : n;
[soln(i) ,dummy]=rand(dummy);
f(i) = 0.0e0;
end; i = fix(n+1);
for k = 1 : nelt;
f(ia(k)) = f(ia(k)) + a(k).*soln(ja(k));
if( isym~=0 && ia(k)~=ja(k) )
f(ja(k)) = f(ja(k)) + a(k).*soln(ia(k));
end;
end; k = fix(nelt+1);
end %subroutine srmgen