function mp = mpapi(t,f,tol,center) %MPAPI multivariate polynomial interpolation %v5function [mp,mpfactor] = mpapi(t,f,tol,center) % % MP = MPAPI(T,F[,TOL[,CENTER]]) %v4 %v5% [MP,MPFACTOR] = MPAPI(T,F[,TOL[,CENTER]]) % % returns the least polynomial interpolant to the data % % (T(:,i),F(:,i)), i=1,...,n, % %v5% and, optionally, the cell array MPFACTOR containing all the information %v5% generated during factorization for the construction of an interpolant %v5% at the same data sites to given data F, using the subsequent call %v5% %v5% MP = MPAPI(MPFACTOR,F) % % The program is a close copy of the informal program in [de Boor and Ron, % Computational aspects of multivariate polynomial interpolation, Math.Comp. % 58 (198); 1992; 705--727], but is explicit about the linear ordering of the % sets A_jd := {a in Z_+^d : norm(a)_1 = j} and the resulting computation of % each block of the blocked Vandermonde from its preceding block. The same % linear ordering is used in MPVAL to evaluate the interpolating polynomial % by Horner's method. Also, the coefficients are returned in that ordering. % cb: 28jul97, 9aug97; revert to v4 (ugghh!) 7-13sep97 % Copyright (c) 1997 by C. de Boor %v5if iscell(t), % no need to repeat the factorization, just extract all the %v5 % needed information from the cell array T = MPFACTOR %v5 % generated earlier by the statement %v5 % mpfactor = {center LU p K dimpjd mm gstring}; %v5 center = t{1}; d = length(center); %v5 LU = t{2}; p = t{3}; n = length(p); K = t{4}; k = K(end); %v5 dimpjd = t{5}; mm = t{6}; gstring = t{7}; %v5else, % need to carry out factorization [d,n] = size(t); if nargin<4, center = (t*ones(n,1))/n; end %v4 %v5if nargin<4, center = sum(t,2)/n; end % subtract the CENTER from the data sites ... nones = ones(n,1); t = t - center(:,nones); % and make the point closest to the origin the first point to be considered % (to avoid a zero row later in the interesting part of Vk for k>0 ) % by setting P(1) equal to its index, where P is the n-sequence % keeping track of pivoting. Precisely, at the end of the j-th elimination % step, row P(i) is known to be pivot row for step i=1:j, while rows P(i), % i=j+1:n, are not yet used as pivot rows. p = 1:n; [ignored,i] = min(ones(1,d)*abs(t)); p(1)=i; p(i)=1; %v4 %v5p = 1:n; [ignored,i] = min(sum(abs(t),1)); p(1)=i; p(i)=1; % Use the mechanism for generating Vkp1 from Vk to generate the weight % vector for the scalar product, by adjoining the point e := ones(d,1) : t = [ones(d,1) t]; if nargin<3, tol = 1.e-14; end % Elimination needs a tolerance for % recognizing when a prospective pivot element is zero. % initialize Gauss elimination items: LU = zeros(n,n); % to contain the interesting parts of L and U ; K = zeros(1,n); % to record the pivot degree for each unknown; % the vector P for recording interchanges was set earlier % initialize block Vandermonde: % Keep around only the current block of the Vandermonde, but accumulate the % the various (g_j)_least, in GSTRING. Vk = ones(n+1,1); np=2:n+1; Wk = Vk(np,:); gstring = []; % also compute the initial size of each entry of Vk: iVk = Vk(np,:).^2./Vk(nones,:); % initialize indexing scheme for A_jd = {a in Z_+^d : |a| = j} [mm,blaises,m,pm] = mpnext(d) ; dimpjd = [0 1]; k = 0; kfact = 1; for j=1:n while 1 % determine the next pivot row, generating the next block of % the Vandermonde if necessary. [pivot,ipivot] = ... max(((Wk(p(j:n),:).^2./Vk(ones(n+1-j,1),:))*ones(blaises(d),1))... ./iVk(p(j:n))); %v5 max(sum(Wk(p(j:n),:).^2./Vk(ones(n+1-j,1),:),2)./iVk(p(j:n))); % [pivot,ipivot] = max(sum(abs(Wk(p(j:n),:)),2)./iVk(p(j:n))); if pivot>tol, break, end k = k+1; kfact= kfact*k; % generate the next block % first, get next piece of MM and % update the other combinatorial information, specifically nakd = blaises(d); % = (k-1 + d-1 choose d-1) [mk,blaises,m,pm] = mpnext(blaises,m,pm); mm = [mm mk]; % now BLAISES = ( (k + h-1 choose h-1) : h=1:d+1 ) dimpjd = [dimpjd blaises(d+1)]; % = ( dim Pi_h(R^d) : h=-1:k ) % Initialize the next V-block to be zero ... Vkp1 = zeros(n+1,blaises(d)); % ... then add to it a properly ordered t(i,:)*Vk, i=1:d for i=1:d ti(:) = t(i,:); %v4 Vkp1(:,mk(i,:)) = Vkp1(:,mk(i,:)) + ... ti(:,ones(1,nakd)).*Vk; %v5 reshape(t(i,:,ones(1,nakd)),n+1,nakd).*Vk; end Vk = Vkp1; %v5 iVk = sum(Vk(np,:).^2./Vk(nones,:),2); iVk = (Vk(np,:).^2./Vk(nones,:))*ones(blaises(d),1); % iVk = sum(abs(Vk(np,:)),2); Wk = Vk(np,:); % apply the previous elimination steps to the new block: % for jj=1:j-1 % Wk(p(jj+1:n),:) = Wk(p(jj+1:n),:) - LU(jj+1:n,jj)*Wk(p(jj),:); % end Wk(p,:) = (tril(LU,-1)+eye(n,n))\Wk(p,:); end if ipivot>1 % exchange the relevant row indices: p(j-1+[ipivot,1]) = p(j-1+[1,ipivot]); LU(j-1+[ipivot,1],:) = LU(j-1+[1,ipivot],:); end K(j) = k; % record the degree of g_jleast i=1:n; (Wk(p(j),:)./Vk(1,:))/kfact; % = coefficients of g_jleast %v5 LU(i,j) = sum(Wk(p(i),:).*ans(nones,:),2); LU(i,j) = (Wk(p(i),:).*ans(nones,:))*ones(blaises(d),1); i=j+1:n; LU(i,j) = LU(i,j)/LU(j,j); Wk(p(i),:) = Wk(p(i),:) - LU(i,j)*Wk(p(j),:); gstring = [gstring ans]; end mm = [mm m(ones(d,1),:)]; % the most recent M is stored (d-fold, in order % to conform to MM's size) in the last part of MM, for % possible use, in degree-raising, integration, etc. %v5if nargout>1 % store the information obtained during factorization %v5 % in the cell array MPFACTOR for later use. %v5 mpfactor = {center LU p K dimpjd mm gstring}; %v5end %v5end [df,nf] = size(f); if nf~=n error('The number of data values does not equal the number of data sites.') end % for j=1:n-1 % i=j+1:n; % a(i) = a(i) - LU(i,j)*a(j); % end % a(n) = a(n)/LU(n,n); % for i=n-1:-1:1 % a(i) = (a(i) - LU(i,i+1:n)*a(i+1:n))/LU(i,i); % end a = (triu(LU)\((tril(LU,-1)+eye(n,n))\f(:,p).')).'; % a = diag(LU).*a; % having suppressed division of Wk(p(j),:) by LU(j,j) % earlier, no need now to multiply a by diag(LU) . coefs = zeros(df,dimpjd(k+2)); gend = 0; %v5coefs = zeros(df,dimpjd(end)); gend = 0; for j=1:n rangej = dimpjd(1+K(j))+1:dimpjd(2+K(j)); gbeg = gend+1; gend = gend + length(rangej); coefs(:,rangej) = coefs(:,rangej) + a(:,j)*gstring(gbeg:gend); end mp = mpmak(coefs,d,k,center,dimpjd,mm);