SUBROUTINE CDSVFT ( * * inputs * : X, Y, SIG, X0, XS, IPOWER, WORK, NDATA, NX, NXP, : MA, MP, NP, B, * * outputs * : A, U, V, W, CHISQ, STATUS) * * Module number: * * Module name: * * Keyphrase: * ---------- * Linear least squares fitting of MA functions of NX variables using * singular value decomposition. * * Description: * ------------ * adapted from Numerical Recipes by William Press et. al., p. 518, * subroutine SVDFIT * Double precision version * * Given a set of NDATA points X(NX,I),Y(I) with individual standard * deviations SIG(I), use chi-square minimization to determine the MA * coefficients A of the fitting function: * * MA * y(x) =Sum A AFUNC (x). * i=1 i i * * Here we solve the fitting equations using singular value decomposition * of the NDATA by MA matrix, as in section 2.9 of Numerical Recipes. * Arrays U,V,W provide workspace on input, on output they define the singular * value decomposition, and can be used to obtain the covariance matrix. * MP,NP are the physical dimensions of the matrices * U,V,W, as indicated below. It is necessary that MP>=NDATA and * NP>=MA. The program returns values for the MA fit parameters A, and * chi-square, CHISQ. The user supplies a subroutine * that returns the MA basis functions evaluated at X in the array AFUNC. * * MAMAX is maximum expected number of fitting functions MA. * * FORTRAN name: CDSVFT.FOR * * Keywords of accessed files and tables: * -------------------------------------- * Name I/O Description / Comments * * Subroutines Called: * ------------------- * CDBS: * CDSVDC, CDSVBK * SDAS: * UMSPUT * Others: * None * * History: * -------- * Version Date Author Description * 1 12-11-88 J.-C. HSU coding *------------------------------------------------------------------------------- * *== input: * --number of input data points INTEGER NDATA, * --number of independent variables : NX, * --physical dimension of X array : NXP, * --number of fitting functions : MA, * --polynomial power in each dimension * --of X : IPOWER(1), * --physical dimensions of U, V, W : MP, NP * --input independent variables DOUBLE PRECISION X(NXP, 1), * --zero point offset of X : X0(1), * --scaling factor of X : XS(1), * --input dependent variable and its * --dtandard deviation : Y(1), SIG(1), * --work area (>= NDATA) : B(1), * --work area (>= MA) : WORK(1) * *== output: * --coefficients of fitting functions DOUBLE PRECISION A(1), * --SVD component matrices * --zero values of W indicate that the * --corresponding functions were excluded * --from the fit by the SVD. : U(MP, NP), V(NP, NP), W(NP), * --chi-squared of the fit : CHISQ * --error status INTEGER STATUS * *== local: * --error message CHARACTER*130 CONTXT, MESS CHARACTER*3 CHAR3 INTEGER I, J, MAMAX, NXMAX, STATOK PARAMETER (MAMAX = 50) PARAMETER (NXMAX = 7) DOUBLE PRECISION TMP, WMAX, SUM, THRESH, DUMMY(NXMAX), * --tolerance : TOL, * --basis functions evaluated at X : AFUNC(MAMAX) PARAMETER (TOL = 1.D-12) * *=========================begin hsp.inc========================================= * --status return code INTEGER OK, ERRNUM(20) INTEGER DEST, PRIO DATA OK /0/ DATA ERRNUM /701, 702, 703, 704, 705, 706, 707, 708, 709, 710, : 711, 712, 713, 714, 715, 716, 717, 718, 719, 720/ * --message destination and priority DATA DEST, PRIO /1, 0/ *=========================end hsp.inc=========================================== *------------------------------------------------------------------------------ * * check that the number of X's dimension * IF (NX .GT. NXMAX) THEN STATUS = ERRNUM(1) WRITE (CHAR3, '(I3)') NXMAX CONTXT = 'dimension of input data is greater than ' : // CHAR3 GO TO 999 END IF * * check that the number of basis functions is within range * IF (MA .GT. MAMAX) THEN STATUS = ERRNUM(1) WRITE (CHAR3, '(I3)') MAMAX CONTXT = 'number of basis functions is greater than ' : // CHAR3 GO TO 999 END IF * * Accumulate coefficients of the fitting matrix * DO 120 I = 1, NDATA DO 100 J = 1, NX DUMMY(J) = X(J, I) 100 CONTINUE * CALL VPNFN2 (DUMMY, NX, X0, XS, IPOWER, AFUNC) * TMP = 1.D0 / SIG(I) DO 110 J = 1, MA U(I,J) = AFUNC(J) * TMP 110 CONTINUE B(I) = Y(I) * TMP 120 CONTINUE * * singular value decomposition * CALL CDSVDC (NDATA, MA, MP, NP, WORK, U, W, V, STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'cannot perform singular value decomposition' GO TO 999 END IF * * Edit the singular values, given TOL from the parameter statement. * WMAX = 0.D0 DO 130 J = 1, MA IF (W(J) .GT. WMAX) WMAX = W(J) 130 CONTINUE THRESH = TOL * WMAX DO 140 J = 1, MA IF(W(J) .LT. THRESH) W(J) = 0.D0 140 CONTINUE * * Solve U * A = B, where U is the matrix before SVD. * CALL CDSVBK (U, W, V, NDATA, MA, MP, NP, B, WORK, A) * * Evaluate chi-square * CHISQ = 0.D0 DO 170 I = 1, NDATA DO 150 J = 1, NX DUMMY(J) = X(J, I) 150 CONTINUE * CALL VPNFN2 (DUMMY, NX, X0, XS, IPOWER, AFUNC) * SUM = 0.D0 DO 160 J = 1, MA SUM = SUM + A(J) * AFUNC(J) 160 CONTINUE CHISQ = CHISQ + ((Y(I)-SUM)/SIG(I))**2 170 CONTINUE * GO TO 1000 * 999 MESS = 'CDSVFT: ' // CONTXT CALL UMSPUT (MESS, DEST, PRIO, STATOK) * 1000 RETURN END