1*DECK FCMN
2      SUBROUTINE FCMN (NDATA, XDATA, YDATA, SDDATA, NORD, NBKPT, BKPTIN,
3     +   NCONST, XCONST, YCONST, NDERIV, MODE, COEFF, BF, XTEMP, PTEMP,
4     +   BKPT, G, MDG, W, MDW, WORK, IWORK)
5C***BEGIN PROLOGUE  FCMN
6C***SUBSIDIARY
7C***PURPOSE  Subsidiary to FC
8C***LIBRARY   SLATEC
9C***TYPE      SINGLE PRECISION (FCMN-S, DFCMN-D)
10C***AUTHOR  (UNKNOWN)
11C***DESCRIPTION
12C
13C     This is a companion subprogram to FC( ).
14C     The documentation for FC( ) has complete usage instructions.
15C
16C***SEE ALSO  FC
17C***ROUTINES CALLED  BNDACC, BNDSOL, BSPLVD, BSPLVN, LSEI, SAXPY, SCOPY,
18C                    SSCAL, SSORT, XERMSG
19C***REVISION HISTORY  (YYMMDD)
20C   780801  DATE WRITTEN
21C   890531  Changed all specific intrinsics to generic.  (WRB)
22C   890618  Completely restructured and extensively revised (WRB & RWC)
23C   891214  Prologue converted to Version 4.0 format.  (BAB)
24C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
25C   900328  Added TYPE section.  (WRB)
26C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
27C***END PROLOGUE  FCMN
28      INTEGER IWORK(*), MDG, MDW, MODE, NBKPT, NCONST, NDATA, NDERIV(*),
29     *   NORD
30      REAL             BF(NORD,*), BKPT(*), BKPTIN(*), COEFF(*),
31     *   G(MDG,*), PTEMP(*), SDDATA(*), W(MDW,*), WORK(*),
32     *   XCONST(*), XDATA(*), XTEMP(*), YCONST(*), YDATA(*)
33C
34      EXTERNAL BNDACC, BNDSOL, BSPLVD, BSPLVN, LSEI, SAXPY, SCOPY,
35     *    SSCAL, SSORT, XERMSG
36C
37      REAL             DUMMY, PRGOPT(10), RNORM, RNORME, RNORML, XMAX,
38     *   XMIN, XVAL, YVAL
39      INTEGER I, IDATA, IDERIV, ILEFT, INTRVL, INTW1, IP, IR, IROW,
40     *   ITYPE, IW1, IW2, L, LW, MT, N, NB, NEQCON, NINCON, NORDM1,
41     *   NORDP1, NP1
42      LOGICAL BAND, NEW, VAR
43      CHARACTER*8 XERN1
44C
45C***FIRST EXECUTABLE STATEMENT  FCMN
46C
47C     Analyze input.
48C
49      IF (NORD.LT.1 .OR. NORD.GT.20) THEN
50         CALL XERMSG ('SLATEC', 'FCMN',
51     +      'IN FC, THE ORDER OF THE B-SPLINE MUST BE 1 THRU 20.',
52     +      2, 1)
53         MODE = -1
54         RETURN
55C
56      ELSEIF (NBKPT.LT.2*NORD) THEN
57         CALL XERMSG ('SLATEC', 'FCMN',
58     +      'IN FC, THE NUMBER OF KNOTS MUST BE AT LEAST TWICE ' //
59     +      'THE B-SPLINE ORDER.', 2, 1)
60         MODE = -1
61         RETURN
62      ENDIF
63C
64      IF (NDATA.LT.0) THEN
65         CALL XERMSG ('SLATEC', 'FCMN',
66     +      'IN FC, THE NUMBER OF DATA POINTS MUST BE NONNEGATIVE.',
67     +      2, 1)
68         MODE = -1
69         RETURN
70      ENDIF
71C
72C     Amount of storage allocated for W(*), IW(*).
73C
74      IW1 = IWORK(1)
75      IW2 = IWORK(2)
76      NB = (NBKPT-NORD+3)*(NORD+1) + 2*MAX(NDATA,NBKPT) + NBKPT +
77     +     NORD**2
78C
79C     See if sufficient storage has been allocated.
80C
81      IF (IW1.LT.NB) THEN
82         WRITE (XERN1, '(I8)') NB
83         CALL XERMSG ('SLATEC', 'FCMN',
84     *      'IN FC, INSUFFICIENT STORAGE FOR W(*).  CHECK NB = ' //
85     *      XERN1, 2, 1)
86         MODE = -1
87         RETURN
88      ENDIF
89C
90      IF (MODE.EQ.1) THEN
91         BAND = .TRUE.
92         VAR = .FALSE.
93         NEW = .TRUE.
94      ELSEIF (MODE.EQ.2) THEN
95         BAND = .FALSE.
96         VAR = .TRUE.
97         NEW = .TRUE.
98      ELSEIF (MODE.EQ.3) THEN
99         BAND = .TRUE.
100         VAR = .FALSE.
101         NEW = .FALSE.
102      ELSEIF (MODE.EQ.4) THEN
103         BAND = .FALSE.
104         VAR = .TRUE.
105         NEW = .FALSE.
106      ELSE
107         CALL XERMSG ('SLATEC', 'FCMN',
108     +      'IN FC, INPUT VALUE OF MODE MUST BE 1-4.', 2, 1)
109         MODE = -1
110         RETURN
111      ENDIF
112      MODE = 0
113C
114C     Sort the breakpoints.
115C
116      CALL SCOPY (NBKPT, BKPTIN, 1, BKPT, 1)
117      CALL SSORT (BKPT, DUMMY, NBKPT, 1)
118C
119C     Initialize variables.
120C
121      NEQCON = 0
122      NINCON = 0
123      DO 100 I = 1,NCONST
124         L = NDERIV(I)
125         ITYPE = MOD(L,4)
126         IF (ITYPE.LT.2) THEN
127            NINCON = NINCON + 1
128         ELSE
129            NEQCON = NEQCON + 1
130         ENDIF
131  100 CONTINUE
132C
133C     Compute the number of variables.
134C
135      N = NBKPT - NORD
136      NP1 = N + 1
137      LW = NB + (NP1+NCONST)*NP1 + 2*(NEQCON+NP1) + (NINCON+NP1) +
138     +     (NINCON+2)*(NP1+6)
139      INTW1 = NINCON + 2*NP1
140C
141C     Save interval containing knots.
142C
143      XMIN = BKPT(NORD)
144      XMAX = BKPT(NP1)
145C
146C     Find the smallest referenced independent variable value in any
147C     constraint.
148C
149      DO 110 I = 1,NCONST
150         XMIN = MIN(XMIN,XCONST(I))
151         XMAX = MAX(XMAX,XCONST(I))
152  110 CONTINUE
153      NORDM1 = NORD - 1
154      NORDP1 = NORD + 1
155C
156C     Define the option vector PRGOPT(1-10) for use in LSEI( ).
157C
158      PRGOPT(1) = 4
159C
160C     Set the covariance matrix computation flag.
161C
162      PRGOPT(2) = 1
163      IF (VAR) THEN
164         PRGOPT(3) = 1
165      ELSE
166         PRGOPT(3) = 0
167      ENDIF
168C
169C     Increase the rank determination tolerances for both equality
170C     constraint equations and least squares equations.
171C
172      PRGOPT(4) = 7
173      PRGOPT(5) = 4
174      PRGOPT(6) = 1.E-4
175C
176      PRGOPT(7) = 10
177      PRGOPT(8) = 5
178      PRGOPT(9) = 1.E-4
179C
180      PRGOPT(10) = 1
181C
182C     Turn off work array length checking in LSEI( ).
183C
184      IWORK(1) = 0
185      IWORK(2) = 0
186C
187C     Initialize variables and analyze input.
188C
189      IF (NEW) THEN
190C
191C        To process least squares equations sort data and an array of
192C        pointers.
193C
194         CALL SCOPY (NDATA, XDATA, 1, XTEMP, 1)
195         DO 120 I = 1,NDATA
196            PTEMP(I) = I
197  120    CONTINUE
198C
199         IF (NDATA.GT.0) THEN
200            CALL SSORT (XTEMP, PTEMP, NDATA, 2)
201            XMIN = MIN(XMIN,XTEMP(1))
202            XMAX = MAX(XMAX,XTEMP(NDATA))
203         ENDIF
204C
205C        Fix breakpoint array if needed.
206C
207         DO 130 I = 1,NORD
208            BKPT(I) = MIN(BKPT(I),XMIN)
209  130    CONTINUE
210C
211         DO 140 I = NP1,NBKPT
212            BKPT(I) = MAX(BKPT(I),XMAX)
213  140    CONTINUE
214C
215C        Initialize parameters of banded matrix processor, BNDACC( ).
216C
217         MT = 0
218         IP = 1
219         IR = 1
220         ILEFT = NORD
221         DO 160 IDATA = 1,NDATA
222C
223C           Sorted indices are in PTEMP(*).
224C
225            L = PTEMP(IDATA)
226            XVAL = XDATA(L)
227C
228C           When interval changes, process equations in the last block.
229C
230            IF (XVAL.GE.BKPT(ILEFT+1)) THEN
231               CALL BNDACC (G, MDG, NORD, IP, IR, MT, ILEFT-NORDM1)
232               MT = 0
233C
234C              Move pointer up to have BKPT(ILEFT).LE.XVAL,
235C                 ILEFT.LT.NP1.
236C
237  150          IF (XVAL.GE.BKPT(ILEFT+1) .AND. ILEFT.LT.N) THEN
238                  ILEFT = ILEFT + 1
239                  GO TO 150
240               ENDIF
241            ENDIF
242C
243C           Obtain B-spline function value.
244C
245            CALL BSPLVN (BKPT, NORD, 1, XVAL, ILEFT, BF)
246C
247C           Move row into place.
248C
249            IROW = IR + MT
250            MT = MT + 1
251            CALL SCOPY (NORD, BF, 1, G(IROW,1), MDG)
252            G(IROW,NORDP1) = YDATA(L)
253C
254C           Scale data if uncertainty is nonzero.
255C
256            IF (SDDATA(L).NE.0.E0) CALL SSCAL (NORDP1, 1.E0/SDDATA(L),
257     +                                  G(IROW,1), MDG)
258C
259C           When staging work area is exhausted, process rows.
260C
261            IF (IROW.EQ.MDG-1) THEN
262               CALL BNDACC (G, MDG, NORD, IP, IR, MT, ILEFT-NORDM1)
263               MT = 0
264            ENDIF
265  160    CONTINUE
266C
267C        Process last block of equations.
268C
269         CALL BNDACC (G, MDG, NORD, IP, IR, MT, ILEFT-NORDM1)
270C
271C        Last call to adjust block positioning.
272C
273         CALL SCOPY (NORDP1, 0.E0, 0, G(IR,1), MDG)
274         CALL BNDACC (G, MDG, NORD, IP, IR, 1, NP1)
275      ENDIF
276C
277      BAND = BAND .AND. NCONST.EQ.0
278      DO 170 I = 1,N
279         BAND = BAND .AND. G(I,1).NE.0.E0
280  170 CONTINUE
281C
282C     Process banded least squares equations.
283C
284      IF (BAND) THEN
285         CALL BNDSOL (1, G, MDG, NORD, IP, IR, COEFF, N, RNORM)
286         RETURN
287      ENDIF
288C
289C     Check further for sufficient storage in working arrays.
290C
291      IF (IW1.LT.LW) THEN
292         WRITE (XERN1, '(I8)') LW
293         CALL XERMSG ('SLATEC', 'FCMN',
294     *      'IN FC, INSUFFICIENT STORAGE FOR W(*).  CHECK LW = ' //
295     *      XERN1, 2, 1)
296         MODE = -1
297         RETURN
298      ENDIF
299C
300      IF (IW2.LT.INTW1) THEN
301         WRITE (XERN1, '(I8)') INTW1
302         CALL XERMSG ('SLATEC', 'FCMN',
303     *      'IN FC, INSUFFICIENT STORAGE FOR IW(*).  CHECK IW1 = ' //
304     *      XERN1, 2, 1)
305         MODE = -1
306         RETURN
307      ENDIF
308C
309C     Write equality constraints.
310C     Analyze constraint indicators for an equality constraint.
311C
312      NEQCON = 0
313      DO 220 IDATA = 1,NCONST
314         L = NDERIV(IDATA)
315         ITYPE = MOD(L,4)
316         IF (ITYPE.GT.1) THEN
317            IDERIV = L/4
318            NEQCON = NEQCON + 1
319            ILEFT = NORD
320            XVAL = XCONST(IDATA)
321C
322  180       IF (XVAL.LT.BKPT(ILEFT+1) .OR. ILEFT.GE.N) GO TO 190
323            ILEFT = ILEFT + 1
324            GO TO 180
325C
326  190       CALL BSPLVD (BKPT, NORD, XVAL, ILEFT, BF, IDERIV+1)
327            CALL SCOPY (NP1, 0.E0, 0, W(NEQCON,1), MDW)
328            CALL SCOPY (NORD, BF(1,IDERIV+1), 1, W(NEQCON,ILEFT-NORDM1),
329     +                  MDW)
330C
331            IF (ITYPE.EQ.2) THEN
332               W(NEQCON,NP1) = YCONST(IDATA)
333            ELSE
334               ILEFT = NORD
335               YVAL = YCONST(IDATA)
336C
337  200          IF (YVAL.LT.BKPT(ILEFT+1) .OR. ILEFT.GE.N) GO TO 210
338               ILEFT = ILEFT + 1
339               GO TO 200
340C
341  210          CALL BSPLVD (BKPT, NORD, YVAL, ILEFT, BF, IDERIV+1)
342               CALL SAXPY (NORD, -1.E0, BF(1, IDERIV+1), 1,
343     +                     W(NEQCON, ILEFT-NORDM1), MDW)
344            ENDIF
345         ENDIF
346  220 CONTINUE
347C
348C     Transfer least squares data.
349C
350      DO 230 I = 1,NP1
351         IROW = I + NEQCON
352         CALL SCOPY (N, 0.E0, 0, W(IROW,1), MDW)
353         CALL SCOPY (MIN(NP1-I, NORD), G(I,1), MDG, W(IROW,I), MDW)
354         W(IROW,NP1) = G(I,NORDP1)
355  230 CONTINUE
356C
357C     Write inequality constraints.
358C     Analyze constraint indicators for inequality constraints.
359C
360      NINCON = 0
361      DO 260 IDATA = 1,NCONST
362         L = NDERIV(IDATA)
363         ITYPE = MOD(L,4)
364         IF (ITYPE.LT.2) THEN
365            IDERIV = L/4
366            NINCON = NINCON + 1
367            ILEFT = NORD
368            XVAL = XCONST(IDATA)
369C
370  240       IF (XVAL.LT.BKPT(ILEFT+1) .OR. ILEFT.GE.N) GO TO 250
371            ILEFT = ILEFT + 1
372            GO TO 240
373C
374  250       CALL BSPLVD (BKPT, NORD, XVAL, ILEFT, BF, IDERIV+1)
375            IROW = NEQCON + NP1 + NINCON
376            CALL SCOPY (N, 0.E0, 0, W(IROW,1), MDW)
377            INTRVL = ILEFT - NORDM1
378            CALL SCOPY (NORD, BF(1, IDERIV+1), 1, W(IROW, INTRVL), MDW)
379C
380            IF (ITYPE.EQ.1) THEN
381               W(IROW,NP1) = YCONST(IDATA)
382            ELSE
383               W(IROW,NP1) = -YCONST(IDATA)
384               CALL SSCAL (NORD, -1.E0, W(IROW, INTRVL), MDW)
385            ENDIF
386         ENDIF
387  260 CONTINUE
388C
389C     Solve constrained least squares equations.
390C
391      CALL LSEI(W, MDW, NEQCON, NP1, NINCON, N, PRGOPT, COEFF, RNORME,
392     +          RNORML, MODE, WORK, IWORK)
393      RETURN
394      END
395