* * $Id: dsumsq.F,v 1.1.1.1 1996/04/01 15:02:20 mclareni Exp $ * * $Log: dsumsq.F,v $ * Revision 1.1.1.1 1996/04/01 15:02:20 mclareni * Mathlib gen * * #include "gen/pilot.h" SUBROUTINE DSUMSQ(SUB,M,N,NC,A,AL,AU,MODE,EPS,MAXIT,IPRT, + MFR,IAFR,PHI,DPHI,COV,STD,W,NERROR) #include "gen/imp64.inc" DIMENSION A(*),AL(*),AU(*),DPHI(*),IAFR(*) DIMENSION COV(NC,*),STD(*) DIMENSION W(*) EXTERNAL SUB *********************************************************************** * LEAMAX, VERSION: 15.03.1993 ************************************************************************ * * DSUMSQ IS THE STEERING ROUTINE FOR MINIMIZING A SUM OF SQUARES * * SUBROUTINE CALLED: D501L1 * * * THE CONSTANTS, VARIABLES AND ARRAYS HAVE THE FOLLOWING MEANING. * * SUB NAME OF USER-SUPPLIED SUBROUTINE SUBPROGRAM, DECLARED * EXTERNAL IN THE CALLING PROGRAM. THIS SUBPROGRAM MUST PROVIDE * THE VALUES OF THE FUNCTION AND, IF MODE=1, THE VALUES OF THE * DERIVATIVES (SEE EXAMPLE) . * M (INTEGER) NUMBER OF NONLINEAR FUNCTIONS. * N (INTEGER) NUMBER OF UNKNOWN PARAMETERS A. * NC (INTEGER) DECLARED FIRST DIMENSION OF ARRAY COV IN THE * CALLING PROGRAM, WITH NC .GE. N . * A (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH N . * ON ENTRY, A MUST CONTAIN THE STARTING VALUES OF THE UNKNOWN * PARAMETERS FOR THE LEVENBERG-MARQUARDT ALGORITHM. * ON EXIT, A CONTAINS AN APPROXIMATION OF THE MINIMUM POINT. * AL (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH N . * ON ENTRY, AL MUST CONTAIN THE LOWER BOUNDS OF A . * AU (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH N . * ON ENTRY, AU MUST CONTAIN THE UPPER BOUNDS OF A . * MODE (INTEGER) * = 0: THE JACOBIAN IS COMPUTED NUMERICALLY * = 1: THE JACOBIAN HAS TO BE EVALUATED IN SUBPROGRAM SUB . * EPS (DOUBLE PRECISION) USER-SUPPLIED TOLERANCE USED TO CONTROL * THE TERMINATION CRITERION. EPS SHOULD BE CHOSEN ACCORDING * TO THE ACCURACY REQUIRED BY THE UNDERLYING PROBLEM AND TO * THE MACHINE ACCURACY ALSO (RECOMMENDED VALUE ON ENTRY: * 1D-6 ... 1D-12 ). * MAXIT (INTEGER) MAXIMUM PERMITTED NUMBER OF ITERATIONS. * IPRT (INTEGER) PRINTING CONTROL. * = 0 : NO PRINTING OF INTERMEDIATE RESULTS * = +/- L : PRINTING OF INTERMEDIATE RESULTS AT EVERY ABS(L)-TH * ITERATION; IF IPRT < 0, PRINTING OF ALL INPUT * PARAMETERS OF DSUMSQ IN ADDITION. * MFR (INTEGER) ON EXIT, MFR CONTAINS THE NUMBER OF FREE VARIABLES * AT THE SOLUTION POINT. * IAFR (INTEGER) ONE-DIMENSIONAL ARRAY OF LENGTH 2 * N , USED AS * WORKING SPACE. ON EXIT, THE FIRST MFR ELEMENTS OF IAFR * CONTAIN THE INDICES OF THE FREE VARIABLES AT THE SOLUTION * POINT. * PHI (DOUBLE PRECISION) ON EXIT, PHI CONTAINS THE VALUE OF THE * OBJECTIVE FUNCTION AT THE MINIMUM POINT. * DPHI (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH N . * ON EXIT, DPHI CONTAINS THE DERIVATIVES OF THE OBJECTIVE * FUNCTION WITH RESPECT TO A (THE GRADIENT) AT THE LAST * ITERATION POINT. * COV (DOUBLE PRECISION) TWO-DIMENSIONAL ARRAY OF DIMENSION (NC,N). * ON EXIT, COV CONTAINS AN APPROXIMATION TO THE COVARIANCE * MATRIX. * STD (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH N . * ON EXIT, STD CONTAINS APPROXIMATIONS TO THE STANDARD * DEVIATIONS OF THE MODEL PARAMETER ESTIMATORS. * W (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY OF LENGTH * 9*N+4*M+2*M*N+3*N*N , USED AS WORKING SPACE. * NERROR (INTEGER) ERROR INDICATOR. ON EXIT: * = 0: NO ERROR OR WARNING DETECTED. * = 1: AT LEAST ONE OF THE CONSTANTS M, N, NC, MAXIT IS * ILLEGAL OR AT LEAST FOR ONE J THE RELATION * AL(J) .LE. AU(J) IS NOT TRUE. * = 2: THE MAXIMUM NUMBER MAXIT OF ITERATIONS HAS BEEN * REACHED. * = 3: THE OBJECTIVE FUNCTION PHI OR ITS DERIVATIVE IS NOT * DEFINED FOR THE CURRENT VALUES OF THE UNKNOWN * PARAMETER VECTOR A. * = 4: THE ROUTINES DGEQPF , DORMQR , DTRTRS OF THE LINEAR * ALGEBRA PACKAGE LAPACK (F001) WERE UNABLE TO SOLVE * THE LINEAR LEAST SQUARES PROBLEMS * OR THE ROUTINE DSINV (F012) WAS UNABLE TO COMPUTE THE * COVARIANCE MATRIX . * ************************************************************************* * * THE FOLLOWING SUBROUTINE IS A SIMPLE EXAMPLE FOR SUB. * * SUBROUTINE SUB(N,A,M,F,DF,MODE,NERROR) * IMPLICIT DOUBLE PRECISION (A-H,O-Z) * PARAMETER (Z0 = 0) * DIMENSION A(*),F(*),DF(M,*) * NERROR=0 * F(1)=A(1)-1D6 * F(2)=A(2)-2D-6 * F(3)=A(1)*A(2)-2 * IF(MODE .EQ. 0) RETURN * CALL DMSET(M,N,Z0,DF(1,1),DF(1,2),DF(2,1)) * DF(1,1)=1 * DF(2,2)=1 * DF(3,1)=A(2) * DF(3,2)=A(1) * RETURN * END * ************************************************************************* M1=1 M2=M1+N M3=M2+N M4=M3+N M5=M4+N M6=M5+2*M M7=M6+3*N M8=M7+N M9=M8+M+N MA=M9+(N+M)*N MB=MA+N*N MC=MB+N*N MD=MC+M CALL D501L1('DSUMSQ',SUB,1,M,X,1,Y,SY,MODE,EPS,MAXIT, 1 IPRT,N,A,AL,AU,PHI,DPHI,IAFR,MFR,COV,NC,STD, 2 W(M1),W(M2),W(M3),W(M4),W(M5),W(M6),W(M7),W(M8), 3 W(M9),W(MA),W(MB),W(MC),W(MD),IAFR(N+1),NERROR) RETURN END