MINICERN/mathlib/gen/d/dfunft.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:20  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE DFUNFT(SUB,K,M,N,NX,NC,X,Y,SY,A,AL,AU,MODE,EPS,MAXIT,
  11      +                  IPRT,MFR,IAFR,PHI,DPHI,COV,STD,W,NERROR)
  12
  13 #include "gen/imp64.inc"
  14
  15       DIMENSION X(*),Y(*),SY(*)
  16       DIMENSION A(*),AL(*),AU(*),DPHI(*),IAFR(*)
  17       DIMENSION COV(NC,*),STD(*)
  18       DIMENSION W(*)
  19
  20       EXTERNAL SUB
  21
  22 ***********************************************************************
  23 *   LEAMAX, VERSION: 15.03.1993
  24 ***********************************************************************
  25 *
  26 *   DFUNFT IS THE STEERING ROUTINE FOR NONLINEAR LEAST-SQUARES FITTING.
  27 *
  28 *   SUBROUTINE CALLED:     D501L1
  29 *
  30 *
  31 *   THE CONSTANTS, VARIABLES AND ARRAYS HAVE THE FOLLOWING MEANING.
  32 *
  33 *   SUB    NAME OF USER-SUPPLIED SUBROUTINE SUBPROGRAM, DECLARED
  34 *          EXTERNAL IN THE CALLING PROGRAM. THIS SUBPROGRAM MUST PROVIDE
  35 *          THE VALUES OF THE FUNCTION AND, IF MODE=1, THE VALUES OF THE
  36 *          DERIVATIVES (SEE EXAMPLE) .
  37 *   K      (INTEGER) DIMENSION OF A SINGLE DATA POINT (OBSERVATION) X .
  38 *   M      (INTEGER) NUMBER OF DATA POINTS (OBSERVATIONS) .
  39 *   N      (INTEGER) NUMBER OF UNKNOWN PARAMETERS A.
  40 *   NX     (INTEGER) DECLARED FIRST DIMENSION OF ARRAY  X  IN THE
  41 *          CALLING PROGRAM, WITH  NX .GE. K .
  42 *   NC     (INTEGER) DECLARED FIRST DIMENSION OF ARRAY  COV  IN THE
  43 *          CALLING PROGRAM, WITH  NC .GE. N .
  44 *   X      (DOUBLE PRECISION) TWO-DIMENSIONAL ARRAY OF DIMENSION (NX,M).
  45 *          ON ENTRY, X MUST CONTAIN THE DATA SET  (X(I)) (THE I-TH
  46 *          COLUMN OF X BELONGS TO THE DATA POINT X(I), I=1,...,M).
  47 *   Y      (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH  M ,
  48 *          CONTAINS THE DATA SET  (Y(I)) ON ENTRY.
  49 *   SY     (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH  M ,
  50 *          CONTAINS THE WEIGTHS (SIGMA(I))  OF THE DATA POINTS ON ENTRY.
  51 *   A      (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH  N .
  52 *          ON ENTRY, A  MUST CONTAIN THE STARTING VALUES OF THE UNKNOWN
  53 *          PARAMETERS FOR THE LEVENBERG-MARQUARDT ALGORITHM.
  54 *          ON EXIT, A  CONTAINS AN APPROXIMATION OF THE MINIMUM POINT.
  55 *   AL     (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH  N .
  56 *          ON ENTRY, AL  MUST CONTAIN THE LOWER BOUNDS OF  A .
  57 *   AU     (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH  N .
  58 *          ON ENTRY, AU  MUST CONTAIN THE UPPER BOUNDS OF  A .
  59 *   MODE   (INTEGER)
  60 *          = 0: THE DERIVATIVE IS COMPUTED NUMERICALLY.
  61 *          = 1: THE DERIVATIVE HAS TO BE EVALUATED IN SUBPROGRAM  SUB .
  62 *   EPS    (DOUBLE PRECISION) USER-SUPPLIED TOLERANCE USED TO CONTROL
  63 *          THE TERMINATION CRITERION. EPS SHOULD BE CHOSEN ACCORDING
  64 *          TO THE ACCURACY REQUIRED BY THE UNDERLYING PROBLEM AND TO
  65 *          THE MACHINE ACCURACY ALSO (RECOMMENDED VALUE ON ENTRY:
  66 *          1D-6 ... 1D-12 ).
  67 *   MAXIT  (INTEGER) MAXIMUM PERMITTED NUMBER OF ITERATIONS.
  68 *   IPRT   (INTEGER) PRINTING CONTROL.                                  G
  69 *          = 0     : NO PRINTING OF INTERMEDIATE RESULTS
  70 *          = +/- L : PRINTING OF INTERMEDIATE RESULTS AT EVERY ABS(L)-TH
  71 *                    ITERATION; IF  IPRT < 0, PRINTING OF ALL INPUT
  72 *                    PARAMETERS OF DSUMSQ IN ADDITION.
  73 *   MFR    (INTEGER) ON EXIT, MFR CONTAINS THE NUMBER OF FREE VARIABLES
  74 *          AT THE SOLUTION POINT.
  75 *   IAFR   (INTEGER) ONE-DIMENSIONAL ARRAY OF LENGTH  2 * N , USED AS
  76 *          WORKING SPACE. ON EXIT, THE FIRST  MFR  ELEMENTS OF  IAFR
  77 *          CONTAIN THE INDICES OF THE FREE VARIABLES AT THE SOLUTION
  78 *          POINT.
  79 *   PHI    (DOUBLE PRECISION) ON EXIT, PHI  CONTAINS THE VALUE OF THE
  80 *          OBJECTIVE FUNCTION AT THE MINIMUM POINT.
  81 *   DPHI   (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH  N .
  82 *          ON EXIT, DPHI  CONTAINS THE DERIVATIVES OF THE OBJECTIVE
  83 *          FUNCTION WITH RESPECT TO A (THE GRADIENT) AT THE LAST
  84 *          ITERATION POINT.
  85 *   COV    (DOUBLE PRECISION) TWO-DIMENSIONAL ARRAY OF DIMENSION (NC,N).
  86 *          ON EXIT, COV CONTAINS AN APPROXIMATION TO THE COVARIANCE
  87 *          MATRIX.
  88 *   STD    (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH  N .
  89 *          ON EXIT, STD  CONTAINS APPROXIMATIONS TO THE STANDARD
  90 *          DEVIATIONS OF THE MODEL PARAMETER ESTIMATORS.
  91 *   W      (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH
  92 *          9*N+4*M+2*M*N+3*N*N , USED AS WORKING SPACE.
  93 *   NERROR (INTEGER) ERROR INDICATOR. ON EXIT:
  94 *           = 0: NO ERROR OR WARNING DETECTED.
  95 *           = 1: AT LEAST ONE OF THE CONSTANTS K, M, N, NX, NC, MAXIT IS
  96 *                ILLEGAL OR AT LEAST FOR ONE J THE RELATION
  97 *                AL(J) .LE. AU(J)  IS NOT TRUE.
  98 *           = 2: THE MAXIMUM NUMBER  MAXIT  OF ITERATIONS HAS BEEN
  99 *                REACHED.
 100 *           = 3: THE OBJECTIVE FUNCTION  PHI  OR ITS DERIVATIVE IS NOT
 101 *                DEFINED FOR THE CURRENT VALUES OF THE UNKNOWN
 102 *                PARAMETER VECTOR  A.
 103 *           = 4: THE ROUTINES  DGEQPF , DORMQR , DTRTRS  OF THE LINEAR
 104 *                ALGEBRA PACKAGE  LAPACK (F001)  WERE UNABLE TO SOLVE
 105 *                THE LINEAR LEAST SQUARES PROBLEMS
 106 *                OR THE ROUTINE  DSINV (F012)  WAS UNABLE TO COMPUTE THE
 107 *                COVARIANCE MATRIX .
 108 *
 109 *************************************************************************
 110 *
 111 *   THE FOLLOWING SUBROUTINE IS A SIMPLE EXAMPLE FOR SUB.
 112 *
 113 *     SUBROUTINE SUB (K,X,N,A,F,DF,MODE,NERROR)
 114 *     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 115 *     DIMENSION A(*),X(*),DF(*)
 116 *
 117 *     NERROR=0
 118 *     T=X(2)*A(2)+X(3)*A(3)
 119 *     IF (T .EQ. 0) THEN
 120 *      NERROR=1
 121 *      RETURN
 122 *     ENDIF
 123 *
 124 *     F=A(1)+X(1)/T
 125 *
 126 *     IF(MODE .EQ. 0) RETURN
 127 *     DF(1)=1
 128 *     DF(2)=-X(1)*X(2)/T**2
 129 *     DF(3)=-X(1)*X(3)/T**2
 130 *     RETURN
 131 *     END
 132 *
 133 *************************************************************************
 134
 135       M1=1
 136       M2=M1+N
 137       M3=M2+N
 138       M4=M3+N
 139       M5=M4+N
 140       M6=M5+2*M
 141       M7=M6+3*N
 142       M8=M7+N
 143       M9=M8+M+N
 144       MA=M9+(N+M)*N
 145       MB=MA+N*N
 146       MC=MB+N*N
 147       MD=MC+M
 148
 149       CALL D501L1('DFUNFT',SUB,K,M,X,NX,Y,SY,MODE,EPS,MAXIT,
 150      1            IPRT,N,A,AL,AU,PHI,DPHI,IAFR,MFR,COV,NC,STD,
 151      2            W(M1),W(M2),W(M3),W(M4),W(M5),W(M6),W(M7),W(M8),
 152      3            W(M9),W(MA),W(MB),W(MC),W(MD),IAFR(N+1),NERROR)
 153
 154       RETURN
 155
 156       END
 157
 158
 159