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