 
------------------------------------------------------------------
                        ODRPACK 1.3  MANUAL
------------------------------------------------------------------
 
1
                               Reference Guide
 
                                     for
 
                                   ODRPACK
 
            software for weighted orthogonal distance regression
 
 
                                version  1.3
                                  02-04-87
 
 
 
                                Paul T. Boggs
              Optimization Group/Scientific Computing Division
            National Bureau of Standards, Gaithersburg, MD 20899
 
                               Richard H. Byrd
                       Department of Computer Science
                  University of Colorado, Boulder, CO 80309
 
                             Janet R. Donaldson
              Optimization Group/Scientific Computing Division
            National Bureau of Standards, Boulder, CO 80303-3328
 
                             Robert B. Schnabel
                       Department of Computer Science
                  University of Colorado, Boulder, CO 80309
                                     and
              Optimization Group/Scientific Computing Division
            National Bureau of Standards, Boulder, CO 80303-3328
 
 
 
 
                               ***keywords***
                       orthogonal distance regression
                           nonlinear least squares
                             errors in variables
 
 
 
 
                              ***categories***
                                  G2E,I1B1
 
 
 
1    TABLE OF CONTENTS
 
 
       I.  Introduction
 
      II.  Background
 
     III.  Starting Values for BETA and DELTA
 
      IV.  Default Values and Structured Arguments
 
       V.  Subroutine Declaration and Call Statements
 
      VI.  Subroutine Argument Descriptions
           A.  Synopsis
           B.  Detailed Descriptions of
               ODRPACK User-Callable Subroutine Arguments
 
     VII.  Examples
           A.  DODR Example Program, Data and ODRPACK Generated Report
           B.  DODRC Example Program, Data and ODRPACK Generated Report
 
    VIII.  Scaling Algorithms
           A.  Beta Scaling
           B.  Delta Scaling
 
      IX.  Extracting Information from the Work Vectors
           A.  Extracting Information from Vector WORK
           B.  Extracting Information from Vector IWORK
 
       X.  Acknowledgments
 
      XI.  References
 
 
 
1   I.  INTRODUCTION
        ------------
 
    ODRPACK is a portable collection of ANSI '77 Fortran subroutines for
    fitting a model to data.  It is designed primarily for instances
    when the independent as well as the dependent variables have
    significant errors, implementing a highly efficient algorithm for
    solving the weighted orthogonal distance regression problem, i.e.,
    for minimizing the sum of the squares of the weighted orthogonal
    distances between each data point and the curve described by the
    model equation.  It can also be used to solve the ordinary least
    squares problem where all of the errors are attributed to the
    observations of the dependent variable.
 
    ODRPACK is designed to accommodate many levels of user sophistication
    and problem difficulty.
 
    * It is easy to use, providing two levels of user-control of the
      computations, extensive error handling facilities, comprehensive
      printed reports and no size restrictions other than effective
      machine size.
 
    * The necessary derivatives (Jacobian matrices) are approximated
      numerically if they are not supplied by the user.
 
    * The correctness of user-supplied derivatives can be verified by
      the derivative checking procedure provided.
 
    * Both weighted and unweighted analysis can be performed.
 
    * Subsets of the unknowns can be treated as constants with their
      values held fixed at their input values, allowing the user to
      examine the results obtained by estimating subsets of the
      unknowns of a general model without rewriting the model
      subroutine.
 
    * The ODRPACK scaling algorithm automatically accommodates poorly
      scaled problems, in which the model parameters and/or unknown
      errors in the independent variables vary widely in magnitude.
 
    * The trust region Levenberg-Marquardt algorithm implemented by
      ODRPACK has a computational effort per step which is of the same
      order as that required for ordinary least squares, even though the
      number of unknowns estimated in the orthogonal distance regression
      problem is the number of unknown model parameters plus the number
      of independent variables, while the number of unknowns estimated
      in the ordinary least squares problem is simply the number of
      unknown model parameters.
 
    * ODRPACK is portable and is easily used with other Fortran
      subroutine libraries.
 
    The following sections describe ODRPACK in greater detail.  Users
    are directed to section II for a brief description of the orthogonal
    distance regression algorithm.  This section introduces notation and
    provides background material for understanding the remainder of the
    documentation.  Section III describes the need for starting values
    for BETA and DELTA, and section IV describes two features of ODRPACK
    which simplify the user interface with the package.  The information
    in these two sections will be especially important to first time
    users of ODRPACK.  The subroutine declaration and call statements
    are given in section V and the subroutine arguments are defined in
    section VI.  The sample programs shown in section VII can be used as
    templates for creating the user's own program.  The information
    provided in section VIII describes the scaling algorithm and section
    IX describes how the user can extract computed results from the work
    vectors.  The information in these two sections is generally not
    needed by first time users of ODRPACK.
 
 
 
1   II.  BACKGROUND
         ----------
 
    Let
 
    Y(I) = FN(X(I,*)+DELTA(I,*);BETA) - EPSILON(I)                (eq.1)
 
    for I=1,...,N, where
 
    N                                 is the number of observations
                                      (see subroutine argument N);
 
    Y(I), I=1,...,N                   are the observed values of the
                                      dependent variable, where Y(I)
                                      depends on X(I,J), J=1,...,M (see
                                      subroutine argument Y);
 
    FN                                is the function used to predict
                                      values of the dependent variable
                                      (see subroutine argument FUN);
 
    X(I,J), I=1,...,N & J=1,...,M     are the observed values of the
                                      independent variable (see
                                      subroutine argument X);
 
    DELTA(I,J), I=1,...,N & J=1,...,M are the unknown errors in X(I,J)
                                      which are to be estimated (see
                                      subroutine arguments JOB and
                                      WORK);
 
    BETA(K), K=1,...,NP               are the function parameters which
                                      are to be estimated (see
                                      subroutine argument BETA);
 
    EPSILON(I), I=1,...,N             are the unknown errors in Y(I)
                                      which are to be estimated (see
                                      subroutine argument WORK).
 
    The square of the weighted orthogonal distance from the point
    (X(I,*),Y(I)) to the point FN(X(I,*)+DELTA(I,*);BETA) on the curve
    described by the model equation, i.e., the square of the
    observational errors, is given by
 
    R(I)**2 = [FN(X(I,*)+DELTA(I,*);BETA) - Y(I)]**2 +
                                                                  (eq.2)
               M
              SUM [D(I,J)*DELTA(I,J)]**2
              J=1
 
    for I = 1,...,N, where
 
    D(I,J), I=1,...,N & J=1,...,M     are the DELTA weights, which can
                                      be used to compensate for
                                      instances when the precision of
                                      the X observations is different
                                      from that of the Y observations
                                      (see subroutine argument RHO).
 
    The least squares orthogonal distance solution is then that which
    minimizes with respect to BETA and DELTA the sum of the squared
    weighted observational errors,
 
     N
    SUM [ ( W(I) * R(I) ) **2 ]                                   (eq.3)
    I=1
 
    where
 
    W(I), I=1,...,N                   are the observational error
                                      weights, which can be used to
                                      compensate for unequal precision
                                      in the observational errors (see
                                      subroutine argument W).
 
    The solution is found using a trust region Levenberg-Marquardt
    method [Boggs, Byrd and Schnabel (1985)], with scaling used to
    accommodate problems in which estimated values have widely varying
    magnitudes.  The Jacobian matrices, i.e., the matrices of first
    partial derivatives of FN with respect to each BETA and each X, are
    computed at every iteration either by finite differences or by a
    user-supplied subroutine, as specified by subroutine argument JOB
    (see section VI.B).  The iterations are stopped when any one of
    three stopping criteria are met.  Two of these indicate the
    iterations have converged to a solution.  These are "sum-of-squares
    convergence", which indicates that the change in the sum of the
    squared weighted observational errors is sufficiently small, and
    "parameter convergence", which indicates the change in the
    parameters is sufficiently small.  The third stopping criteria is a
    limit on the number of iterations.
 
 
 
1   III.  STARTING VALUES FOR BETA AND DELTA
          ----------------------------------
 
    Starting values for BETA must be provided by the user.  Users
    familiar with the ordinary nonlinear least squares problem are
    generally aware of the importance of obtaining good starting values
    for the estimated function parameters.  It is equally important to
    obtain good starting values for the parameters when using the
    orthogonal distance regression technique.  Good starting values can
    significantly decrease the number of iterations required to find a
    solution; a poor starting value may even prevent the solution from
    being found at all.  Reasonable starting values are often available
    from previous analysis or experiments.  When good starting values
    are not readily available, the user may have to do some preliminary
    analysis to obtain them.  Himmelblau (1970) offers several
    suggestions for obtaining starting values when they are not
    available from other sources.
 
    When using the technique of orthogonal distance regression it is
    also important to have good starting values for the estimated
    errors, DELTA, in the independent variables.  The ODRPACK default is
    to initialize DELTA to zero, which is the most obvious initial value
    for the DELTA's.  (Note that zero starting values for DELTA do not
    cause the scaling problems discussed in section VI that zero
    starting values for BETA cause.) However initializing the DELTA's to
    zero is equivalent to initially assigning all of the errors to the
    dependent variable as is done for ordinary least squares.  While
    initializing the DELTA's to zero is quite adequate in many cases, in
    others it is not.  A plot of the curve described by the model
    function and observed data for the initial parameters may indicate
    whether or not zero starting values for DELTA are reasonable.  Often
    it is visually possible to determine better starting values for the
    DELTA's, especially when an asymptote is involved.  For example, in
    the case of an asymptote, the user may need to initialize some of
    the DELTA's to the horizontal distance to the curve, while leaving
    the other DELTA's initialized to zero in order to obtain a
    reasonable solution.  Proper initialization of DELTA can mean the
    difference between solving a difficult problem and not solving it.
 
 
 
1   IV.  DEFAULT VALUES AND STRUCTURED ARGUMENTS
         ---------------------------------------
 
    ODRPACK uses default values and structured arguments to simplify the
    user interface.  The availability of default values in ODRPACK means
    that the user does not have to be concerned with determining values
    for many of the ODRPACK arguments unless the problem being solved
    requires the use of non-default values.  Structured arguments, which
    exploit the possibly symmetric structure of the independent variable
    data, reduce the amount of storage space required for arguments and
    reduces the work required by the user to initialize those arguments.
 
    DEFAULT VALUES.  Default values have been specified for ODRPACK
    subroutine arguments wherever feasible.  These default values are
    invoked by setting the argument to any negative value.  Arrays with
    default values are invoked by setting the first element of the array
    to a negative value, in which case only the first value of the array
    will ever be used.  This allows a scalar to be used to invoke the
    default values of arrays, thus saving space and the need to declare
    such arrays.
 
    Users are encouraged to invoke the default values of arguments
    wherever possible.  The default values have been found to be
    reasonable for a wide class of problems.  Their use will greatly
    simplify the initial use of ODRPACK for a given problem.  Fine
    tuning of these arguments can then be done later if it is found
    necessary.
 
    STRUCTURED ARGUMENTS.  Structured arguments are included in ODRPACK
    because the properties of the individual elements of the possibly
    multi-column independent variable data are often constant throughout
    a given column of the independent variable or even throughout the
    whole independent variable matrix.  For example, section II
    introduces the DELTA weights, specified by subroutine argument RHO,
    that indicate how DELTA and EPSILON of each observation (X,Y) are to
    be weighted in the weighted orthogonal distance.  If each row of the
    independent variable indicates an hourly temperature reading and
    each column a different day on which the temperature readings were
    taken, then the user would probably want to weight each of the
    DELTA's equally.  If one column of the independent variable
    contained hourly temperature readings and the other hourly humidity
    readings, then the user might want to weight each of the DELTA's in
    the first column the same, and to weight each of the DELTA's in the
    second column the same, but not necessarily want to weight the two
    columns equally.  Of course, in other cases, the user might want to
    weight each of the DELTA's differently.
 
    ODRPACK structured arguments exploit this possible symmetry as
    follows.  If each of the N by M elements of an array describing some
    property of the independent variable are identically equal, then a
    single value can be used to specify all N by M elements.  If the
    values of such an array only vary between columns, then each column
    of the array can be specified by a single value.  Thus, it is only
    necessary to supply all N by M elements of the structured argument
    array when the elements of one or more of the columns must be
    individually specified.
 
    The use of ODRPACK structured arguments is summarized as follows.
 
    Structure of      Encoding of      Accessed      Resulting
    property P to     structured       elements of   assignment of
    be specified      argument A and   structured    property P
    by structured     its leading      argument A
    argument A        dimension LDA
    ---------------   --------------   -----------   -----------------
    Property P        A(1,1) < zero    A(1,1)        P(I,J) = -A(1,1),
    constant          with                           I=1,...,N &
    throughout        LDA = 1 or                     J=1,...,M
    independent       LDA >= N
    variable matrix
 
    Property P        A(1,1) >= zero   A(1,J),       P(I,J) = A(1,J),
    varies only       with             J=1,...,M     I=1,...,N &
    between columns   LDA = 1                        J=1,...,M
    of independent
    variable matrix
 
    Property P        A(1,1) >= zero   A(I,J),       P(I,J) = A(I,J),
    varies between    with             I=1,...,N &   I=1,...,N &
    and within        LDA >= N         J=1,...,M     J=1,...,M
    columns of
    independent
    variable matrix
 
    If the first element of the structured argument is negative, then
    each of the N by M elements described by the argument is set to the
    absolute value of the first element.  In this case, only the first
    element of the structured argument is ever referenced, allowing the
    user to set the N by M elements by inputting only a scalar.  (Note
    that in this case, setting the first element to a negative value
    does not necessarily invoke a default value.)  This feature thus
    saves space and the need to declare the structured argument as an
    array.
 
    If the first element of the structured argument is positive, then
    the way the structured argument will be used to designate the N by M
    values specified by it will depend its leading dimension.  The
    leading dimension of the structured argument can be either exactly
    equal to one, or greater than or equal to N.  When the leading
    dimension is exactly equal to one, the structured argument must be
    passed to ODRPACK as a one by M row vector containing the M values
    used to set each of the M columns.  When the leading dimension is
    greater than or equal to N, the structured argument passed to
    ODRPACK must contain an N by M array of values.
 
 
 
1   V.  SUBROUTINE DECLARATION AND CALL STATEMENTS
        ------------------------------------------
 
    The declaration and call statements for ODRPACK's user-callable
    routines, SODR, SODRC, DODR and DODRC, are given below.  SODR and
    SODRC invoke the single precision version of the code and DODR and
    DODRC invoke the double precision version.
 
    SODR and DODR preset many arguments to their default values and
    therefore have shorter call statements than SODRC and DODRC.  SODRC
    and DODRC have expanded call statements that give the user greater
    control in solving the orthogonal distance regression problem.
 
    The information in this section is provided primarily for reference.
    Users are directed to section VII for example programs using DODR
    and DODRC.  These examples, which use Fortran PARAMETER statements
    to dimension ODRPACK arrays, provide a recommended format for
    creating an ODRPACK driver which will allow future changes to be
    made easily.
 
1   SODR:   Single precision subroutine to compute the weighted
            orthogonal distance regression or ordinary linear or
            nonlinear least squares solution.  Derivatives are either
            supplied by the user or numerically approximated by ODRPACK.
            Control values are preset, and a two-part report of the
            results is automatically generated.
 
                PROGRAM MAIN
                .
                .
                .
                EXTERNAL
               +   FUN,JAC
                INTEGER
               +   N,M,NP,
               +   LDX,
               +   LDRHO,
               +   JOB,
               +   LWORK,IWORK(LIWORK),LIWORK,
               +   INFO
                REAL
               +   X(LDX,M),
               +   Y(N),
               +   BETA(NP),
               +   RHO(LDRHO,M),
               +   WORK(LWORK)
                .
                .
                .
                CALL SODR
               +   (FUN,JAC,
               +   N,M,NP,
               +   X,LDX,
               +   Y,
               +   BETA,
               +   RHO,LDRHO,
               +   JOB,
               +   WORK,LWORK,IWORK,LIWORK,
               +   INFO)
                .
                .
                .
                END
 
1   SODRC:  Single precision subroutine to compute the weighted
            orthogonal distance regression or ordinary linear or
            nonlinear least squares solution.  Derivatives are either
            supplied by the user or numerically approximated by ODRPACK.
            Control values are supplied by the user, and a three-part
            report of the results is optionally generated.
 
                PROGRAM MAIN
                .
                .
                .
                EXTERNAL
               +   FUN,JAC
                INTEGER
               +   N,M,NP,
               +   LDX,IFIXX(LDIFX,M),LDIFX,LDSCLD,
               +   IFIXB(NP),
               +   LDRHO,
               +   JOB,
               +   MAXIT,
               +   IPRINT,LUNRPT,LUNERR,
               +   LWORK,IWORK(LIWORK),LIWORK,
               +   INFO
                REAL
               +   X(LDX,M),SCLD(LDSCLD,M),
               +   Y(N),
               +   BETA(NP),SCLB(NP),
               +   RHO(LDRHO,M),W(N),
               +   SSTOL,PARTOL,
               +   WORK(LWORK)
                .
                .
                .
                CALL SODRC
               +   (FUN,JAC,
               +   N,M,NP,
               +   X,LDX,IFIXX,LDIFX,SCLD,LDSCLD,
               +   Y,
               +   BETA,IFIXB,SCLB,
               +   RHO,LDRHO,W,
               +   JOB,TAUFAC,
               +   SSTOL,PARTOL,MAXIT,
               +   IPRINT,LUNERR,LUNRPT,
               +   WORK,LWORK,IWORK,LIWORK,
               +   INFO)
                .
                .
                .
                END
 
1   DODR:   Double precision subroutine to compute the weighted
            orthogonal distance regression or ordinary linear or
            nonlinear least squares solution.  Derivatives are either
            supplied by the user or numerically approximated by ODRPACK.
            Control values are preset, and a two-part report of the
            results is automatically generated.
 
                PROGRAM MAIN
                .
                .
                .
                EXTERNAL
               +   FUN,JAC
                INTEGER
               +   N,M,NP,
               +   LDX,
               +   LDRHO,
               +   JOB,
               +   LWORK,IWORK(LIWORK),LIWORK,
               +   INFO
                DOUBLE PRECISION
               +   X(LDX,M),
               +   Y(N),
               +   BETA(NP),
               +   RHO(LDRHO,M),
               +   WORK(LWORK)
                .
                .
                .
                CALL DODR
               +   (FUN,JAC,
               +   N,M,NP,
               +   X,LDX,
               +   Y,
               +   BETA,
               +   RHO,LDRHO,
               +   JOB,
               +   WORK,LWORK,IWORK,LIWORK,
               +   INFO)
                .
                .
                .
                END
 
1   DODRC:  Double precision subroutine to compute the weighted
            orthogonal distance regression or ordinary linear or
            nonlinear least squares solution.  Derivatives are either
            supplied by the user or numerically approximated by ODRPACK.
            Control values are supplied by the user, and a three-part
            report of the results is optionally generated.
 
                PROGRAM MAIN
                .
                .
                .
                EXTERNAL
               +   FUN,JAC
                INTEGER
               +   N,M,NP,
               +   LDX,IFIXX(LDIFX,M),LDIFX,LDSCLD,
               +   IFIXB(NP),
               +   LDRHO,
               +   JOB,
               +   MAXIT,
               +   IPRINT,LUNRPT,LUNERR,
               +   LWORK,IWORK(LIWORK),LIWORK,
               +   INFO
                DOUBLE PRECISION
               +   X(LDX,M),SCLD(LDSCLD,M),
               +   Y(N),
               +   BETA(NP),SCLB(NP),
               +   RHO(LDRHO,M),W(N),
               +   SSTOL,PARTOL,
               +   WORK(LWORK)
                .
                .
                .
                CALL DODRC
               +   (FUN,JAC,
               +   N,M,NP,
               +   X,LDX,IFIXX,LDIFX,SCLD,LDSCLD,
               +   Y,
               +   BETA,IFIXB,SCLB,
               +   RHO,LDRHO,W,
               +   JOB,TAUFAC,
               +   SSTOL,PARTOL,MAXIT,
               +   IPRINT,LUNERR,LUNRPT,
               +   WORK,LWORK,IWORK,LIWORK,
               +   INFO)
                .
                .
                .
                END
 
 
 
1   VI.   SUBROUTINE ARGUMENT DESCRIPTIONS
          --------------------------------
 
    VI.A  Synopsis
 
    The arguments of the ODRPACK user-callable subroutines are logically
    grouped as shown below.  Arguments shown in parenthesis (...) are
    not included in the SODR and DODR call statements; SODR and DODR
    automatically preset these variables to the default values given in
    section VI.B.  All other arguments are common to all ODRPACK user-
    callable subroutines.
 
    Argument
      Number   Arguments                            Group Description
    --------   ---------                            -----------------
      1 to 2   FUN,JAC                              Names of user-
                                                    supplied subroutines
                                                    for function and
                                                    Jacobian matrix
                                                    computation
 
      3 to 5   N,M,NP                               Problem size
                                                    specification
 
     6 to 11   X,LDX,(IFIXX,LDIFX,SCLD,LDSCLD)      Independent variable
                                                    information
 
          12   Y                                    Dependent variable
 
    13 to 15   BETA,(IFIXB,SCLB)                    Function parameter
                                                    information
 
    16 to 18   RHO,LDRHO,(W)                        Weights
 
    19 to 20   JOB,(TAUFAC)                         Computation and
                                                    initialization
                                                    control
 
    21 to 23   (SSTOL,PARTOL,MAXIT)                 Stopping criteria
 
    24 to 26   (IPRINT,LUNERR,LUNRPT)               Print control
 
    27 to 30   WORK,LWORK,IWORK,LIWORK              Work vectors and
                                                    returned results
 
          31   INFO                                 Stopping condition
 
 
 
1   VI.B  Detailed Descriptions of
          ODRPACK User-Callable Subroutine Arguments
 
    The arguments of ODRPACK's user-callable subroutines are described
    below in order of their occurance in the call statements.
    Appropriate declaration statements for each argument are shown in
    brackets [...] following the argument name; the character string
    <real> denotes REAL when using single precision subroutines SODR and
    SODRC, and denotes DOUBLE PRECISION when using double precision
    subroutines DODR and DODRC.  Each argument is numbered as shown in
    section VI.A, allowing the user to easily find the definition of a
    specific argument.  In addition, three common characteristics of
    ODRPACK subroutine arguments are flagged in the left margin by the
    argument number.  The flags are:
 
    C  which indicates the argument is only included in the call
       statements for SODRC and DODRC (SODR and DODR will preset these
       variables to their default values);
 
    D  which indicates the argument has a default value which can be
       invoked by setting the argument to any negative value; and
 
    S  which indicates the argument exploits possible symmetry in the
       properties of the independent variables as described in
       section IV.
 
 
                                    NOTE
 
            Substitute REAL for <real> when using SODR and SODRC.
 
      Substitute DOUBLE PRECISION for <real> when using DODR and DODRC.
 
                                    ****
 
 
 
1       1. FUN [EXTERNAL FUN]
 
           The name of the user-supplied subroutine that computes the
           predicted values, F, of the dependent variable given the
           current values of the independent variable, XPLUSD=X+DELTA,
           and the function parameters, BETA.  The subroutine argument
           list and declaration statements must be exactly as shown
           below.
 
                 SUBROUTINE FUN(N,NP,M,BETA,XPLUSD,LDXPD,F,IFLAG)
           C
           C  INPUT ARGUMENTS
           C  (WHICH MUST NOT BE CHANGED BY THIS ROUTINE)
           C
                 INTEGER N,NP,M,LDXPD
                 <real> BETA(NP),XPLUSD(LDXPD,M)
           C
           C  OUTPUT ARGUMENTS
           C
                 <real> F(N)
                 INTEGER IFLAG
 
                 < computations for F(I)=FN(XPLUSD;BETA), I=1,...,N >
 
                 RETURN
                 END
 
           where
 
           INTEGER N
              is the number of observations, i.e., the number of
              points (X,Y).
           INTEGER NP
              is the number of function parameters, i.e., the number
              of BETA's.
           INTEGER M
              is the number of columns of data in the independent
              variable matrix XPLUSD.
           <real> BETA(NP)
              is the singly subscripted array that contains the
              current values of the NP function parameters.
           <real> XPLUSD(LDXPD,M)
              is the doubly subscripted array that contains the
              current value of the N by M matrix of the
              independent variables, XPLUSD = X + DELTA.
           INTEGER LDXPD
              is the leading dimension of array XPLUSD.
           <real> F(N)
              is the singly subscripted array that contains
              the N predicted values of the function given the
              current values of the function parameters and the
              independent variables, F = FN(XPLUSD;BETA).
           INTEGER IFLAG
              is a control value which can be used to stop the
              computations if the user so desires.
              If IFLAG < 0 then
                 upon return from the subroutine INFO will be set
                 to IFLAG (see argument INFO) and the computations
                 will be stopped
              else
                 the computations will continue.
 
 
         2. JAC [EXTERNAL JAC]
 
           The name of the user-supplied subroutine that computes the
           Jacobian matrices, i.e., the matrices of first partial
           derivatives of FN with respect to each BETA and each X.  This
           subroutine must be supplied only when digit C of JOB is
           nonzero (see subroutine argument JOB) although the external
           statement must always be used; when digit C of JOB is zero
           the necessary Jacobian matrices will be computed using
           finite differences.
 
           Note that the logical argument ISODR, which is passed to
           subroutine JAC by ODRPACK, can be used to avoid computing the
           Jacobian matrix with repsect to X.  ISODR will be "false"
           when the fit is by the method of ordinary least squares and
           the derivative with respect to X is not needed.
 
           The subroutine argument list and dimension statements must be
           exactly as shown below.
 
1                SUBROUTINE JAC(N,NP,M,BETA,XPLUSD,LDXPD,
                +               FJACB,LDFJB,ISODR,FJACX,LDFJX,IFLAG)
           C
           C  INPUT ARGUMENTS
           C  (WHICH MUST NOT BE CHANGED BY THIS ROUTINE)
           C
                 INTEGER N,NP,M,LDXPD
                 <real> BETA(NP),XPLUSD(LDXPD,M)
                 LOGICAL ISODR
           C
           C  OUTPUT ARGUMENTS
           C
                 <real> FJACB(LDFJB,NP),FJACX(LDFJX,M)
                 INTEGER IFLAG
 
                 < computations for FJACB(I,K)=first partial derivative
                                               of FN with respect to
                                               BETA(K), K=1,...,NP,
                                               at each observation
                                               I=1,...,N >
 
                 IF (ISODR) THEN
 
                 < computations for FJACX(I,J)=first partial derivative
                                               of FN with respect to
                                               X(I,J), J=1,...,M
                                               at each observation
                                               I=1,...,N >
 
                 END IF
 
                 RETURN
                 END
 
           where
 
           INTEGER N
              is the number of observations, i.e., the number of
              points (X,Y).
           INTEGER NP
              is the number of function parameters, i.e., the number
              of BETA's.
           INTEGER M
              is the number of columns of data in the independent
              variable matrix XPLUSD.
           <real> BETA(NP)
              is the singly subscripted array that contains the
              current values of the NP function parameters.
           <real> XPLUSD(LDXPD,M)
              is the doubly subscripted array that contains the
              current value of the N by M matrix of the
              independent variables, XPLUSD = X + DELTA.
           INTEGER LDXPD
              is the leading dimension of array XPLUSD.
           <real> FJACB(LDFJB,NP)
              is the doubly subscripted array that contains
              the N by NP matrix of derivatives with respect to
              BETA at the current values of the function
              parameters and the independent variables.
           INTEGER LDFJB
              is the leading dimension of array FJACB.
           LOGICAL ISODR
              is a control value which can be used to inhibit the
              computation of the derivatives with respect to X
              when the solution is being computed by ordinary
              least squares and the derivatives with respect to
              X are not needed.
              If ISODR is true then
                 the solution is being computed by ODR and the
                 derivatives with respect to X must be computed
              else
                 the solution is being computed by OLS and the
                 derivatives with respect to X are not needed.
           <real> FJACX(LDFJX,M)
              is the doubly subscripted array that contains
              the N by M matrix of derivatives with respect to X
              at the current values of the function parameters
              and the independent variables, needed only when
              ISODR = true.
           INTEGER LDFJX
              is the leading dimension of array FJACX.
           INTEGER IFLAG
              is a control value which can be used to stop the
              computations if the user so desires.
              If IFLAG < 0 then
                 upon return from the subroutine INFO will be set
                 to IFLAG (see argument INFO) and the computations
                 will be stopped
              else
                 the computations will continue.
 
 
        3. N [INTEGER N]
 
           The number of observations, i.e., the number of points (X,Y).
           (See subroutine arguments X and Y.)
 
 
        4. M [INTEGER M]
 
           The number of columns of data in the independent variable
           matrix X (see subroutine argument X).
 
 
1       5. NP [INTEGER NP]
 
           The number of function parameters, i.e., the number of
           BETA's (see subroutine argument BETA).
 
 
        6. X [<real> X(LDX,M)]
 
           The doubly subscripted array that contains the observed
           values of the N by M matrix of independent variables.
 
 
        7. LDX [INTEGER LDX]
 
           The leading dimension of array X.
 
           LDX must equal or exceed N; values of LDX less than N will
           be treated as an input error.
 
 
    CDS 8. IFIXX [INTEGER IFIXX(LDIFX,M)]
 
           The doubly subscripted array that contains the indicator
           values used to designate whether element X(I,J), I=1,...,N &
           J=1,...,M, of the independent variable matrix is to treated
           as without error, i.e., DELTA(I,J) is to be fixed at zero,
           or whether the error DELTA(I,J) in that observation of the
           independent variable is to be estimated.
 
           By default, all of the independent variables are treated as
           "unfixed", i.e.  the errors DELTA(I,J) are estimated for all
           I=1,...,N & J=1,...,M.  The default value is invoked when
           IFIXX(1,1) is set to any negative value.  Other options for
           specifying IFIXX are described below.
 
           If IFIXX(1,1) >= 0 then
 
              the way IFIXX is used depends on the value of the leading
              dimension of IFIXX, i.e., on LDIFX.
 
              If LDIFX = 1 then
 
                 IFIXX must contain a 1 by M matrix of values, where
                 for J=1,...,M
 
                    if IFIXX(1,J) = 0 then
 
                       X(I,J), I=1,...,N, is treated as exact and
                       DELTA(I,J), I=1,...,N, is fixed at zero
 
                    else
 
                       X(I,J), I=1,...,N, is treated as approximate and
                       DELTA(I,J), I=1,...,N, is estimated.
 
              If LDIFX >= N then
 
                 IFIXX must contain an N by M matrix of values, where
                 for I=1,...,N & J=1,...,M
 
                    if IFIXX(I,J) = 0 then
 
                       X(I,J) is treated as exact and DELTA(I,J) is
                       fixed at zero
 
                    else
 
                       X(I,J) is treated as approximate and DELTA(I,J)
                       is estimated.
 
           If IFIXX(1,1) < 0 then
 
              the default option is invoked, i.e., each observation of
              the independent variable, X(I,J), is treated as being
              measured with error DELTA(I,J) which is estimated as
              described above in section II.  In this case, only the
              first element of IFIXX is ever referenced and IFIXX can
              be a scalar.
 
 
      C 9. LDIFX [INTEGER LDIFX]
 
           The leading dimension of array IFIXX.
 
           LDIFX must exactly equal one or must equal or exceed N;
           values of LDIFX less than one or between two and N-1,
           inclusive, will be treated as an input error.
 
           See subroutine argument IFIXX for further details.
 
 
   CDS 10. SCLD [<real> SCLD(LDSCLD,M)]
 
           The doubly subscripted array that contains the scale values,
           i.e., the expected magnitudes or typical sizes, for
           DELTA(I,J), I=1,...,N & J=1,...,M.  The scale values
           specified for each DELTA must be greater than zero; values
           less than zero will be treated as an input error.
 
           By default, the scale values will be set using the algorithm
           given in section VIII.B.  The default values are invoked when
           SCLD(1,1) is set to any negative value.  Other options for
           specifying SCLD are described below.
 
           If SCLD(1,1) > 0 then
 
              each value of SCLD must be greater than zero and the way
              SCLD is used depends on the value of the leading dimension
              of SCLD, i.e., on LDSCLD.
 
              If LDSCLD = 1 then
 
                 SCLD must contain a 1 by M matrix of values, and the
                 scale of DELTA(I,J), I=1,...,N, is set to SCLD(1,J) for
                 J=1,...,M.
 
              If LDSCLD >= N then
 
                 SCLD must contain an N by M matrix of values, and the
                 scale of DELTA(I,J) is set to SCLD(I,J) for I=1,...,N &
                 J=1,...,M.
 
           If SCLD(1,1) < 0 then
 
              the default option is invoked and each DELTA(I,J) is
              scaled as described in section VIII.B.  In this case, only
              the first element of SCLD is ever referenced and SCLD can
              be a scalar.
 
 
      C 9. LDSCLD [INTEGER LDSCLD]
 
           The leading dimension of array SCLD.
 
           LDSCLD must exactly equal one or must equal or exceed N;
           values of LDSCLD less than one or between two and N-1,
           inclusive, will be treated as an input error.
 
           See subroutine argument SCLD for further details.
 
 
       12. Y [<real> Y(N)]
 
           The singly subscripted array that contains the N observed
           values of the dependent variable.
 
 
       13. BETA [<real> BETA(NP)]
 
           The singly subscripted array that contains the (current)
           values of the NP function parameters.
 
           On input:   BETA must contain initial approximations for the
                       function parameters.  Initial approximations
                       should be chosen with care since poor initial
                       approximations can significantly increase the
                       number of iterations required to find a solution
                       and possibly prevent the solution from being
                       found at all.
 
                       Users who do not provide scale information are
                       strongly encouraged not to use zero as an initial
                       approximation since a zero value can result in
                       incorrect scale value selection by the scaling
                       algorithm (see section VIII).  Setting the
                       initial approximation to the largest magnitude
                       which, for the user's problem, is effectively
                       zero rather than the actual value zero will
                       eliminate scaling problems, possibly producing
                       faster convergence.  For example, if BETA(1)
                       represents change in cost in millions of dollars,
                       then the value 10.0 might be considered
                       "effectively zero", while if BETA(1) represents
                       the change in cost in tens of dollars, then the
                       value 0.01 might be considered "effectively
                       zero."
 
           On return:  BETA contains the "best" estimate of the solution
                       at the time the computations stopped.
 
 
    CD 14. IFIXB [INTEGER IFIXB(NP)]
 
           The singly subscripted array that contains the indicator
           values used to designate whether the corresponding value in
           BETA is to be treated as a fixed constant or is to be
           estimated.
 
           By default, all of the function parameters, BETA, are
           treated as "unfixed", i.e.  each of the BETA(K), K=1,...,NP,
           is estimated.  The default value is invoked when IFIXB(1) is
           set to any negative value.  Other options for specifying
           IFIXB are described below.
 
           If IFIXB(1) >= 0 then
 
              IFIXB must contain a vector of NP values, where
              for K=1,...,NP
 
                 if IFIXB(K) = 0 then
 
                    BETA(K) will be held fixed at its input value
 
                 else
 
                    BETA(K) will be estimated as described above.
 
           If IFIXB(1) < 0 then
 
              the default option is invoked, i.e., all BETA(K),
              K=1,...,NP, will be estimated as described above in
              section II.  In this case, only the first element of
              IFIXB is ever referenced and IFIXB can be a scalar.
 
 
    CD 15. SCLB [<real> SCLB(NP)]
 
           The singly subscripted array that contains the scale values,
           i.e., the expected magnitudes or typical sizes, for BETA(K),
           K=1,...,NP.  The scale values specified for each BETA must
           be greater than zero; values less than zero will be treated
           as an input error.
 
           By default, the scale values will be set using the algorithm
           given in section VIII.A.  The default values are invoked when
           SCLB(1) is set to any negative value.  If SCLB(1) >= 0 then
           SCLB must contain a vector of NP values each greater than
           zero and the scale of BETA(K) is set to SCLB(K) for
           K=1,...,NP.
 
 
     S 16. RHO [<real> RHO(LDRHO,M)]
 
           The doubly subscripted array that contains the values that
           specify the DELTA weights, D, which indicate how the
           DELTA's and EPSILON's of the observation (X,Y) are to be
           weighted in the weighted orthogonal distance, R (see eq.2).
           For example, RHO(I,J) might be the the ratio of the precision
           of the X(I,J) observation to that of the Y(I) observation.
 
           All elements of RHO must be nonzero.
 
           If RHO(1,1) < zero then
 
              only the first element of RHO is ever referenced (in this
              case, RHO can be a scalar) and
 
                 D(I,J) = ABS(RHO(1,1)) for I=1,...,N & J=1,...,M,
 
              i.e., D(I,J) is constant and every DELTA is weighted
              equally with respect to each of the EPSILON's.  When
              ABS(RHO(1,1)) = 1, the DELTA's and EPSILON's are both
              weighted equally, possibly indicating the X and Y
              observations are equally precise.
 
           If RHO(1,1) > zero then
 
              then all elements of RHO must be greater than zero and
              the way RHO is used to specify D depends on the value of
              the leading dimension of RHO, i.e., on LDRHO.
 
              If LDRHO = 1 then
 
                 RHO must contain a 1 by M matrix of values, where
                 for J=1,...,M
 
                    D(I,J) = RHO(1,J), I=1,...,N,
 
                 i.e., each column of D is constant.  In this case,
                 each column of DELTA is weighted equally with respect
                 to EPSILON, possibly reflecting that each observation
                 within a given column of X is equally precise, but that
                 the precision between columns varies.
 
              If LDRHO >= N then
 
                 RHO must contain an N by M matrix of values, where
 
                    D(I,J) = RHO(I,J) for I=1,...,N & J=1,...,M,
 
                 i.e., each element of D is individually specified,
                 possibly indicating that the individual observations
                 of X vary significantly in precision both from each
                 other and from the corresponding observations of Y.
 
 
       17. LDRHO [INTEGER LDRHO]
 
           The leading dimension of array RHO.
 
           LDRHO must exactly equal one or must equal or exceed N;
           values of LDRHO less than one or between two and N-1,
           inclusive, will be treated as an input error.
 
           See subroutine argument RHO for further details.
 
 
    CD 18. W [<real> W(N)]
 
           The singly subscripted array that contains the values that
           specify the observational error weights, which can be used
           to compensate for unequal precision in the observational
           errors (see eq.3).
 
           By default, the observational errors are unweighted, i.e.,
           all of the weights are assumed to be identically equal to
           one.  The default value is invoked when W(1) is set to any
           negative value.  Other options for specifying W are
           described below.
 
           If W(1) >= zero then
 
              W must contain a vector of N values, where all elements
              of W must be greater than or equal to zero, and W(I),
              I=1,...,N, specifies the weight for the observational
              error R(I).  Zero weights eliminate the corresponding
              observation from the analysis.
 
           If W(1) < zero then
 
              the default option is invoked, i.e., the observational
              errors are unweighted.  In this case, only the first
              element of W is ever referenced and W can be a scalar.
 
 
     D 19. JOB [INTEGER JOB]
 
           The value used to specify problem initialization and
           computational methods.  The user has the option of specifying
           four different aspects of the problem specification:
              - whether the fit is to be by orthogonal distance
                regression (ODR) or by ordinary least squares (OLS);
              - whether the user has supplied subroutine JAC to compute
                the necessary Jacobian matrices and whether the user-
                supplied Jacobian matrices should be checked;
              - whether the DELTA's have been initialized by the user;
                and
              - whether the fit is a restart.
 
           By default:
              - the solution will be found by ODR;
              - the derivatives will be computed by finite differences;
              - the DELTA's will be initialized to zero; and
              - the fit will not be a restart.
           The default value is invoked by setting JOB to any value less
           than zero.
 
           Setting JOB = 1 will have the same consequence as JOB = -1
           except that the solution will be found by OLS.
 
           If JOB > 0 then
 
              JOB is assumed to be a 4-digit INTEGER with decimal
              expansion ABCD, where each digit controls a different
              aspect of the problem specification.
 
              Digit A indicates whether the fit is a restart.
 
                      A = 0 indicates fit is not a restart.
 
                      A > 0 indicates fit is a restart - computations
                            will continue from where they left off for
                            another 10 iterations.  If the fit is a
                            restart then the elements of vector WORK
                            must be exactly as returned from a previous
                            call to ODRPACK.
 
              Digit B indicates whether the DELTA's have been
                      initialized by the user.
 
                      B = 0 indicates DELTA's have not been initialized
                            by user - DELTA's will be initialized to
                            zero.
 
                      B > 0 indicates DELTA's have been initialized by
                            user (see subroutine argument WORK).
 
              Digit C indicates whether the user has supplied subroutine
                      JAC to compute the necessary Jacobian matrices and
                      whether the user-supplied Jacobian matrices should
                      be checked.
 
                      C = 0 indicates that the Jacobian matrices are to
                            be computed by finite differences and that
                            subroutine JAC will not be used.
 
                      C = 1 indicates that the user has supplied
                            subroutine JAC to compute the necessary
                            Jacobian matrices (see subroutine argument
                            JAC), and that the results of the user-
                            supplied routine will be checked for
                            correctness.
 
                      C > 1 indicates that the user has supplied
                            subroutine JAC to compute the necessary
                            Jacobian matrices (see subroutine argument
                            JAC), and that the results of the user-
                            supplied routine will not be checked for
                            correctness.
 
              Digit D indicates whether the fit is to be by orthogonal
                      distance regression (ODR) or by ordinary least
                      squares (OLS).
 
                      D = 0 indicates an ODR fit.
 
                      D > 0 indicates an OLS fit.
 
           If JOB < 0 then
 
              the "default" value will be used.
 
 
    CD 20. TAUFAC [<real> TAUFAC]
 
           The value used to specify the initializing factor for the
           trust region radius.  The trust region is the region in
           which the local approximation to the user's function is
           considered to be reliable.  The diameter of this region is
           adaptively chosen at each iteration based on information
           from the previous iteration.  At the first iteration, the
           initial diameter is set to the initializing factor times the
           length of the full Gauss-Newton step at the initial
           estimates.
 
           By default, the initialization factor for the trust region
           radius is one, thus allowing the full Gauss-Newton step to
           be taken at the first iteration if it does, in fact, reduce
           the weighted sum-of-squares.  The default value is invoked
           when TAUFAC is set to any value less than or equal to zero.
 
           A value of TAUFAC greater than zero but less than one may be
           appropriate if, at the first iteration, the computed results
           overflow, or the function parameters, BETA, leave the region
           of interest in parameter space.
 
 
    CD 21. SSTOL [<real> SSTOL]
 
           The value used to specify the stopping tolerance for the
           convergence test based on relative change in the sum of the
           squared weighted observational errors (eq.3).
 
           The "default" sum-of-squares convergence stopping tolerance
           is the square root of machine precision, where machine
           precision is defined as the smallest value e such that 1+e>1
           on the computer being used.  The default value is invoked
           when the user-supplied value for SSTOL is outside the
           interval [e,1).
 
 
    CD 22. PARTOL [<real> PARTOL]
 
           The value used to specify the stopping tolerance for the
           convergence test based on relative change in the estimated
           parameters BETA and DELTA.
 
           By default, the stopping tolerance for parameter convergence
           is (machine precision)**(2/3), where machine precision is
           defined as the smallest value e such that 1+e>1 on the
           computer being used.  The default value is invoked when the
           user-supplied value for PARTOL is outside the interval [e,1).
 
 
    CD 23. MAXIT [INTEGER MAXIT]
 
           The value used to specify the maximum number of iterations
           allowed.
 
           By default, the maximum number of iterations is 50.  The
           default value is invoked when the user-supplied value for
           MAXIT is less than or equal to zero.
 
 
    CD 24. IPRINT [INTEGER IPRINT]
 
           The value used to control the generated computation reports,
           which are divided into three sections:
              - the initial summary
              - the iteration summary and
              - the final summary.
           The choice of content for each of these sections is described
           below.
 
           By default, the computation reports include
              - a "long" initial summary
              - no iteration summary and
              - a "short" final summary
           The default value is invoked when the user-supplied value for
           IPRINT is less than or equal to zero.
 
           If IPRINT > 0 then
 
              IPRINT is assumed to be a 4-digit INTEGER with decimal
              expansion ABCD, where each digit controls a different
              part of the generated reports.
 
              Digit A indicates whether the initial summary will be
                      generated.
 
                      A = 0 indicates the initial summary will not be
                            generated.
 
                      A = 1 indicates a "short" initial summary will be
                            generated which will include
                            - the values N, M and NP, which indicate the
                              problem size, and the number of BETA's
                              actually being estimated;
                            - the control values JOB, TAUFAC, SSTOL,
                              PARTOL, and MAXIT; and
                            - the sum of the squared weighted observa-
                              tional errors, the sum of the squared
                              weighted DELTA's and the sum of the
                              squared weighted EPSILON's at the initial
                              values of BETA and DELTA.
 
                      A > 1 indicates a "long" initial summary will be
                            generated, which includes all the
                            information found in the "short" initial
                            summary and, in addition, includes
                            - a summary of the independent variable
                              data, organized by column;
                            - the first and last observation of the
                              dependent variable and the first and last
                              observational error weight; and
                            - for each function parameter BETA, the
                              initial value, whether or not the
                              parameter is treated as fixed or not,
                              and the scale value.
 
              Digit B indicates whether the iteration summary will be
                      generated.
 
                      B = 0 indicates no iteration summary will be
                            generated.
 
                      B = 1 indicates a "short" 1-line, 68-column
                            iteration summary will be generated every
                            Cth iteration beginning with iteration
                            one.  This summary will list
                            - the number of function evaluations;
                            - the sum of the squared weighted observa-
                              tional errors at the current point;
                            - the actual relative reduction in the
                              sum of the squared weighted observa-
                              tional errors due to the most recently
                              tried step (used to check for sum-of-
                              squares convergence);
                            - the predicted relative reduction in the
                              sum of the squared weighted observa-
                              tional errors due to the most recently
                              tried step (used to check for sum-of-
                              squares convergence);
                            - the ratio of the trust region radius to
                              the norm of the BETA's and DELTA's, which
                              is an upper bound on the relative change
                              in the estimated values possible at the
                              next step (used to check for parameter
                              convergence); and
                            - whether the step was a Gauss-Newton step.
 
                      B > 1 indicates an [NP/3]-line, 125-column
                            iteration summary will be generated every
                            Cth iteration beginning with iteration 1.
                            This summary lists all of the information
                            found in the "short" iteration summary and,
                            in addition, includes
                            - current values of the BETA's.  (Note that,
                              at the last iteration, the values listed
                              for BETA will be those that produced the
                              actual and predicted relative reductions
                              shown only if the most recently tried step
                              did in fact make the fit better.  If not,
                              then the values of BETA are those which
                              produced the best fit.
 
              Digit C indicates the frequency of the iteration summary.
 
                      C = 0 indicates no iteration summary will be
                            generated, even if the value of digit B is
                            nonzero.
 
                      C > 0 indicates an iteration summary will be
                            generated every Cth iteration beginning
                            with iteration one.
 
              Digit D indicates whether the final summary will be
                      generated.
 
                      D = 0 indicates the final summary will not be
                            generated.
 
                      D = 1 indicates a "short" final summary will be
                            generated, which includes
                            - the stopping condition;
                            - the number of iterations and function
                              evaluations at the time the computations
                              stopped;
                            - the condition number of the problem at the
                              time the computations stopped;
                            - the rank deficiency of the model at the
                              time the computations stopped;
                            - the sum of the squared weighted observa-
                              tional errors, the sum of the squared
                              weighted DELTA's and the sum of the
                              squared weighted EPSILON's at the final
                              values of BETA and DELTA; and
                            - the final values of BETA, the first 32
                              values of EPSILON, and the first 32 values
                              of each column of DELTA.
 
                      D > 1 indicates a "long" final summary will be
                            generated, which includes the same
                            information as the "short" final summary
                            except that
                            - the values of all of the EPSILON's and
                              DELTA's are listed.
 
           If IPRINT < 0 then
 
              the default reports will be generated.
 
 
    CD 25. LUNERR [INTEGER LUNERR]
 
           The value used to specify the logical unit number of the
           file to be used for error messages.
 
           By default, the error messages are generated on unit 6.
           The default value is invoked when the user-supplied value for
           LUNERR is less than or equal to zero.
 
           If LUNERR > 0  the error messages will be generated on
                          unit LUNERR.
 
           If LUNERR = 0  no error messages will be generated.
 
           If LUNERR < 0  the "default" unit number will be used.
 
 
    CD 26. LUNRPT [INTEGER LUNRPT]
 
           The value used to specify the logical unit number of the
           file to be used for computation reports.
 
           By default, the computation reports are generated on unit 6.
           The default value is invoked when the user-supplied value for
           LUNRPT is less than or equal to zero.
 
           If LUNRPT > 0  the computation reports will be generated on
                          unit LUNRPT.
 
           If LUNRPT = 0  no computation reports will be generated.
 
           If LUNRPT < 0  the "default" unit number will be used.
 
 
       27. WORK [<real> WORK(LWORK)]
 
           The singly subscripted array used for <real> work space and
           the array in which various computed values are returned.  The
           smallest acceptable dimension of WORK is given below in the
           definition of subroutine argument LWORK.
 
           The work area does not need to be initialized by the user
           unless the user wishes to initialize the DELTA's.
 
           The first N*M locations of WORK contain the values for DELTA.
           An easy way to access these values, either for initialization
           (as indicated by digit B of subroutine argument JOB) or for
           analysis upon return from ODRPACK, is to include in the
           user's program the declaration statements
 
              <real> DELTA(<N>,<M>)
              EQUIVALENCE (WORK(1), DELTA(1,1))
 
           where <N> indicates the first dimension of the array DELTA
                     must be exactly the number of observations, N; and
                 <M> indicates the second dimension of the array DELTA
                     must be exactly the number of columns, M, of the
                     independent variable, X.
 
           This allows the error associated with observation X(I,J) of
           the independent variable matrix to be accessed as DELTA(I,J)
           rather than as WORK(I+(J-1)*N).  The input values of array
           DELTA will be over written by the final estimates of the
           errors in the independent variable matrix when this
           equivalencing method is used.
 
           Other values returned in array WORK can be accessed as
           described below in section IX.
 
           N.B., if the fit is a restart, i.e., if digit A of subroutine
           argument JOB is nonzero, then the elements of vector WORK
           must be exactly as returned from a previous call to ODRPACK.
 
 
       28. LWORK [INTEGER LWORK]
 
           The length of array WORK.  LWORK must equal or exceed
 
           9 + 8*N + 10*N*M + 2*N*NP + 8*NP .
 
           Values of LWORK less than this value will be treated as an
           input error.
 
 
       29. IWORK [INTEGER IWORK(LIWORK)]
 
           The singly subscripted array used for INTEGER work space and
           the array in which various computed values are returned.  The
           smallest acceptable dimension of IWORK is given below in the
           definition of subroutine argument LIWORK.
 
           Certain values returned in array IWORK are of general
           interest and can be accessed as described below in section
           IX.
 
 
       30. LIWORK [INTEGER LIWORK]
 
           The length of array IWORK.  LIWORK must equal or exceed
 
           2*NP + M + 10 .
 
           Values of LIWORK less than this value will be treated as an
           input error.
 
 
       31. INFO [INTEGER INFO]
 
           The argument used to indicate why the computations stopped.
 
           If INFO < 0 then
 
              the computations were stopped by the user in user-supplied
              subroutine FUN or JAC.
 
           If 1 <= INFO <= 3 then
 
              the program converged satisfactorily.  The convergence
              condition met is indicated by the value of INFO as
              follows.
              INFO = 1 :  sum-of-squares convergence
              INFO = 2 :  parameter convergence
              INFO = 3 :  sum-of-squares convergence and
                          parameter convergence
 
           If INFO = 4 then
 
              the program reached the maximun number of iterations
              allowed without meeting one of the convergence
              conditions.
 
           If INFO > 4 and INFO < 10000 then
 
              the results from ODRPACK are questionable.  In this case,
              INFO is a 3-digit INTEGER with decimal expansion ABC,
              where digit C indicates the actual stopping condition,
              and digits A and B are used to indicate what questionable
              conditions were found.
 
              Digit A = 1 indicates
 
                 the ODRPACK Jacobian matrix checking procedure
                 determined that the correctness of the user-supplied
                 Jacobian matrices is questionable.  Errors in user-
                 supplied subroutine JAC could be adversely affecting
                 the least squares solution.  Users are encouraged to
                 examine subroutine JAC and the appropriate ODRPACK
                 generated reports to insure that there is not an error
                 in the user-supplied derivatives.
 
              Digit B = 1 indicates
 
                 the results of user-supplied subroutines FUN and/or JAC
                 are unaffected by changes in the unfixed function
                 parameters (BETA), indicating a probable error in these
                 user-supplied subroutines.
 
              Digit C > 0 indicates
 
                 the actual stopping condition, where
                 if C=1 the program met the sum-of-squares convergence
                        criteria
                 if C=2 the program met the parameter convergence
                        criteria
                 if C=3 the program met the sum-of-squares convergence
                        criteria and the parameter convergence criteria
                 if C=4 the program reached the maximun number of
                        iterations allowed without meeting one of the
                        convergence conditions.
 
           If INFO > 9999 then
 
              fatal errors were detected in the user's input.  In this
              case, INFO is a 5-digit INTEGER with decimal expansion
              ABCDE, where each digit indicates a different error
              condition.
 
              Digit A = 1 indicates an error was detected in the
                          arguments used to specify the problem size.
 
                 When digit A = 1 then
 
                    digit B = 1 indicates N < 1
 
                    digit C = 1 indicates M < 1
 
                    digit D = 1 indicates NP < 1 or NP > N
 
              Digit A = 2 indicates an error was detected in the
                          arguments used to specify array dimensions.
 
                 When digit A = 2 then
 
                    digit B = 1 indicates LDX < N
 
                    digit C > 0 indicates LDIFX, LDSCLD and/or LDRHO are
                                unacceptable (see definitions of LDIFX,
                                LDSCLD and LDRHO for acceptable values)
                                where
                                if C=1 LDIFX is bad
                                if C=2 LDSCLD is bad
                                if C=3 LDIFX & LDSCLD are bad
                                if C=4 LDRHO is bad
                                if C=5 LDIFX & LDRHO are bad
                                if C=6 LDSCLD & LDRHO are bad
                                if C=7 LDIFX, LDSCLD & LDRHO are bad
 
                    digit D = 1 indicates LWORK is too small (see
                                definition of LWORK for smallest
                                acceptable value)
 
                    digit E = 1 indicates LIWORK is too small (see
                                definition of LIWORK for smallest
                                acceptable value)
 
              Digit A = 3 indicates an error was detected in the
                          arguments used to specify scaling and/or the
                          in the arguments used to specify the weights.
 
                 When digit A = 3 then
 
                    digit B = 1 indicates an error in SCLD (see defini-
                                tion of SCLD for reasonable values)
 
                    digit C = 1 indicates an error in SCLB (see defini-
                                tion of SCLB for reasonable values)
 
                    digit D > 1 indicates an error in W, where
                                if D=1 one or more of the elements of
                                       W are invalid (see definition
                                       of W for reasonable values)
                                if D=2 the number of non zero values
                                       in W is less than NP
 
                    digit E = 1 indicates an error in RHO (see defini-
                                tion of RHO for reasonable values)
 
              Digit A = 4 indicates an error was detected in the
                          user-supplied Jacobian matrices.
 
                 When digit A = 4 then
 
                    digit B = 1 indicates an error in the Jacobian
                                matrix with respect to BETA (see the
                                generated error reports, or section IX
                                for locations in IWORK that indicate
                                which derivatives are in error)
 
                    digit C = 1 indicates an error in the Jacobian
                                matrix with respect to X (see the
                                generated error reports, or section IX
                                for locations in IWORK that indicate
                                which derivatives are in error)
 
 
 
1   VII.   EXAMPLES
           --------
 
    The following sample programs use DODR and DODRC to solve exercise I
    on page 521 and 522 of Draper and Smith (1981).  The program calling
    DODR uses the default option of computing the derivatives by finite
    differences, while the program calling DODRC uses analytic
    derivatives.  Note that the results of these two examples are not
    identical, primarily because the DODRC example has "fixed" one
    column of the independent variable.  Finite difference derivatives
    generally cause very little change in the results from those
    obtained using analytic derivatives.
 
    Users are encouraged to extract these examples from the online
    ODRPACK documentation, and to then modify them as necessary to form
    their own ODRPACK drivers.  (Single precision sample programs can be
    easily generated from these two programs by changing all REAL
    variables to DOUBLE PRECISION, and substituting DODR for SODR and
    DODRC for SODRC.) Note especially that by using MAXN, MAXM and MAXNP
    to specify the largest problem the program can solve without
    modification, and by specifying LWORK and LIWORK exactly as shown,
    the user greatly reduces the number of changes which must be made to
    the program in order to solve a larger problem.
 
 
 
1   VII.A  DODR Example Program, Data and ODRPACK Generated Report
 
    User-supplied code for DODR example:
 
          PROGRAM DDRV1
    C
    C  SET PARAMETERS FOR MAXIMUM PROBLEM SIZE HANDLED BY THIS DRIVER
    C  WHERE  MAXN IS THE MAXIMUM NUMBER OF OBSERVATIONS ALLOWED,
    C         MAXM IS THE MAXIMUM NUMBER OF COLUMNS IN THE
    C              INDEPENDENT VARIABLE ALLOWED, AND
    C        MAXNP IS THE MAXIMUM NUMBER OF FUNCTION PARAMETERS
    C              ALLOWED.
    C
          INTEGER MAXN,MAXM,MAXNP
          PARAMETER (MAXN=15,MAXM=5,MAXNP=5)
    C
    C  SET PARAMETERS TO SPECIFY DIMENSIONS OF ARRAYS USED BY DODR
    C
          INTEGER LDX,LDRHO,LWORK,LIWORK
          PARAMETER
         +   (LDX=MAXN,
         +   LDRHO=1,
         +   LWORK = 9 + 8*MAXN + 10*MAXN*MAXM + 2*MAXN*MAXNP + 8*MAXNP,
         +   LIWORK = 2*MAXNP + MAXM + 10)
    C
    C  DECLARE USER-SUPPLIED SUBROUTINES AND
    C  ALL OTHER NECESSARY VARIABLES AND ARRAYS
    C
          EXTERNAL
         +   FUN,JAC
          INTEGER
         +   N,M,NP,
         +   JOB,
         +   IWORK(LIWORK),
         +   INFO
          DOUBLE PRECISION
         +   X(LDX,MAXM),
         +   Y(LDX),
         +   BETA(MAXNP),
         +   RHO(LDRHO,MAXM),
         +   WORK(LWORK)
    C
          OPEN(UNIT=5,FILE='DATA1')
          OPEN(UNIT=6,FILE='REPORT')
    C
    C  READ NUMBER OF OBSERVATIONS
    C       NUMBER OF COLUMNS OF DATA IN THE INDEPENDENT VARIABLE
    C       NUMBER OF PARAMETERS
    C       OBSERVED VALUES OF INDEPENDENT AND DEPENDENT VARIABLES
    C       STARTING VALUES OF FUNCTION PARAMETERS
    C
          READ (5,*) N,M,NP
          READ (5,*) ((X(I,J),I=1,N),J=1,M)
          READ (5,*) (Y(I),I=1,N)
          READ (5,*) (BETA(I),I=1,NP)
    C
    C  SPECIFY DELTA WEIGHTS
    C
          RHO(1,1) = 3.0D0
          RHO(1,2) = 5.0D0
    C
    C  SET CONTROL VALUE TO INVOKE DEFAULT SETTING
    C
          JOB = -1
    C
    C  COMPUTE ODR SOLUTION USING FINITE-DIFFERENCE DERIVATIVES
    C
          CALL DODR
         +   (FUN,JAC,
         +   N,M,NP,
         +   X,LDX,
         +   Y,
         +   BETA,
         +   RHO,LDRHO,
         +   JOB,
         +   WORK,LWORK,IWORK,LIWORK,
         +   INFO)
    C
          END
          SUBROUTINE FUN(N,NP,M,BETA,XPLUSD,LDXPD,F,IFLAG)
    C
    C  INPUT ARGUMENTS
    C  (WHICH MUST NOT BE CHANGED BY THIS ROUTINE)
    C
          INTEGER N,NP,M,LDXPD
          DOUBLE PRECISION BETA(NP),XPLUSD(LDXPD,M)
    C
    C  OUTPUT ARGUMENTS
    C
          DOUBLE PRECISION F(N)
          INTEGER IFLAG
    C
          DO 10 I = 1, N
             F(I) = EXP(-BETA(1)*XPLUSD(I,1)*
         +          EXP(-BETA(2)*
         +               (1.0D0/XPLUSD(I,2) - 1.0D0/620.0D0)))
       10 CONTINUE
    C
          RETURN
          END
 
 
1   User-supplied data (file DATA1):
 
        8    2    2
      109.0   65.0 1180.0   66.0
     1270.0   69.0 1230.0   68.0
      600.0  640.0  600.0  640.0
      600.0  640.0  600.0  640.0
      0.912  0.382  0.397  0.376
      0.342  0.358  0.348  0.376
    0.01155 5000.0
 
 
1   Report generated by DODR example program, run using a Cyber 180/855
    under NOS 2.4.2.
 
 
 
 
    ******************************************************
    * ODRPACK VERSION 1.3 OF 02-04-87 (DOUBLE PRECISION) *
    ******************************************************
 
 
 
 
 
    INITIAL SUMMARY FOR FIT BY METHOD OF ODR
    ========================================
 
 
 
    PROBLEM SIZE:
    -------------
 
    NUMBER OF OBSERVATIONS                                8
    NUMBER OF COLUMNS OF DATA IN INDEPENDENT VARIABLE     2
    NUMBER OF FUNCTION PARAMETERS                         2
    NUMBER OF UNFIXED FUNCTION PARAMETERS                 2
 
 
 
    INDEPENDENT VARIABLE AND DELTA WEIGHT SUMMARY:
    ----------------------------------------------
 
                                  COLUMN   1                COLUMN   2
                              OBS 1        OBS N        OBS 1        OBS N
                  X -    .10900D+03   .68000D+02   .60000D+03   .64000D+03
              FIXED -            NO           NO           NO           NO
      INITIAL DELTA -    .00000D+00   .00000D+00   .00000D+00   .00000D+00
        DELTA SCALE -    .78740D-03   .78740D-03   .15625D-02   .15625D-02
      DELTA WEIGHTS -    .30000D+01   .30000D+01   .50000D+01   .50000D+01
 
 
 
    DEPENDENT VARIABLE AND OBSERVATIONAL ERROR WEIGHT SUMMARY:
    ----------------------------------------------------------
 
                              OBS 1        OBS N
                  Y -    .91200D+00   .37600D+00
    OBS. ERROR WTS. -    .10000D+01   .10000D+01
 
 
 
    FUNCTION PARAMETER SUMMARY:
    ---------------------------
 
           INDEX -                1               2
    INITIAL BETA -    .11550000D-01   .50000000D+04
           FIXED -               NO              NO
      BETA SCALE -    .86580087D+02   .20000000D-03
 
 
 
    CONTROL VALUES AND STOPPING CRITERIA:
    --------------------------------------
 
          *
       JOB    TAUFAC     SSTOL    PARTOL  MAXIT
     -0001   .10D+01   .50D-14   .86D-19     50
 
    *
     A.  FIT IS NOT A RESTART.
     B.  DELTAS ARE INITIALIZED TO ZERO.
     C.  DERIVATIVES ARE COMPUTED BY FINITE DIFFERENCES.
     D.  FIT IS BY METHOD OF ORTHOGONAL DISTANCE REGRESSION.
 
 
 
    INITIAL SUMS OF SQUARES:
    ------------------------
 
    SUM OF SQUARED WEIGHTED OBSERVATIONAL ERRORS     .67662011D+00
    SUM OF SQUARED WEIGHTED DELTAS                   .00000000D+00
    SUM OF SQUARED WEIGHTED EPSILONS                 .67662011D+00
 
 
 
 
 
 
 
    FINAL SUMMARY FOR FIT BY METHOD OF ODR
    ======================================
 
 
 
    STOPPING CONDITION (INFO =      1):
    -----------------------------------
 
    THE RELATIVE CHANGE IN THE SUM OF THE SQUARED
    WEIGHTED OBSERVATIONAL ERRORS IS LESS THAN SSTOL
 
                                CONDITION
          NUMBER OF  NUMBER OF     NUMBER        RANK
         ITERATIONS   FN EVALS  (INVERSE)  DEFICIENCY
                  6         38  .1888D-06           0
 
 
 
    FINAL SUMS OF SQUARES:
    ----------------------
 
    SUM OF SQUARED WEIGHTED OBSERVATIONAL ERRORS     .75382323D-03
    SUM OF SQUARED WEIGHTED DELTAS                   .23542099D-07
    SUM OF SQUARED WEIGHTED EPSILONS                 .75379969D-03
 
 
 
    ESTIMATED BETA(I), I = 1, ..., NP:
    ----------------------------------
 
            INDEX            VALUE -------------->
        1 TO    2    .36579727D-02   .27627327D+05
 
 
 
    ESTIMATED EPSILON(I) AND DELTA(I,*), I = 1, ..., N:
    ---------------------------------------------------
 
        I      EPSILON(I)      DELTA(I,1)      DELTA(I,2)
        1   .16752445D-02   .14086172D-06   .42418798D-06
        2   .20434718D-02   .12838222D-05   .20262810D-05
        3  -.20690085D-01  -.71652291D-06  -.23358824D-04
        4   .24305832D-02   .15047092D-05   .24114481D-05
        5   .72777482D-02   .23393281D-06   .82079313D-05
        6   .40793264D-02   .24162846D-05   .40483552D-05
        7   .13043071D-01   .43337344D-06   .14726727D-04
        8  -.85499649D-02  -.51394680D-05  -.84861063D-05
 
 
 
1   VII.B  DODRC Example Program, Data and ODRPACK Generated Report
 
    User-supplied code for DODRC example:
 
          PROGRAM DDRV2
    C
    C  SET PARAMETERS FOR MAXIMUM PROBLEM SIZE HANDLED BY THIS DRIVER
    C  WHERE  MAXN IS THE MAXIMUM NUMBER OF OBSERVATIONS ALLOWED,
    C         MAXM IS THE MAXIMUM NUMBER OF COLUMNS IN THE
    C              INDEPENDENT VARIABLE ALLOWED, AND
    C        MAXNP IS THE MAXIMUM NUMBER OF FUNCTION PARAMETERS
    C              ALLOWED.
    C
          INTEGER MAXN,MAXM,MAXNP
          PARAMETER (MAXN=15,MAXM=5,MAXNP=5)
    C
    C  SET PARAMETERS TO SPECIFY DIMENSIONS OF ARRAYS USED BY DODR
    C
          INTEGER LDX,LDSCLD,LDIFX,LDRHO,LWORK,LIWORK
          PARAMETER
         +   (LDX=MAXN,
         +   LDSCLD = 1,
         +   LDIFX = 1,
         +   LDRHO=1,
         +   LWORK = 9 + 8*MAXN + 10*MAXN*MAXM + 2*MAXN*MAXNP + 8*MAXNP,
         +   LIWORK = 2*MAXNP + MAXM + 10)
    C
    C  DECLARE USER-SUPPLIED SUBROUTINES AND
    C  ALL OTHER NECESSARY VARIABLES AND ARRAYS
    C
          EXTERNAL
         +   FUN,JAC
          INTEGER
         +   N,M,NP,
         +   IFIXX(LDIFX,MAXM),
         +   IFIXB(MAXNP),
         +   JOB,
         +   MAXIT,
         +   IPRINT,LUNRPT,LUNERR,
         +   IWORK(LIWORK),
         +   INFO
          DOUBLE PRECISION
         +   X(LDX,MAXM),SCLD(LDSCLD,MAXM),
         +   Y(LDX),
         +   BETA(MAXNP),SCLB(MAXNP),
         +   RHO(LDRHO,MAXM),W(LDX),
         +   SSTOL,PARTOL,
         +   TAUFAC,
         +   WORK(LWORK)
    C
          OPEN(UNIT=5,FILE='DATA1')
          OPEN(UNIT=6,FILE='REPORT')
    C
    C  READ NUMBER OF OBSERVATIONS
    C       NUMBER OF COLUMNS OF DATA IN THE INDEPENDENT VARIABLE
    C       NUMBER OF PARAMETERS
    C       OBSERVED VALUES OF INDEPENDENT AND DEPENDENT VARIABLES
    C       STARTING VALUES OF FUNCTION PARAMETERS
    C
          READ (5,*) N,M,NP
          READ (5,*) ((X(I,J),I=1,N),J=1,M)
          READ (5,*) (Y(I),I=1,N)
          READ (5,*) (BETA(I),I=1,NP)
    C
    C  FIX SECOND COLUMN OF INDEPENDENT VARIABLE AT OBSERVED VALUES
    C
          IFIXX(1,1) = 1
          IFIXX(1,2) = 0
    C
    C  SPECIFY USE OF DEFAULT SCALING
    C
          SCLD(1,1) = -1.0D0
          SCLB(1) = -1.0D0
    C
    C  INDICATE ALL BETA'S ARE TO BE ESTIMATED
    C
          IFIXB(1) = -1
    C
    C  SPECIFY WEIGHTS
    C
          RHO(1,1) = 3.0D0
          RHO(1,2) = 5.0D0
          W(1) = -1.0D0
    C
    C  SET CONTROL VALUES AND STOPPING CRITERIA
    C
          JOB = 10
          TAUFAC = -1.0D0
          SSTOL = -1.0D0
          PARTOL = -1.0D0
          MAXIT = -1
          IPRINT = 1111
          LUNERR = -1
          LUNRPT = -1
    C
    C  COMPUTE ODR SOLUTION USING USER-SUPPLIED ANALYTIC DERIVATIVES
    C
          CALL DODRC
         +   (FUN,JAC,
         +   N,M,NP,
         +   X,LDX,IFIXX,LDIFX,SCLD,LDSCLD,
         +   Y,
         +   BETA,IFIXB,SCLB,
         +   RHO,LDRHO,W,
         +   JOB,TAUFAC,
         +   SSTOL,PARTOL,MAXIT,
         +   IPRINT,LUNERR,LUNRPT,
         +   WORK,LWORK,IWORK,LIWORK,
         +   INFO)
    C
          END
          SUBROUTINE FUN(N,NP,M,BETA,XPLUSD,LDXPD,F,IFLAG)
    C
    C  INPUT ARGUMENTS
    C  (WHICH MUST NOT BE CHANGED BY THIS ROUTINE)
    C
          INTEGER N,NP,M,LDXPD
          DOUBLE PRECISION BETA(NP),XPLUSD(LDXPD,M)
    C
    C  OUTPUT ARGUMENTS
    C
          DOUBLE PRECISION F(N)
          INTEGER IFLAG
    C
          DO 10 I = 1, N
             F(I) = EXP(-BETA(1)*XPLUSD(I,1)*
         +          EXP(-BETA(2)*
         +               (1.0D0/XPLUSD(I,2) - 1.0D0/620.0D0)))
       10 CONTINUE
    C
          RETURN
          END
          SUBROUTINE JAC(N,NP,M,BETA,XPLUSD,LDXPD,
         +               FJACB,LDFJB,ISODR,FJACX,LDFJX,IFLAG)
    C
    C  INPUT ARGUMENTS
    C  (WHICH MUST NOT BE CHANGED BY THIS ROUTINE)
    C
          INTEGER N,NP,M,LDXPD
          DOUBLE PRECISION BETA(NP),XPLUSD(LDXPD,M)
          LOGICAL ISODR
    C
    C  OUTPUT ARGUMENTS
    C
          DOUBLE PRECISION FJACB(LDFJB,NP),FJACX(LDFJX,M)
          INTEGER IFLAG
    C
    C  INTERNAL VARIABLES
    C
          DOUBLE PRECISION FAC1,FAC2,FAC3,FAC4
    C
          DO 10 I=1,N
             FAC1 = 1.0D0/XPLUSD(I,2) - 1.0D0/620.0D0
             FAC2 = EXP(-BETA(2)*FAC1)
             FAC3 = BETA(1)*XPLUSD(I,1)
             FAC4 = EXP(-FAC3*FAC2)
    C
             FJACB(I,1) = -FAC4*XPLUSD(I,1)*FAC2
             FJACB(I,2) = FAC4*FAC3*FAC2*FAC1
    C
             IF (ISODR) THEN
                FJACX(I,1) = -FAC4*BETA(1)*FAC2
                FJACX(I,2) = -FAC4*FAC3*FAC2*BETA(2)/XPLUSD(I,2)**2
             END IF
       10 CONTINUE
    C
          RETURN
          END
 
 
1   User-supplied data (file DATA1):
 
        8    2    2
      109.0   65.0 1180.0   66.0
     1270.0   69.0 1230.0   68.0
      600.0  640.0  600.0  640.0
      600.0  640.0  600.0  640.0
      0.912  0.382  0.397  0.376
      0.342  0.358  0.348  0.376
    0.01155 5000.0
 
 
1   Report generated by DODRC example program, run using a Cyber 180/855
    under NOS 2.4.2:
 
 
 
 
    ******************************************************
    * ODRPACK VERSION 1.3 OF 02-04-87 (DOUBLE PRECISION) *
    ******************************************************
 
 
 
 
 
    INITIAL SUMMARY FOR FIT BY METHOD OF ODR
    ========================================
 
 
 
    PROBLEM SIZE:
    -------------
 
    NUMBER OF OBSERVATIONS                                8
    NUMBER OF COLUMNS OF DATA IN INDEPENDENT VARIABLE     2
    NUMBER OF FUNCTION PARAMETERS                         2
    NUMBER OF UNFIXED FUNCTION PARAMETERS                 2
 
 
 
    CONTROL VALUES AND STOPPING CRITERIA:
    --------------------------------------
 
          *
       JOB    TAUFAC     SSTOL    PARTOL  MAXIT
      0010   .10D+01   .50D-14   .86D-19     50
 
    *
     A.  FIT IS NOT A RESTART.
     B.  DELTAS ARE INITIALIZED TO ZERO.
     C.  DERIVATIVES ARE SUPPLIED BY USER.
         USER-SUPPLIED DERIVATIVES WERE CHECKED.
         THE DERIVATIVES APPEAR TO BE CORRECT.
     D.  FIT IS BY METHOD OF ORTHOGONAL DISTANCE REGRESSION.
 
 
 
    INITIAL SUMS OF SQUARES:
    ------------------------
 
    SUM OF SQUARED WEIGHTED OBSERVATIONAL ERRORS     .67662011D+00
    SUM OF SQUARED WEIGHTED DELTAS                   .00000000D+00
    SUM OF SQUARED WEIGHTED EPSILONS                 .67662011D+00
 
 
 
 
 
    ITERATION REPORTS FOR FIT BY METHOD OF ODR
    ==========================================
 
 
            CUM.                 ACT. REL.   PRED. REL.
     IT.  NO. FN     WEIGHTED   SUM-OF-SQS   SUM-OF-SQS              G-N
    NUM.   EVALS   SUM-OF-SQS    REDUCTION    REDUCTION  TAU/PNORM  STEP
    ----  ------  -----------  -----------  -----------  ---------  ----
 
       1       2   .19694D+00    .7089D+00    .4162D+00   .151D+01   YES
       2       3   .18660D-02    .9905D+00    .9957D+00   .671D+00   YES
       3       4   .75385D-03    .5960D+00    .5961D+00   .463D-01   YES
       4       5   .75385D-03    .3659D-06    .3659D-06   .224D-04   YES
       5       6   .75385D-03    .3889D-13    .3892D-13   .482D-08   YES
       6       7   .75385D-03    .1675D-19    .1676D-19   .392D-11   YES
 
 
 
 
 
 
 
    FINAL SUMMARY FOR FIT BY METHOD OF ODR
    ======================================
 
 
 
    STOPPING CONDITION (INFO =      1):
    -----------------------------------
 
    THE RELATIVE CHANGE IN THE SUM OF THE SQUARED
    WEIGHTED OBSERVATIONAL ERRORS IS LESS THAN SSTOL
 
                                CONDITION
          NUMBER OF  NUMBER OF     NUMBER        RANK
         ITERATIONS   FN EVALS  (INVERSE)  DEFICIENCY
                  6          8  .1888D-06           0
 
 
 
    FINAL SUMS OF SQUARES:
    ----------------------
 
    SUM OF SQUARED WEIGHTED OBSERVATIONAL ERRORS     .75384644D-03
    SUM OF SQUARED WEIGHTED DELTAS                   .33248273D-09
    SUM OF SQUARED WEIGHTED EPSILONS                 .75384611D-03
 
 
 
    ESTIMATED BETA(I), I = 1, ..., NP:
    ----------------------------------
 
            INDEX            VALUE -------------->
        1 TO    2    .36579727D-02   .27627326D+05
 
 
 
    ESTIMATED EPSILON(I) AND DELTA(I,*), I = 1, ..., N:
    ---------------------------------------------------
 
        I      EPSILON(I)      DELTA(I,1)      DELTA(I,2)
        1   .16752465D-02   .14086189D-06   .00000000D+00
        2   .20435276D-02   .12838572D-05   .00000000D+00
        3  -.20690747D-01  -.71654588D-06   .00000000D+00
        4   .24306485D-02   .15047496D-05   .00000000D+00
        5   .72779764D-02   .23394016D-06   .00000000D+00
        6   .40794324D-02   .24163474D-05   .00000000D+00
        7   .13043483D-01   .43338715D-06   .00000000D+00
        8  -.85501699D-02  -.51395912D-05   .00000000D+00
 
 
 
1   VIII.   SCALING ALGORITHMS
            ------------------
 
    Poorly scaled problems, i.e., problems in which the unknowns BETA
    and DELTA vary over several orders of magnitude, can cause least
    squares procedures difficulty.  ODRPACK's scaling algorithms
    (discussed below) attempt to overcome these difficulties
    automatically , although it is preferable for the user to choose the
    units of the variable space so that the estimated parameters will
    have roughly the same magnitude [Dennis and Schnabel (1983)].  When
    the variables have roughly the same magnitude, the ODRPACK scaling
    algorithm will select scale values which are roughly equal, and the
    resulting computations will be the same (except for the effect of
    finite precision arithmetic) as an unscaled analysis, i.e., an
    analysis in which all of the scale values are set to one.  If the
    user does not do this, the ODRPACK scaling algorithm will select
    varying scale values.  This will not change the optimal solution,
    but it may affect the number of iterations required, or, in some
    cases, whether the algorithm is or is not successful.
 
    Users may substitute their own scaling values using subroutine
    arguments SCLD and SCLB (see section VI).
 
 
 
1   VIII.A  BETA Scaling
 
    ODRPACK chooses the scale values for the estimated BETA's as
    follows.
 
    If some of the starting values of BETA are nonzero then
       let BETA_max = the largest absolute value
                      of the nonzero starting values of BETA,
           BETA_min = the smallest absolute value
                      of the nonzero starting values of BETA, and
             epsmac = the value of machine precision, where machine
                      precision is defined as the smallest value e such
                      that 1 + e > 1 on the computer being used.
 
       For K = 1 to NP do
          if BETA(K) = zero then
             scale_BETA(K) = one/(SQRT(epsmac)*BETA_min)
          else
             if LOG10(BETA_max)-LOG10(BETA_min) > two then
                scale_BETA(K) = one/ABS(BETA(K))
             else
                scale_BETA(K) = one/BETA_max.
 
    If all of the starting values of BETA are zero then
       for K = 1 to NP do
           scale_BETA(K) = one.
 
    Users may substitute their own BETA scaling values via subroutine
    argument SCLB.
 
 
 
1   VIII.B  DELTA Scaling
 
    ODRPACK chooses scale values for the estimated errors in the
    independent variables, i.e., for the DELTA's, as follows.
 
    For J = 1 to M do
       let  X_max = the largest nonzero absolute value in the J-th
                    column of array X, and
            X_min = the smallest nonzero absolute value in the J-th
                    column of array X.
 
       For I = 1 to N do
          if X(I,J) = zero then
             scale_X(I,J) = one
          else
             if LOG10(X_max)-LOG10(X_min) < two then
                scale_X(I,J) = one/X_max
             else
                scale_X(I,J) = one/ABS(X(I,J)).
 
    Users may substitute their own DELTA scaling values via subroutine
    argument SCLD.
 
 
 
1   IX.   EXTRACTING INFORMATION FROM THE WORK VECTORS
          --------------------------------------------
 
    IX.A  Extracting Information from Vector WORK
 
    Upon return from a call to ODRPACK, array WORK contains various
    values, some of which may be of interest to the user.
 
    To extract information from WORK, the following declaration
    statement must be added to the user's program:
 
       INTEGER
      +   N,M,NP,
      +   DELTAI,EPSI,RNORMI,PARTLI,SSTOLI,TAUFCI,EPSMAI,OLMAVI,
      +   FJACBI,FJACXI,XPLUSI,BETACI,BETASI,BETANI,
      +   DELTSI,DELTNI,DDELTI,FSI,FNI,SI,SSSI,SSI,SSFI,TI,TTI,
      +   TAUI,ALPHAI,TFJACI,OMEGAI,YTI,UI,QRAUXI,WRK1I,WRK2I,STPI,
      +   RCONDI
 
    where the variables DELTAI through RCONDI indicate the starting
    locations within WORK of the stored values.  The appropriate values
    of DELTAI through RCONDI are obtained by invoking subroutine SWINF
    when using either of the single precision ODRPACK subroutines, SODR
    or SODRC, and by invoking DWINF when using either of the double
    precision subroutines, DODR or DODRC.  The call statements for SWINF
    and DWINF are as follows.
 
    When using single precision subroutines SODR and SODRC:
 
       CALL SWINF
      +   (N,M,NP,
      +   DELTAI,EPSI,RNORMI,PARTLI,SSTOLI,TAUFCI,EPSMAI,OLMAVI,
      +   FJACBI,FJACXI,XPLUSI,BETACI,BETASI,BETANI,
      +   DELTSI,DELTNI,DDELTI,FSI,FNI,SI,SSSI,SSI,SSFI,TI,TTI,
      +   TAUI,ALPHAI,TFJACI,OMEGAI,YTI,UI,QRAUXI,WRK1I,WRK2I,STPI,
      +   RCONDI)
 
    When using double precision subroutines DODR and DODRC:
 
       CALL DWINF
      +   (N,M,NP,
      +   DELTAI,EPSI,RNORMI,PARTLI,SSTOLI,TAUFCI,EPSMAI,OLMAVI,
      +   FJACBI,FJACXI,XPLUSI,BETACI,BETASI,BETANI,
      +   DELTSI,DELTNI,DDELTI,FSI,FNI,SI,SSSI,SSI,SSFI,TI,TTI,
      +   TAUI,ALPHAI,TFJACI,OMEGAI,YTI,UI,QRAUXI,WRK1I,WRK2I,STPI,
      +   RCONDI)
 
    Users should extract these declaration and call statements from
    online ODRPACK documentation, if possible, to avoid typographical
    errors.
 
    In the following descriptions of the information returned in WORK,
    (*) indicates values which are likely to be of greatest interest.
 
    (*) WORK(DELTAI)  is the first element of the N by M matrix, DELTA,
                      containing the estimated errors in the independent
                      variables,
                         DELTA(I,J) = WORK(DELTAI-1+I+(J-1)*N)
                      for I=1,...,N & J=1,...,M.
 
    (*) WORK(EPSI)    is the first element of the N vector, EPSILON,
                      containing the estimated errors in the dependent
                      variable,
                         EPSILON(I) = WORK(EPSI-1+I)
                      for I=1,...,N.
 
    (*) WORK(RNORMI)  is the square root of the sum of the squared
                      weighted observational errors (eq.3) at the
                      time the computations stopped.
 
        WORK(PARTLI)  is the value of the stopping tolerance used to
                      detect parameter convergence.
 
        WORK(SSTOLI)  is the value of the stopping tolerance used to
                      detect sum-of-squares convergence.
 
        WORK(TAUFCI)  is the value of the factor used to compute the
                      initial trust region radius.
 
        WORK(EPSMAI)  is the value of machine precision, which is the
                      smallest value e such that 1+e>1.
 
        WORK(OLMAVI)  is the average number of Levenberg-Marquardt steps
                      per iteration.
 
        WORK(FJACBI)  is the first element of the N by NP matrix, FJACB,
                      containing the weighted Jacobian with respect to
                      BETA,
                         FJACB(I,J) = WORK(FJACBI-1+I+(J-1)*N)
                      for I=1,...,N & J=1,...,NP.
 
        WORK(FJACXI)  is the first element of the N by M matrix, FJACX,
                      containing the weighted Jacobian with respect to
                      X,
                         FJACX(I,J) = WORK(FJACXI-1+I+(J-1)*N)
                      for I=1,...,N & J=1,...,M.
 
    (*) WORK(XPLUSI)  is the first element of the N by M matrix, XPLUSD,
                      containing the final values of X + DELTA,
                         XPLUSD(I,J) = WORK(XPLUSI-1+I+(J-1)*N)
                      for I=1,...,N & J=1,...,M.
 
        WORK(BETACI)  is the first element of the NP vector, BETAC,
                      containing the current estimates of the unknown
                      function parameters,
                         BETAC(I) = WORK(BETACI-1+I)
                      for I=1,...,NP.
 
        WORK(BETASI)  is the first element of the NP vector, BETAS,
                      containing the saved estimates of the unknown
                      function parameters,
                         BETAS(I) = WORK(BETASI-1+I)
                      for I=1,...,NP.
 
        WORK(BETANI)  is the first element of the NP vector, BETAN,
                      containing the new estimates of the unknown
                      function parameters,
                         BETAN(I) = WORK(BETANI-1+I)
                      for I=1,...,NP.
 
        WORK(DELTSI)  is the first element of the N by M matrix, DELTAS,
                      containing the saved estimated errors in the
                      independent variables,
                         DELTAS(I,J) = WORK(DELTASI-1+I+(J-1)*N)
                      for I=1,...,N & J=1,...,M.
 
        WORK(DELTNI)  is the first element of the N by M matrix, DELTAN,
                      containing the new estimated errors in the
                      independent variables,
                         DELTAN(I,J) = WORK(DELTANI-1+I+(J-1)*N)
                      for I=1,...,N & J=1,...,M.
 
        WORK(DDELTI)  is the first element of the N by M matrix, DDELTA,
                      containing the weighted estimated errors in the
                      independent variables (W*D)**2 * DELTA,
                         DDELTA(I,J) = WORK(DDELTI-1+I+(J-1)*N)
                      for I=1,...,N & J=1,...,M.
 
        WORK(FSI)     is the first element of the N vector, FS,
                      containing the saved weighted estimated errors in
                      the dependent variable,
                         FS(I) = WORK(FSI-1+I)
                      for I=1,...,N.
 
        WORK(FNI)     is the first element of the N vector, FN,
                      containing the new weighted estimated errors in
                      the dependent variable,
                         FN(I) = WORK(FNI-1+I)
                      for I=1,...,N.
 
        WORK(SI)      is the first element of the NP vector, S,
                      containing the step in the estimated function
                      parameters,
                         S(I) = WORK(SI-1+I)
                      for I=1,...,NP.
 
        WORK(SSSI)    is the first element of the N + N*M vector, SSS,
                      used for computing sums of squares of various
                      quantities,
                         SSS(I) = WORK(SSSI-1+I)
                      for I=1,...,N + N*M.
 
        WORK(SSI)     is the first element of the NP vector, SS,
                      containing the scale of the estimated function
                      parameters,
                         SS(I) = WORK(SSI-1+I)
                      for I=1,...,NP.
 
        WORK(SSFI)    is the first element of the NP vector, SSF,
                      containing the scale of each of the function
                      parameters,
                         SSF(I) = WORK(SSFI-1+I)
                      for I=1,...,NP.
 
        WORK(TI)      is the first element of the N by M array, T,
                      containing the step in the estimated errors in
                      the independent variable,
                         T(I,J) = WORK(TI-1+I+(J-1)*N)
                      for I=1,...,N & J=1,...,M.
 
        WORK(TTI)     is the first element of the N by M array, TT,
                      containing the scale of each the estimated errors
                      in the independent variable,
                         TT(I,J) = WORK(TTI-1+I+(J-1)*N)
                      for I=1,...,N & J=1,...,M.
 
        WORK(TAUI)    is the trust region radius at the time the
                      computations stopped.
 
        WORK(ALPHAI)  is the Levenburg-Marquardt parameter at the time
                      the computations stopped.
 
        WORK(TFJACI)  is the first element of the N by M array, TFJACB,
                      containing scaled Jacobian matrix
                      inv(diag(sqrt(OMEGA(I),I=1,...,N)) * FJACB,
                         TFJACB(I,J) = WORK(TFJACB-1+I+(J-1)*N)
                      for I=1,...,N & J=1,...,M.
 
        WORK(OMEGAI)  is the first element of the N vector, OMEGA,
                      containing the diagonal matrix
                      (I-JFACX*inv(P)*trans(FJACX)**(-1/2) where
                      P=trans(FJACX)*FJACX + D**2 + ALPHA*TT**2,
                         OMEGA(I) = WORK(OMEGAI-1+I)
                      for I=1,...,N.
 
        WORK(YTI)     is the first element of the N vector, YTI,
                      containing the array
                      -diag(sqrt(OMEGA(I),I=1,...,N)*(G1-V*inv(E)*D*G2)
                         YT(I) = WORK(YTI-1+I)
                      for I=1,...,N.
 
        WORK(UI)      is the first element of the N vector, UI,
                      containing the approximate null vector for TFJACB,
                         U(I) = WORK(UI-1+I)
                      for I=1,...,N.
 
        WORK(QRAUXI)  is the first element of the NP vector, QRAUX,
                      required to recover the Q-R decomposition of
                      TFJACB,
                         QRAUX(I) = WORK(QRAUXI-1+I)
                      for I=1,...,NP.
 
        WORK(WRK1I)   is the first element of the N by M matrix, WRK1,
                      required for work space,
                         WRK1(I,J) = WORK(WRK1I-1+I+(J-1)*N)
                      for I=1,...,N & J=1,...,M.
 
        WORK(WRK2I)   is the first element of the NP vector, WRK2,
                      required for work space,
                         WRK2(I) = WORK(WRK2I-1+I)
                      for I=1,...,NP.
 
        WORK(STPI)    is the first element of the N vector, STP,
                      containing the step used for computing finite
                      difference derivatives,
                         STP(I) = WORK(STPI-1+I)
                      for I=1,...,N.
 
    (*) WORK(RCONDI)  is the reciprocal condition number at the time
                      the computations stopped.
 
 
 
1   IX.B  Extracting Information from Vector IWORK
 
    Upon return from a call to ODRPACK, array IWORK contains various
    values, some of which may be of interest to the user.
 
    To extract information from IWORK, the following declaration
    statement must be added to the user's program
 
       INTEGER
      +   M,NP,
      +   MSGB,MSGX,JPVTI,MAXITI,IPRINI,NFEVI,NJEVI,INT2I,IRANKI,
      +   LUNERI,LUNRPI
 
    where the variables MSGB through LUNRPI indicate the starting
    locations within IWORK of the stored values.  The appropriate values
    of MSGB through LUNRPI are obtained by invoking subroutine SIWINF
    when using either of the single precision ODRPACK subroutines, SODR
    or SODRC, and by invoking DIWINF when using either of the double
    precision subroutines, DODR or DODRC.  The call statements for
    SIWINF and DIWINF are as follows.
 
    When using single precision subroutines SODR and SODRC:
 
       CALL SIWINF
      +   (M,NP,
      +   MSGB,MSGX,JPVTI,MAXITI,IPRINI,NFEVI,NJEVI,INT2I,IRANKI,
      +   LUNERI,LUNRPI)
 
    When using double precision subroutines DODR and DODRC:
 
       CALL DIWINF
      +   (M,NP,
      +   MSGB,MSGX,JPVTI,MAXITI,IPRINI,NFEVI,NJEVI,INT2I,IRANKI,
      +   LUNERI,LUNRPI)
 
    Users should extract these declaration and call statements from
    online ODRPACK documentation, if possible, to avoid typographical
    errors.
 
    In the following descriptions of the information returned in IWORK,
    (*) indicates values which are likely to be of greatest interest.
 
    (*) IWORK(MSGB)   is the first element of the NP+1 vector, MSGB,
                      used to indicate the results of checking the
                      partial derivatives with respect to BETA.
 
                      The value of IWORK(MSGB) summarizes the results
                      over all of the BETA's.
 
                      If IWORK(MSGB) < 0, the partial derivatives with
                                          respect to each of the BETA's
                                          were not checked.
 
                      If IWORK(MSGB) = 0, the partial derivatives with
                                          respect to each of the BETA's
                                          appear to be correct.
 
                      If IWORK(MSGB) = 1, the partial derivative with
                                          respect to at least one of the
                                          BETA's appears to be
                                          incorrect.
 
                      If IWORK(MSGB) = 2, the partial derivative with
                                          respect to at least one of the
                                          BETA's is questionable.
 
                      The value of IWORK(MSGB+K), K=1,...,NP, indicates
                      the individual results for the partial derivative
                      with respect to BETA(K), K=1,...,NP.
 
                      If IWORK(MSGB+K) = 0, the partial derivative with
                                            respect to BETA(K) appears
                                            to be correct.
 
                      If IWORK(MSGB+K) = 1, the partial derivative with
                                            respect to BETA(K) appears
                                            to be incorrect, i.e., the
                                            user-supplied derivative and
                                            the finite difference value
                                            it is checked against do not
                                            agree to within the required
                                            tolerance and there is no
                                            reason to question the
                                            results.
 
                      If IWORK(MSGB+K) = 2, the partial derivative with
                                            respect to BETA(K) appears
                                            to be questionable because
                                            the user-supplied derivative
                                            and the finite difference
                                            value it is checked against
                                            are both zero.
 
                      If IWORK(MSGB+K) = 3, the partial derivative with
                                            respect to BETA(K) appears
                                            to be questionable because
                                            the user-supplied derivative
                                            is exactly zero and the
                                            finite difference value it
                                            is checked against is only
                                            approximately zero.
 
                      If IWORK(MSGB+K) = 4, the partial derivative with
                                            respect to BETA(K) appears
                                            to be questionable because
                                            the user-supplied derivative
                                            is exactly zero and the
                                            finite difference value it
                                            is checked against is not
                                            even approximately zero.
 
                      If IWORK(MSGB+K) = 5, the partial derivative with
                                            respect to BETA(K) appears
                                            to be questionable because
                                            the finite difference value
                                            it is being checked against
                                            is questionable due to a
                                            high ratio of relative
                                            curvature to relative slope
                                            or to an incorrect scale
                                            value.
 
                      If IWORK(MSGB+K) = 6, the partial derivative with
                                            respect to BETA(K) appears
                                            to be questionable because
                                            the finite difference value
                                            it is being checked against
                                            is questionable due to a
                                            high ratio of relative
                                            curvature to relative slope.
 
    (*) IWORK(MSGX)   is the first element of the M+1 vector, MSGX,
                      used to indicate the results of checking the
                      partial derivatives with respect to X.
 
                      The value of IWORK(MSGX) summarizes the results
                      over all of the X's.
 
                      If IWORK(MSGX) < 0, the partial derivatives with
                                          respect to each of the X's
                                          were not checked.
 
                      If IWORK(MSGX) = 0, the partial derivatives with
                                          respect to each of the X's
                                          appear to be correct.
 
                      If IWORK(MSGX) = 1, the partial derivative with
                                          respect to at least one of the
                                          X's appears to be incorrect.
 
                      If IWORK(MSGX) = 2, the partial derivative with
                                          respect to at least one of the
                                          X's is questionable.
 
                      The value of IWORK(MSGX+J), J=1,...,M, indicates
                      the individual results for the partial derivative
                      with respect to the Jth column of X, J=1,...,M.
 
                      If IWORK(MSGX+J) = 0, the partial derivative with
                                            respect to the J-th column
                                            of X appears to be correct.
 
                      If IWORK(MSGX+J) = 1, the partial derivative with
                                            respect to the J-th column
                                            of X to be incorrect, i.e.,
                                            the user-supplied derivative
                                            and the finite difference
                                            value it is checked against
                                            do not agree to within the
                                            required tolerance and there
                                            is no reason to question the
                                            results.
 
                      If IWORK(MSGX+J) = 2, the partial derivative with
                                            respect to the J-th column
                                            of X appears to be
                                            questionable because the
                                            user-supplied derivative and
                                            the finite difference value
                                            it is checked against are
                                            both zero.
 
                      If IWORK(MSGX+J) = 3, the partial derivative with
                                            respect to the J-th column
                                            of X appears to be
                                            questionable because the
                                            user-supplied derivative is
                                            exactly zero and the finite
                                            difference value it is
                                            checked against is only
                                            approximately zero.
 
                      If IWORK(MSGX+J) = 4, the partial derivative with
                                            respect to the J-th column
                                            of X appears to be
                                            questionable because the
                                            user-supplied derivative is
                                            exactly zero and the finite
                                            difference value it is
                                            checked against is not even
                                            approximately zero.
 
                      If IWORK(MSGX+J) = 5, the partial derivative with
                                            respect to the J-th column
                                            of X appears to be
                                            questionable because the
                                            finite difference value it
                                            is being checked against is
                                            questionable due to a high
                                            ratio of relative curvature
                                            to relative slope or to an
                                            incorrect scale value.
 
                      If IWORK(MSGX+J) = 6, the partial derivative with
                                            respect to the J-th column
                                            of X appears to be
                                            questionable because the
                                            finite difference value it
                                            is being checked against is
                                            questionable due to a high
                                            ratio of relative curvature
                                            to relative slope.
 
        IWORK(JPVTI)  is the first element of the NP vector, JPVT,
                      containing the pivot vector,
                         JPVT(I) = WORK(JPVTI-1+I)
                      for I=1,...,NP.
 
        IWORK(MAXITI) is the maximum number of iterations allowed.
 
        IWORK(IPRINI) is the print control value used.
 
    (*) IWORK(NFEVI)  is the number of function evaluations at the time
                      the computations stopped.
 
    (*) IWORK(NJEVI)  is the number of Jacobian matrix evaluations, and,
                      equivalently, the number of iterations, at the
                      time the computations stopped.
 
        IWORK(INT2I)  is the number of internal doubling steps taken at
                      the time the computations stopped.
 
    (*) IWORK(IRANKI) is the rank deficiency at the solution.
 
        IWORK(LUNERI) is the logical unit number used for error reports.
 
        IWORK(LUNRPI) is the logical unit number used for computation
                      reports.
 
 
 
1   X.  ACKNOWLEDGMENTS
        ---------------
 
    Code for ODRPACK has been developed at the National Bureau of
    Standards, Gaithersburg, MD.  The subroutine that supplies the value
    of machine precision has been modeled after subroutines R1MACH and
    D1MACH from the Bell Laboratories "Framework for a Portable Library"
    [Fox et al.  1978].  We have also used subroutines from LINPACK
    [Dongarra et al.  1979] and from the "Basic Linear Algebra
    Subprograms for Fortran Usage (BLAS)" [Lawson et al.  1979].  The
    code used to check user-supplied derivatives was adapted from
    STARPAC [Donaldson and Tryon (1986)] using algorithms developed by
    Schnabel (1982).
 
 
 
1   XI.  REFERENCES
         ----------
 
    Boggs, P. T., R. H. Byrd, and R. B. Schnabel (1985),
      "A stable and efficient algorithm for nonlinear orthogonal
      distance regression," University of Colorado Department of
      Computer Science Technical Report Number CU-CS-317-85.
 
    Dennis, J. E., and R. B. Schnabel (1983),
      Numerical Methods for Unconstrained Optimization and Nonlinear
      Equations, Prentice-Hall, Englewood Cliffs, NJ.
 
    Donaldson, J. R., and P. V. Tryon (1986),
      "STARPAC - the standards time series and regression package,"
      National Bureau of Standards (U.S.) Interim Report 86-3448.
 
    Dongarra, J. J., C. B. Moler, J. R. Bunch, and G. W. Stewart (1979),
      Linpack Users' Guide, Siam, Philadelphia, PA.
 
    Draper, N. R., and H. Smith (1981),
      Applied Regression Analysis, Second Edition, John Wiley and Sons,
      New York, NY.
 
    Fox, P. A., A. D. Hall and N. L. Schryer (1978),
      "Algorithm 528:  framework for a portable library [z]," ACM Trans.
      Math.  Software, 4(2):177-188.
 
    Himmelblau, D. M. (1970),
      Process Analysis by Statistical Methods, John Wiley and Sons, New
      York, NY.
 
    Lawson, C., R. Hanson, D. Kincaid, and F. Krogh (1979),
      "Basic linear algebra subprograms for fortran usage", ACM Trans.
      Math.  Software, 5(3):308-371.
 
    Schnabel, R. B. (1982),
      "Finite difference derivatives - theory and practice",
      (unpublished, available from author).
 
 
