% FIGURE 2.4.1
      function YOUT = DLSQSP(X,N,XD,YD,M,XOUT,NOUT)
%
%  FUNCTION DLSQSP CALCULATES A NATURAL CUBIC SPLINE WITH KNOTS AT
%    X(1),...,X(N) WHICH IS THE LEAST SQUARES FIT TO THE DATA POINTS
%    (XD(I),YD(I)), I=1,...,M, AND EVALUATES THIS SPLINE AT THE OUTPUT
%    POINTS XOUT(1),...,XOUT(NOUT).
%
%  ARGUMENTS
%
%             ON INPUT                          ON OUTPUT
%             --------                          ---------
%
%    X      - A VECTOR OF LENGTH N CONTAINING
%             THE SPLINE KNOTS.
%
%    N      - THE NUMBER OF KNOTS.
%             (N >= 3).
%
%    XD     - A VECTOR OF LENGTH M CONTAINING
%             THE X-COORDINATES OF THE DATA
%             POINTS.
%
%    YD     - A VECTOR OF LENGTH M CONTAINING
%             THE Y-COORDINATES OF THE DATA
%             POINTS.
%
%    M      - THE NUMBER OF DATA POINTS.
%
%    XOUT   - A VECTOR OF LENGTH NOUT CONTAINING
%             THE X-COORDINATES AT WHICH THE
%             CUBIC SPLINE IS EVALUATED.  THE
%             ELEMENTS OF XOUT MUST BE IN
%             ASCENDING ORDER.
%
%    YOUT   -                                   A VECTOR OF LENGTH NOUT.
%                                               YOUT(I) CONTAINS THE
%                                               VALUE OF THE SPLINE
%                                               AT XOUT(I).
%
%    NOUT   - THE NUMBER OF OUTPUT POINTS.
%
%-----------------------------------------------------------------------
      ZERO = 0.0;
      for J=1:N
         Y(J) = 0.0;
      end
      for J=1:N
         Y(J) = 1.0;
%                             CALCULATE PHI(J,X), NATURAL CUBIC SPLINE
%                             WHICH IS EQUAL TO ONE AT KNOT X(J) AND
%                             ZERO AT OTHER KNOTS.  THEN SET
%                                 A(I,J) = PHI(J,XD(I)), I=1,...,M
         AJ = DSPLN(X,Y,N,ZERO,ZERO,XD,M);
         for I=1:M
            A(I,J) = AJ(I);
         end
         Y(J) = 0.0;
      end
%                             CALL DLLSQR TO MINIMIZE NORM OF A*Y-YD
      Y = DLLSQR(A,M,N,YD);
%                             LEAST SQUARES SPLINE IS
%                              Y(1)*PHI(1,X) + ... + Y(N)*PHI(N,X).
%                             EVALUATE SPLINE AT XOUT(1),...,XOUT(NOUT)
      YOUT = DSPLN(X,Y,N,ZERO,ZERO,XOUT,NOUT);