% MCDIS Multiple-component discretization procedure.
%
% This routine performs a sequence of discretizations of the
% given weight function (or measure), each discretization being
% followed by an application of the Stieltjes, or Lanczos,
% procedure to produce approximations to the desired recurrence
% coefficients. The fineness of the discretization is
% characterized by a discretization parameter N. The support of
% the continuous part of the weight function is decomposed into
% a given number mc of subintervals (some or all of which may
% be identical). The routine then applies to each subinterval
% an N-point quadrature rule to discretize the weight function
% on that subinterval. The discrete part of the weight function
% (if there is any) is added on to the discretized continuous
% weight function. The sequence of discretizations, if chosen
% judiciously, leads to convergence of the recurrence
% coefficients for the discretized measures to those of the
% given measure. If convergence to within a prescribed accuracy
% eps0 occurs before N reaches its maximum allowed value Nmax,
% then the value of N that yields convergence is output as
% Ncap, and so is the number of iterations, kount. If there is
% no convergence, the routine displays the message "Ncap
% exceeds Nmax in mcdis" prior to exiting.
%
% The choice between the Stieltjes and the Lanczos procedure is
% made by setting the parameter irout equal to 1 in the former,
% and different from 1, in the latter case.
%
% The details of the discretization are to be specified prior
% to calling the procedure. They are embodied in the following
% global parameters:
%
% mc = the number of component intervals
% mp = the number of points in the discrete part of the
% measure (mp=0 if there is none)
% iq = a parameter to be set equal to 1, if the user
% provides his or her own quadrature routine, and
% different from 1 otherwise
% idelta = a parameter whose default value is 1, but is
% preferably set equal to 2, if iq=1 and the user
% provides Gauss-type quadrature routines
%
% The component intervals have to be specified (in the order
% left to right) by a global mcx2 array AB=[[a1 b1];[a2 b2];
% ...;[amc bmc]], where for infinite extreme intervals a1=-Inf
% resp. bmc=Inf. The discrete spectrum (if mp>0) is similarly
% specified by a global mpx2 array DM=[[x1 y1];[x2 y2];...;
% ;[xmp ymp]] containing the abscissae and jumps.
%
% If the user provides his or her own quadrature routine
% "quadown", the routine mcdis must be called with the input
% parameter "quad" replaced by "@quadown", otherwise with
% "quad" replaced by "@quadgp", a general-purpose routine
% provided in the package. The quadrature routine must have
% the form
%
% function xw=quad(N,i)
%
% where N is the number of nodes and i identifies the interval
% to which the routine is to be applied.
%
% The routine mcdis also applies to measures given originally
% in multi-component form.
%
function [ab,Ncap,kount]=mcdis(N,eps0,quad,Nmax)
global mc mp iq idelta irout DM uv AB
f='Ncap exceeds Nmax in mcdis with irout=%2.0f\n';
if N<1, error('Input variable N out of range'), end
iNcap=1; kount=-1;
ab(:,2)=zeros(N,1); b=ones(N,1);
Ncap=floor((2*N-1)/idelta);
while any(abs(ab(:,2)-b)>eps0*abs(ab(:,2)))
b=ab(:,2);
kount=kount+1;
if kount>1, iNcap=2^(floor(kount/5))*N; end
Ncap=Ncap+iNcap;
if Ncap>Nmax
fprintf(f,irout)
return
end
mtNcap=mc*Ncap;
if iq~=1, uv=fejer(Ncap); end
for i=1:mc
im1tn=(i-1)*Ncap;
xw=feval(quad,Ncap,i);
xwm(im1tn+1:im1tn+Ncap,:)=xw;
end
if mp~=0, xwm(mtNcap+1:mtNcap+mp,:)=DM; end
if irout==1
ab=stieltjes(N,xwm);
else
ab=lanczos(N,xwm);
end
end
% Ncap, kount