1c     program DRDNLSFB
2c>> 2001-05-24 DRDNLSFB Krogh Minor change for making .f90 version.
3c>> 1994-11-02 DRDNLSFB Krogh  Changes to use M77CON
4c>> 1994-09-14 DRDNLSFB CLL Set IV(OUTLEV) = 0 for comparing output.
5c>> 1992-02-03 CLL @ JPL
6c>> 1990-07-02 CLL @ JPL
7c>> 1990-06-27 CLL @ JPL
8c>> 1990-06-14 CLL @ JPL
9c>> 1990-04-05 CLL @ JPL
10c>> 1990-03-29 CLL @ JPL
11c     Demo driver for DNLSFB. A variant of the nonlinear LS code NL2SOL.
12c     DNLSFB solves the "separable" problem.
13c     DNLSFB handles bounds on the nonlinear variables.
14c     DNLSFB requires function values only.
15c     Note:  The expressions below set LIV and LV larger than
16c     necessary because the precise formulas cannot be written in a
17c     Fortran parameter statement.
18c     ------------------------------------------------------------------
19c--D replaces "?": DR?NLSFB, ?NLSFB, ?CALCA, ?IVSET
20c     ------------------------------------------------------------------
21      external DCALCA
22      integer  J,ITERM,IVAR,LIV,LV,MA,MB,MDATA,MDIR,NA,NB,NDATA
23      parameter(MDATA = 30, MA = 2, MB = 5, MDIR = 4)
24      parameter(LIV = 122 + 2*MDIR + 7*MA + 2*MB + MB+1 + 6*MA)
25      parameter(LV = 105 + MDATA*(2*MB + 6 + MA) +
26     *               (MB*(MB+3))/2 + MA*(2*MA + 22))
27      integer IND(MB+1,MA), IV(LIV)
28      integer F, OUTLEV, SOLPRT, STATPR, X0PRT
29      parameter(F=10, OUTLEV=19, SOLPRT=22, STATPR=23, X0PRT=24)
30      double precision ALF(MA), BET(MB), BND(2,MA)
31      double precision DOF, V(LV), YDATA(MDATA)
32      data (BND(1,J),J=1,2) /  2*5.0d0 /
33      data (BND(2,J),J=1,2) / 2*10.0d0 /
34      data ((IND(ITERM,IVAR),IVAR = 1,2),ITERM = 1,6)/
35     *     0, 0,
36     *     1, 0,
37     *     1, 0,
38     *     0, 1,
39     *     0, 1,
40     *     0, 0/
41      data YDATA /
42     *     1.700641d0, 1.793512d0, 1.838309d0, 1.838416d0, 1.792204d0,
43     *     1.700501d0, 1.579804d0, 1.426268d0, 1.260724d0, 1.084901d0,
44     *     0.917094d0, 0.761920d0, 0.627304d0, 0.522146d0, 0.446645d0,
45     *     0.404920d0, 0.392033d0, 0.409622d0, 0.453045d0, 0.510765d0,
46     *     0.584554d0, 0.663109d0, 0.747613d0, 0.829439d0, 0.908496d0,
47     *     0.983178d0, 1.051046d0, 1.114072d0, 1.171746d0, 1.227823d0/
48c     ------------------------------------------------------------------
49      NDATA = MDATA
50      NA = MA
51      NB = MB
52      ALF(1) =    5.0d0
53      ALF(2) =   10.0d0
54      IV(1) = 0
55
56      print'(1x,a)',
57     * 'Program DRDNLSFB.. Demo driver for DNLSFB.',
58     * '   A variant of NL2SOL.',
59     * '   DNLSFB handles the Separable problem.',
60     * '   DNLSFB requires function values but not the Jacobian.',
61     * '   DNLSFB handles bounds on the nonlinear variables.',
62     * ' ',
63     * 'Sample problem is a nonlinear curve fit to data.',
64     * 'Model function is B1 + B2 * cos(A1*t) + B3 * sin(A1*t) +',
65     * '                       B4 * cos(A2*t) + B5 * sin(A2*t) + Noise',
66     * 'Data generated using',
67     * '(A1, A2, B1, ..., B5) = (6, 9, 1, 0.5, 0.4, 0.2, 0.1)',
68     * 'and Gaussian noise with mean 0 and',
69     * 'sample standard deviation 0.001',
70     * ' '
71
72      call DIVSET(1, IV, LIV, LV, V)
73      IV( X0PRT) = 1
74      IV(OUTLEV) = 0
75      IV(STATPR) = 1
76      IV(SOLPRT) = 1
77
78      call DNLSFB(NDATA, NA, NB, ALF, BND, BET, YDATA, DCALCA,
79     *         IND, NB+1, IV, LIV, LV, V)
80
81      DOF = max(NDATA - NA - NB, 1)
82      print'(1x/1x,a,g12.4)',
83     *    'SIGFAC: sqrt((2 * V(F))/DOF) =',
84     *             sqrt(2.0d0 * V(F)/DOF)
85      stop
86      end
87c     ==================================================================
88      subroutine DCALCA(NDATA, NA, NB, ALF, NCOUNT, PHI)
89c     Test case for separable nonlinear least squares computation.
90c     Computes MDATA x NB matrix PHI as a function of the
91c     nonlinear parameters ALF().
92c     For J .le. NB the (I,J) term of PHI is the coefficient of the
93c     linear coefficient B(J) in row I of the model.
94c     In this example the model does not have a term that is not
95c     multiplied by a linear coefficient.  If such a term is present
96c     then PHI must have an (NB+1)st column to hold this term.
97c     ------------------------------------------------------------------
98      integer I, MDATA, NA, NB, NCOUNT, NDATA
99      parameter(MDATA = 30)
100      double precision ALF(NA)
101      double precision DEL, PHI(NDATA,NB), T
102c     ------------------------------------------------------------------
103      T = 0.0D0
104      DEL = 1.0D0 / 29.0D0
105      do 10 I = 1,NDATA
106         PHI(I,1) = 1.0D0
107         PHI(I,2) = cos(ALF(1)*T)
108         PHI(I,3) = sin(ALF(1)*T)
109         PHI(I,4) = cos(ALF(2)*T)
110         PHI(I,5) = sin(ALF(2)*T)
111         T = T + DEL
112   10 continue
113      return
114      end
115