SUBROUTINE CHKDER(M,N,X,FVEC,FJAC,LDFJAC,XP,FVECP,MODE,ERR) INTEGER M,N,LDFJAC,MODE REAL X(N),FVEC(M),FJAC(LDFJAC,N),XP(N),FVECP(M),ERR(M) C ********** C C SUBROUTINE CHKDER C C THIS SUBROUTINE CHECKS THE GRADIENTS OF M NONLINEAR FUNCTIONS C IN N VARIABLES, EVALUATED AT A POINT X, FOR CONSISTENCY WITH C THE FUNCTIONS THEMSELVES. THE USER MUST CALL CHKDER TWICE, C FIRST WITH MODE = 1 AND THEN WITH MODE = 2. C C MODE = 1. ON INPUT, X MUST CONTAIN THE POINT OF EVALUATION. C ON OUTPUT, XP IS SET TO A NEIGHBORING POINT. C C MODE = 2. ON INPUT, FVEC MUST CONTAIN THE FUNCTIONS AND THE C ROWS OF FJAC MUST CONTAIN THE GRADIENTS C OF THE RESPECTIVE FUNCTIONS EACH EVALUATED C AT X, AND FVECP MUST CONTAIN THE FUNCTIONS C EVALUATED AT XP. C ON OUTPUT, ERR CONTAINS MEASURES OF CORRECTNESS OF C THE RESPECTIVE GRADIENTS. C C THE SUBROUTINE DOES NOT PERFORM RELIABLY IF CANCELLATION OR C ROUNDING ERRORS CAUSE A SEVERE LOSS OF SIGNIFICANCE IN THE C EVALUATION OF A FUNCTION. THEREFORE, NONE OF THE COMPONENTS C OF X SHOULD BE UNUSUALLY SMALL (IN PARTICULAR, ZERO) OR ANY C OTHER VALUE WHICH MAY CAUSE LOSS OF SIGNIFICANCE. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE CHKDER(M,N,X,FVEC,FJAC,LDFJAC,XP,FVECP,MODE,ERR) C C WHERE C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VARIABLES. C C X IS AN INPUT ARRAY OF LENGTH N. C C FVEC IS AN ARRAY OF LENGTH M. ON INPUT WHEN MODE = 2, C FVEC MUST CONTAIN THE FUNCTIONS EVALUATED AT X. C C FJAC IS AN M BY N ARRAY. ON INPUT WHEN MODE = 2, C THE ROWS OF FJAC MUST CONTAIN THE GRADIENTS OF C THE RESPECTIVE FUNCTIONS EVALUATED AT X. C C LDFJAC IS A POSITIVE INTEGER INPUT PARAMETER NOT LESS THAN M C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. C C XP IS AN ARRAY OF LENGTH N. ON OUTPUT WHEN MODE = 1, C XP IS SET TO A NEIGHBORING POINT OF X. C C FVECP IS AN ARRAY OF LENGTH M. ON INPUT WHEN MODE = 2, C FVECP MUST CONTAIN THE FUNCTIONS EVALUATED AT XP. C C MODE IS AN INTEGER INPUT VARIABLE SET TO 1 ON THE FIRST CALL C AND 2 ON THE SECOND. OTHER VALUES OF MODE ARE EQUIVALENT C TO MODE = 1. C C ERR IS AN ARRAY OF LENGTH M. ON OUTPUT WHEN MODE = 2, C ERR CONTAINS MEASURES OF CORRECTNESS OF THE RESPECTIVE C GRADIENTS. IF THERE IS NO SEVERE LOSS OF SIGNIFICANCE, C THEN IF ERR(I) IS 1.0 THE I-TH GRADIENT IS CORRECT, C WHILE IF ERR(I) IS 0.0 THE I-TH GRADIENT IS INCORRECT. C FOR VALUES OF ERR BETWEEN 0.0 AND 1.0, THE CATEGORIZATION C IS LESS CERTAIN. IN GENERAL, A VALUE OF ERR(I) GREATER C THAN 0.5 INDICATES THAT THE I-TH GRADIENT IS PROBABLY C CORRECT, WHILE A VALUE OF ERR(I) LESS THAN 0.5 INDICATES C THAT THE I-TH GRADIENT IS PROBABLY INCORRECT. C C SUBPROGRAMS CALLED C C MINPACK SUPPLIED ... SPMPAR C C FORTRAN SUPPLIED ... ABS,ALOG10,SQRT C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,J REAL EPS,EPSF,EPSLOG,EPSMCH,FACTOR,ONE,TEMP,ZERO REAL SPMPAR DATA FACTOR,ONE,ZERO /1.0E2,1.0E0,0.0E0/ C C EPSMCH IS THE MACHINE PRECISION. C EPSMCH = SPMPAR(1) C EPS = SQRT(EPSMCH) C IF (MODE .EQ. 2) GO TO 20 C C MODE = 1. C DO 10 J = 1, N TEMP = EPS*ABS(X(J)) IF (TEMP .EQ. ZERO) TEMP = EPS XP(J) = X(J) + TEMP 10 CONTINUE GO TO 70 20 CONTINUE C C MODE = 2. C EPSF = FACTOR*EPSMCH EPSLOG = ALOG10(EPS) DO 30 I = 1, M ERR(I) = ZERO 30 CONTINUE DO 50 J = 1, N TEMP = ABS(X(J)) IF (TEMP .EQ. ZERO) TEMP = ONE DO 40 I = 1, M ERR(I) = ERR(I) + TEMP*FJAC(I,J) 40 CONTINUE 50 CONTINUE DO 60 I = 1, M TEMP = ONE IF (FVEC(I) .NE. ZERO .AND. FVECP(I) .NE. ZERO * .AND. ABS(FVECP(I)-FVEC(I)) .GE. EPSF*ABS(FVEC(I))) * TEMP = EPS*ABS((FVECP(I)-FVEC(I))/EPS-ERR(I)) * /(ABS(FVEC(I)) + ABS(FVECP(I))) ERR(I) = ONE IF (TEMP .GT. EPSMCH .AND. TEMP .LT. EPS) * ERR(I) = (ALOG10(TEMP) - EPSLOG)/EPSLOG IF (TEMP .GE. EPS) ERR(I) = ZERO 60 CONTINUE 70 CONTINUE C RETURN C C LAST CARD OF SUBROUTINE CHKDER. C END SUBROUTINE DOGLEG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2) INTEGER N,LR REAL DELTA REAL R(LR),DIAG(N),QTB(N),X(N),WA1(N),WA2(N) C ********** C C SUBROUTINE DOGLEG C C GIVEN AN M BY N MATRIX A, AN N BY N NONSINGULAR DIAGONAL C MATRIX D, AN M-VECTOR B, AND A POSITIVE NUMBER DELTA, THE C PROBLEM IS TO DETERMINE THE CONVEX COMBINATION X OF THE C GAUSS-NEWTON AND SCALED GRADIENT DIRECTIONS THAT MINIMIZES C (A*X - B) IN THE LEAST SQUARES SENSE, SUBJECT TO THE C RESTRICTION THAT THE EUCLIDEAN NORM OF D*X BE AT MOST DELTA. C C THIS SUBROUTINE COMPLETES THE SOLUTION OF THE PROBLEM C IF IT IS PROVIDED WITH THE NECESSARY INFORMATION FROM THE C QR FACTORIZATION OF A. THAT IS, IF A = Q*R, WHERE Q HAS C ORTHOGONAL COLUMNS AND R IS AN UPPER TRIANGULAR MATRIX, C THEN DOGLEG EXPECTS THE FULL UPPER TRIANGLE OF R AND C THE FIRST N COMPONENTS OF (Q TRANSPOSE)*B. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE DOGLEG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER OF R. C C R IS AN INPUT ARRAY OF LENGTH LR WHICH MUST CONTAIN THE UPPER C TRIANGULAR MATRIX R STORED BY ROWS. C C LR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN C (N*(N+1))/2. C C DIAG IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE C DIAGONAL ELEMENTS OF THE MATRIX D. C C QTB IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE FIRST C N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*B. C C DELTA IS A POSITIVE INPUT VARIABLE WHICH SPECIFIES AN UPPER C BOUND ON THE EUCLIDEAN NORM OF D*X. C C X IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE DESIRED C CONVEX COMBINATION OF THE GAUSS-NEWTON DIRECTION AND THE C SCALED GRADIENT DIRECTION. C C WA1 AND WA2 ARE WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... SPMPAR,ENORM C C FORTRAN-SUPPLIED ... ABS,AMAX1,AMIN1,SQRT C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,J,JJ,JP1,K,L REAL ALPHA,BNORM,EPSMCH,GNORM,ONE,QNORM,SGNORM,SUM,TEMP,ZERO REAL SPMPAR,ENORM DATA ONE,ZERO /1.0E0,0.0E0/ C C EPSMCH IS THE MACHINE PRECISION. C EPSMCH = SPMPAR(1) C C FIRST, CALCULATE THE GAUSS-NEWTON DIRECTION. C JJ = (N*(N + 1))/2 + 1 DO 50 K = 1, N J = N - K + 1 JP1 = J + 1 JJ = JJ - K L = JJ + 1 SUM = ZERO IF (N .LT. JP1) GO TO 20 DO 10 I = JP1, N SUM = SUM + R(L)*X(I) L = L + 1 10 CONTINUE 20 CONTINUE TEMP = R(JJ) IF (TEMP .NE. ZERO) GO TO 40 L = J DO 30 I = 1, J TEMP = AMAX1(TEMP,ABS(R(L))) L = L + N - I 30 CONTINUE TEMP = EPSMCH*TEMP IF (TEMP .EQ. ZERO) TEMP = EPSMCH 40 CONTINUE X(J) = (QTB(J) - SUM)/TEMP 50 CONTINUE C C TEST WHETHER THE GAUSS-NEWTON DIRECTION IS ACCEPTABLE. C DO 60 J = 1, N WA1(J) = ZERO WA2(J) = DIAG(J)*X(J) 60 CONTINUE QNORM = ENORM(N,WA2) IF (QNORM .LE. DELTA) GO TO 140 C C THE GAUSS-NEWTON DIRECTION IS NOT ACCEPTABLE. C NEXT, CALCULATE THE SCALED GRADIENT DIRECTION. C L = 1 DO 80 J = 1, N TEMP = QTB(J) DO 70 I = J, N WA1(I) = WA1(I) + R(L)*TEMP L = L + 1 70 CONTINUE WA1(J) = WA1(J)/DIAG(J) 80 CONTINUE C C CALCULATE THE NORM OF THE SCALED GRADIENT AND TEST FOR C THE SPECIAL CASE IN WHICH THE SCALED GRADIENT IS ZERO. C GNORM = ENORM(N,WA1) SGNORM = ZERO ALPHA = DELTA/QNORM IF (GNORM .EQ. ZERO) GO TO 120 C C CALCULATE THE POINT ALONG THE SCALED GRADIENT C AT WHICH THE QUADRATIC IS MINIMIZED. C DO 90 J = 1, N WA1(J) = (WA1(J)/GNORM)/DIAG(J) 90 CONTINUE L = 1 DO 110 J = 1, N SUM = ZERO DO 100 I = J, N SUM = SUM + R(L)*WA1(I) L = L + 1 100 CONTINUE WA2(J) = SUM 110 CONTINUE TEMP = ENORM(N,WA2) SGNORM = (GNORM/TEMP)/TEMP C C TEST WHETHER THE SCALED GRADIENT DIRECTION IS ACCEPTABLE. C ALPHA = ZERO IF (SGNORM .GE. DELTA) GO TO 120 C C THE SCALED GRADIENT DIRECTION IS NOT ACCEPTABLE. C FINALLY, CALCULATE THE POINT ALONG THE DOGLEG C AT WHICH THE QUADRATIC IS MINIMIZED. C BNORM = ENORM(N,QTB) TEMP = (BNORM/GNORM)*(BNORM/QNORM)*(SGNORM/DELTA) TEMP = TEMP - (DELTA/QNORM)*(SGNORM/DELTA)**2 * + SQRT((TEMP-(DELTA/QNORM))**2 * +(ONE-(DELTA/QNORM)**2)*(ONE-(SGNORM/DELTA)**2)) ALPHA = ((DELTA/QNORM)*(ONE - (SGNORM/DELTA)**2))/TEMP 120 CONTINUE C C FORM APPROPRIATE CONVEX COMBINATION OF THE GAUSS-NEWTON C DIRECTION AND THE SCALED GRADIENT DIRECTION. C TEMP = (ONE - ALPHA)*AMIN1(SGNORM,DELTA) DO 130 J = 1, N X(J) = TEMP*WA1(J) + ALPHA*X(J) 130 CONTINUE 140 CONTINUE RETURN C C LAST CARD OF SUBROUTINE DOGLEG. C END REAL FUNCTION ENORM(N,X) INTEGER N REAL X(N) C ********** C C FUNCTION ENORM C C GIVEN AN N-VECTOR X, THIS FUNCTION CALCULATES THE C EUCLIDEAN NORM OF X. C C THE EUCLIDEAN NORM IS COMPUTED BY ACCUMULATING THE SUM OF C SQUARES IN THREE DIFFERENT SUMS. THE SUMS OF SQUARES FOR THE C SMALL AND LARGE COMPONENTS ARE SCALED SO THAT NO OVERFLOWS C OCCUR. NON-DESTRUCTIVE UNDERFLOWS ARE PERMITTED. UNDERFLOWS C AND OVERFLOWS DO NOT OCCUR IN THE COMPUTATION OF THE UNSCALED C SUM OF SQUARES FOR THE INTERMEDIATE COMPONENTS. C THE DEFINITIONS OF SMALL, INTERMEDIATE AND LARGE COMPONENTS C DEPEND ON TWO CONSTANTS, RDWARF AND RGIANT. THE MAIN C RESTRICTIONS ON THESE CONSTANTS ARE THAT RDWARF**2 NOT C UNDERFLOW AND RGIANT**2 NOT OVERFLOW. THE CONSTANTS C GIVEN HERE ARE SUITABLE FOR EVERY KNOWN COMPUTER. C C THE FUNCTION STATEMENT IS C C REAL FUNCTION ENORM(N,X) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE. C C X IS AN INPUT ARRAY OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... ABS,SQRT C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I REAL AGIANT,FLOATN,ONE,RDWARF,RGIANT,S1,S2,S3,XABS,X1MAX,X3MAX, * ZERO DATA ONE,ZERO,RDWARF,RGIANT /1.0E0,0.0E0,3.834E-20,1.304E19/ S1 = ZERO S2 = ZERO S3 = ZERO X1MAX = ZERO X3MAX = ZERO FLOATN = N AGIANT = RGIANT/FLOATN DO 90 I = 1, N XABS = ABS(X(I)) IF (XABS .GT. RDWARF .AND. XABS .LT. AGIANT) GO TO 70 IF (XABS .LE. RDWARF) GO TO 30 C C SUM FOR LARGE COMPONENTS. C IF (XABS .LE. X1MAX) GO TO 10 S1 = ONE + S1*(X1MAX/XABS)**2 X1MAX = XABS GO TO 20 10 CONTINUE S1 = S1 + (XABS/X1MAX)**2 20 CONTINUE GO TO 60 30 CONTINUE C C SUM FOR SMALL COMPONENTS. C IF (XABS .LE. X3MAX) GO TO 40 S3 = ONE + S3*(X3MAX/XABS)**2 X3MAX = XABS GO TO 50 40 CONTINUE IF (XABS .NE. ZERO) S3 = S3 + (XABS/X3MAX)**2 50 CONTINUE 60 CONTINUE GO TO 80 70 CONTINUE C C SUM FOR INTERMEDIATE COMPONENTS. C S2 = S2 + XABS**2 80 CONTINUE 90 CONTINUE C C CALCULATION OF NORM. C IF (S1 .EQ. ZERO) GO TO 100 ENORM = X1MAX*SQRT(S1+(S2/X1MAX)/X1MAX) GO TO 130 100 CONTINUE IF (S2 .EQ. ZERO) GO TO 110 IF (S2 .GE. X3MAX) * ENORM = SQRT(S2*(ONE+(X3MAX/S2)*(X3MAX*S3))) IF (S2 .LT. X3MAX) * ENORM = SQRT(X3MAX*((S2/X3MAX)+(X3MAX*S3))) GO TO 120 110 CONTINUE ENORM = X3MAX*SQRT(S3) 120 CONTINUE 130 CONTINUE RETURN C C LAST CARD OF FUNCTION ENORM. C END SUBROUTINE FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN, * WA1,WA2) INTEGER N,LDFJAC,IFLAG,ML,MU REAL EPSFCN REAL X(N),FVEC(N),FJAC(LDFJAC,N),WA1(N),WA2(N) C ********** C C SUBROUTINE FDJAC1 C C THIS SUBROUTINE COMPUTES A FORWARD-DIFFERENCE APPROXIMATION C TO THE N BY N JACOBIAN MATRIX ASSOCIATED WITH A SPECIFIED C PROBLEM OF N FUNCTIONS IN N VARIABLES. IF THE JACOBIAN HAS C A BANDED FORM, THEN FUNCTION EVALUATIONS ARE SAVED BY ONLY C APPROXIMATING THE NONZERO TERMS. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN, C WA1,WA2) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED C IN AN EXTERNAL STATEMENT IN THE USER CALLING C PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(N,X,FVEC,IFLAG) C INTEGER N,IFLAG C REAL X(N),FVEC(N) C ---------- C CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. C ---------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF FDJAC1. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS AND VARIABLES. C C X IS AN INPUT ARRAY OF LENGTH N. C C FVEC IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE C FUNCTIONS EVALUATED AT X. C C FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE C APPROXIMATION TO THE JACOBIAN MATRIX EVALUATED AT X. C C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. C C IFLAG IS AN INTEGER VARIABLE WHICH CAN BE USED TO TERMINATE C THE EXECUTION OF FDJAC1. SEE DESCRIPTION OF FCN. C C ML IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES C THE NUMBER OF SUBDIAGONALS WITHIN THE BAND OF THE C JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET C ML TO AT LEAST N - 1. C C EPSFCN IS AN INPUT VARIABLE USED IN DETERMINING A SUITABLE C STEP LENGTH FOR THE FORWARD-DIFFERENCE APPROXIMATION. THIS C APPROXIMATION ASSUMES THAT THE RELATIVE ERRORS IN THE C FUNCTIONS ARE OF THE ORDER OF EPSFCN. IF EPSFCN IS LESS C THAN THE MACHINE PRECISION, IT IS ASSUMED THAT THE RELATIVE C ERRORS IN THE FUNCTIONS ARE OF THE ORDER OF THE MACHINE C PRECISION. C C MU IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES C THE NUMBER OF SUPERDIAGONALS WITHIN THE BAND OF THE C JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET C MU TO AT LEAST N - 1. C C WA1 AND WA2 ARE WORK ARRAYS OF LENGTH N. IF ML + MU + 1 IS AT C LEAST N, THEN THE JACOBIAN IS CONSIDERED DENSE, AND WA2 IS C NOT REFERENCED. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... SPMPAR C C FORTRAN-SUPPLIED ... ABS,AMAX1,SQRT C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,J,K,MSUM REAL EPS,EPSMCH,H,TEMP,ZERO REAL SPMPAR DATA ZERO /0.0E0/ C C EPSMCH IS THE MACHINE PRECISION. C EPSMCH = SPMPAR(1) C EPS = SQRT(AMAX1(EPSFCN,EPSMCH)) MSUM = ML + MU + 1 IF (MSUM .LT. N) GO TO 40 C C COMPUTATION OF DENSE APPROXIMATE JACOBIAN. C DO 20 J = 1, N TEMP = X(J) H = EPS*ABS(TEMP) IF (H .EQ. ZERO) H = EPS X(J) = TEMP + H CALL FCN(N,X,WA1,IFLAG) IF (IFLAG .LT. 0) GO TO 30 X(J) = TEMP DO 10 I = 1, N FJAC(I,J) = (WA1(I) - FVEC(I))/H 10 CONTINUE 20 CONTINUE 30 CONTINUE GO TO 110 40 CONTINUE C C COMPUTATION OF BANDED APPROXIMATE JACOBIAN. C DO 90 K = 1, MSUM DO 60 J = K, N, MSUM WA2(J) = X(J) H = EPS*ABS(WA2(J)) IF (H .EQ. ZERO) H = EPS X(J) = WA2(J) + H 60 CONTINUE CALL FCN(N,X,WA1,IFLAG) IF (IFLAG .LT. 0) GO TO 100 DO 80 J = K, N, MSUM X(J) = WA2(J) H = EPS*ABS(WA2(J)) IF (H .EQ. ZERO) H = EPS DO 70 I = 1, N FJAC(I,J) = ZERO IF (I .GE. J - MU .AND. I .LE. J + ML) * FJAC(I,J) = (WA1(I) - FVEC(I))/H 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE 110 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FDJAC1. C END SUBROUTINE FDJAC2(FCN,M,N,X,FVEC,FJAC,LDFJAC,IFLAG,EPSFCN,WA) INTEGER M,N,LDFJAC,IFLAG REAL EPSFCN REAL X(N),FVEC(M),FJAC(LDFJAC,N),WA(M) C ********** C C SUBROUTINE FDJAC2 C C THIS SUBROUTINE COMPUTES A FORWARD-DIFFERENCE APPROXIMATION C TO THE M BY N JACOBIAN MATRIX ASSOCIATED WITH A SPECIFIED C PROBLEM OF M FUNCTIONS IN N VARIABLES. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE FDJAC2(FCN,M,N,X,FVEC,FJAC,LDFJAC,IFLAG,EPSFCN,WA) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED C IN AN EXTERNAL STATEMENT IN THE USER CALLING C PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(M,N,X,FVEC,IFLAG) C INTEGER M,N,IFLAG C REAL X(N),FVEC(M) C ---------- C CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. C ---------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF FDJAC2. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VARIABLES. N MUST NOT EXCEED M. C C X IS AN INPUT ARRAY OF LENGTH N. C C FVEC IS AN INPUT ARRAY OF LENGTH M WHICH MUST CONTAIN THE C FUNCTIONS EVALUATED AT X. C C FJAC IS AN OUTPUT M BY N ARRAY WHICH CONTAINS THE C APPROXIMATION TO THE JACOBIAN MATRIX EVALUATED AT X. C C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. C C IFLAG IS AN INTEGER VARIABLE WHICH CAN BE USED TO TERMINATE C THE EXECUTION OF FDJAC2. SEE DESCRIPTION OF FCN. C C EPSFCN IS AN INPUT VARIABLE USED IN DETERMINING A SUITABLE C STEP LENGTH FOR THE FORWARD-DIFFERENCE APPROXIMATION. THIS C APPROXIMATION ASSUMES THAT THE RELATIVE ERRORS IN THE C FUNCTIONS ARE OF THE ORDER OF EPSFCN. IF EPSFCN IS LESS C THAN THE MACHINE PRECISION, IT IS ASSUMED THAT THE RELATIVE C ERRORS IN THE FUNCTIONS ARE OF THE ORDER OF THE MACHINE C PRECISION. C C WA IS A WORK ARRAY OF LENGTH M. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ...... FCN C C MINPACK-SUPPLIED ... SPMPAR C C FORTRAN-SUPPLIED ... ABS,AMAX1,SQRT C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,J REAL EPS,EPSMCH,H,TEMP,ZERO REAL SPMPAR DATA ZERO /0.0E0/ C C EPSMCH IS THE MACHINE PRECISION. C EPSMCH = SPMPAR(1) C EPS = SQRT(AMAX1(EPSFCN,EPSMCH)) DO 20 J = 1, N TEMP = X(J) H = EPS*ABS(TEMP) IF (H .EQ. ZERO) H = EPS X(J) = TEMP + H CALL FCN(M,N,X,WA,IFLAG) IF (IFLAG .LT. 0) GO TO 30 X(J) = TEMP DO 10 I = 1, M FJAC(I,J) = (WA(I) - FVEC(I))/H 10 CONTINUE 20 CONTINUE 30 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FDJAC2. C END SUBROUTINE HYBRD(FCN,N,X,FVEC,XTOL,MAXFEV,ML,MU,EPSFCN,DIAG, * MODE,FACTOR,NPRINT,INFO,NFEV,FJAC,LDFJAC,R,LR, * QTF,WA1,WA2,WA3,WA4) INTEGER N,MAXFEV,ML,MU,MODE,NPRINT,INFO,NFEV,LDFJAC,LR REAL XTOL,EPSFCN,FACTOR REAL X(N),FVEC(N),DIAG(N),FJAC(LDFJAC,N),R(LR),QTF(N),WA1(N), * WA2(N),WA3(N),WA4(N) EXTERNAL FCN C ********** C C SUBROUTINE HYBRD C C THE PURPOSE OF HYBRD IS TO FIND A ZERO OF A SYSTEM OF C N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION C OF THE POWELL HYBRID METHOD. THE USER MUST PROVIDE A C SUBROUTINE WHICH CALCULATES THE FUNCTIONS. THE JACOBIAN IS C THEN CALCULATED BY A FORWARD-DIFFERENCE APPROXIMATION. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE HYBRD(FCN,N,X,FVEC,XTOL,MAXFEV,ML,MU,EPSFCN, C DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,FJAC, C LDFJAC,R,LR,QTF,WA1,WA2,WA3,WA4) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED C IN AN EXTERNAL STATEMENT IN THE USER CALLING C PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(N,X,FVEC,IFLAG) C INTEGER N,IFLAG C REAL X(N),FVEC(N) C ---------- C CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. C --------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF HYBRD. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS AND VARIABLES. C C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. C C FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS C THE FUNCTIONS EVALUATED AT THE OUTPUT X. C C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION C OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE C ITERATES IS AT MOST XTOL. C C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION C OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST MAXFEV C BY THE END OF AN ITERATION. C C ML IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES C THE NUMBER OF SUBDIAGONALS WITHIN THE BAND OF THE C JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET C ML TO AT LEAST N - 1. C C MU IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES C THE NUMBER OF SUPERDIAGONALS WITHIN THE BAND OF THE C JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET C MU TO AT LEAST N - 1. C C EPSFCN IS AN INPUT VARIABLE USED IN DETERMINING A SUITABLE C STEP LENGTH FOR THE FORWARD-DIFFERENCE APPROXIMATION. THIS C APPROXIMATION ASSUMES THAT THE RELATIVE ERRORS IN THE C FUNCTIONS ARE OF THE ORDER OF EPSFCN. IF EPSFCN IS LESS C THAN THE MACHINE PRECISION, IT IS ASSUMED THAT THE RELATIVE C ERRORS IN THE FUNCTIONS ARE OF THE ORDER OF THE MACHINE C PRECISION. C C DIAG IS AN ARRAY OF LENGTH N. IF MODE = 1 (SEE C BELOW), DIAG IS INTERNALLY SET. IF MODE = 2, DIAG C MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS C MULTIPLICATIVE SCALE FACTORS FOR THE VARIABLES. C C MODE IS AN INTEGER INPUT VARIABLE. IF MODE = 1, THE C VARIABLES WILL BE SCALED INTERNALLY. IF MODE = 2, C THE SCALING IS SPECIFIED BY THE INPUT DIAG. OTHER C VALUES OF MODE ARE EQUIVALENT TO MODE = 1. C C FACTOR IS A POSITIVE INPUT VARIABLE USED IN DETERMINING THE C INITIAL STEP BOUND. THIS BOUND IS SET TO THE PRODUCT OF C FACTOR AND THE EUCLIDEAN NORM OF DIAG*X IF NONZERO, OR ELSE C TO FACTOR ITSELF. IN MOST CASES FACTOR SHOULD LIE IN THE C INTERVAL (.1,100.). 100. IS A GENERALLY RECOMMENDED VALUE. C C NPRINT IS AN INTEGER INPUT VARIABLE THAT ENABLES CONTROLLED C PRINTING OF ITERATES IF IT IS POSITIVE. IN THIS CASE, C FCN IS CALLED WITH IFLAG = 0 AT THE BEGINNING OF THE FIRST C ITERATION AND EVERY NPRINT ITERATIONS THEREAFTER AND C IMMEDIATELY PRIOR TO RETURN, WITH X AND FVEC AVAILABLE C FOR PRINTING. IF NPRINT IS NOT POSITIVE, NO SPECIAL CALLS C OF FCN WITH IFLAG = 0 ARE MADE. C C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE) C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE, C INFO IS SET AS FOLLOWS. C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO = 1 RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES C IS AT MOST XTOL. C C INFO = 2 NUMBER OF CALLS TO FCN HAS REACHED OR EXCEEDED C MAXFEV. C C INFO = 3 XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN C THE APPROXIMATE SOLUTION X IS POSSIBLE. C C INFO = 4 ITERATION IS NOT MAKING GOOD PROGRESS, AS C MEASURED BY THE IMPROVEMENT FROM THE LAST C FIVE JACOBIAN EVALUATIONS. C C INFO = 5 ITERATION IS NOT MAKING GOOD PROGRESS, AS C MEASURED BY THE IMPROVEMENT FROM THE LAST C TEN ITERATIONS. C C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF C CALLS TO FCN. C C FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE C ORTHOGONAL MATRIX Q PRODUCED BY THE QR FACTORIZATION C OF THE FINAL APPROXIMATE JACOBIAN. C C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. C C R IS AN OUTPUT ARRAY OF LENGTH LR WHICH CONTAINS THE C UPPER TRIANGULAR MATRIX PRODUCED BY THE QR FACTORIZATION C OF THE FINAL APPROXIMATE JACOBIAN, STORED ROWWISE. C C LR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN C (N*(N+1))/2. C C QTF IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS C THE VECTOR (Q TRANSPOSE)*FVEC. C C WA1, WA2, WA3, AND WA4 ARE WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ...... FCN C C MINPACK-SUPPLIED ... DOGLEG,SPMPAR,ENORM,FDJAC1, C QFORM,QRFAC,R1MPYQ,R1UPDT C C FORTRAN-SUPPLIED ... ABS,AMAX1,AMIN1,MIN0,MOD C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,IFLAG,ITER,J,JM1,L,MSUM,NCFAIL,NCSUC,NSLOW1,NSLOW2 INTEGER IWA(1) LOGICAL JEVAL,SING REAL ACTRED,DELTA,EPSMCH,FNORM,FNORM1,ONE,PNORM,PRERED,P1,P5, * P001,P0001,RATIO,SUM,TEMP,XNORM,ZERO REAL SPMPAR,ENORM DATA ONE,P1,P5,P001,P0001,ZERO * /1.0E0,1.0E-1,5.0E-1,1.0E-3,1.0E-4,0.0E0/ C C EPSMCH IS THE MACHINE PRECISION. C EPSMCH = SPMPAR(1) C INFO = 0 IFLAG = 0 NFEV = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. XTOL .LT. ZERO .OR. MAXFEV .LE. 0 * .OR. ML .LT. 0 .OR. MU .LT. 0 .OR. FACTOR .LE. ZERO * .OR. LDFJAC .LT. N .OR. LR .LT. (N*(N + 1))/2) GO TO 300 IF (MODE .NE. 2) GO TO 20 DO 10 J = 1, N IF (DIAG(J) .LE. ZERO) GO TO 300 10 CONTINUE 20 CONTINUE C C EVALUATE THE FUNCTION AT THE STARTING POINT C AND CALCULATE ITS NORM. C IFLAG = 1 CALL FCN(N,X,FVEC,IFLAG) NFEV = 1 IF (IFLAG .LT. 0) GO TO 300 FNORM = ENORM(N,FVEC) C C DETERMINE THE NUMBER OF CALLS TO FCN NEEDED TO COMPUTE C THE JACOBIAN MATRIX. C MSUM = MIN0(ML+MU+1,N) C C INITIALIZE ITERATION COUNTER AND MONITORS. C ITER = 1 NCSUC = 0 NCFAIL = 0 NSLOW1 = 0 NSLOW2 = 0 C C BEGINNING OF THE OUTER LOOP. C 30 CONTINUE JEVAL = .TRUE. C C CALCULATE THE JACOBIAN MATRIX. C IFLAG = 2 CALL FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN,WA1, * WA2) NFEV = NFEV + MSUM IF (IFLAG .LT. 0) GO TO 300 C C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN. C CALL QRFAC(N,N,FJAC,LDFJAC,.FALSE.,IWA,1,WA1,WA2,WA3) C C ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING C TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN. C IF (ITER .NE. 1) GO TO 70 IF (MODE .EQ. 2) GO TO 50 DO 40 J = 1, N DIAG(J) = WA2(J) IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE 40 CONTINUE 50 CONTINUE C C ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X C AND INITIALIZE THE STEP BOUND DELTA. C DO 60 J = 1, N WA3(J) = DIAG(J)*X(J) 60 CONTINUE XNORM = ENORM(N,WA3) DELTA = FACTOR*XNORM IF (DELTA .EQ. ZERO) DELTA = FACTOR 70 CONTINUE C C FORM (Q TRANSPOSE)*FVEC AND STORE IN QTF. C DO 80 I = 1, N QTF(I) = FVEC(I) 80 CONTINUE DO 120 J = 1, N IF (FJAC(J,J) .EQ. ZERO) GO TO 110 SUM = ZERO DO 90 I = J, N SUM = SUM + FJAC(I,J)*QTF(I) 90 CONTINUE TEMP = -SUM/FJAC(J,J) DO 100 I = J, N QTF(I) = QTF(I) + FJAC(I,J)*TEMP 100 CONTINUE 110 CONTINUE 120 CONTINUE C C COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R. C SING = .FALSE. DO 150 J = 1, N L = J JM1 = J - 1 IF (JM1 .LT. 1) GO TO 140 DO 130 I = 1, JM1 R(L) = FJAC(I,J) L = L + N - I 130 CONTINUE 140 CONTINUE R(L) = WA1(J) IF (WA1(J) .EQ. ZERO) SING = .TRUE. 150 CONTINUE C C ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC. C CALL QFORM(N,N,FJAC,LDFJAC,WA1) C C RESCALE IF NECESSARY. C IF (MODE .EQ. 2) GO TO 170 DO 160 J = 1, N DIAG(J) = AMAX1(DIAG(J),WA2(J)) 160 CONTINUE 170 CONTINUE C C BEGINNING OF THE INNER LOOP. C 180 CONTINUE C C IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES. C IF (NPRINT .LE. 0) GO TO 190 IFLAG = 0 IF (MOD(ITER-1,NPRINT) .EQ. 0) CALL FCN(N,X,FVEC,IFLAG) IF (IFLAG .LT. 0) GO TO 300 190 CONTINUE C C DETERMINE THE DIRECTION P. C CALL DOGLEG(N,R,LR,DIAG,QTF,DELTA,WA1,WA2,WA3) C C STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P. C DO 200 J = 1, N WA1(J) = -WA1(J) WA2(J) = X(J) + WA1(J) WA3(J) = DIAG(J)*WA1(J) 200 CONTINUE PNORM = ENORM(N,WA3) C C ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND. C IF (ITER .EQ. 1) DELTA = AMIN1(DELTA,PNORM) C C EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM. C IFLAG = 1 CALL FCN(N,WA2,WA4,IFLAG) NFEV = NFEV + 1 IF (IFLAG .LT. 0) GO TO 300 FNORM1 = ENORM(N,WA4) C C COMPUTE THE SCALED ACTUAL REDUCTION. C ACTRED = -ONE IF (FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2 C C COMPUTE THE SCALED PREDICTED REDUCTION. C L = 1 DO 220 I = 1, N SUM = ZERO DO 210 J = I, N SUM = SUM + R(L)*WA1(J) L = L + 1 210 CONTINUE WA3(I) = QTF(I) + SUM 220 CONTINUE TEMP = ENORM(N,WA3) PRERED = ZERO IF (TEMP .LT. FNORM) PRERED = ONE - (TEMP/FNORM)**2 C C COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED C REDUCTION. C RATIO = ZERO IF (PRERED .GT. ZERO) RATIO = ACTRED/PRERED C C UPDATE THE STEP BOUND. C IF (RATIO .GE. P1) GO TO 230 NCSUC = 0 NCFAIL = NCFAIL + 1 DELTA = P5*DELTA GO TO 240 230 CONTINUE NCFAIL = 0 NCSUC = NCSUC + 1 IF (RATIO .GE. P5 .OR. NCSUC .GT. 1) * DELTA = AMAX1(DELTA,PNORM/P5) IF (ABS(RATIO-ONE) .LE. P1) DELTA = PNORM/P5 240 CONTINUE C C TEST FOR SUCCESSFUL ITERATION. C IF (RATIO .LT. P0001) GO TO 260 C C SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS. C DO 250 J = 1, N X(J) = WA2(J) WA2(J) = DIAG(J)*X(J) FVEC(J) = WA4(J) 250 CONTINUE XNORM = ENORM(N,WA2) FNORM = FNORM1 ITER = ITER + 1 260 CONTINUE C C DETERMINE THE PROGRESS OF THE ITERATION. C NSLOW1 = NSLOW1 + 1 IF (ACTRED .GE. P001) NSLOW1 = 0 IF (JEVAL) NSLOW2 = NSLOW2 + 1 IF (ACTRED .GE. P1) NSLOW2 = 0 C C TEST FOR CONVERGENCE. C IF (DELTA .LE. XTOL*XNORM .OR. FNORM .EQ. ZERO) INFO = 1 IF (INFO .NE. 0) GO TO 300 C C TESTS FOR TERMINATION AND STRINGENT TOLERANCES. C IF (NFEV .GE. MAXFEV) INFO = 2 IF (P1*AMAX1(P1*DELTA,PNORM) .LE. EPSMCH*XNORM) INFO = 3 IF (NSLOW2 .EQ. 5) INFO = 4 IF (NSLOW1 .EQ. 10) INFO = 5 IF (INFO .NE. 0) GO TO 300 C C CRITERION FOR RECALCULATING JACOBIAN APPROXIMATION C BY FORWARD DIFFERENCES. C IF (NCFAIL .EQ. 2) GO TO 290 C C CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN C AND UPDATE QTF IF NECESSARY. C DO 280 J = 1, N SUM = ZERO DO 270 I = 1, N SUM = SUM + FJAC(I,J)*WA4(I) 270 CONTINUE WA2(J) = (SUM - WA3(J))/PNORM WA1(J) = DIAG(J)*((DIAG(J)*WA1(J))/PNORM) IF (RATIO .GE. P0001) QTF(J) = SUM 280 CONTINUE C C COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN. C CALL R1UPDT(N,N,R,LR,WA1,WA2,WA3,SING) CALL R1MPYQ(N,N,FJAC,LDFJAC,WA2,WA3) CALL R1MPYQ(1,N,QTF,1,WA2,WA3) C C END OF THE INNER LOOP. C JEVAL = .FALSE. GO TO 180 290 CONTINUE C C END OF THE OUTER LOOP. C GO TO 30 300 CONTINUE C C TERMINATION, EITHER NORMAL OR USER IMPOSED. C IF (IFLAG .LT. 0) INFO = IFLAG IFLAG = 0 IF (NPRINT .GT. 0) CALL FCN(N,X,FVEC,IFLAG) RETURN C C LAST CARD OF SUBROUTINE HYBRD. C END SUBROUTINE HYBRD1(FCN,N,X,FVEC,TOL,INFO,WA,LWA) INTEGER N,INFO,LWA REAL TOL REAL X(N),FVEC(N),WA(LWA) EXTERNAL FCN C ********** C C SUBROUTINE HYBRD1 C C THE PURPOSE OF HYBRD1 IS TO FIND A ZERO OF A SYSTEM OF C N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION C OF THE POWELL HYBRID METHOD. THIS IS DONE BY USING THE C MORE GENERAL NONLINEAR EQUATION SOLVER HYBRD. THE USER C MUST PROVIDE A SUBROUTINE WHICH CALCULATES THE FUNCTIONS. C THE JACOBIAN IS THEN CALCULATED BY A FORWARD-DIFFERENCE C APPROXIMATION. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE HYBRD1(FCN,N,X,FVEC,TOL,INFO,WA,LWA) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED C IN AN EXTERNAL STATEMENT IN THE USER CALLING C PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(N,X,FVEC,IFLAG) C INTEGER N,IFLAG C REAL X(N),FVEC(N) C ---------- C CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. C --------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF HYBRD1. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS AND VARIABLES. C C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. C C FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS C THE FUNCTIONS EVALUATED AT THE OUTPUT X. C C TOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS C WHEN THE ALGORITHM ESTIMATES THAT THE RELATIVE ERROR C BETWEEN X AND THE SOLUTION IS AT MOST TOL. C C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE) C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE, C INFO IS SET AS FOLLOWS. C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO = 1 ALGORITHM ESTIMATES THAT THE RELATIVE ERROR C BETWEEN X AND THE SOLUTION IS AT MOST TOL. C C INFO = 2 NUMBER OF CALLS TO FCN HAS REACHED OR EXCEEDED C 200*(N+1). C C INFO = 3 TOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN C THE APPROXIMATE SOLUTION X IS POSSIBLE. C C INFO = 4 ITERATION IS NOT MAKING GOOD PROGRESS. C C WA IS A WORK ARRAY OF LENGTH LWA. C C LWA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN C (N*(3*N+13))/2. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ...... FCN C C MINPACK-SUPPLIED ... HYBRD C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER INDEX,J,LR,MAXFEV,ML,MODE,MU,NFEV,NPRINT REAL EPSFCN,FACTOR,ONE,XTOL,ZERO DATA FACTOR,ONE,ZERO /1.0E2,1.0E0,0.0E0/ INFO = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. TOL .LT. ZERO .OR. LWA .LT. (N*(3*N + 13))/2) * GO TO 20 C C CALL HYBRD. C MAXFEV = 200*(N + 1) XTOL = TOL ML = N - 1 MU = N - 1 EPSFCN = ZERO MODE = 2 DO 10 J = 1, N WA(J) = ONE 10 CONTINUE NPRINT = 0 LR = (N*(N + 1))/2 INDEX = 6*N + LR CALL HYBRD(FCN,N,X,FVEC,XTOL,MAXFEV,ML,MU,EPSFCN,WA(1),MODE, * FACTOR,NPRINT,INFO,NFEV,WA(INDEX+1),N,WA(6*N+1),LR, * WA(N+1),WA(2*N+1),WA(3*N+1),WA(4*N+1),WA(5*N+1)) IF (INFO .EQ. 5) INFO = 4 20 CONTINUE RETURN C C LAST CARD OF SUBROUTINE HYBRD1. C END SUBROUTINE HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,DIAG,MODE, * FACTOR,NPRINT,INFO,NFEV,NJEV,R,LR,QTF,WA1,WA2, * WA3,WA4) INTEGER N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,LR REAL XTOL,FACTOR REAL X(N),FVEC(N),FJAC(LDFJAC,N),DIAG(N),R(LR),QTF(N),WA1(N), * WA2(N),WA3(N),WA4(N) C ********** C C SUBROUTINE HYBRJ C C THE PURPOSE OF HYBRJ IS TO FIND A ZERO OF A SYSTEM OF C N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION C OF THE POWELL HYBRID METHOD. THE USER MUST PROVIDE A C SUBROUTINE WHICH CALCULATES THE FUNCTIONS AND THE JACOBIAN. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,DIAG, C MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,R,LR,QTF, C WA1,WA2,WA3,WA4) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS AND THE JACOBIAN. FCN MUST C BE DECLARED IN AN EXTERNAL STATEMENT IN THE USER C CALLING PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG) C INTEGER N,LDFJAC,IFLAG C REAL X(N),FVEC(N),FJAC(LDFJAC,N) C ---------- C IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC. C IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND C RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC. C --------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF HYBRJ. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS AND VARIABLES. C C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. C C FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS C THE FUNCTIONS EVALUATED AT THE OUTPUT X. C C FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE C ORTHOGONAL MATRIX Q PRODUCED BY THE QR FACTORIZATION C OF THE FINAL APPROXIMATE JACOBIAN. C C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. C C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION C OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE C ITERATES IS AT MOST XTOL. C C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION C OCCURS WHEN THE NUMBER OF CALLS TO FCN WITH IFLAG = 1 C HAS REACHED MAXFEV. C C DIAG IS AN ARRAY OF LENGTH N. IF MODE = 1 (SEE C BELOW), DIAG IS INTERNALLY SET. IF MODE = 2, DIAG C MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS C MULTIPLICATIVE SCALE FACTORS FOR THE VARIABLES. C C MODE IS AN INTEGER INPUT VARIABLE. IF MODE = 1, THE C VARIABLES WILL BE SCALED INTERNALLY. IF MODE = 2, C THE SCALING IS SPECIFIED BY THE INPUT DIAG. OTHER C VALUES OF MODE ARE EQUIVALENT TO MODE = 1. C C FACTOR IS A POSITIVE INPUT VARIABLE USED IN DETERMINING THE C INITIAL STEP BOUND. THIS BOUND IS SET TO THE PRODUCT OF C FACTOR AND THE EUCLIDEAN NORM OF DIAG*X IF NONZERO, OR ELSE C TO FACTOR ITSELF. IN MOST CASES FACTOR SHOULD LIE IN THE C INTERVAL (.1,100.). 100. IS A GENERALLY RECOMMENDED VALUE. C C NPRINT IS AN INTEGER INPUT VARIABLE THAT ENABLES CONTROLLED C PRINTING OF ITERATES IF IT IS POSITIVE. IN THIS CASE, C FCN IS CALLED WITH IFLAG = 0 AT THE BEGINNING OF THE FIRST C ITERATION AND EVERY NPRINT ITERATIONS THEREAFTER AND C IMMEDIATELY PRIOR TO RETURN, WITH X AND FVEC AVAILABLE C FOR PRINTING. FVEC AND FJAC SHOULD NOT BE ALTERED. C IF NPRINT IS NOT POSITIVE, NO SPECIAL CALLS OF FCN C WITH IFLAG = 0 ARE MADE. C C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE) C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE, C INFO IS SET AS FOLLOWS. C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO = 1 RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES C IS AT MOST XTOL. C C INFO = 2 NUMBER OF CALLS TO FCN WITH IFLAG = 1 HAS C REACHED MAXFEV. C C INFO = 3 XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN C THE APPROXIMATE SOLUTION X IS POSSIBLE. C C INFO = 4 ITERATION IS NOT MAKING GOOD PROGRESS, AS C MEASURED BY THE IMPROVEMENT FROM THE LAST C FIVE JACOBIAN EVALUATIONS. C C INFO = 5 ITERATION IS NOT MAKING GOOD PROGRESS, AS C MEASURED BY THE IMPROVEMENT FROM THE LAST C TEN ITERATIONS. C C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF C CALLS TO FCN WITH IFLAG = 1. C C NJEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF C CALLS TO FCN WITH IFLAG = 2. C C R IS AN OUTPUT ARRAY OF LENGTH LR WHICH CONTAINS THE C UPPER TRIANGULAR MATRIX PRODUCED BY THE QR FACTORIZATION C OF THE FINAL APPROXIMATE JACOBIAN, STORED ROWWISE. C C LR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN C (N*(N+1))/2. C C QTF IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS C THE VECTOR (Q TRANSPOSE)*FVEC. C C WA1, WA2, WA3, AND WA4 ARE WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ...... FCN C C MINPACK-SUPPLIED ... DOGLEG,SPMPAR,ENORM, C QFORM,QRFAC,R1MPYQ,R1UPDT C C FORTRAN-SUPPLIED ... ABS,AMAX1,AMIN1,MOD C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,IFLAG,ITER,J,JM1,L,NCFAIL,NCSUC,NSLOW1,NSLOW2 INTEGER IWA(1) LOGICAL JEVAL,SING REAL ACTRED,DELTA,EPSMCH,FNORM,FNORM1,ONE,PNORM,PRERED,P1,P5, * P001,P0001,RATIO,SUM,TEMP,XNORM,ZERO REAL SPMPAR,ENORM DATA ONE,P1,P5,P001,P0001,ZERO * /1.0E0,1.0E-1,5.0E-1,1.0E-3,1.0E-4,0.0E0/ C C EPSMCH IS THE MACHINE PRECISION. C EPSMCH = SPMPAR(1) C INFO = 0 IFLAG = 0 NFEV = 0 NJEV = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. LDFJAC .LT. N .OR. XTOL .LT. ZERO * .OR. MAXFEV .LE. 0 .OR. FACTOR .LE. ZERO * .OR. LR .LT. (N*(N + 1))/2) GO TO 300 IF (MODE .NE. 2) GO TO 20 DO 10 J = 1, N IF (DIAG(J) .LE. ZERO) GO TO 300 10 CONTINUE 20 CONTINUE C C EVALUATE THE FUNCTION AT THE STARTING POINT C AND CALCULATE ITS NORM. C IFLAG = 1 CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG) NFEV = 1 IF (IFLAG .LT. 0) GO TO 300 FNORM = ENORM(N,FVEC) C C INITIALIZE ITERATION COUNTER AND MONITORS. C ITER = 1 NCSUC = 0 NCFAIL = 0 NSLOW1 = 0 NSLOW2 = 0 C C BEGINNING OF THE OUTER LOOP. C 30 CONTINUE JEVAL = .TRUE. C C CALCULATE THE JACOBIAN MATRIX. C IFLAG = 2 CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG) NJEV = NJEV + 1 IF (IFLAG .LT. 0) GO TO 300 C C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN. C CALL QRFAC(N,N,FJAC,LDFJAC,.FALSE.,IWA,1,WA1,WA2,WA3) C C ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING C TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN. C IF (ITER .NE. 1) GO TO 70 IF (MODE .EQ. 2) GO TO 50 DO 40 J = 1, N DIAG(J) = WA2(J) IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE 40 CONTINUE 50 CONTINUE C C ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X C AND INITIALIZE THE STEP BOUND DELTA. C DO 60 J = 1, N WA3(J) = DIAG(J)*X(J) 60 CONTINUE XNORM = ENORM(N,WA3) DELTA = FACTOR*XNORM IF (DELTA .EQ. ZERO) DELTA = FACTOR 70 CONTINUE C C FORM (Q TRANSPOSE)*FVEC AND STORE IN QTF. C DO 80 I = 1, N QTF(I) = FVEC(I) 80 CONTINUE DO 120 J = 1, N IF (FJAC(J,J) .EQ. ZERO) GO TO 110 SUM = ZERO DO 90 I = J, N SUM = SUM + FJAC(I,J)*QTF(I) 90 CONTINUE TEMP = -SUM/FJAC(J,J) DO 100 I = J, N QTF(I) = QTF(I) + FJAC(I,J)*TEMP 100 CONTINUE 110 CONTINUE 120 CONTINUE C C COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R. C SING = .FALSE. DO 150 J = 1, N L = J JM1 = J - 1 IF (JM1 .LT. 1) GO TO 140 DO 130 I = 1, JM1 R(L) = FJAC(I,J) L = L + N - I 130 CONTINUE 140 CONTINUE R(L) = WA1(J) IF (WA1(J) .EQ. ZERO) SING = .TRUE. 150 CONTINUE C C ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC. C CALL QFORM(N,N,FJAC,LDFJAC,WA1) C C RESCALE IF NECESSARY. C IF (MODE .EQ. 2) GO TO 170 DO 160 J = 1, N DIAG(J) = AMAX1(DIAG(J),WA2(J)) 160 CONTINUE 170 CONTINUE C C BEGINNING OF THE INNER LOOP. C 180 CONTINUE C C IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES. C IF (NPRINT .LE. 0) GO TO 190 IFLAG = 0 IF (MOD(ITER-1,NPRINT) .EQ. 0) * CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG) IF (IFLAG .LT. 0) GO TO 300 190 CONTINUE C C DETERMINE THE DIRECTION P. C CALL DOGLEG(N,R,LR,DIAG,QTF,DELTA,WA1,WA2,WA3) C C STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P. C DO 200 J = 1, N WA1(J) = -WA1(J) WA2(J) = X(J) + WA1(J) WA3(J) = DIAG(J)*WA1(J) 200 CONTINUE PNORM = ENORM(N,WA3) C C ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND. C IF (ITER .EQ. 1) DELTA = AMIN1(DELTA,PNORM) C C EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM. C IFLAG = 1 CALL FCN(N,WA2,WA4,FJAC,LDFJAC,IFLAG) NFEV = NFEV + 1 IF (IFLAG .LT. 0) GO TO 300 FNORM1 = ENORM(N,WA4) C C COMPUTE THE SCALED ACTUAL REDUCTION. C ACTRED = -ONE IF (FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2 C C COMPUTE THE SCALED PREDICTED REDUCTION. C L = 1 DO 220 I = 1, N SUM = ZERO DO 210 J = I, N SUM = SUM + R(L)*WA1(J) L = L + 1 210 CONTINUE WA3(I) = QTF(I) + SUM 220 CONTINUE TEMP = ENORM(N,WA3) PRERED = ZERO IF (TEMP .LT. FNORM) PRERED = ONE - (TEMP/FNORM)**2 C C COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED C REDUCTION. C RATIO = ZERO IF (PRERED .GT. ZERO) RATIO = ACTRED/PRERED C C UPDATE THE STEP BOUND. C IF (RATIO .GE. P1) GO TO 230 NCSUC = 0 NCFAIL = NCFAIL + 1 DELTA = P5*DELTA GO TO 240 230 CONTINUE NCFAIL = 0 NCSUC = NCSUC + 1 IF (RATIO .GE. P5 .OR. NCSUC .GT. 1) * DELTA = AMAX1(DELTA,PNORM/P5) IF (ABS(RATIO-ONE) .LE. P1) DELTA = PNORM/P5 240 CONTINUE C C TEST FOR SUCCESSFUL ITERATION. C IF (RATIO .LT. P0001) GO TO 260 C C SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS. C DO 250 J = 1, N X(J) = WA2(J) WA2(J) = DIAG(J)*X(J) FVEC(J) = WA4(J) 250 CONTINUE XNORM = ENORM(N,WA2) FNORM = FNORM1 ITER = ITER + 1 260 CONTINUE C C DETERMINE THE PROGRESS OF THE ITERATION. C NSLOW1 = NSLOW1 + 1 IF (ACTRED .GE. P001) NSLOW1 = 0 IF (JEVAL) NSLOW2 = NSLOW2 + 1 IF (ACTRED .GE. P1) NSLOW2 = 0 C C TEST FOR CONVERGENCE. C IF (DELTA .LE. XTOL*XNORM .OR. FNORM .EQ. ZERO) INFO = 1 IF (INFO .NE. 0) GO TO 300 C C TESTS FOR TERMINATION AND STRINGENT TOLERANCES. C IF (NFEV .GE. MAXFEV) INFO = 2 IF (P1*AMAX1(P1*DELTA,PNORM) .LE. EPSMCH*XNORM) INFO = 3 IF (NSLOW2 .EQ. 5) INFO = 4 IF (NSLOW1 .EQ. 10) INFO = 5 IF (INFO .NE. 0) GO TO 300 C C CRITERION FOR RECALCULATING JACOBIAN. C IF (NCFAIL .EQ. 2) GO TO 290 C C CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN C AND UPDATE QTF IF NECESSARY. C DO 280 J = 1, N SUM = ZERO DO 270 I = 1, N SUM = SUM + FJAC(I,J)*WA4(I) 270 CONTINUE WA2(J) = (SUM - WA3(J))/PNORM WA1(J) = DIAG(J)*((DIAG(J)*WA1(J))/PNORM) IF (RATIO .GE. P0001) QTF(J) = SUM 280 CONTINUE C C COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN. C CALL R1UPDT(N,N,R,LR,WA1,WA2,WA3,SING) CALL R1MPYQ(N,N,FJAC,LDFJAC,WA2,WA3) CALL R1MPYQ(1,N,QTF,1,WA2,WA3) C C END OF THE INNER LOOP. C JEVAL = .FALSE. GO TO 180 290 CONTINUE C C END OF THE OUTER LOOP. C GO TO 30 300 CONTINUE C C TERMINATION, EITHER NORMAL OR USER IMPOSED. C IF (IFLAG .LT. 0) INFO = IFLAG IFLAG = 0 IF (NPRINT .GT. 0) CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG) RETURN C C LAST CARD OF SUBROUTINE HYBRJ. C END SUBROUTINE HYBRJ1(FCN,N,X,FVEC,FJAC,LDFJAC,TOL,INFO,WA,LWA) INTEGER N,LDFJAC,INFO,LWA REAL TOL REAL X(N),FVEC(N),FJAC(LDFJAC,N),WA(LWA) EXTERNAL FCN C ********** C C SUBROUTINE HYBRJ1 C C THE PURPOSE OF HYBRJ1 IS TO FIND A ZERO OF A SYSTEM OF C N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION C OF THE POWELL HYBRID METHOD. THIS IS DONE BY USING THE C MORE GENERAL NONLINEAR EQUATION SOLVER HYBRJ. THE USER C MUST PROVIDE A SUBROUTINE WHICH CALCULATES THE FUNCTIONS C AND THE JACOBIAN. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE HYBRJ1(FCN,N,X,FVEC,FJAC,LDFJAC,TOL,INFO,WA,LWA) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS AND THE JACOBIAN. FCN MUST C BE DECLARED IN AN EXTERNAL STATEMENT IN THE USER C CALLING PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG) C INTEGER N,LDFJAC,IFLAG C REAL X(N),FVEC(N),FJAC(LDFJAC,N) C ---------- C IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC. C IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND C RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC. C --------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF HYBRJ1. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS AND VARIABLES. C C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. C C FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS C THE FUNCTIONS EVALUATED AT THE OUTPUT X. C C FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE C ORTHOGONAL MATRIX Q PRODUCED BY THE QR FACTORIZATION C OF THE FINAL APPROXIMATE JACOBIAN. C C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. C C TOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS C WHEN THE ALGORITHM ESTIMATES THAT THE RELATIVE ERROR C BETWEEN X AND THE SOLUTION IS AT MOST TOL. C C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE) C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE, C INFO IS SET AS FOLLOWS. C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO = 1 ALGORITHM ESTIMATES THAT THE RELATIVE ERROR C BETWEEN X AND THE SOLUTION IS AT MOST TOL. C C INFO = 2 NUMBER OF CALLS TO FCN WITH IFLAG = 1 HAS C REACHED 100*(N+1). C C INFO = 3 TOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN C THE APPROXIMATE SOLUTION X IS POSSIBLE. C C INFO = 4 ITERATION IS NOT MAKING GOOD PROGRESS. C C WA IS A WORK ARRAY OF LENGTH LWA. C C LWA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN C (N*(N+13))/2. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ...... FCN C C MINPACK-SUPPLIED ... HYBRJ C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER J,LR,MAXFEV,MODE,NFEV,NJEV,NPRINT REAL FACTOR,ONE,XTOL,ZERO DATA FACTOR,ONE,ZERO /1.0E2,1.0E0,0.0E0/ INFO = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. LDFJAC .LT. N .OR. TOL .LT. ZERO * .OR. LWA .LT. (N*(N + 13))/2) GO TO 20 C C CALL HYBRJ. C MAXFEV = 100*(N + 1) XTOL = TOL MODE = 2 DO 10 J = 1, N WA(J) = ONE 10 CONTINUE NPRINT = 0 LR = (N*(N + 1))/2 CALL HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,WA(1),MODE, * FACTOR,NPRINT,INFO,NFEV,NJEV,WA(6*N+1),LR,WA(N+1), * WA(2*N+1),WA(3*N+1),WA(4*N+1),WA(5*N+1)) IF (INFO .EQ. 5) INFO = 4 20 CONTINUE RETURN C C LAST CARD OF SUBROUTINE HYBRJ1. C END SUBROUTINE LMDER(FCN,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,GTOL, * MAXFEV,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,NJEV, * IPVT,QTF,WA1,WA2,WA3,WA4) INTEGER M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV INTEGER IPVT(N) REAL FTOL,XTOL,GTOL,FACTOR REAL X(N),FVEC(M),FJAC(LDFJAC,N),DIAG(N),QTF(N),WA1(N),WA2(N), * WA3(N),WA4(M) C ********** C C SUBROUTINE LMDER C C THE PURPOSE OF LMDER IS TO MINIMIZE THE SUM OF THE SQUARES OF C M NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION OF C THE LEVENBERG-MARQUARDT ALGORITHM. THE USER MUST PROVIDE A C SUBROUTINE WHICH CALCULATES THE FUNCTIONS AND THE JACOBIAN. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE LMDER(FCN,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,GTOL, C MAXFEV,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV, C NJEV,IPVT,QTF,WA1,WA2,WA3,WA4) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS AND THE JACOBIAN. FCN MUST C BE DECLARED IN AN EXTERNAL STATEMENT IN THE USER C CALLING PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(M,N,X,FVEC,FJAC,LDFJAC,IFLAG) C INTEGER M,N,LDFJAC,IFLAG C REAL X(N),FVEC(M),FJAC(LDFJAC,N) C ---------- C IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC. C IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND C RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC. C ---------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF LMDER. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VARIABLES. N MUST NOT EXCEED M. C C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. C C FVEC IS AN OUTPUT ARRAY OF LENGTH M WHICH CONTAINS C THE FUNCTIONS EVALUATED AT THE OUTPUT X. C C FJAC IS AN OUTPUT M BY N ARRAY. THE UPPER N BY N SUBMATRIX C OF FJAC CONTAINS AN UPPER TRIANGULAR MATRIX R WITH C DIAGONAL ELEMENTS OF NONINCREASING MAGNITUDE SUCH THAT C C T T T C P *(JAC *JAC)*P = R *R, C C WHERE P IS A PERMUTATION MATRIX AND JAC IS THE FINAL C CALCULATED JACOBIAN. COLUMN J OF P IS COLUMN IPVT(J) C (SEE BELOW) OF THE IDENTITY MATRIX. THE LOWER TRAPEZOIDAL C PART OF FJAC CONTAINS INFORMATION GENERATED DURING C THE COMPUTATION OF R. C C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. C C FTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION C OCCURS WHEN BOTH THE ACTUAL AND PREDICTED RELATIVE C REDUCTIONS IN THE SUM OF SQUARES ARE AT MOST FTOL. C THEREFORE, FTOL MEASURES THE RELATIVE ERROR DESIRED C IN THE SUM OF SQUARES. C C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION C OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE C ITERATES IS AT MOST XTOL. THEREFORE, XTOL MEASURES THE C RELATIVE ERROR DESIRED IN THE APPROXIMATE SOLUTION. C C GTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION C OCCURS WHEN THE COSINE OF THE ANGLE BETWEEN FVEC AND C ANY COLUMN OF THE JACOBIAN IS AT MOST GTOL IN ABSOLUTE C VALUE. THEREFORE, GTOL MEASURES THE ORTHOGONALITY C DESIRED BETWEEN THE FUNCTION VECTOR AND THE COLUMNS C OF THE JACOBIAN. C C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION C OCCURS WHEN THE NUMBER OF CALLS TO FCN WITH IFLAG = 1 C HAS REACHED MAXFEV. C C DIAG IS AN ARRAY OF LENGTH N. IF MODE = 1 (SEE C BELOW), DIAG IS INTERNALLY SET. IF MODE = 2, DIAG C MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS C MULTIPLICATIVE SCALE FACTORS FOR THE VARIABLES. C C MODE IS AN INTEGER INPUT VARIABLE. IF MODE = 1, THE C VARIABLES WILL BE SCALED INTERNALLY. IF MODE = 2, C THE SCALING IS SPECIFIED BY THE INPUT DIAG. OTHER C VALUES OF MODE ARE EQUIVALENT TO MODE = 1. C C FACTOR IS A POSITIVE INPUT VARIABLE USED IN DETERMINING THE C INITIAL STEP BOUND. THIS BOUND IS SET TO THE PRODUCT OF C FACTOR AND THE EUCLIDEAN NORM OF DIAG*X IF NONZERO, OR ELSE C TO FACTOR ITSELF. IN MOST CASES FACTOR SHOULD LIE IN THE C INTERVAL (.1,100.). 100. IS A GENERALLY RECOMMENDED VALUE. C C NPRINT IS AN INTEGER INPUT VARIABLE THAT ENABLES CONTROLLED C PRINTING OF ITERATES IF IT IS POSITIVE. IN THIS CASE, C FCN IS CALLED WITH IFLAG = 0 AT THE BEGINNING OF THE FIRST C ITERATION AND EVERY NPRINT ITERATIONS THEREAFTER AND C IMMEDIATELY PRIOR TO RETURN, WITH X, FVEC, AND FJAC C AVAILABLE FOR PRINTING. FVEC AND FJAC SHOULD NOT BE C ALTERED. IF NPRINT IS NOT POSITIVE, NO SPECIAL CALLS C OF FCN WITH IFLAG = 0 ARE MADE. C C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE) C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE, C INFO IS SET AS FOLLOWS. C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO = 1 BOTH ACTUAL AND PREDICTED RELATIVE REDUCTIONS C IN THE SUM OF SQUARES ARE AT MOST FTOL. C C INFO = 2 RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES C IS AT MOST XTOL. C C INFO = 3 CONDITIONS FOR INFO = 1 AND INFO = 2 BOTH HOLD. C C INFO = 4 THE COSINE OF THE ANGLE BETWEEN FVEC AND ANY C COLUMN OF THE JACOBIAN IS AT MOST GTOL IN C ABSOLUTE VALUE. C C INFO = 5 NUMBER OF CALLS TO FCN WITH IFLAG = 1 HAS C REACHED MAXFEV. C C INFO = 6 FTOL IS TOO SMALL. NO FURTHER REDUCTION IN C THE SUM OF SQUARES IS POSSIBLE. C C INFO = 7 XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN C THE APPROXIMATE SOLUTION X IS POSSIBLE. C C INFO = 8 GTOL IS TOO SMALL. FVEC IS ORTHOGONAL TO THE C COLUMNS OF THE JACOBIAN TO MACHINE PRECISION. C C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF C CALLS TO FCN WITH IFLAG = 1. C C NJEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF C CALLS TO FCN WITH IFLAG = 2. C C IPVT IS AN INTEGER OUTPUT ARRAY OF LENGTH N. IPVT C DEFINES A PERMUTATION MATRIX P SUCH THAT JAC*P = Q*R, C WHERE JAC IS THE FINAL CALCULATED JACOBIAN, Q IS C ORTHOGONAL (NOT STORED), AND R IS UPPER TRIANGULAR C WITH DIAGONAL ELEMENTS OF NONINCREASING MAGNITUDE. C COLUMN J OF P IS COLUMN IPVT(J) OF THE IDENTITY MATRIX. C C QTF IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS C THE FIRST N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*FVEC. C C WA1, WA2, AND WA3 ARE WORK ARRAYS OF LENGTH N. C C WA4 IS A WORK ARRAY OF LENGTH M. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ...... FCN C C MINPACK-SUPPLIED ... SPMPAR,ENORM,LMPAR,QRFAC C C FORTRAN-SUPPLIED ... ABS,AMAX1,AMIN1,SQRT,MOD C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,IFLAG,ITER,J,L REAL ACTRED,DELTA,DIRDER,EPSMCH,FNORM,FNORM1,GNORM,ONE,PAR, * PNORM,PRERED,P1,P5,P25,P75,P0001,RATIO,SUM,TEMP,TEMP1, * TEMP2,XNORM,ZERO REAL SPMPAR,ENORM DATA ONE,P1,P5,P25,P75,P0001,ZERO * /1.0E0,1.0E-1,5.0E-1,2.5E-1,7.5E-1,1.0E-4,0.0E0/ C C EPSMCH IS THE MACHINE PRECISION. C EPSMCH = SPMPAR(1) C INFO = 0 IFLAG = 0 NFEV = 0 NJEV = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. M .LT. N .OR. LDFJAC .LT. M * .OR. FTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. GTOL .LT. ZERO * .OR. MAXFEV .LE. 0 .OR. FACTOR .LE. ZERO) GO TO 300 IF (MODE .NE. 2) GO TO 20 DO 10 J = 1, N IF (DIAG(J) .LE. ZERO) GO TO 300 10 CONTINUE 20 CONTINUE C C EVALUATE THE FUNCTION AT THE STARTING POINT C AND CALCULATE ITS NORM. C IFLAG = 1 CALL FCN(M,N,X,FVEC,FJAC,LDFJAC,IFLAG) NFEV = 1 IF (IFLAG .LT. 0) GO TO 300 FNORM = ENORM(M,FVEC) C C INITIALIZE LEVENBERG-MARQUARDT PARAMETER AND ITERATION COUNTER. C PAR = ZERO ITER = 1 C C BEGINNING OF THE OUTER LOOP. C 30 CONTINUE C C CALCULATE THE JACOBIAN MATRIX. C IFLAG = 2 CALL FCN(M,N,X,FVEC,FJAC,LDFJAC,IFLAG) NJEV = NJEV + 1 IF (IFLAG .LT. 0) GO TO 300 C C IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES. C IF (NPRINT .LE. 0) GO TO 40 IFLAG = 0 IF (MOD(ITER-1,NPRINT) .EQ. 0) * CALL FCN(M,N,X,FVEC,FJAC,LDFJAC,IFLAG) IF (IFLAG .LT. 0) GO TO 300 40 CONTINUE C C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN. C CALL QRFAC(M,N,FJAC,LDFJAC,.TRUE.,IPVT,N,WA1,WA2,WA3) C C ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING C TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN. C IF (ITER .NE. 1) GO TO 80 IF (MODE .EQ. 2) GO TO 60 DO 50 J = 1, N DIAG(J) = WA2(J) IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE 50 CONTINUE 60 CONTINUE C C ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X C AND INITIALIZE THE STEP BOUND DELTA. C DO 70 J = 1, N WA3(J) = DIAG(J)*X(J) 70 CONTINUE XNORM = ENORM(N,WA3) DELTA = FACTOR*XNORM IF (DELTA .EQ. ZERO) DELTA = FACTOR 80 CONTINUE C C FORM (Q TRANSPOSE)*FVEC AND STORE THE FIRST N COMPONENTS IN C QTF. C DO 90 I = 1, M WA4(I) = FVEC(I) 90 CONTINUE DO 130 J = 1, N IF (FJAC(J,J) .EQ. ZERO) GO TO 120 SUM = ZERO DO 100 I = J, M SUM = SUM + FJAC(I,J)*WA4(I) 100 CONTINUE TEMP = -SUM/FJAC(J,J) DO 110 I = J, M WA4(I) = WA4(I) + FJAC(I,J)*TEMP 110 CONTINUE 120 CONTINUE FJAC(J,J) = WA1(J) QTF(J) = WA4(J) 130 CONTINUE C C COMPUTE THE NORM OF THE SCALED GRADIENT. C GNORM = ZERO IF (FNORM .EQ. ZERO) GO TO 170 DO 160 J = 1, N L = IPVT(J) IF (WA2(L) .EQ. ZERO) GO TO 150 SUM = ZERO DO 140 I = 1, J SUM = SUM + FJAC(I,J)*(QTF(I)/FNORM) 140 CONTINUE GNORM = AMAX1(GNORM,ABS(SUM/WA2(L))) 150 CONTINUE 160 CONTINUE 170 CONTINUE C C TEST FOR CONVERGENCE OF THE GRADIENT NORM. C IF (GNORM .LE. GTOL) INFO = 4 IF (INFO .NE. 0) GO TO 300 C C RESCALE IF NECESSARY. C IF (MODE .EQ. 2) GO TO 190 DO 180 J = 1, N DIAG(J) = AMAX1(DIAG(J),WA2(J)) 180 CONTINUE 190 CONTINUE C C BEGINNING OF THE INNER LOOP. C 200 CONTINUE C C DETERMINE THE LEVENBERG-MARQUARDT PARAMETER. C CALL LMPAR(N,FJAC,LDFJAC,IPVT,DIAG,QTF,DELTA,PAR,WA1,WA2, * WA3,WA4) C C STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P. C DO 210 J = 1, N WA1(J) = -WA1(J) WA2(J) = X(J) + WA1(J) WA3(J) = DIAG(J)*WA1(J) 210 CONTINUE PNORM = ENORM(N,WA3) C C ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND. C IF (ITER .EQ. 1) DELTA = AMIN1(DELTA,PNORM) C C EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM. C IFLAG = 1 CALL FCN(M,N,WA2,WA4,FJAC,LDFJAC,IFLAG) NFEV = NFEV + 1 IF (IFLAG .LT. 0) GO TO 300 FNORM1 = ENORM(M,WA4) C C COMPUTE THE SCALED ACTUAL REDUCTION. C ACTRED = -ONE IF (P1*FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2 C C COMPUTE THE SCALED PREDICTED REDUCTION AND C THE SCALED DIRECTIONAL DERIVATIVE. C DO 230 J = 1, N WA3(J) = ZERO L = IPVT(J) TEMP = WA1(L) DO 220 I = 1, J WA3(I) = WA3(I) + FJAC(I,J)*TEMP 220 CONTINUE 230 CONTINUE TEMP1 = ENORM(N,WA3)/FNORM TEMP2 = (SQRT(PAR)*PNORM)/FNORM PRERED = TEMP1**2 + TEMP2**2/P5 DIRDER = -(TEMP1**2 + TEMP2**2) C C COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED C REDUCTION. C RATIO = ZERO IF (PRERED .NE. ZERO) RATIO = ACTRED/PRERED C C UPDATE THE STEP BOUND. C IF (RATIO .GT. P25) GO TO 240 IF (ACTRED .GE. ZERO) TEMP = P5 IF (ACTRED .LT. ZERO) * TEMP = P5*DIRDER/(DIRDER + P5*ACTRED) IF (P1*FNORM1 .GE. FNORM .OR. TEMP .LT. P1) TEMP = P1 DELTA = TEMP*AMIN1(DELTA,PNORM/P1) PAR = PAR/TEMP GO TO 260 240 CONTINUE IF (PAR .NE. ZERO .AND. RATIO .LT. P75) GO TO 250 DELTA = PNORM/P5 PAR = P5*PAR 250 CONTINUE 260 CONTINUE C C TEST FOR SUCCESSFUL ITERATION. C IF (RATIO .LT. P0001) GO TO 290 C C SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS. C DO 270 J = 1, N X(J) = WA2(J) WA2(J) = DIAG(J)*X(J) 270 CONTINUE DO 280 I = 1, M FVEC(I) = WA4(I) 280 CONTINUE XNORM = ENORM(N,WA2) FNORM = FNORM1 ITER = ITER + 1 290 CONTINUE C C TESTS FOR CONVERGENCE. C IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL * .AND. P5*RATIO .LE. ONE) INFO = 1 IF (DELTA .LE. XTOL*XNORM) INFO = 2 IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL * .AND. P5*RATIO .LE. ONE .AND. INFO .EQ. 2) INFO = 3 IF (INFO .NE. 0) GO TO 300 C C TESTS FOR TERMINATION AND STRINGENT TOLERANCES. C IF (NFEV .GE. MAXFEV) INFO = 5 IF (ABS(ACTRED) .LE. EPSMCH .AND. PRERED .LE. EPSMCH * .AND. P5*RATIO .LE. ONE) INFO = 6 IF (DELTA .LE. EPSMCH*XNORM) INFO = 7 IF (GNORM .LE. EPSMCH) INFO = 8 IF (INFO .NE. 0) GO TO 300 C C END OF THE INNER LOOP. REPEAT IF ITERATION UNSUCCESSFUL. C IF (RATIO .LT. P0001) GO TO 200 C C END OF THE OUTER LOOP. C GO TO 30 300 CONTINUE C C TERMINATION, EITHER NORMAL OR USER IMPOSED. C IF (IFLAG .LT. 0) INFO = IFLAG IFLAG = 0 IF (NPRINT .GT. 0) CALL FCN(M,N,X,FVEC,FJAC,LDFJAC,IFLAG) RETURN C C LAST CARD OF SUBROUTINE LMDER. C END SUBROUTINE LMDER1(FCN,M,N,X,FVEC,FJAC,LDFJAC,TOL,INFO,IPVT,WA, * LWA) INTEGER M,N,LDFJAC,INFO,LWA INTEGER IPVT(N) REAL TOL REAL X(N),FVEC(M),FJAC(LDFJAC,N),WA(LWA) EXTERNAL FCN C ********** C C SUBROUTINE LMDER1 C C THE PURPOSE OF LMDER1 IS TO MINIMIZE THE SUM OF THE SQUARES OF C M NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION OF THE C LEVENBERG-MARQUARDT ALGORITHM. THIS IS DONE BY USING THE MORE C GENERAL LEAST-SQUARES SOLVER LMDER. THE USER MUST PROVIDE A C SUBROUTINE WHICH CALCULATES THE FUNCTIONS AND THE JACOBIAN. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE LMDER1(FCN,M,N,X,FVEC,FJAC,LDFJAC,TOL,INFO, C IPVT,WA,LWA) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS AND THE JACOBIAN. FCN MUST C BE DECLARED IN AN EXTERNAL STATEMENT IN THE USER C CALLING PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(M,N,X,FVEC,FJAC,LDFJAC,IFLAG) C INTEGER M,N,LDFJAC,IFLAG C REAL X(N),FVEC(M),FJAC(LDFJAC,N) C ---------- C IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC. C IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND C RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC. C ---------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF LMDER1. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VARIABLES. N MUST NOT EXCEED M. C C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. C C FVEC IS AN OUTPUT ARRAY OF LENGTH M WHICH CONTAINS C THE FUNCTIONS EVALUATED AT THE OUTPUT X. C C FJAC IS AN OUTPUT M BY N ARRAY. THE UPPER N BY N SUBMATRIX C OF FJAC CONTAINS AN UPPER TRIANGULAR MATRIX R WITH C DIAGONAL ELEMENTS OF NONINCREASING MAGNITUDE SUCH THAT C C T T T C P *(JAC *JAC)*P = R *R, C C WHERE P IS A PERMUTATION MATRIX AND JAC IS THE FINAL C CALCULATED JACOBIAN. COLUMN J OF P IS COLUMN IPVT(J) C (SEE BELOW) OF THE IDENTITY MATRIX. THE LOWER TRAPEZOIDAL C PART OF FJAC CONTAINS INFORMATION GENERATED DURING C THE COMPUTATION OF R. C C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. C C TOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS C WHEN THE ALGORITHM ESTIMATES EITHER THAT THE RELATIVE C ERROR IN THE SUM OF SQUARES IS AT MOST TOL OR THAT C THE RELATIVE ERROR BETWEEN X AND THE SOLUTION IS AT C MOST TOL. C C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE) C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE, C INFO IS SET AS FOLLOWS. C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO = 1 ALGORITHM ESTIMATES THAT THE RELATIVE ERROR C IN THE SUM OF SQUARES IS AT MOST TOL. C C INFO = 2 ALGORITHM ESTIMATES THAT THE RELATIVE ERROR C BETWEEN X AND THE SOLUTION IS AT MOST TOL. C C INFO = 3 CONDITIONS FOR INFO = 1 AND INFO = 2 BOTH HOLD. C C INFO = 4 FVEC IS ORTHOGONAL TO THE COLUMNS OF THE C JACOBIAN TO MACHINE PRECISION. C C INFO = 5 NUMBER OF CALLS TO FCN WITH IFLAG = 1 HAS C REACHED 100*(N+1). C C INFO = 6 TOL IS TOO SMALL. NO FURTHER REDUCTION IN C THE SUM OF SQUARES IS POSSIBLE. C C INFO = 7 TOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN C THE APPROXIMATE SOLUTION X IS POSSIBLE. C C IPVT IS AN INTEGER OUTPUT ARRAY OF LENGTH N. IPVT C DEFINES A PERMUTATION MATRIX P SUCH THAT JAC*P = Q*R, C WHERE JAC IS THE FINAL CALCULATED JACOBIAN, Q IS C ORTHOGONAL (NOT STORED), AND R IS UPPER TRIANGULAR C WITH DIAGONAL ELEMENTS OF NONINCREASING MAGNITUDE. C COLUMN J OF P IS COLUMN IPVT(J) OF THE IDENTITY MATRIX. C C WA IS A WORK ARRAY OF LENGTH LWA. C C LWA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN 5*N+M. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ...... FCN C C MINPACK-SUPPLIED ... LMDER C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER MAXFEV,MODE,NFEV,NJEV,NPRINT REAL FACTOR,FTOL,GTOL,XTOL,ZERO DATA FACTOR,ZERO /1.0E2,0.0E0/ INFO = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. M .LT. N .OR. LDFJAC .LT. M .OR. TOL .LT. ZERO * .OR. LWA .LT. 5*N + M) GO TO 10 C C CALL LMDER. C MAXFEV = 100*(N + 1) FTOL = TOL XTOL = TOL GTOL = ZERO MODE = 1 NPRINT = 0 CALL LMDER(FCN,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,GTOL,MAXFEV, * WA(1),MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,IPVT,WA(N+1), * WA(2*N+1),WA(3*N+1),WA(4*N+1),WA(5*N+1)) IF (INFO .EQ. 8) INFO = 4 10 CONTINUE RETURN C C LAST CARD OF SUBROUTINE LMDER1. C END SUBROUTINE LMDIF(FCN,M,N,X,FVEC,FTOL,XTOL,GTOL,MAXFEV,EPSFCN, * DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,FJAC,LDFJAC, * IPVT,QTF,WA1,WA2,WA3,WA4) INTEGER M,N,MAXFEV,MODE,NPRINT,INFO,NFEV,LDFJAC INTEGER IPVT(N) REAL FTOL,XTOL,GTOL,EPSFCN,FACTOR REAL X(N),FVEC(M),DIAG(N),FJAC(LDFJAC,N),QTF(N),WA1(N),WA2(N), * WA3(N),WA4(M) EXTERNAL FCN C ********** C C SUBROUTINE LMDIF C C THE PURPOSE OF LMDIF IS TO MINIMIZE THE SUM OF THE SQUARES OF C M NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION OF C THE LEVENBERG-MARQUARDT ALGORITHM. THE USER MUST PROVIDE A C SUBROUTINE WHICH CALCULATES THE FUNCTIONS. THE JACOBIAN IS C THEN CALCULATED BY A FORWARD-DIFFERENCE APPROXIMATION. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE LMDIF(FCN,M,N,X,FVEC,FTOL,XTOL,GTOL,MAXFEV,EPSFCN, C DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,FJAC, C LDFJAC,IPVT,QTF,WA1,WA2,WA3,WA4) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED C IN AN EXTERNAL STATEMENT IN THE USER CALLING C PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(M,N,X,FVEC,IFLAG) C INTEGER M,N,IFLAG C REAL X(N),FVEC(M) C ---------- C CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. C ---------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF LMDIF. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VARIABLES. N MUST NOT EXCEED M. C C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. C C FVEC IS AN OUTPUT ARRAY OF LENGTH M WHICH CONTAINS C THE FUNCTIONS EVALUATED AT THE OUTPUT X. C C FTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION C OCCURS WHEN BOTH THE ACTUAL AND PREDICTED RELATIVE C REDUCTIONS IN THE SUM OF SQUARES ARE AT MOST FTOL. C THEREFORE, FTOL MEASURES THE RELATIVE ERROR DESIRED C IN THE SUM OF SQUARES. C C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION C OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE C ITERATES IS AT MOST XTOL. THEREFORE, XTOL MEASURES THE C RELATIVE ERROR DESIRED IN THE APPROXIMATE SOLUTION. C C GTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION C OCCURS WHEN THE COSINE OF THE ANGLE BETWEEN FVEC AND C ANY COLUMN OF THE JACOBIAN IS AT MOST GTOL IN ABSOLUTE C VALUE. THEREFORE, GTOL MEASURES THE ORTHOGONALITY C DESIRED BETWEEN THE FUNCTION VECTOR AND THE COLUMNS C OF THE JACOBIAN. C C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION C OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST C MAXFEV BY THE END OF AN ITERATION. C C EPSFCN IS AN INPUT VARIABLE USED IN DETERMINING A SUITABLE C STEP LENGTH FOR THE FORWARD-DIFFERENCE APPROXIMATION. THIS C APPROXIMATION ASSUMES THAT THE RELATIVE ERRORS IN THE C FUNCTIONS ARE OF THE ORDER OF EPSFCN. IF EPSFCN IS LESS C THAN THE MACHINE PRECISION, IT IS ASSUMED THAT THE RELATIVE C ERRORS IN THE FUNCTIONS ARE OF THE ORDER OF THE MACHINE C PRECISION. C C DIAG IS AN ARRAY OF LENGTH N. IF MODE = 1 (SEE C BELOW), DIAG IS INTERNALLY SET. IF MODE = 2, DIAG C MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS C MULTIPLICATIVE SCALE FACTORS FOR THE VARIABLES. C C MODE IS AN INTEGER INPUT VARIABLE. IF MODE = 1, THE C VARIABLES WILL BE SCALED INTERNALLY. IF MODE = 2, C THE SCALING IS SPECIFIED BY THE INPUT DIAG. OTHER C VALUES OF MODE ARE EQUIVALENT TO MODE = 1. C C FACTOR IS A POSITIVE INPUT VARIABLE USED IN DETERMINING THE C INITIAL STEP BOUND. THIS BOUND IS SET TO THE PRODUCT OF C FACTOR AND THE EUCLIDEAN NORM OF DIAG*X IF NONZERO, OR ELSE C TO FACTOR ITSELF. IN MOST CASES FACTOR SHOULD LIE IN THE C INTERVAL (.1,100.). 100. IS A GENERALLY RECOMMENDED VALUE. C C NPRINT IS AN INTEGER INPUT VARIABLE THAT ENABLES CONTROLLED C PRINTING OF ITERATES IF IT IS POSITIVE. IN THIS CASE, C FCN IS CALLED WITH IFLAG = 0 AT THE BEGINNING OF THE FIRST C ITERATION AND EVERY NPRINT ITERATIONS THEREAFTER AND C IMMEDIATELY PRIOR TO RETURN, WITH X AND FVEC AVAILABLE C FOR PRINTING. IF NPRINT IS NOT POSITIVE, NO SPECIAL CALLS C OF FCN WITH IFLAG = 0 ARE MADE. C C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE) C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE, C INFO IS SET AS FOLLOWS. C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO = 1 BOTH ACTUAL AND PREDICTED RELATIVE REDUCTIONS C IN THE SUM OF SQUARES ARE AT MOST FTOL. C C INFO = 2 RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES C IS AT MOST XTOL. C C INFO = 3 CONDITIONS FOR INFO = 1 AND INFO = 2 BOTH HOLD. C C INFO = 4 THE COSINE OF THE ANGLE BETWEEN FVEC AND ANY C COLUMN OF THE JACOBIAN IS AT MOST GTOL IN C ABSOLUTE VALUE. C C INFO = 5 NUMBER OF CALLS TO FCN HAS REACHED OR C EXCEEDED MAXFEV. C C INFO = 6 FTOL IS TOO SMALL. NO FURTHER REDUCTION IN C THE SUM OF SQUARES IS POSSIBLE. C C INFO = 7 XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN C THE APPROXIMATE SOLUTION X IS POSSIBLE. C C INFO = 8 GTOL IS TOO SMALL. FVEC IS ORTHOGONAL TO THE C COLUMNS OF THE JACOBIAN TO MACHINE PRECISION. C C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF C CALLS TO FCN. C C FJAC IS AN OUTPUT M BY N ARRAY. THE UPPER N BY N SUBMATRIX C OF FJAC CONTAINS AN UPPER TRIANGULAR MATRIX R WITH C DIAGONAL ELEMENTS OF NONINCREASING MAGNITUDE SUCH THAT C C T T T C P *(JAC *JAC)*P = R *R, C C WHERE P IS A PERMUTATION MATRIX AND JAC IS THE FINAL C CALCULATED JACOBIAN. COLUMN J OF P IS COLUMN IPVT(J) C (SEE BELOW) OF THE IDENTITY MATRIX. THE LOWER TRAPEZOIDAL C PART OF FJAC CONTAINS INFORMATION GENERATED DURING C THE COMPUTATION OF R. C C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. C C IPVT IS AN INTEGER OUTPUT ARRAY OF LENGTH N. IPVT C DEFINES A PERMUTATION MATRIX P SUCH THAT JAC*P = Q*R, C WHERE JAC IS THE FINAL CALCULATED JACOBIAN, Q IS C ORTHOGONAL (NOT STORED), AND R IS UPPER TRIANGULAR C WITH DIAGONAL ELEMENTS OF NONINCREASING MAGNITUDE. C COLUMN J OF P IS COLUMN IPVT(J) OF THE IDENTITY MATRIX. C C QTF IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS C THE FIRST N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*FVEC. C C WA1, WA2, AND WA3 ARE WORK ARRAYS OF LENGTH N. C C WA4 IS A WORK ARRAY OF LENGTH M. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ...... FCN C C MINPACK-SUPPLIED ... SPMPAR,ENORM,FDJAC2,LMPAR,QRFAC C C FORTRAN-SUPPLIED ... ABS,AMAX1,AMIN1,SQRT,MOD C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,IFLAG,ITER,J,L REAL ACTRED,DELTA,DIRDER,EPSMCH,FNORM,FNORM1,GNORM,ONE,PAR, * PNORM,PRERED,P1,P5,P25,P75,P0001,RATIO,SUM,TEMP,TEMP1, * TEMP2,XNORM,ZERO REAL SPMPAR,ENORM DATA ONE,P1,P5,P25,P75,P0001,ZERO * /1.0E0,1.0E-1,5.0E-1,2.5E-1,7.5E-1,1.0E-4,0.0E0/ C C EPSMCH IS THE MACHINE PRECISION. C EPSMCH = SPMPAR(1) C INFO = 0 IFLAG = 0 NFEV = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. M .LT. N .OR. LDFJAC .LT. M * .OR. FTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. GTOL .LT. ZERO * .OR. MAXFEV .LE. 0 .OR. FACTOR .LE. ZERO) GO TO 300 IF (MODE .NE. 2) GO TO 20 DO 10 J = 1, N IF (DIAG(J) .LE. ZERO) GO TO 300 10 CONTINUE 20 CONTINUE C C EVALUATE THE FUNCTION AT THE STARTING POINT C AND CALCULATE ITS NORM. C IFLAG = 1 CALL FCN(M,N,X,FVEC,IFLAG) NFEV = 1 IF (IFLAG .LT. 0) GO TO 300 FNORM = ENORM(M,FVEC) C C INITIALIZE LEVENBERG-MARQUARDT PARAMETER AND ITERATION COUNTER. C PAR = ZERO ITER = 1 C C BEGINNING OF THE OUTER LOOP. C 30 CONTINUE C C CALCULATE THE JACOBIAN MATRIX. C IFLAG = 2 CALL FDJAC2(FCN,M,N,X,FVEC,FJAC,LDFJAC,IFLAG,EPSFCN,WA4) NFEV = NFEV + N IF (IFLAG .LT. 0) GO TO 300 C C IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES. C IF (NPRINT .LE. 0) GO TO 40 IFLAG = 0 IF (MOD(ITER-1,NPRINT) .EQ. 0) CALL FCN(M,N,X,FVEC,IFLAG) IF (IFLAG .LT. 0) GO TO 300 40 CONTINUE C C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN. C CALL QRFAC(M,N,FJAC,LDFJAC,.TRUE.,IPVT,N,WA1,WA2,WA3) C C ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING C TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN. C IF (ITER .NE. 1) GO TO 80 IF (MODE .EQ. 2) GO TO 60 DO 50 J = 1, N DIAG(J) = WA2(J) IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE 50 CONTINUE 60 CONTINUE C C ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X C AND INITIALIZE THE STEP BOUND DELTA. C DO 70 J = 1, N WA3(J) = DIAG(J)*X(J) 70 CONTINUE XNORM = ENORM(N,WA3) DELTA = FACTOR*XNORM IF (DELTA .EQ. ZERO) DELTA = FACTOR 80 CONTINUE C C FORM (Q TRANSPOSE)*FVEC AND STORE THE FIRST N COMPONENTS IN C QTF. C DO 90 I = 1, M WA4(I) = FVEC(I) 90 CONTINUE DO 130 J = 1, N IF (FJAC(J,J) .EQ. ZERO) GO TO 120 SUM = ZERO DO 100 I = J, M SUM = SUM + FJAC(I,J)*WA4(I) 100 CONTINUE TEMP = -SUM/FJAC(J,J) DO 110 I = J, M WA4(I) = WA4(I) + FJAC(I,J)*TEMP 110 CONTINUE 120 CONTINUE FJAC(J,J) = WA1(J) QTF(J) = WA4(J) 130 CONTINUE C C COMPUTE THE NORM OF THE SCALED GRADIENT. C GNORM = ZERO IF (FNORM .EQ. ZERO) GO TO 170 DO 160 J = 1, N L = IPVT(J) IF (WA2(L) .EQ. ZERO) GO TO 150 SUM = ZERO DO 140 I = 1, J SUM = SUM + FJAC(I,J)*(QTF(I)/FNORM) 140 CONTINUE GNORM = AMAX1(GNORM,ABS(SUM/WA2(L))) 150 CONTINUE 160 CONTINUE 170 CONTINUE C C TEST FOR CONVERGENCE OF THE GRADIENT NORM. C IF (GNORM .LE. GTOL) INFO = 4 IF (INFO .NE. 0) GO TO 300 C C RESCALE IF NECESSARY. C IF (MODE .EQ. 2) GO TO 190 DO 180 J = 1, N DIAG(J) = AMAX1(DIAG(J),WA2(J)) 180 CONTINUE 190 CONTINUE C C BEGINNING OF THE INNER LOOP. C 200 CONTINUE C C DETERMINE THE LEVENBERG-MARQUARDT PARAMETER. C CALL LMPAR(N,FJAC,LDFJAC,IPVT,DIAG,QTF,DELTA,PAR,WA1,WA2, * WA3,WA4) C C STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P. C DO 210 J = 1, N WA1(J) = -WA1(J) WA2(J) = X(J) + WA1(J) WA3(J) = DIAG(J)*WA1(J) 210 CONTINUE PNORM = ENORM(N,WA3) C C ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND. C IF (ITER .EQ. 1) DELTA = AMIN1(DELTA,PNORM) C C EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM. C IFLAG = 1 CALL FCN(M,N,WA2,WA4,IFLAG) NFEV = NFEV + 1 IF (IFLAG .LT. 0) GO TO 300 FNORM1 = ENORM(M,WA4) C C COMPUTE THE SCALED ACTUAL REDUCTION. C ACTRED = -ONE IF (P1*FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2 C C COMPUTE THE SCALED PREDICTED REDUCTION AND C THE SCALED DIRECTIONAL DERIVATIVE. C DO 230 J = 1, N WA3(J) = ZERO L = IPVT(J) TEMP = WA1(L) DO 220 I = 1, J WA3(I) = WA3(I) + FJAC(I,J)*TEMP 220 CONTINUE 230 CONTINUE TEMP1 = ENORM(N,WA3)/FNORM TEMP2 = (SQRT(PAR)*PNORM)/FNORM PRERED = TEMP1**2 + TEMP2**2/P5 DIRDER = -(TEMP1**2 + TEMP2**2) C C COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED C REDUCTION. C RATIO = ZERO IF (PRERED .NE. ZERO) RATIO = ACTRED/PRERED C C UPDATE THE STEP BOUND. C IF (RATIO .GT. P25) GO TO 240 IF (ACTRED .GE. ZERO) TEMP = P5 IF (ACTRED .LT. ZERO) * TEMP = P5*DIRDER/(DIRDER + P5*ACTRED) IF (P1*FNORM1 .GE. FNORM .OR. TEMP .LT. P1) TEMP = P1 DELTA = TEMP*AMIN1(DELTA,PNORM/P1) PAR = PAR/TEMP GO TO 260 240 CONTINUE IF (PAR .NE. ZERO .AND. RATIO .LT. P75) GO TO 250 DELTA = PNORM/P5 PAR = P5*PAR 250 CONTINUE 260 CONTINUE C C TEST FOR SUCCESSFUL ITERATION. C IF (RATIO .LT. P0001) GO TO 290 C C SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS. C DO 270 J = 1, N X(J) = WA2(J) WA2(J) = DIAG(J)*X(J) 270 CONTINUE DO 280 I = 1, M FVEC(I) = WA4(I) 280 CONTINUE XNORM = ENORM(N,WA2) FNORM = FNORM1 ITER = ITER + 1 290 CONTINUE C C TESTS FOR CONVERGENCE. C IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL * .AND. P5*RATIO .LE. ONE) INFO = 1 IF (DELTA .LE. XTOL*XNORM) INFO = 2 IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL * .AND. P5*RATIO .LE. ONE .AND. INFO .EQ. 2) INFO = 3 IF (INFO .NE. 0) GO TO 300 C C TESTS FOR TERMINATION AND STRINGENT TOLERANCES. C IF (NFEV .GE. MAXFEV) INFO = 5 IF (ABS(ACTRED) .LE. EPSMCH .AND. PRERED .LE. EPSMCH * .AND. P5*RATIO .LE. ONE) INFO = 6 IF (DELTA .LE. EPSMCH*XNORM) INFO = 7 IF (GNORM .LE. EPSMCH) INFO = 8 IF (INFO .NE. 0) GO TO 300 C C END OF THE INNER LOOP. REPEAT IF ITERATION UNSUCCESSFUL. C IF (RATIO .LT. P0001) GO TO 200 C C END OF THE OUTER LOOP. C GO TO 30 300 CONTINUE C C TERMINATION, EITHER NORMAL OR USER IMPOSED. C IF (IFLAG .LT. 0) INFO = IFLAG IFLAG = 0 IF (NPRINT .GT. 0) CALL FCN(M,N,X,FVEC,IFLAG) RETURN C C LAST CARD OF SUBROUTINE LMDIF. C END SUBROUTINE LMDIF1(FCN,M,N,X,FVEC,TOL,INFO,IWA,WA,LWA) INTEGER M,N,INFO,LWA INTEGER IWA(N) REAL TOL REAL X(N),FVEC(M),WA(LWA) EXTERNAL FCN C ********** C C SUBROUTINE LMDIF1 C C THE PURPOSE OF LMDIF1 IS TO MINIMIZE THE SUM OF THE SQUARES OF C M NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION OF THE C LEVENBERG-MARQUARDT ALGORITHM. THIS IS DONE BY USING THE MORE C GENERAL LEAST-SQUARES SOLVER LMDIF. THE USER MUST PROVIDE A C SUBROUTINE WHICH CALCULATES THE FUNCTIONS. THE JACOBIAN IS C THEN CALCULATED BY A FORWARD-DIFFERENCE APPROXIMATION. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE LMDIF1(FCN,M,N,X,FVEC,TOL,INFO,IWA,WA,LWA) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED C IN AN EXTERNAL STATEMENT IN THE USER CALLING C PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(M,N,X,FVEC,IFLAG) C INTEGER M,N,IFLAG C REAL X(N),FVEC(M) C ---------- C CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. C ---------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF LMDIF1. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VARIABLES. N MUST NOT EXCEED M. C C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. C C FVEC IS AN OUTPUT ARRAY OF LENGTH M WHICH CONTAINS C THE FUNCTIONS EVALUATED AT THE OUTPUT X. C C TOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS C WHEN THE ALGORITHM ESTIMATES EITHER THAT THE RELATIVE C ERROR IN THE SUM OF SQUARES IS AT MOST TOL OR THAT C THE RELATIVE ERROR BETWEEN X AND THE SOLUTION IS AT C MOST TOL. C C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE) C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE, C INFO IS SET AS FOLLOWS. C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO = 1 ALGORITHM ESTIMATES THAT THE RELATIVE ERROR C IN THE SUM OF SQUARES IS AT MOST TOL. C C INFO = 2 ALGORITHM ESTIMATES THAT THE RELATIVE ERROR C BETWEEN X AND THE SOLUTION IS AT MOST TOL. C C INFO = 3 CONDITIONS FOR INFO = 1 AND INFO = 2 BOTH HOLD. C C INFO = 4 FVEC IS ORTHOGONAL TO THE COLUMNS OF THE C JACOBIAN TO MACHINE PRECISION. C C INFO = 5 NUMBER OF CALLS TO FCN HAS REACHED OR C EXCEEDED 200*(N+1). C C INFO = 6 TOL IS TOO SMALL. NO FURTHER REDUCTION IN C THE SUM OF SQUARES IS POSSIBLE. C C INFO = 7 TOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN C THE APPROXIMATE SOLUTION X IS POSSIBLE. C C IWA IS AN INTEGER WORK ARRAY OF LENGTH N. C C WA IS A WORK ARRAY OF LENGTH LWA. C C LWA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN C M*N+5*N+M. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ...... FCN C C MINPACK-SUPPLIED ... LMDIF C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER MAXFEV,MODE,MP5N,NFEV,NPRINT REAL EPSFCN,FACTOR,FTOL,GTOL,XTOL,ZERO DATA FACTOR,ZERO /1.0E2,0.0E0/ INFO = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. M .LT. N .OR. TOL .LT. ZERO * .OR. LWA .LT. M*N + 5*N + M) GO TO 10 C C CALL LMDIF. C MAXFEV = 200*(N + 1) FTOL = TOL XTOL = TOL GTOL = ZERO EPSFCN = ZERO MODE = 1 NPRINT = 0 MP5N = M + 5*N CALL LMDIF(FCN,M,N,X,FVEC,FTOL,XTOL,GTOL,MAXFEV,EPSFCN,WA(1), * MODE,FACTOR,NPRINT,INFO,NFEV,WA(MP5N+1),M,IWA, * WA(N+1),WA(2*N+1),WA(3*N+1),WA(4*N+1),WA(5*N+1)) IF (INFO .EQ. 8) INFO = 4 10 CONTINUE RETURN C C LAST CARD OF SUBROUTINE LMDIF1. C END SUBROUTINE LMPAR(N,R,LDR,IPVT,DIAG,QTB,DELTA,PAR,X,SDIAG,WA1, * WA2) INTEGER N,LDR INTEGER IPVT(N) REAL DELTA,PAR REAL R(LDR,N),DIAG(N),QTB(N),X(N),SDIAG(N),WA1(N),WA2(N) C ********** C C SUBROUTINE LMPAR C C GIVEN AN M BY N MATRIX A, AN N BY N NONSINGULAR DIAGONAL C MATRIX D, AN M-VECTOR B, AND A POSITIVE NUMBER DELTA, C THE PROBLEM IS TO DETERMINE A VALUE FOR THE PARAMETER C PAR SUCH THAT IF X SOLVES THE SYSTEM C C A*X = B , SQRT(PAR)*D*X = 0 , C C IN THE LEAST SQUARES SENSE, AND DXNORM IS THE EUCLIDEAN C NORM OF D*X, THEN EITHER PAR IS ZERO AND C C (DXNORM-DELTA) .LE. 0.1*DELTA , C C OR PAR IS POSITIVE AND C C ABS(DXNORM-DELTA) .LE. 0.1*DELTA . C C THIS SUBROUTINE COMPLETES THE SOLUTION OF THE PROBLEM C IF IT IS PROVIDED WITH THE NECESSARY INFORMATION FROM THE C QR FACTORIZATION, WITH COLUMN PIVOTING, OF A. THAT IS, IF C A*P = Q*R, WHERE P IS A PERMUTATION MATRIX, Q HAS ORTHOGONAL C COLUMNS, AND R IS AN UPPER TRIANGULAR MATRIX WITH DIAGONAL C ELEMENTS OF NONINCREASING MAGNITUDE, THEN LMPAR EXPECTS C THE FULL UPPER TRIANGLE OF R, THE PERMUTATION MATRIX P, C AND THE FIRST N COMPONENTS OF (Q TRANSPOSE)*B. ON OUTPUT C LMPAR ALSO PROVIDES AN UPPER TRIANGULAR MATRIX S SUCH THAT C C T T T C P *(A *A + PAR*D*D)*P = S *S . C C S IS EMPLOYED WITHIN LMPAR AND MAY BE OF SEPARATE INTEREST. C C ONLY A FEW ITERATIONS ARE GENERALLY NEEDED FOR CONVERGENCE C OF THE ALGORITHM. IF, HOWEVER, THE LIMIT OF 10 ITERATIONS C IS REACHED, THEN THE OUTPUT PAR WILL CONTAIN THE BEST C VALUE OBTAINED SO FAR. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE LMPAR(N,R,LDR,IPVT,DIAG,QTB,DELTA,PAR,X,SDIAG, C WA1,WA2) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER OF R. C C R IS AN N BY N ARRAY. ON INPUT THE FULL UPPER TRIANGLE C MUST CONTAIN THE FULL UPPER TRIANGLE OF THE MATRIX R. C ON OUTPUT THE FULL UPPER TRIANGLE IS UNALTERED, AND THE C STRICT LOWER TRIANGLE CONTAINS THE STRICT UPPER TRIANGLE C (TRANSPOSED) OF THE UPPER TRIANGULAR MATRIX S. C C LDR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY R. C C IPVT IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH DEFINES THE C PERMUTATION MATRIX P SUCH THAT A*P = Q*R. COLUMN J OF P C IS COLUMN IPVT(J) OF THE IDENTITY MATRIX. C C DIAG IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE C DIAGONAL ELEMENTS OF THE MATRIX D. C C QTB IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE FIRST C N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*B. C C DELTA IS A POSITIVE INPUT VARIABLE WHICH SPECIFIES AN UPPER C BOUND ON THE EUCLIDEAN NORM OF D*X. C C PAR IS A NONNEGATIVE VARIABLE. ON INPUT PAR CONTAINS AN C INITIAL ESTIMATE OF THE LEVENBERG-MARQUARDT PARAMETER. C ON OUTPUT PAR CONTAINS THE FINAL ESTIMATE. C C X IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE LEAST C SQUARES SOLUTION OF THE SYSTEM A*X = B, SQRT(PAR)*D*X = 0, C FOR THE OUTPUT PAR. C C SDIAG IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE C DIAGONAL ELEMENTS OF THE UPPER TRIANGULAR MATRIX S. C C WA1 AND WA2 ARE WORK ARRAYS OF LENGTH N. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... SPMPAR,ENORM,QRSOLV C C FORTRAN-SUPPLIED ... ABS,AMAX1,AMIN1,SQRT C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,ITER,J,JM1,JP1,K,L,NSING REAL DXNORM,DWARF,FP,GNORM,PARC,PARL,PARU,P1,P001,SUM,TEMP,ZERO REAL SPMPAR,ENORM DATA P1,P001,ZERO /1.0E-1,1.0E-3,0.0E0/ C C DWARF IS THE SMALLEST POSITIVE MAGNITUDE. C DWARF = SPMPAR(2) C C COMPUTE AND STORE IN X THE GAUSS-NEWTON DIRECTION. IF THE C JACOBIAN IS RANK-DEFICIENT, OBTAIN A LEAST SQUARES SOLUTION. C NSING = N DO 10 J = 1, N WA1(J) = QTB(J) IF (R(J,J) .EQ. ZERO .AND. NSING .EQ. N) NSING = J - 1 IF (NSING .LT. N) WA1(J) = ZERO 10 CONTINUE IF (NSING .LT. 1) GO TO 50 DO 40 K = 1, NSING J = NSING - K + 1 WA1(J) = WA1(J)/R(J,J) TEMP = WA1(J) JM1 = J - 1 IF (JM1 .LT. 1) GO TO 30 DO 20 I = 1, JM1 WA1(I) = WA1(I) - R(I,J)*TEMP 20 CONTINUE 30 CONTINUE 40 CONTINUE 50 CONTINUE DO 60 J = 1, N L = IPVT(J) X(L) = WA1(J) 60 CONTINUE C C INITIALIZE THE ITERATION COUNTER. C EVALUATE THE FUNCTION AT THE ORIGIN, AND TEST C FOR ACCEPTANCE OF THE GAUSS-NEWTON DIRECTION. C ITER = 0 DO 70 J = 1, N WA2(J) = DIAG(J)*X(J) 70 CONTINUE DXNORM = ENORM(N,WA2) FP = DXNORM - DELTA IF (FP .LE. P1*DELTA) GO TO 220 C C IF THE JACOBIAN IS NOT RANK DEFICIENT, THE NEWTON C STEP PROVIDES A LOWER BOUND, PARL, FOR THE ZERO OF C THE FUNCTION. OTHERWISE SET THIS BOUND TO ZERO. C PARL = ZERO IF (NSING .LT. N) GO TO 120 DO 80 J = 1, N L = IPVT(J) WA1(J) = DIAG(L)*(WA2(L)/DXNORM) 80 CONTINUE DO 110 J = 1, N SUM = ZERO JM1 = J - 1 IF (JM1 .LT. 1) GO TO 100 DO 90 I = 1, JM1 SUM = SUM + R(I,J)*WA1(I) 90 CONTINUE 100 CONTINUE WA1(J) = (WA1(J) - SUM)/R(J,J) 110 CONTINUE TEMP = ENORM(N,WA1) PARL = ((FP/DELTA)/TEMP)/TEMP 120 CONTINUE C C CALCULATE AN UPPER BOUND, PARU, FOR THE ZERO OF THE FUNCTION. C DO 140 J = 1, N SUM = ZERO DO 130 I = 1, J SUM = SUM + R(I,J)*QTB(I) 130 CONTINUE L = IPVT(J) WA1(J) = SUM/DIAG(L) 140 CONTINUE GNORM = ENORM(N,WA1) PARU = GNORM/DELTA IF (PARU .EQ. ZERO) PARU = DWARF/AMIN1(DELTA,P1) C C IF THE INPUT PAR LIES OUTSIDE OF THE INTERVAL (PARL,PARU), C SET PAR TO THE CLOSER ENDPOINT. C PAR = AMAX1(PAR,PARL) PAR = AMIN1(PAR,PARU) IF (PAR .EQ. ZERO) PAR = GNORM/DXNORM C C BEGINNING OF AN ITERATION. C 150 CONTINUE ITER = ITER + 1 C C EVALUATE THE FUNCTION AT THE CURRENT VALUE OF PAR. C IF (PAR .EQ. ZERO) PAR = AMAX1(DWARF,P001*PARU) TEMP = SQRT(PAR) DO 160 J = 1, N WA1(J) = TEMP*DIAG(J) 160 CONTINUE CALL QRSOLV(N,R,LDR,IPVT,WA1,QTB,X,SDIAG,WA2) DO 170 J = 1, N WA2(J) = DIAG(J)*X(J) 170 CONTINUE DXNORM = ENORM(N,WA2) TEMP = FP FP = DXNORM - DELTA C C IF THE FUNCTION IS SMALL ENOUGH, ACCEPT THE CURRENT VALUE C OF PAR. ALSO TEST FOR THE EXCEPTIONAL CASES WHERE PARL C IS ZERO OR THE NUMBER OF ITERATIONS HAS REACHED 10. C IF (ABS(FP) .LE. P1*DELTA * .OR. PARL .EQ. ZERO .AND. FP .LE. TEMP * .AND. TEMP .LT. ZERO .OR. ITER .EQ. 10) GO TO 220 C C COMPUTE THE NEWTON CORRECTION. C DO 180 J = 1, N L = IPVT(J) WA1(J) = DIAG(L)*(WA2(L)/DXNORM) 180 CONTINUE DO 210 J = 1, N WA1(J) = WA1(J)/SDIAG(J) TEMP = WA1(J) JP1 = J + 1 IF (N .LT. JP1) GO TO 200 DO 190 I = JP1, N WA1(I) = WA1(I) - R(I,J)*TEMP 190 CONTINUE 200 CONTINUE 210 CONTINUE TEMP = ENORM(N,WA1) PARC = ((FP/DELTA)/TEMP)/TEMP C C DEPENDING ON THE SIGN OF THE FUNCTION, UPDATE PARL OR PARU. C IF (FP .GT. ZERO) PARL = AMAX1(PARL,PAR) IF (FP .LT. ZERO) PARU = AMIN1(PARU,PAR) C C COMPUTE AN IMPROVED ESTIMATE FOR PAR. C PAR = AMAX1(PARL,PAR+PARC) C C END OF AN ITERATION. C GO TO 150 220 CONTINUE C C TERMINATION. C IF (ITER .EQ. 0) PAR = ZERO RETURN C C LAST CARD OF SUBROUTINE LMPAR. C END SUBROUTINE LMSTR(FCN,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,GTOL, * MAXFEV,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,NJEV, * IPVT,QTF,WA1,WA2,WA3,WA4) INTEGER M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV INTEGER IPVT(N) LOGICAL SING REAL FTOL,XTOL,GTOL,FACTOR REAL X(N),FVEC(M),FJAC(LDFJAC,N),DIAG(N),QTF(N),WA1(N),WA2(N), * WA3(N),WA4(M) C ********** C C SUBROUTINE LMSTR C C THE PURPOSE OF LMSTR IS TO MINIMIZE THE SUM OF THE SQUARES OF C M NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION OF C THE LEVENBERG-MARQUARDT ALGORITHM WHICH USES MINIMAL STORAGE. C THE USER MUST PROVIDE A SUBROUTINE WHICH CALCULATES THE C FUNCTIONS AND THE ROWS OF THE JACOBIAN. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE LMSTR(FCN,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,GTOL, C MAXFEV,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV, C NJEV,IPVT,QTF,WA1,WA2,WA3,WA4) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS AND THE ROWS OF THE JACOBIAN. C FCN MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE C USER CALLING PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(M,N,X,FVEC,FJROW,IFLAG) C INTEGER M,N,IFLAG C REAL X(N),FVEC(M),FJROW(N) C ---------- C IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. C IF IFLAG = I CALCULATE THE (I-1)-ST ROW OF THE C JACOBIAN AT X AND RETURN THIS VECTOR IN FJROW. C ---------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF LMSTR. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VARIABLES. N MUST NOT EXCEED M. C C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. C C FVEC IS AN OUTPUT ARRAY OF LENGTH M WHICH CONTAINS C THE FUNCTIONS EVALUATED AT THE OUTPUT X. C C FJAC IS AN OUTPUT N BY N ARRAY. THE UPPER TRIANGLE OF FJAC C CONTAINS AN UPPER TRIANGULAR MATRIX R SUCH THAT C C T T T C P *(JAC *JAC)*P = R *R, C C WHERE P IS A PERMUTATION MATRIX AND JAC IS THE FINAL C CALCULATED JACOBIAN. COLUMN J OF P IS COLUMN IPVT(J) C (SEE BELOW) OF THE IDENTITY MATRIX. THE LOWER TRIANGULAR C PART OF FJAC CONTAINS INFORMATION GENERATED DURING C THE COMPUTATION OF R. C C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. C C FTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION C OCCURS WHEN BOTH THE ACTUAL AND PREDICTED RELATIVE C REDUCTIONS IN THE SUM OF SQUARES ARE AT MOST FTOL. C THEREFORE, FTOL MEASURES THE RELATIVE ERROR DESIRED C IN THE SUM OF SQUARES. C C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION C OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE C ITERATES IS AT MOST XTOL. THEREFORE, XTOL MEASURES THE C RELATIVE ERROR DESIRED IN THE APPROXIMATE SOLUTION. C C GTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION C OCCURS WHEN THE COSINE OF THE ANGLE BETWEEN FVEC AND C ANY COLUMN OF THE JACOBIAN IS AT MOST GTOL IN ABSOLUTE C VALUE. THEREFORE, GTOL MEASURES THE ORTHOGONALITY C DESIRED BETWEEN THE FUNCTION VECTOR AND THE COLUMNS C OF THE JACOBIAN. C C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION C OCCURS WHEN THE NUMBER OF CALLS TO FCN WITH IFLAG = 1 C HAS REACHED MAXFEV. C C DIAG IS AN ARRAY OF LENGTH N. IF MODE = 1 (SEE C BELOW), DIAG IS INTERNALLY SET. IF MODE = 2, DIAG C MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS C MULTIPLICATIVE SCALE FACTORS FOR THE VARIABLES. C C MODE IS AN INTEGER INPUT VARIABLE. IF MODE = 1, THE C VARIABLES WILL BE SCALED INTERNALLY. IF MODE = 2, C THE SCALING IS SPECIFIED BY THE INPUT DIAG. OTHER C VALUES OF MODE ARE EQUIVALENT TO MODE = 1. C C FACTOR IS A POSITIVE INPUT VARIABLE USED IN DETERMINING THE C INITIAL STEP BOUND. THIS BOUND IS SET TO THE PRODUCT OF C FACTOR AND THE EUCLIDEAN NORM OF DIAG*X IF NONZERO, OR ELSE C TO FACTOR ITSELF. IN MOST CASES FACTOR SHOULD LIE IN THE C INTERVAL (.1,100.). 100. IS A GENERALLY RECOMMENDED VALUE. C C NPRINT IS AN INTEGER INPUT VARIABLE THAT ENABLES CONTROLLED C PRINTING OF ITERATES IF IT IS POSITIVE. IN THIS CASE, C FCN IS CALLED WITH IFLAG = 0 AT THE BEGINNING OF THE FIRST C ITERATION AND EVERY NPRINT ITERATIONS THEREAFTER AND C IMMEDIATELY PRIOR TO RETURN, WITH X AND FVEC AVAILABLE C FOR PRINTING. IF NPRINT IS NOT POSITIVE, NO SPECIAL CALLS C OF FCN WITH IFLAG = 0 ARE MADE. C C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE) C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE, C INFO IS SET AS FOLLOWS. C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO = 1 BOTH ACTUAL AND PREDICTED RELATIVE REDUCTIONS C IN THE SUM OF SQUARES ARE AT MOST FTOL. C C INFO = 2 RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES C IS AT MOST XTOL. C C INFO = 3 CONDITIONS FOR INFO = 1 AND INFO = 2 BOTH HOLD. C C INFO = 4 THE COSINE OF THE ANGLE BETWEEN FVEC AND ANY C COLUMN OF THE JACOBIAN IS AT MOST GTOL IN C ABSOLUTE VALUE. C C INFO = 5 NUMBER OF CALLS TO FCN WITH IFLAG = 1 HAS C REACHED MAXFEV. C C INFO = 6 FTOL IS TOO SMALL. NO FURTHER REDUCTION IN C THE SUM OF SQUARES IS POSSIBLE. C C INFO = 7 XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN C THE APPROXIMATE SOLUTION X IS POSSIBLE. C C INFO = 8 GTOL IS TOO SMALL. FVEC IS ORTHOGONAL TO THE C COLUMNS OF THE JACOBIAN TO MACHINE PRECISION. C C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF C CALLS TO FCN WITH IFLAG = 1. C C NJEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF C CALLS TO FCN WITH IFLAG = 2. C C IPVT IS AN INTEGER OUTPUT ARRAY OF LENGTH N. IPVT C DEFINES A PERMUTATION MATRIX P SUCH THAT JAC*P = Q*R, C WHERE JAC IS THE FINAL CALCULATED JACOBIAN, Q IS C ORTHOGONAL (NOT STORED), AND R IS UPPER TRIANGULAR. C COLUMN J OF P IS COLUMN IPVT(J) OF THE IDENTITY MATRIX. C C QTF IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS C THE FIRST N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*FVEC. C C WA1, WA2, AND WA3 ARE WORK ARRAYS OF LENGTH N. C C WA4 IS A WORK ARRAY OF LENGTH M. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ...... FCN C C MINPACK-SUPPLIED ... SPMPAR,ENORM,LMPAR,QRFAC,RWUPDT C C FORTRAN-SUPPLIED ... ABS,AMAX1,AMIN1,SQRT,MOD C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, DUDLEY V. GOETSCHEL, KENNETH E. HILLSTROM, C JORGE J. MORE C C ********** INTEGER I,IFLAG,ITER,J,L REAL ACTRED,DELTA,DIRDER,EPSMCH,FNORM,FNORM1,GNORM,ONE,PAR, * PNORM,PRERED,P1,P5,P25,P75,P0001,RATIO,SUM,TEMP,TEMP1, * TEMP2,XNORM,ZERO REAL SPMPAR,ENORM DATA ONE,P1,P5,P25,P75,P0001,ZERO * /1.0E0,1.0E-1,5.0E-1,2.5E-1,7.5E-1,1.0E-4,0.0E0/ C C EPSMCH IS THE MACHINE PRECISION. C EPSMCH = SPMPAR(1) C INFO = 0 IFLAG = 0 NFEV = 0 NJEV = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. M .LT. N .OR. LDFJAC .LT. N * .OR. FTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. GTOL .LT. ZERO * .OR. MAXFEV .LE. 0 .OR. FACTOR .LE. ZERO) GO TO 340 IF (MODE .NE. 2) GO TO 20 DO 10 J = 1, N IF (DIAG(J) .LE. ZERO) GO TO 340 10 CONTINUE 20 CONTINUE C C EVALUATE THE FUNCTION AT THE STARTING POINT C AND CALCULATE ITS NORM. C IFLAG = 1 CALL FCN(M,N,X,FVEC,WA3,IFLAG) NFEV = 1 IF (IFLAG .LT. 0) GO TO 340 FNORM = ENORM(M,FVEC) C C INITIALIZE LEVENBERG-MARQUARDT PARAMETER AND ITERATION COUNTER. C PAR = ZERO ITER = 1 C C BEGINNING OF THE OUTER LOOP. C 30 CONTINUE C C IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES. C IF (NPRINT .LE. 0) GO TO 40 IFLAG = 0 IF (MOD(ITER-1,NPRINT) .EQ. 0) CALL FCN(M,N,X,FVEC,WA3,IFLAG) IF (IFLAG .LT. 0) GO TO 340 40 CONTINUE C C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN MATRIX C CALCULATED ONE ROW AT A TIME, WHILE SIMULTANEOUSLY C FORMING (Q TRANSPOSE)*FVEC AND STORING THE FIRST C N COMPONENTS IN QTF. C DO 60 J = 1, N QTF(J) = ZERO DO 50 I = 1, N FJAC(I,J) = ZERO 50 CONTINUE 60 CONTINUE IFLAG = 2 DO 70 I = 1, M CALL FCN(M,N,X,FVEC,WA3,IFLAG) IF (IFLAG .LT. 0) GO TO 340 TEMP = FVEC(I) CALL RWUPDT(N,FJAC,LDFJAC,WA3,QTF,TEMP,WA1,WA2) IFLAG = IFLAG + 1 70 CONTINUE NJEV = NJEV + 1 C C IF THE JACOBIAN IS RANK DEFICIENT, CALL QRFAC TO C REORDER ITS COLUMNS AND UPDATE THE COMPONENTS OF QTF. C SING = .FALSE. DO 80 J = 1, N IF (FJAC(J,J) .EQ. ZERO) SING = .TRUE. IPVT(J) = J WA2(J) = ENORM(J,FJAC(1,J)) 80 CONTINUE IF (.NOT.SING) GO TO 130 CALL QRFAC(N,N,FJAC,LDFJAC,.TRUE.,IPVT,N,WA1,WA2,WA3) DO 120 J = 1, N IF (FJAC(J,J) .EQ. ZERO) GO TO 110 SUM = ZERO DO 90 I = J, N SUM = SUM + FJAC(I,J)*QTF(I) 90 CONTINUE TEMP = -SUM/FJAC(J,J) DO 100 I = J, N QTF(I) = QTF(I) + FJAC(I,J)*TEMP 100 CONTINUE 110 CONTINUE FJAC(J,J) = WA1(J) 120 CONTINUE 130 CONTINUE C C ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING C TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN. C IF (ITER .NE. 1) GO TO 170 IF (MODE .EQ. 2) GO TO 150 DO 140 J = 1, N DIAG(J) = WA2(J) IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE 140 CONTINUE 150 CONTINUE C C ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X C AND INITIALIZE THE STEP BOUND DELTA. C DO 160 J = 1, N WA3(J) = DIAG(J)*X(J) 160 CONTINUE XNORM = ENORM(N,WA3) DELTA = FACTOR*XNORM IF (DELTA .EQ. ZERO) DELTA = FACTOR 170 CONTINUE C C COMPUTE THE NORM OF THE SCALED GRADIENT. C GNORM = ZERO IF (FNORM .EQ. ZERO) GO TO 210 DO 200 J = 1, N L = IPVT(J) IF (WA2(L) .EQ. ZERO) GO TO 190 SUM = ZERO DO 180 I = 1, J SUM = SUM + FJAC(I,J)*(QTF(I)/FNORM) 180 CONTINUE GNORM = AMAX1(GNORM,ABS(SUM/WA2(L))) 190 CONTINUE 200 CONTINUE 210 CONTINUE C C TEST FOR CONVERGENCE OF THE GRADIENT NORM. C IF (GNORM .LE. GTOL) INFO = 4 IF (INFO .NE. 0) GO TO 340 C C RESCALE IF NECESSARY. C IF (MODE .EQ. 2) GO TO 230 DO 220 J = 1, N DIAG(J) = AMAX1(DIAG(J),WA2(J)) 220 CONTINUE 230 CONTINUE C C BEGINNING OF THE INNER LOOP. C 240 CONTINUE C C DETERMINE THE LEVENBERG-MARQUARDT PARAMETER. C CALL LMPAR(N,FJAC,LDFJAC,IPVT,DIAG,QTF,DELTA,PAR,WA1,WA2, * WA3,WA4) C C STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P. C DO 250 J = 1, N WA1(J) = -WA1(J) WA2(J) = X(J) + WA1(J) WA3(J) = DIAG(J)*WA1(J) 250 CONTINUE PNORM = ENORM(N,WA3) C C ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND. C IF (ITER .EQ. 1) DELTA = AMIN1(DELTA,PNORM) C C EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM. C IFLAG = 1 CALL FCN(M,N,WA2,WA4,WA3,IFLAG) NFEV = NFEV + 1 IF (IFLAG .LT. 0) GO TO 340 FNORM1 = ENORM(M,WA4) C C COMPUTE THE SCALED ACTUAL REDUCTION. C ACTRED = -ONE IF (P1*FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2 C C COMPUTE THE SCALED PREDICTED REDUCTION AND C THE SCALED DIRECTIONAL DERIVATIVE. C DO 270 J = 1, N WA3(J) = ZERO L = IPVT(J) TEMP = WA1(L) DO 260 I = 1, J WA3(I) = WA3(I) + FJAC(I,J)*TEMP 260 CONTINUE 270 CONTINUE TEMP1 = ENORM(N,WA3)/FNORM TEMP2 = (SQRT(PAR)*PNORM)/FNORM PRERED = TEMP1**2 + TEMP2**2/P5 DIRDER = -(TEMP1**2 + TEMP2**2) C C COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED C REDUCTION. C RATIO = ZERO IF (PRERED .NE. ZERO) RATIO = ACTRED/PRERED C C UPDATE THE STEP BOUND. C IF (RATIO .GT. P25) GO TO 280 IF (ACTRED .GE. ZERO) TEMP = P5 IF (ACTRED .LT. ZERO) * TEMP = P5*DIRDER/(DIRDER + P5*ACTRED) IF (P1*FNORM1 .GE. FNORM .OR. TEMP .LT. P1) TEMP = P1 DELTA = TEMP*AMIN1(DELTA,PNORM/P1) PAR = PAR/TEMP GO TO 300 280 CONTINUE IF (PAR .NE. ZERO .AND. RATIO .LT. P75) GO TO 290 DELTA = PNORM/P5 PAR = P5*PAR 290 CONTINUE 300 CONTINUE C C TEST FOR SUCCESSFUL ITERATION. C IF (RATIO .LT. P0001) GO TO 330 C C SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS. C DO 310 J = 1, N X(J) = WA2(J) WA2(J) = DIAG(J)*X(J) 310 CONTINUE DO 320 I = 1, M FVEC(I) = WA4(I) 320 CONTINUE XNORM = ENORM(N,WA2) FNORM = FNORM1 ITER = ITER + 1 330 CONTINUE C C TESTS FOR CONVERGENCE. C IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL * .AND. P5*RATIO .LE. ONE) INFO = 1 IF (DELTA .LE. XTOL*XNORM) INFO = 2 IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL * .AND. P5*RATIO .LE. ONE .AND. INFO .EQ. 2) INFO = 3 IF (INFO .NE. 0) GO TO 340 C C TESTS FOR TERMINATION AND STRINGENT TOLERANCES. C IF (NFEV .GE. MAXFEV) INFO = 5 IF (ABS(ACTRED) .LE. EPSMCH .AND. PRERED .LE. EPSMCH * .AND. P5*RATIO .LE. ONE) INFO = 6 IF (DELTA .LE. EPSMCH*XNORM) INFO = 7 IF (GNORM .LE. EPSMCH) INFO = 8 IF (INFO .NE. 0) GO TO 340 C C END OF THE INNER LOOP. REPEAT IF ITERATION UNSUCCESSFUL. C IF (RATIO .LT. P0001) GO TO 240 C C END OF THE OUTER LOOP. C GO TO 30 340 CONTINUE C C TERMINATION, EITHER NORMAL OR USER IMPOSED. C IF (IFLAG .LT. 0) INFO = IFLAG IFLAG = 0 IF (NPRINT .GT. 0) CALL FCN(M,N,X,FVEC,WA3,IFLAG) RETURN C C LAST CARD OF SUBROUTINE LMSTR. C END SUBROUTINE LMSTR1(FCN,M,N,X,FVEC,FJAC,LDFJAC,TOL,INFO,IPVT,WA, * LWA) INTEGER M,N,LDFJAC,INFO,LWA INTEGER IPVT(N) REAL TOL REAL X(N),FVEC(M),FJAC(LDFJAC,N),WA(LWA) EXTERNAL FCN C ********** C C SUBROUTINE LMSTR1 C C THE PURPOSE OF LMSTR1 IS TO MINIMIZE THE SUM OF THE SQUARES OF C M NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION OF C THE LEVENBERG-MARQUARDT ALGORITHM WHICH USES MINIMAL STORAGE. C THIS IS DONE BY USING THE MORE GENERAL LEAST-SQUARES SOLVER C LMSTR. THE USER MUST PROVIDE A SUBROUTINE WHICH CALCULATES C THE FUNCTIONS AND THE ROWS OF THE JACOBIAN. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE LMSTR1(FCN,M,N,X,FVEC,FJAC,LDFJAC,TOL,INFO, C IPVT,WA,LWA) C C WHERE C C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH C CALCULATES THE FUNCTIONS AND THE ROWS OF THE JACOBIAN. C FCN MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE C USER CALLING PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS. C C SUBROUTINE FCN(M,N,X,FVEC,FJROW,IFLAG) C INTEGER M,N,IFLAG C REAL X(N),FVEC(M),FJROW(N) C ---------- C IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND C RETURN THIS VECTOR IN FVEC. C IF IFLAG = I CALCULATE THE (I-1)-ST ROW OF THE C JACOBIAN AT X AND RETURN THIS VECTOR IN FJROW. C ---------- C RETURN C END C C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS C THE USER WANTS TO TERMINATE EXECUTION OF LMSTR1. C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF FUNCTIONS. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VARIABLES. N MUST NOT EXCEED M. C C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. C C FVEC IS AN OUTPUT ARRAY OF LENGTH M WHICH CONTAINS C THE FUNCTIONS EVALUATED AT THE OUTPUT X. C C FJAC IS AN OUTPUT N BY N ARRAY. THE UPPER TRIANGLE OF FJAC C CONTAINS AN UPPER TRIANGULAR MATRIX R SUCH THAT C C T T T C P *(JAC *JAC)*P = R *R, C C WHERE P IS A PERMUTATION MATRIX AND JAC IS THE FINAL C CALCULATED JACOBIAN. COLUMN J OF P IS COLUMN IPVT(J) C (SEE BELOW) OF THE IDENTITY MATRIX. THE LOWER TRIANGULAR C PART OF FJAC CONTAINS INFORMATION GENERATED DURING C THE COMPUTATION OF R. C C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. C C TOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS C WHEN THE ALGORITHM ESTIMATES EITHER THAT THE RELATIVE C ERROR IN THE SUM OF SQUARES IS AT MOST TOL OR THAT C THE RELATIVE ERROR BETWEEN X AND THE SOLUTION IS AT C MOST TOL. C C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE) C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE, C INFO IS SET AS FOLLOWS. C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO = 1 ALGORITHM ESTIMATES THAT THE RELATIVE ERROR C IN THE SUM OF SQUARES IS AT MOST TOL. C C INFO = 2 ALGORITHM ESTIMATES THAT THE RELATIVE ERROR C BETWEEN X AND THE SOLUTION IS AT MOST TOL. C C INFO = 3 CONDITIONS FOR INFO = 1 AND INFO = 2 BOTH HOLD. C C INFO = 4 FVEC IS ORTHOGONAL TO THE COLUMNS OF THE C JACOBIAN TO MACHINE PRECISION. C C INFO = 5 NUMBER OF CALLS TO FCN WITH IFLAG = 1 HAS C REACHED 100*(N+1). C C INFO = 6 TOL IS TOO SMALL. NO FURTHER REDUCTION IN C THE SUM OF SQUARES IS POSSIBLE. C C INFO = 7 TOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN C THE APPROXIMATE SOLUTION X IS POSSIBLE. C C IPVT IS AN INTEGER OUTPUT ARRAY OF LENGTH N. IPVT C DEFINES A PERMUTATION MATRIX P SUCH THAT JAC*P = Q*R, C WHERE JAC IS THE FINAL CALCULATED JACOBIAN, Q IS C ORTHOGONAL (NOT STORED), AND R IS UPPER TRIANGULAR. C COLUMN J OF P IS COLUMN IPVT(J) OF THE IDENTITY MATRIX. C C WA IS A WORK ARRAY OF LENGTH LWA. C C LWA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN 5*N+M. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ...... FCN C C MINPACK-SUPPLIED ... LMSTR C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, DUDLEY V. GOETSCHEL, KENNETH E. HILLSTROM, C JORGE J. MORE C C ********** INTEGER MAXFEV,MODE,NFEV,NJEV,NPRINT REAL FACTOR,FTOL,GTOL,XTOL,ZERO DATA FACTOR,ZERO /1.0E2,0.0E0/ INFO = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. M .LT. N .OR. LDFJAC .LT. N .OR. TOL .LT. ZERO * .OR. LWA .LT. 5*N + M) GO TO 10 C C CALL LMSTR. C MAXFEV = 100*(N + 1) FTOL = TOL XTOL = TOL GTOL = ZERO MODE = 1 NPRINT = 0 CALL LMSTR(FCN,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,GTOL,MAXFEV, * WA(1),MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,IPVT,WA(N+1), * WA(2*N+1),WA(3*N+1),WA(4*N+1),WA(5*N+1)) IF (INFO .EQ. 8) INFO = 4 10 CONTINUE RETURN C C LAST CARD OF SUBROUTINE LMSTR1. C END SUBROUTINE QFORM(M,N,Q,LDQ,WA) INTEGER M,N,LDQ REAL Q(LDQ,M),WA(M) C ********** C C SUBROUTINE QFORM C C THIS SUBROUTINE PROCEEDS FROM THE COMPUTED QR FACTORIZATION OF C AN M BY N MATRIX A TO ACCUMULATE THE M BY M ORTHOGONAL MATRIX C Q FROM ITS FACTORED FORM. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE QFORM(M,N,Q,LDQ,WA) C C WHERE C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF ROWS OF A AND THE ORDER OF Q. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C Q IS AN M BY M ARRAY. ON INPUT THE FULL LOWER TRAPEZOID IN C THE FIRST MIN(M,N) COLUMNS OF Q CONTAINS THE FACTORED FORM. C ON OUTPUT Q HAS BEEN ACCUMULATED INTO A SQUARE MATRIX. C C LDQ IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY Q. C C WA IS A WORK ARRAY OF LENGTH M. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... MIN0 C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,J,JM1,K,L,MINMN,NP1 REAL ONE,SUM,TEMP,ZERO DATA ONE,ZERO /1.0E0,0.0E0/ C C ZERO OUT UPPER TRIANGLE OF Q IN THE FIRST MIN(M,N) COLUMNS. C MINMN = MIN0(M,N) IF (MINMN .LT. 2) GO TO 30 DO 20 J = 2, MINMN JM1 = J - 1 DO 10 I = 1, JM1 Q(I,J) = ZERO 10 CONTINUE 20 CONTINUE 30 CONTINUE C C INITIALIZE REMAINING COLUMNS TO THOSE OF THE IDENTITY MATRIX. C NP1 = N + 1 IF (M .LT. NP1) GO TO 60 DO 50 J = NP1, M DO 40 I = 1, M Q(I,J) = ZERO 40 CONTINUE Q(J,J) = ONE 50 CONTINUE 60 CONTINUE C C ACCUMULATE Q FROM ITS FACTORED FORM. C DO 120 L = 1, MINMN K = MINMN - L + 1 DO 70 I = K, M WA(I) = Q(I,K) Q(I,K) = ZERO 70 CONTINUE Q(K,K) = ONE IF (WA(K) .EQ. ZERO) GO TO 110 DO 100 J = K, M SUM = ZERO DO 80 I = K, M SUM = SUM + Q(I,J)*WA(I) 80 CONTINUE TEMP = SUM/WA(K) DO 90 I = K, M Q(I,J) = Q(I,J) - TEMP*WA(I) 90 CONTINUE 100 CONTINUE 110 CONTINUE 120 CONTINUE RETURN C C LAST CARD OF SUBROUTINE QFORM. C END SUBROUTINE QRFAC(M,N,A,LDA,PIVOT,IPVT,LIPVT,RDIAG,ACNORM,WA) INTEGER M,N,LDA,LIPVT INTEGER IPVT(LIPVT) LOGICAL PIVOT REAL A(LDA,N),RDIAG(N),ACNORM(N),WA(N) C ********** C C SUBROUTINE QRFAC C C THIS SUBROUTINE USES HOUSEHOLDER TRANSFORMATIONS WITH COLUMN C PIVOTING (OPTIONAL) TO COMPUTE A QR FACTORIZATION OF THE C M BY N MATRIX A. THAT IS, QRFAC DETERMINES AN ORTHOGONAL C MATRIX Q, A PERMUTATION MATRIX P, AND AN UPPER TRAPEZOIDAL C MATRIX R WITH DIAGONAL ELEMENTS OF NONINCREASING MAGNITUDE, C SUCH THAT A*P = Q*R. THE HOUSEHOLDER TRANSFORMATION FOR C COLUMN K, K = 1,2,...,MIN(M,N), IS OF THE FORM C C T C I - (1/U(K))*U*U C C WHERE U HAS ZEROS IN THE FIRST K-1 POSITIONS. THE FORM OF C THIS TRANSFORMATION AND THE METHOD OF PIVOTING FIRST C APPEARED IN THE CORRESPONDING LINPACK SUBROUTINE. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE QRFAC(M,N,A,LDA,PIVOT,IPVT,LIPVT,RDIAG,ACNORM,WA) C C WHERE C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF ROWS OF A. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C A IS AN M BY N ARRAY. ON INPUT A CONTAINS THE MATRIX FOR C WHICH THE QR FACTORIZATION IS TO BE COMPUTED. ON OUTPUT C THE STRICT UPPER TRAPEZOIDAL PART OF A CONTAINS THE STRICT C UPPER TRAPEZOIDAL PART OF R, AND THE LOWER TRAPEZOIDAL C PART OF A CONTAINS A FACTORED FORM OF Q (THE NON-TRIVIAL C ELEMENTS OF THE U VECTORS DESCRIBED ABOVE). C C LDA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A. C C PIVOT IS A LOGICAL INPUT VARIABLE. IF PIVOT IS SET TRUE, C THEN COLUMN PIVOTING IS ENFORCED. IF PIVOT IS SET FALSE, C THEN NO COLUMN PIVOTING IS DONE. C C IPVT IS AN INTEGER OUTPUT ARRAY OF LENGTH LIPVT. IPVT C DEFINES THE PERMUTATION MATRIX P SUCH THAT A*P = Q*R. C COLUMN J OF P IS COLUMN IPVT(J) OF THE IDENTITY MATRIX. C IF PIVOT IS FALSE, IPVT IS NOT REFERENCED. C C LIPVT IS A POSITIVE INTEGER INPUT VARIABLE. IF PIVOT IS FALSE, C THEN LIPVT MAY BE AS SMALL AS 1. IF PIVOT IS TRUE, THEN C LIPVT MUST BE AT LEAST N. C C RDIAG IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE C DIAGONAL ELEMENTS OF R. C C ACNORM IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE C NORMS OF THE CORRESPONDING COLUMNS OF THE INPUT MATRIX A. C IF THIS INFORMATION IS NOT NEEDED, THEN ACNORM CAN COINCIDE C WITH RDIAG. C C WA IS A WORK ARRAY OF LENGTH N. IF PIVOT IS FALSE, THEN WA C CAN COINCIDE WITH RDIAG. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... SPMPAR,ENORM C C FORTRAN-SUPPLIED ... AMAX1,SQRT,MIN0 C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,J,JP1,K,KMAX,MINMN REAL AJNORM,EPSMCH,ONE,P05,SUM,TEMP,ZERO REAL SPMPAR,ENORM DATA ONE,P05,ZERO /1.0E0,5.0E-2,0.0E0/ C C EPSMCH IS THE MACHINE PRECISION. C EPSMCH = SPMPAR(1) C C COMPUTE THE INITIAL COLUMN NORMS AND INITIALIZE SEVERAL ARRAYS. C DO 10 J = 1, N ACNORM(J) = ENORM(M,A(1,J)) RDIAG(J) = ACNORM(J) WA(J) = RDIAG(J) IF (PIVOT) IPVT(J) = J 10 CONTINUE C C REDUCE A TO R WITH HOUSEHOLDER TRANSFORMATIONS. C MINMN = MIN0(M,N) DO 110 J = 1, MINMN IF (.NOT.PIVOT) GO TO 40 C C BRING THE COLUMN OF LARGEST NORM INTO THE PIVOT POSITION. C KMAX = J DO 20 K = J, N IF (RDIAG(K) .GT. RDIAG(KMAX)) KMAX = K 20 CONTINUE IF (KMAX .EQ. J) GO TO 40 DO 30 I = 1, M TEMP = A(I,J) A(I,J) = A(I,KMAX) A(I,KMAX) = TEMP 30 CONTINUE RDIAG(KMAX) = RDIAG(J) WA(KMAX) = WA(J) K = IPVT(J) IPVT(J) = IPVT(KMAX) IPVT(KMAX) = K 40 CONTINUE C C COMPUTE THE HOUSEHOLDER TRANSFORMATION TO REDUCE THE C J-TH COLUMN OF A TO A MULTIPLE OF THE J-TH UNIT VECTOR. C AJNORM = ENORM(M-J+1,A(J,J)) IF (AJNORM .EQ. ZERO) GO TO 100 IF (A(J,J) .LT. ZERO) AJNORM = -AJNORM DO 50 I = J, M A(I,J) = A(I,J)/AJNORM 50 CONTINUE A(J,J) = A(J,J) + ONE C C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS C AND UPDATE THE NORMS. C JP1 = J + 1 IF (N .LT. JP1) GO TO 100 DO 90 K = JP1, N SUM = ZERO DO 60 I = J, M SUM = SUM + A(I,J)*A(I,K) 60 CONTINUE TEMP = SUM/A(J,J) DO 70 I = J, M A(I,K) = A(I,K) - TEMP*A(I,J) 70 CONTINUE IF (.NOT.PIVOT .OR. RDIAG(K) .EQ. ZERO) GO TO 80 TEMP = A(J,K)/RDIAG(K) RDIAG(K) = RDIAG(K)*SQRT(AMAX1(ZERO,ONE-TEMP**2)) IF (P05*(RDIAG(K)/WA(K))**2 .GT. EPSMCH) GO TO 80 RDIAG(K) = ENORM(M-J,A(JP1,K)) WA(K) = RDIAG(K) 80 CONTINUE 90 CONTINUE 100 CONTINUE RDIAG(J) = -AJNORM 110 CONTINUE RETURN C C LAST CARD OF SUBROUTINE QRFAC. C END SUBROUTINE QRSOLV(N,R,LDR,IPVT,DIAG,QTB,X,SDIAG,WA) INTEGER N,LDR INTEGER IPVT(N) REAL R(LDR,N),DIAG(N),QTB(N),X(N),SDIAG(N),WA(N) C ********** C C SUBROUTINE QRSOLV C C GIVEN AN M BY N MATRIX A, AN N BY N DIAGONAL MATRIX D, C AND AN M-VECTOR B, THE PROBLEM IS TO DETERMINE AN X WHICH C SOLVES THE SYSTEM C C A*X = B , D*X = 0 , C C IN THE LEAST SQUARES SENSE. C C THIS SUBROUTINE COMPLETES THE SOLUTION OF THE PROBLEM C IF IT IS PROVIDED WITH THE NECESSARY INFORMATION FROM THE C QR FACTORIZATION, WITH COLUMN PIVOTING, OF A. THAT IS, IF C A*P = Q*R, WHERE P IS A PERMUTATION MATRIX, Q HAS ORTHOGONAL C COLUMNS, AND R IS AN UPPER TRIANGULAR MATRIX WITH DIAGONAL C ELEMENTS OF NONINCREASING MAGNITUDE, THEN QRSOLV EXPECTS C THE FULL UPPER TRIANGLE OF R, THE PERMUTATION MATRIX P, C AND THE FIRST N COMPONENTS OF (Q TRANSPOSE)*B. THE SYSTEM C A*X = B, D*X = 0, IS THEN EQUIVALENT TO C C T T C R*Z = Q *B , P *D*P*Z = 0 , C C WHERE X = P*Z. IF THIS SYSTEM DOES NOT HAVE FULL RANK, C THEN A LEAST SQUARES SOLUTION IS OBTAINED. ON OUTPUT QRSOLV C ALSO PROVIDES AN UPPER TRIANGULAR MATRIX S SUCH THAT C C T T T C P *(A *A + D*D)*P = S *S . C C S IS COMPUTED WITHIN QRSOLV AND MAY BE OF SEPARATE INTEREST. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE QRSOLV(N,R,LDR,IPVT,DIAG,QTB,X,SDIAG,WA) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER OF R. C C R IS AN N BY N ARRAY. ON INPUT THE FULL UPPER TRIANGLE C MUST CONTAIN THE FULL UPPER TRIANGLE OF THE MATRIX R. C ON OUTPUT THE FULL UPPER TRIANGLE IS UNALTERED, AND THE C STRICT LOWER TRIANGLE CONTAINS THE STRICT UPPER TRIANGLE C (TRANSPOSED) OF THE UPPER TRIANGULAR MATRIX S. C C LDR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY R. C C IPVT IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH DEFINES THE C PERMUTATION MATRIX P SUCH THAT A*P = Q*R. COLUMN J OF P C IS COLUMN IPVT(J) OF THE IDENTITY MATRIX. C C DIAG IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE C DIAGONAL ELEMENTS OF THE MATRIX D. C C QTB IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE FIRST C N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*B. C C X IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE LEAST C SQUARES SOLUTION OF THE SYSTEM A*X = B, D*X = 0. C C SDIAG IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE C DIAGONAL ELEMENTS OF THE UPPER TRIANGULAR MATRIX S. C C WA IS A WORK ARRAY OF LENGTH N. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... ABS,SQRT C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,J,JP1,K,KP1,L,NSING REAL COS,COTAN,P5,P25,QTBPJ,SIN,SUM,TAN,TEMP,ZERO DATA P5,P25,ZERO /5.0E-1,2.5E-1,0.0E0/ C C COPY R AND (Q TRANSPOSE)*B TO PRESERVE INPUT AND INITIALIZE S. C IN PARTICULAR, SAVE THE DIAGONAL ELEMENTS OF R IN X. C DO 20 J = 1, N DO 10 I = J, N R(I,J) = R(J,I) 10 CONTINUE X(J) = R(J,J) WA(J) = QTB(J) 20 CONTINUE C C ELIMINATE THE DIAGONAL MATRIX D USING A GIVENS ROTATION. C DO 100 J = 1, N C C PREPARE THE ROW OF D TO BE ELIMINATED, LOCATING THE C DIAGONAL ELEMENT USING P FROM THE QR FACTORIZATION. C L = IPVT(J) IF (DIAG(L) .EQ. ZERO) GO TO 90 DO 30 K = J, N SDIAG(K) = ZERO 30 CONTINUE SDIAG(J) = DIAG(L) C C THE TRANSFORMATIONS TO ELIMINATE THE ROW OF D C MODIFY ONLY A SINGLE ELEMENT OF (Q TRANSPOSE)*B C BEYOND THE FIRST N, WHICH IS INITIALLY ZERO. C QTBPJ = ZERO DO 80 K = J, N C C DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE C APPROPRIATE ELEMENT IN THE CURRENT ROW OF D. C IF (SDIAG(K) .EQ. ZERO) GO TO 70 IF (ABS(R(K,K)) .GE. ABS(SDIAG(K))) GO TO 40 COTAN = R(K,K)/SDIAG(K) SIN = P5/SQRT(P25+P25*COTAN**2) COS = SIN*COTAN GO TO 50 40 CONTINUE TAN = SDIAG(K)/R(K,K) COS = P5/SQRT(P25+P25*TAN**2) SIN = COS*TAN 50 CONTINUE C C COMPUTE THE MODIFIED DIAGONAL ELEMENT OF R AND C THE MODIFIED ELEMENT OF ((Q TRANSPOSE)*B,0). C R(K,K) = COS*R(K,K) + SIN*SDIAG(K) TEMP = COS*WA(K) + SIN*QTBPJ QTBPJ = -SIN*WA(K) + COS*QTBPJ WA(K) = TEMP C C ACCUMULATE THE TRANFORMATION IN THE ROW OF S. C KP1 = K + 1 IF (N .LT. KP1) GO TO 70 DO 60 I = KP1, N TEMP = COS*R(I,K) + SIN*SDIAG(I) SDIAG(I) = -SIN*R(I,K) + COS*SDIAG(I) R(I,K) = TEMP 60 CONTINUE 70 CONTINUE 80 CONTINUE 90 CONTINUE C C STORE THE DIAGONAL ELEMENT OF S AND RESTORE C THE CORRESPONDING DIAGONAL ELEMENT OF R. C SDIAG(J) = R(J,J) R(J,J) = X(J) 100 CONTINUE C C SOLVE THE TRIANGULAR SYSTEM FOR Z. IF THE SYSTEM IS C SINGULAR, THEN OBTAIN A LEAST SQUARES SOLUTION. C NSING = N DO 110 J = 1, N IF (SDIAG(J) .EQ. ZERO .AND. NSING .EQ. N) NSING = J - 1 IF (NSING .LT. N) WA(J) = ZERO 110 CONTINUE IF (NSING .LT. 1) GO TO 150 DO 140 K = 1, NSING J = NSING - K + 1 SUM = ZERO JP1 = J + 1 IF (NSING .LT. JP1) GO TO 130 DO 120 I = JP1, NSING SUM = SUM + R(I,J)*WA(I) 120 CONTINUE 130 CONTINUE WA(J) = (WA(J) - SUM)/SDIAG(J) 140 CONTINUE 150 CONTINUE C C PERMUTE THE COMPONENTS OF Z BACK TO COMPONENTS OF X. C DO 160 J = 1, N L = IPVT(J) X(L) = WA(J) 160 CONTINUE RETURN C C LAST CARD OF SUBROUTINE QRSOLV. C END SUBROUTINE RWUPDT(N,R,LDR,W,B,ALPHA,COS,SIN) INTEGER N,LDR REAL ALPHA REAL R(LDR,N),W(N),B(N),COS(N),SIN(N) C ********** C C SUBROUTINE RWUPDT C C GIVEN AN N BY N UPPER TRIANGULAR MATRIX R, THIS SUBROUTINE C COMPUTES THE QR DECOMPOSITION OF THE MATRIX FORMED WHEN A ROW C IS ADDED TO R. IF THE ROW IS SPECIFIED BY THE VECTOR W, THEN C RWUPDT DETERMINES AN ORTHOGONAL MATRIX Q SUCH THAT WHEN THE C N+1 BY N MATRIX COMPOSED OF R AUGMENTED BY W IS PREMULTIPLIED C BY (Q TRANSPOSE), THE RESULTING MATRIX IS UPPER TRAPEZOIDAL. C THE MATRIX (Q TRANSPOSE) IS THE PRODUCT OF N TRANSFORMATIONS C C G(N)*G(N-1)* ... *G(1) C C WHERE G(I) IS A GIVENS ROTATION IN THE (I,N+1) PLANE WHICH C ELIMINATES ELEMENTS IN THE (N+1)-ST PLANE. RWUPDT ALSO C COMPUTES THE PRODUCT (Q TRANSPOSE)*C WHERE C IS THE C (N+1)-VECTOR (B,ALPHA). Q ITSELF IS NOT ACCUMULATED, RATHER C THE INFORMATION TO RECOVER THE G ROTATIONS IS SUPPLIED. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE RWUPDT(N,R,LDR,W,B,ALPHA,COS,SIN) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER OF R. C C R IS AN N BY N ARRAY. ON INPUT THE UPPER TRIANGULAR PART OF C R MUST CONTAIN THE MATRIX TO BE UPDATED. ON OUTPUT R C CONTAINS THE UPDATED TRIANGULAR MATRIX. C C LDR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY R. C C W IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE ROW C VECTOR TO BE ADDED TO R. C C B IS AN ARRAY OF LENGTH N. ON INPUT B MUST CONTAIN THE C FIRST N ELEMENTS OF THE VECTOR C. ON OUTPUT B CONTAINS C THE FIRST N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*C. C C ALPHA IS A VARIABLE. ON INPUT ALPHA MUST CONTAIN THE C (N+1)-ST ELEMENT OF THE VECTOR C. ON OUTPUT ALPHA CONTAINS C THE (N+1)-ST ELEMENT OF THE VECTOR (Q TRANSPOSE)*C. C C COS IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE C COSINES OF THE TRANSFORMING GIVENS ROTATIONS. C C SIN IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE C SINES OF THE TRANSFORMING GIVENS ROTATIONS. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... ABS,SQRT C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, DUDLEY V. GOETSCHEL, KENNETH E. HILLSTROM, C JORGE J. MORE C C ********** INTEGER I,J,JM1 REAL COTAN,ONE,P5,P25,ROWJ,TAN,TEMP,ZERO DATA ONE,P5,P25,ZERO /1.0E0,5.0E-1,2.5E-1,0.0E0/ C DO 60 J = 1, N ROWJ = W(J) JM1 = J - 1 C C APPLY THE PREVIOUS TRANSFORMATIONS TO C R(I,J), I=1,2,...,J-1, AND TO W(J). C IF (JM1 .LT. 1) GO TO 20 DO 10 I = 1, JM1 TEMP = COS(I)*R(I,J) + SIN(I)*ROWJ ROWJ = -SIN(I)*R(I,J) + COS(I)*ROWJ R(I,J) = TEMP 10 CONTINUE 20 CONTINUE C C DETERMINE A GIVENS ROTATION WHICH ELIMINATES W(J). C COS(J) = ONE SIN(J) = ZERO IF (ROWJ .EQ. ZERO) GO TO 50 IF (ABS(R(J,J)) .GE. ABS(ROWJ)) GO TO 30 COTAN = R(J,J)/ROWJ SIN(J) = P5/SQRT(P25+P25*COTAN**2) COS(J) = SIN(J)*COTAN GO TO 40 30 CONTINUE TAN = ROWJ/R(J,J) COS(J) = P5/SQRT(P25+P25*TAN**2) SIN(J) = COS(J)*TAN 40 CONTINUE C C APPLY THE CURRENT TRANSFORMATION TO R(J,J), B(J), AND ALPHA. C R(J,J) = COS(J)*R(J,J) + SIN(J)*ROWJ TEMP = COS(J)*B(J) + SIN(J)*ALPHA ALPHA = -SIN(J)*B(J) + COS(J)*ALPHA B(J) = TEMP 50 CONTINUE 60 CONTINUE RETURN C C LAST CARD OF SUBROUTINE RWUPDT. C END SUBROUTINE R1MPYQ(M,N,A,LDA,V,W) INTEGER M,N,LDA REAL A(LDA,N),V(N),W(N) C ********** C C SUBROUTINE R1MPYQ C C GIVEN AN M BY N MATRIX A, THIS SUBROUTINE COMPUTES A*Q WHERE C Q IS THE PRODUCT OF 2*(N - 1) TRANSFORMATIONS C C GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1) C C AND GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE WHICH C ELIMINATE ELEMENTS IN THE I-TH AND N-TH PLANES, RESPECTIVELY. C Q ITSELF IS NOT GIVEN, RATHER THE INFORMATION TO RECOVER THE C GV, GW ROTATIONS IS SUPPLIED. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE R1MPYQ(M,N,A,LDA,V,W) C C WHERE C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF ROWS OF A. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF A. C C A IS AN M BY N ARRAY. ON INPUT A MUST CONTAIN THE MATRIX C TO BE POSTMULTIPLIED BY THE ORTHOGONAL MATRIX Q C DESCRIBED ABOVE. ON OUTPUT A*Q HAS REPLACED A. C C LDA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A. C C V IS AN INPUT ARRAY OF LENGTH N. V(I) MUST CONTAIN THE C INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GV(I) C DESCRIBED ABOVE. C C W IS AN INPUT ARRAY OF LENGTH N. W(I) MUST CONTAIN THE C INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GW(I) C DESCRIBED ABOVE. C C SUBROUTINES CALLED C C FORTRAN-SUPPLIED ... ABS,SQRT C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER I,J,NMJ,NM1 REAL COS,ONE,SIN,TEMP DATA ONE /1.0E0/ C C APPLY THE FIRST SET OF GIVENS ROTATIONS TO A. C NM1 = N - 1 IF (NM1 .LT. 1) GO TO 50 DO 20 NMJ = 1, NM1 J = N - NMJ IF (ABS(V(J)) .GT. ONE) COS = ONE/V(J) IF (ABS(V(J)) .GT. ONE) SIN = SQRT(ONE-COS**2) IF (ABS(V(J)) .LE. ONE) SIN = V(J) IF (ABS(V(J)) .LE. ONE) COS = SQRT(ONE-SIN**2) DO 10 I = 1, M TEMP = COS*A(I,J) - SIN*A(I,N) A(I,N) = SIN*A(I,J) + COS*A(I,N) A(I,J) = TEMP 10 CONTINUE 20 CONTINUE C C APPLY THE SECOND SET OF GIVENS ROTATIONS TO A. C DO 40 J = 1, NM1 IF (ABS(W(J)) .GT. ONE) COS = ONE/W(J) IF (ABS(W(J)) .GT. ONE) SIN = SQRT(ONE-COS**2) IF (ABS(W(J)) .LE. ONE) SIN = W(J) IF (ABS(W(J)) .LE. ONE) COS = SQRT(ONE-SIN**2) DO 30 I = 1, M TEMP = COS*A(I,J) + SIN*A(I,N) A(I,N) = -SIN*A(I,J) + COS*A(I,N) A(I,J) = TEMP 30 CONTINUE 40 CONTINUE 50 CONTINUE RETURN C C LAST CARD OF SUBROUTINE R1MPYQ. C END SUBROUTINE R1UPDT(M,N,S,LS,U,V,W,SING) INTEGER M,N,LS LOGICAL SING REAL S(LS),U(M),V(N),W(M) C ********** C C SUBROUTINE R1UPDT C C GIVEN AN M BY N LOWER TRAPEZOIDAL MATRIX S, AN M-VECTOR U, C AND AN N-VECTOR V, THE PROBLEM IS TO DETERMINE AN C ORTHOGONAL MATRIX Q SUCH THAT C C T C (S + U*V )*Q C C IS AGAIN LOWER TRAPEZOIDAL. C C THIS SUBROUTINE DETERMINES Q AS THE PRODUCT OF 2*(N - 1) C TRANSFORMATIONS C C GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1) C C WHERE GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE C WHICH ELIMINATE ELEMENTS IN THE I-TH AND N-TH PLANES, C RESPECTIVELY. Q ITSELF IS NOT ACCUMULATED, RATHER THE C INFORMATION TO RECOVER THE GV, GW ROTATIONS IS RETURNED. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE R1UPDT(M,N,S,LS,U,V,W,SING) C C WHERE C C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF ROWS OF S. C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF COLUMNS OF S. N MUST NOT EXCEED M. C C S IS AN ARRAY OF LENGTH LS. ON INPUT S MUST CONTAIN THE LOWER C TRAPEZOIDAL MATRIX S STORED BY COLUMNS. ON OUTPUT S CONTAINS C THE LOWER TRAPEZOIDAL MATRIX PRODUCED AS DESCRIBED ABOVE. C C LS IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN C (N*(2*M-N+1))/2. C C U IS AN INPUT ARRAY OF LENGTH M WHICH MUST CONTAIN THE C VECTOR U. C C V IS AN ARRAY OF LENGTH N. ON INPUT V MUST CONTAIN THE VECTOR C V. ON OUTPUT V(I) CONTAINS THE INFORMATION NECESSARY TO C RECOVER THE GIVENS ROTATION GV(I) DESCRIBED ABOVE. C C W IS AN OUTPUT ARRAY OF LENGTH M. W(I) CONTAINS INFORMATION C NECESSARY TO RECOVER THE GIVENS ROTATION GW(I) DESCRIBED C ABOVE. C C SING IS A LOGICAL OUTPUT VARIABLE. SING IS SET TRUE IF ANY C OF THE DIAGONAL ELEMENTS OF THE OUTPUT S ARE ZERO. OTHERWISE C SING IS SET FALSE. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... SPMPAR C C FORTRAN-SUPPLIED ... ABS,SQRT C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE, C JOHN L. NAZARETH C C ********** INTEGER I,J,JJ,L,NMJ,NM1 REAL COS,COTAN,GIANT,ONE,P5,P25,SIN,TAN,TAU,TEMP,ZERO REAL SPMPAR DATA ONE,P5,P25,ZERO /1.0E0,5.0E-1,2.5E-1,0.0E0/ C C GIANT IS THE LARGEST MAGNITUDE. C GIANT = SPMPAR(3) C C INITIALIZE THE DIAGONAL ELEMENT POINTER. C JJ = (N*(2*M - N + 1))/2 - (M - N) C C MOVE THE NONTRIVIAL PART OF THE LAST COLUMN OF S INTO W. C L = JJ DO 10 I = N, M W(I) = S(L) L = L + 1 10 CONTINUE C C ROTATE THE VECTOR V INTO A MULTIPLE OF THE N-TH UNIT VECTOR C IN SUCH A WAY THAT A SPIKE IS INTRODUCED INTO W. C NM1 = N - 1 IF (NM1 .LT. 1) GO TO 70 DO 60 NMJ = 1, NM1 J = N - NMJ JJ = JJ - (M - J + 1) W(J) = ZERO IF (V(J) .EQ. ZERO) GO TO 50 C C DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE C J-TH ELEMENT OF V. C IF (ABS(V(N)) .GE. ABS(V(J))) GO TO 20 COTAN = V(N)/V(J) SIN = P5/SQRT(P25+P25*COTAN**2) COS = SIN*COTAN TAU = ONE IF (ABS(COS)*GIANT .GT. ONE) TAU = ONE/COS GO TO 30 20 CONTINUE TAN = V(J)/V(N) COS = P5/SQRT(P25+P25*TAN**2) SIN = COS*TAN TAU = SIN 30 CONTINUE C C APPLY THE TRANSFORMATION TO V AND STORE THE INFORMATION C NECESSARY TO RECOVER THE GIVENS ROTATION. C V(N) = SIN*V(J) + COS*V(N) V(J) = TAU C C APPLY THE TRANSFORMATION TO S AND EXTEND THE SPIKE IN W. C L = JJ DO 40 I = J, M TEMP = COS*S(L) - SIN*W(I) W(I) = SIN*S(L) + COS*W(I) S(L) = TEMP L = L + 1 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 CONTINUE C C ADD THE SPIKE FROM THE RANK 1 UPDATE TO W. C DO 80 I = 1, M W(I) = W(I) + V(N)*U(I) 80 CONTINUE C C ELIMINATE THE SPIKE. C SING = .FALSE. IF (NM1 .LT. 1) GO TO 140 DO 130 J = 1, NM1 IF (W(J) .EQ. ZERO) GO TO 120 C C DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE C J-TH ELEMENT OF THE SPIKE. C IF (ABS(S(JJ)) .GE. ABS(W(J))) GO TO 90 COTAN = S(JJ)/W(J) SIN = P5/SQRT(P25+P25*COTAN**2) COS = SIN*COTAN TAU = ONE IF (ABS(COS)*GIANT .GT. ONE) TAU = ONE/COS GO TO 100 90 CONTINUE TAN = W(J)/S(JJ) COS = P5/SQRT(P25+P25*TAN**2) SIN = COS*TAN TAU = SIN 100 CONTINUE C C APPLY THE TRANSFORMATION TO S AND REDUCE THE SPIKE IN W. C L = JJ DO 110 I = J, M TEMP = COS*S(L) + SIN*W(I) W(I) = -SIN*S(L) + COS*W(I) S(L) = TEMP L = L + 1 110 CONTINUE C C STORE THE INFORMATION NECESSARY TO RECOVER THE C GIVENS ROTATION. C W(J) = TAU 120 CONTINUE C C TEST FOR ZERO DIAGONAL ELEMENTS IN THE OUTPUT S. C IF (S(JJ) .EQ. ZERO) SING = .TRUE. JJ = JJ + (M - J + 1) 130 CONTINUE 140 CONTINUE C C MOVE W BACK INTO THE LAST COLUMN OF THE OUTPUT S. C L = JJ DO 150 I = N, M S(L) = W(I) L = L + 1 150 CONTINUE IF (S(JJ) .EQ. ZERO) SING = .TRUE. RETURN C C LAST CARD OF SUBROUTINE R1UPDT. C END .