*DECK DFCMN SUBROUTINE DFCMN (NDATA, XDATA, YDATA, SDDATA, NORD, NBKPT, + BKPTIN, NCONST, XCONST, YCONST, NDERIV, MODE, COEFF, BF, XTEMP, + PTEMP, BKPT, G, MDG, W, MDW, WORK, IWORK) C***BEGIN PROLOGUE DFCMN C***SUBSIDIARY C***PURPOSE Subsidiary to FC C***LIBRARY SLATEC C***TYPE DOUBLE PRECISION (FCMN-S, DFCMN-D) C***AUTHOR (UNKNOWN) C***DESCRIPTION C C This is a companion subprogram to DFC( ). C The documentation for DFC( ) has complete usage instructions. C C***SEE ALSO DFC C***ROUTINES CALLED DAXPY, DBNDAC, DBNDSL, DCOPY, DFSPVD, DFSPVN, C DLSEI, DSCAL, DSORT, XERMSG C***REVISION HISTORY (YYMMDD) C 780801 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890618 Completely restructured and extensively revised (WRB & RWC) C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900328 Added TYPE section. (WRB) C 900510 Convert XERRWV calls to XERMSG calls. (RWC) C 900604 DP version created from SP version. (RWC) C***END PROLOGUE DFCMN INTEGER IWORK(*), MDG, MDW, MODE, NBKPT, NCONST, NDATA, NDERIV(*), * NORD DOUBLE PRECISION BF(NORD,*), BKPT(*), BKPTIN(*), COEFF(*), * G(MDG,*), PTEMP(*), SDDATA(*), W(MDW,*), WORK(*), * XCONST(*), XDATA(*), XTEMP(*), YCONST(*), YDATA(*) C EXTERNAL DAXPY, DBNDAC, DBNDSL, DCOPY, DFSPVD, DFSPVN, DLSEI, * DSCAL, DSORT, XERMSG C DOUBLE PRECISION DUMMY, PRGOPT(10), RNORM, RNORME, RNORML, XMAX, * XMIN, XVAL, YVAL INTEGER I, IDATA, IDERIV, ILEFT, INTRVL, INTW1, IP, IR, IROW, * ITYPE, IW1, IW2, L, LW, MT, N, NB, NEQCON, NINCON, NORDM1, * NORDP1, NP1 LOGICAL BAND, NEW, VAR CHARACTER*8 XERN1 C C***FIRST EXECUTABLE STATEMENT DFCMN C C Analyze input. C IF (NORD.LT.1 .OR. NORD.GT.20) THEN CALL XERMSG ('SLATEC', 'DFCMN', + 'IN DFC, THE ORDER OF THE B-SPLINE MUST BE 1 THRU 20.', + 2, 1) MODE = -1 RETURN C ELSEIF (NBKPT.LT.2*NORD) THEN CALL XERMSG ('SLATEC', 'DFCMN', + 'IN DFC, THE NUMBER OF KNOTS MUST BE AT LEAST TWICE ' // + 'THE B-SPLINE ORDER.', 2, 1) MODE = -1 RETURN ENDIF C IF (NDATA.LT.0) THEN CALL XERMSG ('SLATEC', 'DFCMN', + 'IN DFC, THE NUMBER OF DATA POINTS MUST BE NONNEGATIVE.', + 2, 1) MODE = -1 RETURN ENDIF C C Amount of storage allocated for W(*), IW(*). C IW1 = IWORK(1) IW2 = IWORK(2) NB = (NBKPT-NORD+3)*(NORD+1) + 2*MAX(NDATA,NBKPT) + NBKPT + + NORD**2 C C See if sufficient storage has been allocated. C IF (IW1.LT.NB) THEN WRITE (XERN1, '(I8)') NB CALL XERMSG ('SLATEC', 'DFCMN', * 'IN DFC, INSUFFICIENT STORAGE FOR W(*). CHECK NB = ' // * XERN1, 2, 1) MODE = -1 RETURN ENDIF C IF (MODE.EQ.1) THEN BAND = .TRUE. VAR = .FALSE. NEW = .TRUE. ELSEIF (MODE.EQ.2) THEN BAND = .FALSE. VAR = .TRUE. NEW = .TRUE. ELSEIF (MODE.EQ.3) THEN BAND = .TRUE. VAR = .FALSE. NEW = .FALSE. ELSEIF (MODE.EQ.4) THEN BAND = .FALSE. VAR = .TRUE. NEW = .FALSE. ELSE CALL XERMSG ('SLATEC', 'DFCMN', + 'IN DFC, INPUT VALUE OF MODE MUST BE 1-4.', 2, 1) MODE = -1 RETURN ENDIF MODE = 0 C C Sort the breakpoints. C CALL DCOPY (NBKPT, BKPTIN, 1, BKPT, 1) CALL DSORT (BKPT, DUMMY, NBKPT, 1) C C Initialize variables. C NEQCON = 0 NINCON = 0 DO 100 I = 1,NCONST L = NDERIV(I) ITYPE = MOD(L,4) IF (ITYPE.LT.2) THEN NINCON = NINCON + 1 ELSE NEQCON = NEQCON + 1 ENDIF 100 CONTINUE C C Compute the number of variables. C N = NBKPT - NORD NP1 = N + 1 LW = NB + (NP1+NCONST)*NP1 + 2*(NEQCON+NP1) + (NINCON+NP1) + + (NINCON+2)*(NP1+6) INTW1 = NINCON + 2*NP1 C C Save interval containing knots. C XMIN = BKPT(NORD) XMAX = BKPT(NP1) C C Find the smallest referenced independent variable value in any C constraint. C DO 110 I = 1,NCONST XMIN = MIN(XMIN,XCONST(I)) XMAX = MAX(XMAX,XCONST(I)) 110 CONTINUE NORDM1 = NORD - 1 NORDP1 = NORD + 1 C C Define the option vector PRGOPT(1-10) for use in DLSEI( ). C PRGOPT(1) = 4 C C Set the covariance matrix computation flag. C PRGOPT(2) = 1 IF (VAR) THEN PRGOPT(3) = 1 ELSE PRGOPT(3) = 0 ENDIF C C Increase the rank determination tolerances for both equality C constraint equations and least squares equations. C PRGOPT(4) = 7 PRGOPT(5) = 4 PRGOPT(6) = 1.D-4 C PRGOPT(7) = 10 PRGOPT(8) = 5 PRGOPT(9) = 1.D-4 C PRGOPT(10) = 1 C C Turn off work array length checking in DLSEI( ). C IWORK(1) = 0 IWORK(2) = 0 C C Initialize variables and analyze input. C IF (NEW) THEN C C To process least squares equations sort data and an array of C pointers. C CALL DCOPY (NDATA, XDATA, 1, XTEMP, 1) DO 120 I = 1,NDATA PTEMP(I) = I 120 CONTINUE C IF (NDATA.GT.0) THEN CALL DSORT (XTEMP, PTEMP, NDATA, 2) XMIN = MIN(XMIN,XTEMP(1)) XMAX = MAX(XMAX,XTEMP(NDATA)) ENDIF C C Fix breakpoint array if needed. C DO 130 I = 1,NORD BKPT(I) = MIN(BKPT(I),XMIN) 130 CONTINUE C DO 140 I = NP1,NBKPT BKPT(I) = MAX(BKPT(I),XMAX) 140 CONTINUE C C Initialize parameters of banded matrix processor, DBNDAC( ). C MT = 0 IP = 1 IR = 1 ILEFT = NORD DO 160 IDATA = 1,NDATA C C Sorted indices are in PTEMP(*). C L = PTEMP(IDATA) XVAL = XDATA(L) C C When interval changes, process equations in the last block. C IF (XVAL.GE.BKPT(ILEFT+1)) THEN CALL DBNDAC (G, MDG, NORD, IP, IR, MT, ILEFT-NORDM1) MT = 0 C C Move pointer up to have BKPT(ILEFT).LE.XVAL, C ILEFT.LT.NP1. C 150 IF (XVAL.GE.BKPT(ILEFT+1) .AND. ILEFT.LT.N) THEN ILEFT = ILEFT + 1 GO TO 150 ENDIF ENDIF C C Obtain B-spline function value. C CALL DFSPVN (BKPT, NORD, 1, XVAL, ILEFT, BF) C C Move row into place. C IROW = IR + MT MT = MT + 1 CALL DCOPY (NORD, BF, 1, G(IROW,1), MDG) G(IROW,NORDP1) = YDATA(L) C C Scale data if uncertainty is nonzero. C IF (SDDATA(L).NE.0.D0) CALL DSCAL (NORDP1, 1.D0/SDDATA(L), + G(IROW,1), MDG) C C When staging work area is exhausted, process rows. C IF (IROW.EQ.MDG-1) THEN CALL DBNDAC (G, MDG, NORD, IP, IR, MT, ILEFT-NORDM1) MT = 0 ENDIF 160 CONTINUE C C Process last block of equations. C CALL DBNDAC (G, MDG, NORD, IP, IR, MT, ILEFT-NORDM1) C C Last call to adjust block positioning. C CALL DCOPY (NORDP1, 0.D0, 0, G(IR,1), MDG) CALL DBNDAC (G, MDG, NORD, IP, IR, 1, NP1) ENDIF C BAND = BAND .AND. NCONST.EQ.0 DO 170 I = 1,N BAND = BAND .AND. G(I,1).NE.0.D0 170 CONTINUE C C Process banded least squares equations. C IF (BAND) THEN CALL DBNDSL (1, G, MDG, NORD, IP, IR, COEFF, N, RNORM) RETURN ENDIF C C Check further for sufficient storage in working arrays. C IF (IW1.LT.LW) THEN WRITE (XERN1, '(I8)') LW CALL XERMSG ('SLATEC', 'DFCMN', * 'IN DFC, INSUFFICIENT STORAGE FOR W(*). CHECK LW = ' // * XERN1, 2, 1) MODE = -1 RETURN ENDIF C IF (IW2.LT.INTW1) THEN WRITE (XERN1, '(I8)') INTW1 CALL XERMSG ('SLATEC', 'DFCMN', * 'IN DFC, INSUFFICIENT STORAGE FOR IW(*). CHECK IW1 = ' // * XERN1, 2, 1) MODE = -1 RETURN ENDIF C C Write equality constraints. C Analyze constraint indicators for an equality constraint. C NEQCON = 0 DO 220 IDATA = 1,NCONST L = NDERIV(IDATA) ITYPE = MOD(L,4) IF (ITYPE.GT.1) THEN IDERIV = L/4 NEQCON = NEQCON + 1 ILEFT = NORD XVAL = XCONST(IDATA) C 180 IF (XVAL.LT.BKPT(ILEFT+1) .OR. ILEFT.GE.N) GO TO 190 ILEFT = ILEFT + 1 GO TO 180 C 190 CALL DFSPVD (BKPT, NORD, XVAL, ILEFT, BF, IDERIV+1) CALL DCOPY (NP1, 0.D0, 0, W(NEQCON,1), MDW) CALL DCOPY (NORD, BF(1,IDERIV+1), 1, W(NEQCON,ILEFT-NORDM1), + MDW) C IF (ITYPE.EQ.2) THEN W(NEQCON,NP1) = YCONST(IDATA) ELSE ILEFT = NORD YVAL = YCONST(IDATA) C 200 IF (YVAL.LT.BKPT(ILEFT+1) .OR. ILEFT.GE.N) GO TO 210 ILEFT = ILEFT + 1 GO TO 200 C 210 CALL DFSPVD (BKPT, NORD, YVAL, ILEFT, BF, IDERIV+1) CALL DAXPY (NORD, -1.D0, BF(1, IDERIV+1), 1, + W(NEQCON, ILEFT-NORDM1), MDW) ENDIF ENDIF 220 CONTINUE C C Transfer least squares data. C DO 230 I = 1,NP1 IROW = I + NEQCON CALL DCOPY (N, 0.D0, 0, W(IROW,1), MDW) CALL DCOPY (MIN(NP1-I, NORD), G(I,1), MDG, W(IROW,I), MDW) W(IROW,NP1) = G(I,NORDP1) 230 CONTINUE C C Write inequality constraints. C Analyze constraint indicators for inequality constraints. C NINCON = 0 DO 260 IDATA = 1,NCONST L = NDERIV(IDATA) ITYPE = MOD(L,4) IF (ITYPE.LT.2) THEN IDERIV = L/4 NINCON = NINCON + 1 ILEFT = NORD XVAL = XCONST(IDATA) C 240 IF (XVAL.LT.BKPT(ILEFT+1) .OR. ILEFT.GE.N) GO TO 250 ILEFT = ILEFT + 1 GO TO 240 C 250 CALL DFSPVD (BKPT, NORD, XVAL, ILEFT, BF, IDERIV+1) IROW = NEQCON + NP1 + NINCON CALL DCOPY (N, 0.D0, 0, W(IROW,1), MDW) INTRVL = ILEFT - NORDM1 CALL DCOPY (NORD, BF(1, IDERIV+1), 1, W(IROW, INTRVL), MDW) C IF (ITYPE.EQ.1) THEN W(IROW,NP1) = YCONST(IDATA) ELSE W(IROW,NP1) = -YCONST(IDATA) CALL DSCAL (NORD, -1.D0, W(IROW, INTRVL), MDW) ENDIF ENDIF 260 CONTINUE C C Solve constrained least squares equations. C CALL DLSEI(W, MDW, NEQCON, NP1, NINCON, N, PRGOPT, COEFF, RNORME, + RNORML, MODE, WORK, IWORK) RETURN END