PROGRAM twoexp ! Code converted using TO_F90 by Alan Miller ! Date: 1999-03-15 Time: 11:15:27 ! SAMPLE DRIVER PROGRAM FOR NONLINEAR LEAST SQUARES ROUTINE ! SOLUTION TO SAMPLE PROBLEM (TO 3 SIGNIFICANT DIGITS): ! LINEAR PARAMETERS .375, 1.94, -1.46 ! NONLINEAR PARAMETERS .0129, .0221 USE separable_leastsq IMPLICIT NONE REAL (dp) :: y(50), t(50,1), alf(12), beta(8), a(50,13), w(50), eta INTEGER :: i, ierr, im1, iprint, iv, j, l, lpnl, n, nl, nmax, p INTEGER, PARAMETER :: INPUT = 5, output = 6 INTERFACE SUBROUTINE ada (lp1, nl, n, nmax, lpp2, iv, a, inc, t, alf, isel) USE separable_leastsq IMPLICIT NONE INTEGER, INTENT(IN) :: lp1 INTEGER, INTENT(IN) :: nl INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: nmax INTEGER, INTENT(IN) :: lpp2 INTEGER, INTENT(IN) :: iv REAL (dp), INTENT(OUT) :: a(:,:) INTEGER, INTENT(OUT) :: inc(:,:) REAL (dp), INTENT(IN) :: t(:,:) REAL (dp), INTENT(IN) :: alf(:) INTEGER, INTENT(IN) :: isel END SUBROUTINE ada END INTERFACE ! READ AND LIST THE DATA OPEN (INPUT, FILE='a.dat', STATUS='OLD') OPEN (output, FILE='b.out', STATUS='NEW') nmax = 50 iprint = 1 READ (INPUT,*) n, l, nl, p, iv WRITE (output,110) n, l, nl, p, iv DO i=1,n w(i) = 1.0_dp READ (INPUT,*) t(i,1:iv), y(i) END DO DO i=1,n WRITE (output,210) i, t(i,1:iv), y(i) END DO READ (INPUT,120) alf(1:nl) WRITE (output,130) alf(1:nl) WRITE (output,140) CALL varpro (l, nl, n, nmax, l+p+2, iv, t, y, w, ada, a, iprint, & alf, beta, ierr) WRITE (6,30) ierr 30 FORMAT (t16, 'IERR =', i7) ! PRINT LOWER TRIANGLE OF COVARIANCE MATRIX IF (ierr <= -4) GO TO 90 lpnl = l + nl WRITE (output,160) DO i=1,lpnl WRITE (output,170) a(i,1:i) a(i,i) = SQRT(a(i,i)) END DO ! PRINT STANDARD DEVIATIONS OF PARAMETER ESTIMATES WRITE (output,180) WRITE (output,170) (a(j,j),j=1,lpnl) ! CALCULATE AND PRINT LOWER TRIANG. OF CORRELATION MATRIX IF (lpnl == 1) GO TO 60 DO i=2,lpnl im1 = i-1 DO j=1,im1 a(i,j) = a(i,j) / (a(i,i)*a(j,j)) END DO END DO 60 WRITE (output,190) DO i=1,lpnl a(i,i) = 1.0_dp WRITE (output,170) a(i,1:i) END DO ! PRINT RESIDUALS WRITE (output,200) DO i=1,n eta = y(i) - a(i,lpnl+1) / SQRT(w(i)) WRITE (output,210) i, w(i), t(i,1:iv), y(i), eta, a(i,lpnl+1) END DO 90 STOP 110 FORMAT ('1 NON-LINEAR LEAST SQUARES PROBLEM'// & ' NUMBER OF OBSERVATIONS =',i5, ' NUMBER OF LINEAR PARAMETERS =',i4// & ' NUMBER OF NONLINEAR PARAMETERS =',i4, & ' NUMBER OF NONVANISHING PARTIAL DERIVATIVES =',i4// & ' NUMBER OF INDEPENDENT VARIABLES =',i4// ' I T(I) Y(I)'//) 120 FORMAT (4E20.7) 130 FORMAT ('0 INITIAL NONLINEAR PARAMETERS'//(4E20.7)) 140 FORMAT ('0',50('*')) 160 FORMAT (/' COVARIANCE MATRIX'/) 170 FORMAT (8E15.7) 180 FORMAT ('0 STANDARD DEVIATIONS OF PARAMETER ESTIMATES'/) 190 FORMAT (/' CORRELATION MATRIX'/) 200 FORMAT (/' I W(I) T(I) Y(I)', & ' PREDICTED Y WEIGHTED RESIDUAL'//) 210 FORMAT (i5,7E16.7) END PROGRAM twoexp SUBROUTINE ada (lp1, nl, n, nmax, lpp2, iv, a, inc, t, alf, isel) ! EVALUATE THE FUNCTIONS PHI AND THEIR PARTIAL DERIVATIVES ! D PHI(J)/D ALF(K), AT THE SAMPLE POINTS T(I). ! ISEL = 1 MEANS COMPUTE BOTH FUNCTIONS AND DERIVATIVES AND ! INITIALIZE INC AND CONSTANT PHI'S ! = 2 MEANS COMPUTE ONLY THE NONCONSTANT FUNCTIONS PHI ! = 3 MEANS COMPUTE ONLY THE DERIVATIVES. ! THIS PARTICULAR ROUTINE IS FOR FITTING TWO EXPONENTIALS AND ! ONE CONSTANT TERM: ! ETA(C, ALF; T) = C + C * EXP(-ALF *T) + C * EXP(-ALF * T) ! 1 2 1 3 2 ! WHERE C IS THE VECTOR OF LINEAR PARAMETERS TO BE DETERMINED, ! AND ALF IS THE VECTOR OF NONLINEAR PARAMETERS TO BE DETERMINED. ! HERE N = 33, L = 3, NL = P = 2, NCFUN = 1, PHI(1) = 1, PHI(2) = ! EXP(-ALF(1)*T), AND PHI(3) = EXP(-ALF(2)*T). THIS EXAMPLE AND ! TEST DATA ARE TAKEN FROM REFERENCE 3 (SEE VARPRO). ! .................................................................. USE separable_leastsq IMPLICIT NONE INTEGER, INTENT(IN) :: lp1 INTEGER, INTENT(IN) :: nl INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: nmax INTEGER, INTENT(IN) :: lpp2 INTEGER, INTENT(IN) :: iv REAL (dp), INTENT(OUT) :: a(:,:) ! a(nmax,lpp2) INTEGER, INTENT(OUT) :: inc(:,:) ! inc(12,lp1) REAL (dp), INTENT(IN) :: t(:,:) ! t(nmax) REAL (dp), INTENT(IN) :: alf(:) ! alf(nl) INTEGER, INTENT(IN) :: isel ! Local variable INTEGER :: i SELECT CASE ( isel ) CASE (1) GO TO 10 CASE (2) GO TO 30 CASE (3) GO TO 50 END SELECT ! ISEL = 1, IN THIS CASE THE INCIDENCE MATRIX INC IS: ! 0 1 0 0 <- ALF(1) ! 0 0 1 0 <- ALF(2) 10 inc = 0 inc(1,2) = 1 inc(2,3) = 1 ! ISEL = 1, SET CONSTANT FUNCTION PHI(1) (REMOVE IF NO ! CONSTANT FUNCTIONS) a(1:n,1) = 1.0_dp ! ISEL = 1 OR 2, COMPUTE NONCONSTANT FUNCTIONS PHI 30 DO i=1,n a(i,2) = EXP(-alf(1)*t(i,1)) a(i,3) = EXP(-alf(2)*t(i,1)) END DO ! COLUMN L+1 = 4 IS LEFT FOR WORKSPACE. IF (isel == 2) GO TO 70 ! ISEL = 1 OR 3, COMPUTE DERIVATIVES DPHI(I) / D ALF(J) 50 DO i=1,n a(i,5) = -t(i,1)*EXP(-alf(1)*t(i,1)) a(i,6) = -t(i,1)*EXP(-alf(2)*t(i,1)) END DO 70 RETURN END SUBROUTINE ada