* From research!csnet!UCLA-LOCUS!rjl Tue Jun 25 12:47:53 PDT 1985 * Here is the VARP2 routine which handles multiple right hand sides. * Prof. Randall LeVeque C C ============================================================== C SAMPLE DRIVER PROGRAM FOR NONLINEAR LEAST SQUARES ROUTINE C DATA FOLLOWS SUBROUTINES C ============================================================== C C DOUBLE PRECISION Y(50,2), T(50, 1), ALF(12), BETA(8,2), X A(100,13), W(50), DSQRT, ETA INTEGER S, P, OUTPUT EXTERNAL ADA DATA INPUT /5/, OUTPUT /6/ C C READ AND LIST THE DATA C NMAX = 50 NDIM = 100 LMAX = 8 IPRINT = 1 READ (INPUT, 100) N, L, NL, P, IV, S WRITE (OUTPUT, 101) N, L, NL, P, IV, S DO 12 I = 1, N W(I) = 1.0 12 READ (INPUT, 105) (T(I,J), J = 1, IV), (Y(I,J), J=1,S) DO 15 I = 1, N 15 WRITE (OUTPUT, 213) I, (T(I,J), J = 1, IV), (Y(I,J), J=1,S) C READ (INPUT, 102) (ALF(K), K = 1, NL) WRITE (OUTPUT, 103) (ALF(K), K = 1, NL) WRITE (OUTPUT, 104) C CALL VARP2 (S, L, LMAX, NL, N, NMAX, NDIM, L+P+S+1, IV, T, Y, W, X ADA, A, IPRINT, ALF, BETA, IERR) C DO 60 K=1,S LK = L+K WRITE(OUTPUT,215) K WRITE(OUTPUT,212) DO 60 I=1,N ETA = Y(I,K) - A(I,LK) / DSQRT(W(I)) 60 WRITE(OUTPUT,213) I,W(I), (T(I,J),J=1,IV), Y(I,K), X ETA, A(I,LK) 80 STOP C 100 FORMAT (6I5) 101 FORMAT (36H1 NON-LINEAR LEAST SQUARES PROBLEM // 1 25H NUMBER OF OBSERVATIONS =, I5, 32H NUMBER OF LINEAR PARAMETE 2RS = , I4 // 33H NUMBER OF NONLINEAR PARAMETERS =, I4, 3 45H NUMBER OF NONVANISHING PARTIAL DERIVATIVES =, I4 // 4 34H NUMBER OF INDEPENDENT VARIABLES =, I4 // 5 29H NUMBER OF RIGHT HAND SIDES =,I4 // 6 33H I T(I) Y(I) //) 102 FORMAT (4E20.7) 103 FORMAT (30H0 INITIAL NONLINEAR PARAMETERS //(4E20.7)) 104 FORMAT (1H0, 50(1H*)) 105 FORMAT (5E15.7) 212 FORMAT (/89H I W(I) T(I) Y(I) X PREDICTED Y WEIGHTED RESIDUAL //) 213 FORMAT (I5, 7E16.7) 215 FORMAT (//,' RIGHT HAND SIDE NUMBER',I3,/) END C ============================================================== SUBROUTINE ADA (S, LP1, NL, N, NMAX, NDIM, LPP2, IV, A, INC, T, X ALF, ISEL) C ============================================================== C C EVALUATE THE FUNCTIONS PHI AND THEIR PARTIAL DERIVATIVES C D PHI(J)/D ALF(K), AT THE SAMPLE POINTS T(I). C ISEL = 1 MEANS COMPUTE BOTH FUNCTIONS AND DERIVATIVES AND C INITIALIZE INC AND CONSTANT PHI'S C = 2 MEANS COMPUTE ONLY THE NONCONSTANT FUNCTIONS PHI C = 3 MEANS COMPUTE ONLY THE DERIVATIVES. C C THIS PARTICULAR ROUTINE IS FOR FITTING NL EXPONENTIALS AND C A CONSTANT TERM IF LP1 = NL+2, NO CONSTANT IF LP1 = NL+1 C C NL C ETA (C,ALF; T) = C + SUM C * EXP(ALF *T) C J 0 I=1 I,J I C C WHERE C IS THE MATRIX OF LINEAR PARAMETERS TO BE DETERMINED, C AND ALF IS THE VECTOR OF NONLINEAR PARAMETERS TO BE DETERMINED. C C .................................................................. C DOUBLE PRECISION A(NDIM, LPP2), ALF(NL), T(NMAX), DEXP INTEGER INC(12, LP1),S C C SET KONST TO 1 IF THERE IS A CONSTANT FCN, 0 OTHERWISE: KONST = IABS(LP1-NL-1) GO TO (10, 20, 30), ISEL C C ISEL = 1 C 10 DO 11 I=1,NL 11 INC(I,I+KONST) = 1 C C ISEL = 1, SET CONSTANT FUNCTION PHI(1) IF PRESENT IF (KONST .EQ. 0) GO TO 20 C DO 12 I = 1, N 12 A(I,1) = 1.0 C C ISEL = 1 OR 2, COMPUTE NONCONSTANT FUNCTIONS PHI C 20 DO 22 I = 1, N DO 22 J=1,NL 22 A(I,J+KONST) = DEXP(ALF(J) * T(I)) C COLUMNS L+1 THRU L+S ARE LEFT FOR WORKSPACE. IF (ISEL .EQ. 2) GO TO 35 C C ISEL = 1 OR 3, COMPUTE DERIVATIVES DPHI(I) / D ALF(J) C 30 DO 32 I = 1, N DO 32 J=1,NL 32 A(I,J+LP1-1+S) = T(I) * DEXP(ALF(J)*T(I)) 35 RETURN END C C C C C C ******************************* C * * C * END OF TEST PROGRAM * C * * C * BEGINNING OF VARP2 ROUTINES * C * * C ******************************* C C C C ============================================================== SUBROUTINE VARP2 (S, L, LMAX, NL, N, NMAX, NDIM, LPPS1, IV, T, X Y, W, ADA, A, IPRINT, ALF, BETA, IERR) C ============================================================== C C GIVEN S SETS OF N OBSERVATIONS EACH, Y(1,1), ..., Y(N,S), OF A C DEPENDENT VARIABLE Y, WHERE Y(I,J) CORRESPONDS TO THE IV C INDEPENDENT VARIABLE(S) T(I,1), T(I,2), ..., T(I,IV), VARP2 C ATTEMPTS TO COMPUTE A WEIGHTED LEAST SQUARES FIT TO A FUNCTION C ETA (THE 'MODEL') WHICH IS A LINEAR COMBINATION C C L C ETA (ALF,BETA; T) = SUM BETA * PHI (ALF;T) + PHI (ALF;T) C K J=1 J,K J L+1 C C OF NONLINEAR FUNCTIONS PHI(J) (E.G., A SUM OF EXPONENTIALS AND/ C OR GAUSSIANS). THAT IS, DETERMINE THE LINEAR PARAMETERS C BETA(J,K) FOR J=1,2,...,L, K=1,2,...,S, AND THE VECTOR OF C NONLINEAR PARAMETERS ALF BY MINIMIZING THE FROBENIUS NORM OF C THE MATRIX OF RESIDUALS: C C 2 S N 2 C NORM(RESIDUAL) = SUM SUM W *(Y -ETA (ALF,BETA;T )) . C K=1 I=1 I I,K K I C C THE (L+1)-ST TERM IS OPTIONAL, AND IS USED WHEN IT IS DESIRED C TO FIX ONE OR MORE OF THE BETA'S (RATHER THAN LET THEM BE C DETERMINED). VARP2 REQUIRES FIRST DERIVATIVES OF THE PHI'S. C C NOTES: C C A) FOR S=1, THE PROBLEM IS A NONLINEAR LEAST SQUARES PROBLEM C OF THE TYPE HANDLED IN REFERENCE 1. THE CASE WITH S>1 SOLVES C A SIMILAR PROBLEM WHITH MULTIPLE RIGHT HAND SIDES, EACH ALLOWED C TO HAVE DIFFERENT LINEAR COEFFICIENTS, BUT CONSTRAINED TO HAVE C THE SAME NONLINEAR PARAMETERS. SEE REFERENCE 8. C C B) THE ORIGINAL PROGRAM VARPRO (OF WHICH THIS IS A MODIFI- C CATION) IS AVAILABLE FOR THE SPECIAL CASE S=1. FOR THAT CASE C VARPRO IS EASIER TO USE AND HAS THE ADDED ADVANTAGE THAT IT C RETURNS THE COVARIANCE MATRIX OF THE PARAMETERS AND THE ESTI- C MATED VARIANCE OF THE OBSERVATIONS. C C C) AN ETA OF THE ABOVE FORM IS CALLED 'SEPARABLE'. THE C CASE OF A NONSEPARABLE ETA CAN BE HANDLED BY SETTING L = 0 C AND USING PHI(L+1). C C D) VARP2 MAY ALSO BE USED TO SOLVE LINEAR LEAST SQUARES C PROBLEMS (IN THAT CASE NO ITERATIONS ARE PERFORMED). SET C NL = 0. C C E) THE MAIN ADVANTAGE OF VARP2 OVER OTHER LEAST SQUARES C PROGRAMS IS THAT NO INITIAL GUESSES ARE NEEDED FOR THE LINEAR C PARAMETERS. NOT ONLY DOES THIS MAKE IT EASIER TO USE, BUT IT C OFTEN LEADS TO FASTER CONVERGENCE. C C C DESCRIPTION OF PARAMETERS C C S NUMBER OF RIGHT HAND SIDES. C L NUMBER OF LINEAR PARAMETERS BETA FOR EACH RIGHT C SIDE (MUST BE .GE. 0). C LMAX THE DECLARED ROW DIMENSION OF THE MATRIX BETA. C NL NUMBER OF NONLINEAR PARAMETERS ALF (MUST BE .GE. 0). C N NUMBER OF OBSERVATIONS FOR EACH RIGHT HAND SIDE. C S*N MUST BE GREATER THAN S*L + NL C (I.E., THE NUMBER OF OBSERVATIONS MUST EXCEED THE C NUMBER OF PARAMETERS). C IV NUMBER OF INDEPENDENT VARIABLES T. C T REAL N BY IV MATRIX OF INDEPENDENT VARIABLES. T(I, J) C CONTAINS THE VALUE OF THE I-TH OBSERVATION OF THE J-TH C INDEPENDENT VARIABLE. C Y N BY S MATRIX OF OBSERVATIONS CORRESPONDING TO THE S C RIGHT HAND SIDES, EACH OF WHICH HAS N VALUES, ONE C FOR EACH ROW OF T. C W N-VECTOR OF NONNEGATIVE WEIGHTS. W(I) IS THE WEIGHT C OF THE I'TH OBSERVATION FOR ALL OF THE S RIGHT HAND C SIDES. THERE IS CURRENTLY NO PROVISION FOR GIVING C DIFFERENT WEIGHTS FOR EACH RHS. W SHOULD BE SET TO C 1'S IF WEIGHTS ARE NOT DESIRED. IF VARIANCES OF THE C INDIVIDUAL OBSERVATIONS ARE KNOWN, W(I) SHOULD BE SET C TO 1./VARIANCE(I). C INC NL X (L+1) INTEGER INCIDENCE MATRIX. INC(K, J) = 1 IF C NON-LINEAR PARAMETER ALF(K) APPEARS IN THE J-TH C FUNCTION PHI(J). (THE PROGRAM SETS ALL OTHER INC(K, J) C TO ZERO.) IF PHI(L+1) IS INCLUDED IN THE MODEL, C THE APPROPRIATE ELEMENTS OF THE (L+1)-ST COLUMN SHOULD C BE SET TO 1'S. INC IS NOT NEEDED WHEN L = 0 OR NL = 0. C CAUTION: THE DECLARED ROW DIMENSION OF INC (IN ADA) C MUST CURRENTLY BE SET TO 12. SEE 'RESTRICTIONS' BELOW. C NMAX THE DECLARED ROW DIMENSION OF THE MATRICES Y AND T. C IT MUST BE AT LEAST N. C NDIM THE DECLARED ROW DIMENSION OF THE MATRIX A. IT MUST C BE AT LEAST MAX(N, 2*NL+3, S*N - (S-1)*L). C LPPS1 L+P+S+1, WHERE P IS THE NUMBER OF ONES IN THE MATRIX C INC. THE DECLARED COLUMN DIMENSION OF A MUST BE AT C LEAST LPPS1. (IF L = 0, SET LPPS1 = NL+S+1. IF NL = C 0, SET LPPS1 = L+S+1.) C A REAL MATRIX OF SIZE NDIM BY LPPS1. C IPRINT INPUT INTEGER CONTROLLING PRINTED OUTPUT. IF IPRINT IS C POSITIVE, THE NONLINEAR PARAMETERS, THE NORM OF THE C RESIDUAL, AND THE MARQUARDT PARAMETER WILL BE OUTPUT C EVERY IPRINT-TH ITERATION (AND INITIALLY, AND AT THE C FINAL ITERATION). THE LINEAR PARAMETERS WILL BE C PRINTED AT THE FINAL ITERATION. ANY ERROR MESSAGES C WILL ALSO BE PRINTED. (IPRINT = 1 IS RECOMMENDED AT C FIRST.) IF IPRINT = 0, ONLY THE FINAL QUANTITIES WILL C BE PRINTED, AS WELL AS ANY ERROR MESSAGES. IF IPRINT = C -1, NO PRINTING WILL BE DONE. THE USER IS THEN C RESPONSIBLE FOR CHECKING THE PARAMETER IERR FOR ERRORS. C ALF NL-VECTOR OF ESTIMATES OF NONLINEAR PARAMETERS C (INPUT). ON OUTPUT IT WILL CONTAIN OPTIMAL VALUES OF C THE NONLINEAR PARAMETERS. C BETA L BY S MATRIX OF LINEAR PARAMETERS WITH DECLARED C ROW DIMENSION LMAX. C IERR INTEGER ERROR FLAG (OUTPUT): C .GT. 0 - SUCCESSFUL CONVERGENCE, IERR IS THE NUMBER OF C ITERATIONS TAKEN. C -1 TERMINATED FOR TOO MANY ITERATIONS. C -2 TERMINATED FOR ILL-CONDITIONING (MARQUARDT C PARAMETER TOO LARGE.) ALSO SEE IERR = -8 BELOW. C -4 INPUT ERROR IN PARAMETER N, L, NL, LPPS1, OR NMAX. C -5 INC MATRIX IMPROPERLY SPECIFIED, OR P DISAGREES C WITH LPPS1. C -6 A WEIGHT WAS NEGATIVE. C -7 'CONSTANT' COLUMN WAS COMPUTED MORE THAN ONCE. C -8 CATASTROPHIC FAILURE - A COLUMN OF THE A MATRIX HAS C BECOME ZERO. SEE 'CONVERGENCE FAILURES' BELOW. C C C SUBROUTINES REQUIRED C C NINE SUBROUTINES, DPA, ORFAC1, ORFAC2, BACSUB, POSTPR, COV, C XNORM, INIT, AND VARERR ARE PROVIDED. IN ADDITION, THE USER C MUST PROVIDE A SUBROUTINE (CORRESPONDING TO THE ARGUMENT ADA) C WHICH, GIVEN ALF, WILL EVALUATE THE FUNCTIONS PHI(J) AND THEIR C PARTIAL DERIVATIVES D PHI(J)/D ALF(K), AT THE SAMPLE POINTS C T(I). THIS ROUTINE MUST BE DECLARED 'EXTERNAL' IN THE CALLING C PROGRAM. ITS CALLING SEQUENCE IS C C SUBROUTINE ADA (S, L+1, NL, N, NMAX, NDIM, LPPS1, IV, A, C INC, T, ALF, ISEL) C C THE USER SHOULD MODIFY THE EXAMPLE SUBROUTINE 'ADA' (GIVEN C ELSEWHERE) FOR HIS OWN FUNCTIONS. C C THE VECTOR SAMPLED FUNCTIONS PHI(J) SHOULD BE STORED IN THE C FIRST N ROWS AND FIRST L+1 COLUMNS OF THE MATRIX A, I.E., C A(I, J) SHOULD CONTAIN PHI(J, ALF; T(I,1), T(I,2), ..., C T(I,IV)), I = 1, ..., N; J = 1, ..., L (OR L+1). THE (L+1)-ST C COLUMN OF A CONTAINS PHI(L+1) IF PHI(L+1) IS IN THE MODEL, C OTHERWISE IT IS RESERVED FOR WORKSPACE. IF S>1, COLUMNS C L+2 THROUGH L+S ARE ALSO RESERVED. THE 'CONSTANT' FUNCTIONS C (THESE ARE FUNCTIONS PHI(J) WHICH DO NOT DEPEND UPON ANY C NONLINEAR PARAMETERS ALF, E.G., T(I)**J) (IF ANY) MUST APPEAR C FIRST, STARTING IN COLUMN 1. THE COLUMN N-VECTORS OF NONZERO C PARTIAL DERIVATIVES D PHI(J) / D ALF(K) SHOULD BE STORED C SEQUENTIALLY IN THE MATRIX A IN COLUMNS L+S+1 THROUGH L+S+P. C THE ORDER IS C C D PHI(1) D PHI(2) D PHI(L) D PHI(L+1) D PHI(1) C --------, --------, ..., --------, ----------, --------, C D ALF(1) D ALF(1) D ALF(1) D ALF(1) D ALF(2) C C D PHI(2) D PHI(L+1) D PHI(1) D PHI(L+1) C --------, ..., ----------, ..., ---------, ..., ----------, C D ALF(2) D ALF(2) D ALF(NL) D ALF(NL) C C OMITTING COLUMNS OF DERIVATIVES WHICH ARE ZERO, AND OMITTING C PHI(L+1) COLUMNS IF PHI(L+1) IS NOT IN THE MODEL. NOTE THAT C THE LINEAR PARAMETERS BETA ARE NOT USED IN THE MATRIX A. C COLUMN L+P+S+1 IS RESERVED FOR WORKSPACE. C C THE CODING OF ADA SHOULD BE ARRANGED SO THAT: C C ISEL = 1 (WHICH OCCURS THE FIRST TIME ADA IS CALLED) MEANS: C A. FILL IN THE INCIDENCE MATRIX INC C B. STORE ANY CONSTANT PHI'S IN A. C C. COMPUTE NONCONSTANT PHI'S AND PARTIAL DERIVA- C TIVES. C = 2 MEANS COMPUTE ONLY THE NONCONSTANT FUNCTIONS PHI C = 3 MEANS COMPUTE ONLY THE DERIVATIVES C C (WHEN THE PROBLEM IS LINEAR (NL = 0) ONLY ISEL = 1 IS USED, AND C DERIVATIVES ARE NOT NEEDED.) C C RESTRICTIONS C C THE SUBROUTINES DPA, INIT (AND ADA) CONTAIN THE LOCALLY C DIMENSIONED MATRIX INC, WHOSE DIMENSIONS ARE CURRENTLY SET FOR C MAXIMA OF L+1 = 8, NL = 12. THEY MUST BE CHANGED FOR LARGER C PROBLEMS. DATA PLACED IN ARRAY A IS OVERWRITTEN ('DESTROYED'). C DATA PLACED IN ARRAYS T, Y AND INC IS LEFT INTACT. THE PROGRAM C RUNS IN WATFIV, EXCEPT WHEN L = 0 OR NL = 0. C C IT IS ASSUMED THAT THE MATRIX PHI(J, ALF; T(I)) HAS FULL C COLUMN RANK. THIS MEANS THAT THE FIRST L COLUMNS OF THE MATRIX C A MUST BE LINEARLY INDEPENDENT. C C OPTIONAL NOTE: AS WILL BE NOTED FROM THE SAMPLE SUBPROGRAM C ADA, THE DERIVATIVES D PHI(J)/D ALF(K) (ISEL = 3) MUST BE C COMPUTED INDEPENDENTLY OF THE FUNCTIONS PHI(J) (ISEL = 2), C SINCE THE FUNCTION VALUES ARE OVERWRITTEN AFTER ADA IS CALLED C WITH ISEL = 2. THIS IS DONE TO MINIMIZE STORAGE, AT THE POS- C SIBLE EXPENSE OF SOME RECOMPUTATION (SINCE THE FUNCTIONS AND C DERIVATIVES FREQUENTLY HAVE SOME COMMON SUBEXPRESSIONS). TO C REDUCE THE AMOUNT OF COMPUTATION AT THE EXPENSE OF SOME C STORAGE, CREATE A MATRIX B OF DIMENSION NMAX BY L+1 IN ADA, AND C AFTER THE COMPUTATION OF THE PHI'S (ISEL = 2), COPY THE VALUES C INTO B. THESE VALUES CAN THEN BE USED TO CALCULATE THE DERIV- C ATIVES (ISEL = 3). (THIS MAKES USE OF THE FACT THAT WHEN A C CALL TO ADA WITH ISEL = 3 FOLLOWS A CALL WITH ISEL = 2, THE C ALFS ARE THE SAME.) C C TO CONVERT TO OTHER MACHINES, CHANGE THE OUTPUT UNIT IN THE C DATA STATEMENTS IN VARP2, DPA, POSTPR, AND VARERR. THE C PROGRAM HAS BEEN CHECKED FOR PORTABILITY BY THE BELL LABS PFORT C VERIFIER. FOR MACHINES WITHOUT DOUBLE PRECISION HARDWARE, IT C MAY BE DESIRABLE TO CONVERT TO SINGLE PRECISION. THIS CAN BE C DONE BY CHANGING (A) THE DECLARATIONS 'DOUBLE PRECISION' TO C 'REAL', (B) THE PATTERN '.D' TO '.E' IN THE 'DATA' STATEMENT IN C VARP2, (C) DSIGN, DSQRT AND DABS TO SIGN, SQRT AND ABS, C RESPECTIVELY, AND (D) DEXP TO EXP IN THE SAMPLE PROGRAMS ONLY. C C CONVERGENCE FAILURES C C IF CONVERGENCE FAILURES OCCUR, FIRST CHECK FOR INCORRECT C CODING OF THE SUBROUTINE ADA. CHECK ESPECIALLY THE ACTION OF C ISEL, AND THE COMPUTATION OF THE PARTIAL DERIVATIVES. IF THESE C ARE CORRECT, TRY SEVERAL STARTING GUESSES FOR ALF. IF ADA C IS CODED CORRECTLY, AND IF ERROR RETURNS IERR = -2 OR -8 C PERSISTENTLY OCCUR, THIS IS A SIGN OF ILL-CONDITIONING, WHICH C MAY BE CAUSED BY SEVERAL THINGS. ONE IS POOR SCALING OF THE C PARAMETERS; ANOTHER IS AN UNFORTUNATE INITIAL GUESS FOR THE C PARAMETERS, STILL ANOTHER IS A POOR CHOICE OF THE MODEL. C C ALGORITHM C C THE RESIDUAL R IS MODIFIED TO INCORPORATE, FOR ANY FIXED C ALF, THE OPTIMAL LINEAR PARAMETERS FOR THAT ALF. IT IS THEN C POSSIBLE TO MINIMIZE ONLY ON THE NONLINEAR PARAMETERS. AFTER C THE OPTIMAL VALUES OF THE NONLINEAR PARAMETERS HAVE BEEN DETER- C MINED, THE LINEAR PARAMETERS CAN BE RECOVERED BY LINEAR LEAST C SQUARES TECHNIQUES (SEE REF. 1). C C THE MINIMIZATION IS BY A MODIFICATION OF OSBORNE'S (REF. 3) C MODIFICATION OF THE LEVENBERG-MARQUARDT ALGORITHM. INSTEAD OF C SOLVING THE NORMAL EQUATIONS WITH MATRIX C C T 2 C (J J + NU * D), WHERE J = D(ETA)/D(ALF), C C STABLE ORTHOGONAL (HOUSEHOLDER) REFLECTIONS ARE USED ON A C MODIFICATION OF THE MATRIX C ( J ) C (------) , C ( NU*D ) C C WHERE D IS A DIAGONAL MATRIX CONSISTING OF THE LENGTHS OF THE C COLUMNS OF J. THIS MARQUARDT STABILIZATION ALLOWS THE ROUTINE C TO RECOVER FROM SOME RANK DEFICIENCIES IN THE JACOBIAN. C OSBORNE'S EMPIRICAL STRATEGY FOR CHOOSING THE MARQUARDT PARAM- C ETER HAS PROVEN REASONABLY SUCCESSFUL IN PRACTICE. (GAUSS- C NEWTON WITH STEP CONTROL CAN BE OBTAINED BY MAKING THE CHANGE C INDICATED BEFORE THE INSTRUCTION LABELED 5). A DESCRIPTION CAN C BE FOUND IN REF. (3), AND A FLOW CHART IN (2), P. 22. C C FOR REFERENCE, SEE C C 1. GENE H. GOLUB AND V. PEREYRA, 'THE DIFFERENTIATION OF C PSEUDO-INVERSES AND NONLINEAR LEAST SQUARES PROBLEMS WHOSE C VARIABLES SEPARATE,' SIAM J. NUMER. ANAL. 10, 413-432 C (1973). C 2. ------, SAME TITLE, STANFORD C.S. REPORT 72-261, FEB. 1972. C 3. OSBORNE, MICHAEL R., 'SOME ASPECTS OF NON-LINEAR LEAST C SQUARES CALCULATIONS,' IN LOOTSMA, ED., 'NUMERICAL METHODS C FOR NON-LINEAR OPTIMIZATION,' ACADEMIC PRESS, LONDON, 1972. C 4. KROGH, FRED, 'EFFICIENT IMPLEMENTATION OF A VARIABLE PRO- C JECTION ALGORITHM FOR NONLINEAR LEAST SQUARES PROBLEMS,' C COMM. ACM 17, PP. 167-169 (MARCH, 1974). C 5. KAUFMAN, LINDA, 'A VARIABLE PROJECTION METHOD FOR SOLVING C SEPARABLE NONLINEAR LEAST SQUARES PROBLEMS', B.I.T. 15, C 49-57 (1975). C 6. DRAPER, N., AND SMITH, H., APPLIED REGRESSION ANALYSIS, C WILEY, N.Y., 1966 (FOR STATISTICAL INFORMATION ONLY). C 7. C. LAWSON AND R. HANSON, SOLVING LEAST SQUARES PROBLEMS, C PRENTICE-HALL, ENGLEWOOD CLIFFS, N. J., 1974. C 8. GOLUB, G. AND LEVEQUE, R., EXTENSIONS AND USES OF THE C VARIABLE PROJECTION ALGORITHM FOR SOLVING NONLINEAR LEAST C SQUARES PROBLEMS, PROC. 1979 ARMY NUM. ANAL. AND COMPUTERS C CONF., ARO REPORT 79-3, PP. 1-12. C C VICTOR PEREYRA C ESCUELA DE COMPUTACION C FACULTAD DE CIENCIAS C UNIVERSIDAD CENTRAL DE VENEZUELA C CARACAS, VENEZUELA C C JOHN BOLSTAD C COMPUTER SCIENCE DEPT., SERRA HOUSE C STANFORD UNIVERSITY C JANUARY, 1977 C C RANDY LEVEQUE C COMPUTER SCIENCE DEPT., SERRA HOUSE C STANFORD UNIVERSITY C DECEMBER, 1978 C C .................................................................. C INTEGER B1, OUTPUT, S DOUBLE PRECISION A(NDIM, LPPS1), BETA(LMAX,S), ALF(NL), 2 W(N), Y(NMAX,S), ACUM, EPS1, GNSTEP, NU, PRJRES, R, RNEW, XNORM, 3 T(NMAX,IV) LOGICAL SKIP EXTERNAL ADA DATA EPS1 /1.D-6/, ITMAX /50/, OUTPUT /6/ C C THE FOLLOWING TWO PARAMETERS ARE USED IN THE CONVERGENCE C TEST: EPS1 IS AN ABSOLUTE AND RELATIVE TOLERANCE FOR THE C NORM OF THE PROJECTION OF THE RESIDUAL ONTO THE RANGE OF THE C JACOBIAN OF THE VARIABLE PROJECTION FUNCTIONAL. C ITMAX IS THE MAXIMUM NUMBER OF FUNCTION AND DERIVATIVE C EVALUATIONS ALLOWED. CAUTION: EPS1 MUST NOT BE C SET SMALLER THAN 10 TIMES THE UNIT ROUND-OFF OF THE MACHINE. C----------------------------------------------------------------- IERR = 1 ITER = 0 LP1 = L + 1 LPS = L+S B1 = L + S + 1 LNLS1 = L + NL + S + 1 NLP1 = NL + 1 SKIP = .FALSE. MODIT = IPRINT IF (IPRINT .LE. 0) MODIT = ITMAX + 2 NU = 0. C IF GAUSS-NEWTON IS DESIRED REMOVE THE NEXT STATEMENT. NU = 1. C C BEGIN OUTER ITERATION LOOP TO UPDATE ALF. C CALCULATE THE NORM OF THE RESIDUAL AND THE DERIVATIVE OF C THE MODIFIED RESIDUAL THE FIRST TIME, BUT ONLY THE C DERIVATIVE IN SUBSEQUENT ITERATIONS. C 5 CALL DPA (S, L, LMAX, NL, N, NMAX, NDIM, LPPS1, IV, T, Y, W, ALF, X ADA, IERR, IPRINT, A, BETA, A(1, LP1), R) GNSTEP = 1.0 ITERIN = 0 IF (ITER .GT. 0) GO TO 10 IF (NL .EQ. 0) GO TO 90 IF (IERR .NE. 1) GO TO 99 C IF (IPRINT .LE. 0) GO TO 10 WRITE (OUTPUT, 207) ITERIN, R WRITE (OUTPUT, 200) NU C BEGIN TWO-STAGE ORTHOGONAL FACTORIZATION 10 CALL ORFAC1(S, NLP1, NDIM, N, L, IPRINT, A(1, B1), PRJRES, IERR) IF (IERR .LT. 0) GO TO 99 IERR = 2 IF (NU .EQ. 0.) GO TO 30 C C BEGIN INNER ITERATION LOOP FOR GENERATING NEW ALF AND C TESTING IT FOR ACCEPTANCE. C 25 CALL ORFAC2(NLP1, NDIM, NU, A(1, B1)) C C SOLVE A NL X NL UPPER TRIANGULAR SYSTEM FOR DELTA-ALF. C THE TRANSFORMED RESIDUAL (IN COL. LNLS1 OF A) IS OVER- C WRITTEN BY THE RESULT DELTA-ALF. C 30 CALL BACSUB (NDIM, NL, A(1, B1), A(1, LNLS1)) DO 35 K = 1, NL 35 A(K, B1) = ALF(K) + A(K, LNLS1) C NEW ALF(K) = ALF(K) + DELTA ALF(K) C C STEP TO THE NEW POINT NEW ALF, AND COMPUTE THE NEW C NORM OF RESIDUAL. NEW ALF IS STORED IN COLUMN B1 OF A. C 40 CALL DPA (S, L, LMAX, NL, N, NMAX,NDIM, LPPS1, IV, T, Y, W, X A(1,B1), ADA, IERR, IPRINT, A, BETA, A(1, LP1), RNEW) IF (IERR .NE. 2) GO TO 99 ITER = ITER + 1 ITERIN = ITERIN + 1 SKIP = MOD(ITER, MODIT) .NE. 0 IF (SKIP) GO TO 45 WRITE (OUTPUT, 203) ITER WRITE (OUTPUT, 216) (A(K, B1), K = 1, NL) WRITE (OUTPUT, 207) ITERIN, RNEW C 45 IF (ITER .LT. ITMAX) GO TO 50 IERR = -1 CALL VARERR (IPRINT, IERR, 1) GO TO 95 50 IF (RNEW - R .LT. EPS1*(R + 1.D0)) GO TO 75 C C RETRACT THE STEP JUST TAKEN C IF (NU .NE. 0.) GO TO 60 C GAUSS-NEWTON OPTION ONLY GNSTEP = 0.5*GNSTEP IF (GNSTEP .LT. EPS1) GO TO 95 DO 55 K = 1, NL 55 A(K, B1) = ALF(K) + GNSTEP*A(K, LNLS1) GO TO 40 C ENLARGE THE MARQUARDT PARAMETER 60 NU = 1.5*NU IF (.NOT. SKIP) WRITE (OUTPUT, 206) NU IF (NU .LE. 100.) GO TO 65 IERR = -2 CALL VARERR (IPRINT, IERR, 1) GO TO 95 C RETRIEVE UPPER TRIANGULAR FORM C AND RESIDUAL OF FIRST STAGE. 65 DO 70 K = 1, NL KSUB = LPS + K DO 70 J = K, NLP1 JSUB = LPS + J ISUB = NLP1 + J 70 A(K, JSUB) = A(ISUB, KSUB) GO TO 25 C END OF INNER ITERATION LOOP C ACCEPT THE STEP JUST TAKEN C 75 R = RNEW DO 80 K = 1, NL 80 ALF(K) = A(K, B1) C CALC. NORM(DELTA ALF)/NORM(ALF) ACUM = GNSTEP*XNORM(NL, A(1, LNLS1))/XNORM(NL, ALF) C C IF ITERIN IS GREATER THAN 1, A STEP WAS RETRACTED DURING C THIS OUTER ITERATION. C IF (ITERIN .EQ. 1) NU = 0.5*NU IF (SKIP) GO TO 85 WRITE (OUTPUT, 200) NU WRITE (OUTPUT, 208) ACUM 85 IERR = 3 IF (PRJRES .GT. EPS1*(R + 1.D0)) GO TO 5 C END OF OUTER ITERATION LOOP C C CALCULATE FINAL QUANTITIES -- LINEAR PARAMETERS, RESIDUALS, C 90 IERR = ITER 95 IF (NL .GT. 0) CALL DPA(S, L, LMAX, NL, N, NMAX, NDIM, LPPS1, IV, X T, Y, W, ALF, ADA, 3, IPRINT, A, BETA, A(1, LP1), R) CALL POSTPR(S, L, LMAX, NL, N, NMAX, NDIM, LNLS1, EPS1, R, IPRINT, X ALF, W, A, A(1, LP1), BETA, IERR) 99 RETURN C 200 FORMAT (9H NU =, E15.7) 203 FORMAT (12H0 ITERATION, I4, 24H NONLINEAR PARAMETERS) 206 FORMAT (25H STEP RETRACTED, NU =, E15.7) 207 FORMAT (1H0, I5, 20H NORM OF RESIDUAL =, E15.7) 208 FORMAT (34H NORM(DELTA-ALF) / NORM(ALF) =, E12.3) 216 FORMAT (1H0, 7E15.7) END C C ============================================================== SUBROUTINE ORFAC1(S, NLP1, NDIM, N, L, IPRINT, B, PRJRES, IERR) C ============================================================== C C STAGE 1: HOUSEHOLDER REDUCTION OF C C ( . ) ( DR'. R3 ) NL C ( DR . R2 ) TO (----. -- ), C ( . ) ( 0 . R4 ) N-L-NL C C NL 1 NL 1 C C WHERE DR = -D(Q2)*Y IS THE DERIVATIVE OF THE MODIFIED RESIDUAL C PRODUCED BY DPA, R2 IS THE TRANSFORMED RESIDUAL FROM DPA, AND C DR' IS IN UPPER TRIANGULAR FORM (AS IN REF. (2), P. 18). C DR IS STORED IN ROWS L+1 TO N*S-L*(S-1) AND COLUMNS L+S+1 TO C L+S+NL OF THE MATRIX A (I.E., COLUMNS 1 TO NL OF THE MATRIX B). C R2 IS STORED IN COLUMN L+NL+S+1 OF THE MATRIX A (COLUMN NL + 1 C OF B). FOR K = 1, 2, ..., NL, FIND REFLECTION I - U * U' / C BETA WHICH ZEROES B(I, K), I = L+K+1, ..., N*S - L*(S-1). C C .................................................................. C INTEGER S DOUBLE PRECISION ACUM, ALPHA, B(NDIM, NLP1), BETA, DSIGN, PRJRES, X U, XNORM C NL = NLP1 - 1 NSLS1 = N*S - L*(S-1) NL23 = 2*NL + 3 LP1 = L + 1 C DO 30 K = 1, NL LPK = L + K ALPHA = DSIGN(XNORM(NSLS1+1-LPK, B(LPK, K)), B(LPK, K)) U = B(LPK, K) + ALPHA B(LPK, K) = U BETA = ALPHA * U IF (ALPHA .NE. 0.0) GO TO 13 C COLUMN WAS ZERO IERR = -8 CALL VARERR (IPRINT, IERR, LP1 + K) GO TO 99 C APPLY REFLECTIONS TO REMAINING COLUMNS C OF B AND TO RESIDUAL VECTOR. 13 KP1 = K + 1 DO 25 J = KP1, NLP1 ACUM = 0.0 DO 20 I = LPK, NSLS1 20 ACUM = ACUM + B(I, K) * B(I, J) ACUM = ACUM / BETA DO 25 I = LPK, NSLS1 25 B(I, J) = B(I, J) - B(I, K) * ACUM 30 B(LPK, K) = -ALPHA C PRJRES = XNORM(NL, B(LP1, NLP1)) C C SAVE UPPER TRIANGULAR FORM AND TRANSFORMED RESIDUAL, FOR USE C IN CASE A STEP IS RETRACTED. ALSO COMPUTE COLUMN LENGTHS. C IF (IERR .EQ. 4) GO TO 99 DO 50 K = 1, NL LPK = L + K DO 40 J = K, NLP1 JSUB = NLP1 + J B(K, J) = B(LPK, J) 40 B(JSUB, K) = B(LPK, J) 50 B(NL23, K) = XNORM(K, B(LP1, K)) C 99 RETURN END C C ============================================================== SUBROUTINE ORFAC2(NLP1, NDIM, NU, B) C ============================================================== C C STAGE 2: SPECIAL HOUSEHOLDER REDUCTION OF C C NL ( DR' . R3 ) (DR'' . R5 ) C (-----. -- ) (-----. -- ) C N-L-NL ( 0 . R4 ) TO ( 0 . R4 ) C (-----. -- ) (-----. -- ) C NL (NU*D . 0 ) ( 0 . R6 ) C C NL 1 NL 1 C C WHERE DR', R3, AND R4 ARE AS IN ORFAC1, NU IS THE MARQUARDT C PARAMETER, D IS A DIAGONAL MATRIX CONSISTING OF THE LENGTHS OF C THE COLUMNS OF DR', AND DR'' IS IN UPPER TRIANGULAR FORM. C DETAILS IN (1), PP. 423-424. NOTE THAT THE (N-L-NL) BAND OF C ZEROES, AND R4, ARE OMITTED IN STORAGE. C C .................................................................. C DOUBLE PRECISION ACUM, ALPHA, B(NDIM, NLP1), BETA, DSIGN, NU, U, X XNORM C NL = NLP1 - 1 NL2 = 2*NL NL23 = NL2 + 3 DO 30 K = 1, NL KP1 = K + 1 NLPK = NL + K NLPKM1 = NLPK - 1 B(NLPK, K) = NU * B(NL23, K) B(NL, K) = B(K, K) ALPHA = DSIGN(XNORM(K+1, B(NL, K)), B(K, K)) U = B(K, K) + ALPHA BETA = ALPHA * U B(K, K) = -ALPHA C THE K-TH REFLECTION MODIFIES ONLY ROWS K, C NL+1, NL+2, ..., NL+K, AND COLUMNS K TO NL+1. DO 30 J = KP1, NLP1 B(NLPK, J) = 0. ACUM = U * B(K,J) DO 20 I = NLP1, NLPKM1 20 ACUM = ACUM + B(I,K) * B(I,J) ACUM = ACUM / BETA B(K,J) = B(K,J) - U * ACUM DO 30 I = NLP1, NLPK 30 B(I,J) = B(I,J) - B(I,K) * ACUM C RETURN END C C ============================================================== SUBROUTINE DPA (S, L, LMAX, NL, N, NMAX, NDIM, LPPS1, IV, T, Y, W, X ALF, ADA, ISEL, IPRINT, A, U, R, RNORM) C ============================================================== C C COMPUTE THE NORM OF THE RESIDUAL (IF ISEL = 1 OR 2), OR THE C (N-L) X NL X S DERIVATIVE OF THE MODIFIED RESIDUAL (N-L) BY S C MATRIX Q2*Y (IF ISEL = 1 OR 3). HERE Q * PHI = TRI, I.E., C C L ( Q1 ) ( . . ) (TRI . R1 . F1 ) C (----) ( PHI . Y . D(PHI) ) = (--- . -- . ---- ) C N-L ( Q2 ) ( . . ) ( 0 . R2 . F2 ) C C N L S P L S P C C WHERE Q IS N X N ORTHOGONAL, AND TRI IS L X L UPPER TRIANGULAR. C THE NORM OF THE RESIDUAL = FROBENIUS NORM(R2), AND THE DESIRED C DERIVATIVE ACCORDING TO REF. (5), IS C -1 C D(Q2 * Y) = -Q2 * D(PHI)* TRI * Q1* Y. C C THE THREE-TENSOR DERIVATIVE IS STORED IN COLUMNS L+S+1 THROUGH C L+S+NL AND ROWS L+1 THROUGH S*N - (S-1)*L OF THE MATRIX A. C THE MATRIX SLAB OF THE DERIVATIVE CORRESPONDING TO THE K'TH C RIGHT HAND SIDE (FOR K=1,2,...,S) IS IN ROWS L+(K-1)*(N-L)+1 C THROUGH L+K*(N-L). C C .................................................................. C INTEGER S, FIRSTC, FIRSTR, INC(12, 8) DOUBLE PRECISION A(NDIM, LPPS1), ALF(NL), T(NMAX, IV), W(N), X ACUM, ALPHA, BETA, RNORM, DSIGN, DSQRT, SAVE, R(NDIM,S), U(L), X XNORM, RI, Y(NMAX,S) LOGICAL NOWATE, PHILP1 EXTERNAL ADA C IF (ISEL .NE. 1) GO TO 3 LP1 = L + 1 LPS = L+S LNLS = L+NL+S LNLS1 = LNLS+1 LSP1 = L+S+1 LPPS = LPPS1-1 FIRSTC = 1 LASTC = LPPS FIRSTR = LP1 CALL INIT(S, L, LMAX, NL, N, NMAX, NDIM, LPPS1, IV, T, W, ALF, X ADA, ISEL, IPRINT, A, INC, NCON, NCONP1, PHILP1, NOWATE) IF (ISEL .NE. 1) GO TO 99 GO TO 30 C 3 CALL ADA (S, LP1, NL, N, NMAX, NDIM, LPPS1, IV, A, INC, T, ALF, X MIN0(ISEL,3)) IF (ISEL .EQ. 2) GO TO 6 C ISEL = 3 OR 4 FIRSTC = LSP1 LASTC = LPPS FIRSTR = (4 - ISEL)*L + 1 GO TO 50 C ISEL = 2 6 FIRSTC = NCONP1 LASTC = LPS IF (NCON .EQ. 0) GO TO 30 IF (A(1, NCON) .EQ. SAVE) GO TO 30 ISEL = -7 CALL VARERR (IPRINT, ISEL, NCON) GO TO 99 C ISEL = 1 OR 2 30 IF (PHILP1) GO TO 40 DO 35 I=1,N DO 35 J=1,S 35 R(I,J) = Y(I,J) GO TO 50 40 DO 45 I=1,N RI = R(I,1) DO 45 J=1,S 45 R(I,J) = Y(I,J) - RI C WEIGHT APPROPRIATE COLUMNS 50 IF (NOWATE) GO TO 58 DO 55 I = 1, N ACUM = W(I) DO 55 J = FIRSTC, LASTC 55 A(I, J) = A(I, J) * ACUM C C COMPUTE ORTHOGONAL FACTORIZATIONS BY HOUSEHOLDER C REFLECTIONS. IF ISEL = 1 OR 2, REDUCE PHI (STORED IN THE C FIRST L COLUMNS OF THE MATRIX A) TO UPPER TRIANGULAR FORM, C (Q*PHI = TRI), AND TRANSFORM Y (STORED IN COLUMNS L+1 C THROUGH L+S), GETTING Q*Y = R. IF ISEL = 1, ALSO TRANSFORM C J = D PHI (STORED IN COLUMNS L+S+1 THROUGH L+P+S OF THE C MATRIX A), GETTING Q*J = F. IF ISEL = 3 OR 4, PHI HAS C ALREADY BEEN REDUCED, TRANSFORM ONLY J. TRI, R, AND F C OVERWRITE PHI, Y, AND J, RESPECTIVELY, AND A FACTORED FORM C OF Q IS SAVED IN U AND THE LOWER TRIANGLE OF PHI. C 58 IF (L .EQ. 0) GO TO 75 DO 70 K = 1, L KP1 = K + 1 IF (ISEL .GE. 3 .OR. (ISEL .EQ. 2 .AND. K .LT.NCONP1)) GO TO 66 ALPHA = DSIGN(XNORM(N+1-K, A(K, K)), A(K, K)) U(K) = A(K, K) + ALPHA A(K, K) = -ALPHA FIRSTC = KP1 IF (ALPHA .NE. 0.0) GO TO 66 ISEL = -8 CALL VARERR (IPRINT, ISEL, K) GO TO 99 C APPLY REFLECTIONS TO COLUMNS C FIRSTC TO LASTC. 66 BETA = -A(K, K) * U(K) DO 70 J = FIRSTC, LASTC ACUM = U(K)*A(K, J) DO 68 I = KP1, N 68 ACUM = ACUM + A(I, K)*A(I, J) ACUM = ACUM / BETA A(K,J) = A(K,J) - U(K)*ACUM DO 70 I = KP1, N 70 A(I, J) = A(I, J) - A(I, K)*ACUM C 75 IF (ISEL .GE. 3) GO TO 85 C C COMPUTE THE FROBENIUS NORM OF THE RESIDUAL MATRIX: RNORM = 0. DO 76 J=1,S 76 RNORM = RNORM + XNORM(N-L, R(LP1,J))**2 RNORM = DSQRT(RNORM) C IF (ISEL .EQ. 2) GO TO 99 IF (NCON .GT. 0) SAVE = A(1, NCON) C C F2 IS NOW CONTAINED IN ROWS L+1 TO N AND COLUMNS L+S+1 TO C L+P+S OF THE MATRIX A. NOW SOLVE THE S (L X L) UPPER C TRIANGULAR SYSTEMS TRI*BETA(J) = R1(J) FOR THE LINEAR C PARAMETERS BETA. BETA OVERWRITES R1. C 85 IF (L .EQ. 0) GO TO 87 DO 86 J=1,S 86 CALL BACSUB(NDIM,L,A,R(1,J)) C C MAJOR PART OF KAUFMAN'S SIMPLIFICATION OCCURS HERE. COMPUTE C THE DERIVATIVE OF ETA WITH RESPECT TO THE NONLINEAR C PARAMETERS C C T D ETA T L D PHI(J) D PHI(L+1) C Q * -------- = Q * (SUM BETA(J) -------- + ----------) = F2*BETA C D ALF(K) J=1 D ALF(K) D ALF(K) C C AND STORE THE RESULT IN COLUMNS L+S+1 TO L+NL+S. THE C FIRST L ROWS ARE OMITTED. THIS IS -D(Q2)*Y. THE RESIDUAL C R2 = Q2*Y (IN COLUMNS L+1 TO L+S) IS COPIED TO COLUMN C L+NL+S+1. C 87 DO 95 I = FIRSTR, N DO 95 ISBACK=1,S IS = S - ISBACK + 1 ISUB = (N-L) * (IS-1) + I IF (L .EQ. NCON) GO TO 95 M = LPS DO 90 K=1,NL ACUM = 0. DO 88 J=NCONP1,L IF (INC(K,J) .EQ. 0) GO TO 88 M = M+1 ACUM = ACUM + A(I,M) * R(J,IS) 88 CONTINUE KSUB = LPS+K IF (INC(K,LP1) .EQ. 0) GO TO 90 M = M+1 ACUM = ACUM + A(I,M) 90 A(ISUB,KSUB) = ACUM 95 A(ISUB,LNLS1) = R(I,IS) C 99 RETURN END C C ============================================================== SUBROUTINE INIT(S, L, LMAX, NL, N, NMAX, NDIM, LPPS1, IV, T, W, X ALF, ADA, ISEL, IPRINT, A, INC, NCON, NCONP1, PHILP1, NOWATE) C ============================================================== C C CHECK VALIDITY OF INPUT PARAMETERS, AND DETERMINE NUMBER OF C CONSTANT FUNCTIONS. C C .................................................................. C DOUBLE PRECISION A(NDIM, LPPS1), ALF(NL), T(NMAX, IV), W(N), X DSQRT INTEGER OUTPUT, P, INC(12, 8), S LOGICAL NOWATE, PHILP1 DATA OUTPUT /6/ C LP1 = L + 1 LNLS1 = L + S + NL + 1 C CHECK FOR VALID INPUT IF (L .GE. 0 .AND. NL .GE. 0 .AND. S*L+NL .LT. S*N .AND. LNLS1.LE. X LPPS1 .AND. 2*NL + 3 .LE. NDIM .AND. N .LE. NMAX .AND. N .LE. X NDIM .AND. IV .GT. 0 .AND. .NOT. (NL .EQ. 0 .AND. L .EQ. 0) .AND. X S*N-(S-1)*L .LE. NDIM .AND. S .GT. 0 .AND. L .LE. LMAX) GO TO 1 ISEL = -4 CALL VARERR (IPRINT, ISEL, 1) GO TO 99 C 1 IF (L .EQ. 0 .OR. NL .EQ. 0) GO TO 3 DO 2 J = 1, LP1 DO 2 K = 1, NL 2 INC(K, J) = 0 C 3 CALL ADA (S, LP1, NL, N, NMAX,NDIM, LPPS1, IV, A, INC, T, ALF, X ISEL) C NOWATE = .TRUE. DO 9 I = 1, N NOWATE = NOWATE .AND. (W(I) .EQ. 1.0) IF (W(I) .GE. 0.) GO TO 9 C ERROR IN WEIGHTS ISEL = -6 CALL VARERR (IPRINT, ISEL, I) GO TO 99 9 W(I) = DSQRT(W(I)) C NCON = L NCONP1 = LP1 PHILP1 = L .EQ. 0 IF (PHILP1 .OR. NL .EQ. 0) GO TO 99 C CHECK INC MATRIX FOR VALID INPUT AND C DETERMINE NUMBER OF CONSTANT FCNS. P = 0 DO 11 J = 1, LP1 IF (P .EQ. 0) NCONP1 = J DO 11 K = 1, NL INCKJ = INC(K, J) IF (INCKJ .NE. 0 .AND. INCKJ .NE. 1) GO TO 15 IF (INCKJ .EQ. 1) P = P + 1 11 CONTINUE C NCON = NCONP1 - 1 IF (IPRINT .GE. 0) WRITE (OUTPUT, 210) NCON IF (L+P+S+1 .EQ. LPPS1) GO TO 20 C INPUT ERROR IN INC MATRIX 15 ISEL = -5 CALL VARERR (IPRINT, ISEL, 1) GO TO 99 C DETERMINE IF PHI(L+1) IS IN THE MODEL. 20 DO 25 K = 1, NL 25 IF (INC(K, LP1) .EQ. 1) PHILP1 = .TRUE. C 99 RETURN 210 FORMAT (33H0 NUMBER OF CONSTANT FUNCTIONS =, I4 /) END C ============================================================== SUBROUTINE BACSUB (NDIM, N, A, X) C ============================================================== C C BACKSOLVE THE N X N UPPER TRIANGULAR SYSTEM A*X = B. C THE SOLUTION X OVERWRITES THE RIGHT SIDE B. C DOUBLE PRECISION A(NDIM, N), X(N), ACUM C X(N) = X(N) / A(N, N) IF (N .EQ. 1) GO TO 30 NP1 = N + 1 DO 20 IBACK = 2, N I = NP1 - IBACK C I = N-1, N-2, ..., 2, 1 IP1 = I + 1 ACUM = X(I) DO 10 J = IP1, N 10 ACUM = ACUM - A(I,J)*X(J) 20 X(I) = ACUM / A(I,I) C 30 RETURN END C ============================================================== SUBROUTINE POSTPR(S, L, LMAX, NL, N, NMAX, NDIM, LNLS1, EPS, X RNORM, IPRINT, ALF, W, A, R, U, IERR) C ============================================================== C C CALCULATE RESIDUALS. C ON INPUT, U CONTAINS INFORMATION ABOUT HOUSEHOLDER REFLECTIONS C FROM DPA. ON OUTPUT, IT CONTAINS THE LINEAR PARAMETERS. C INTEGER OUTPUT, S DOUBLE PRECISION A(NDIM, LNLS1), ALF(NL), R(NDIM,S), U(LMAX,S), X W(N), ACUM, EPS, PRJRES, RNORM, SAVE, DABS DATA OUTPUT /6/ C LP1 = L + 1 LPNL = LNLS1 - 2 LNL1 = LPNL + 1 DO 10 I = 1, N 10 W(I) = W(I)**2 C C UNWIND HOUSEHOLDER TRANSFORMATIONS TO GET RESIDUALS, C AND MOVE THE LINEAR PARAMETERS FROM R TO U. C IF (L .EQ. 0) GO TO 30 USAVE = L+S+2 DO 19 I=1,L 19 A(I,USAVE) = U(I,1) DO 25 IS=1,S DO 25 KBACK = 1, L K = LP1 - KBACK KP1 = K + 1 ACUM = 0. DO 20 I = KP1, N 20 ACUM = ACUM + A(I, K) * R(I,IS) U(K,IS) = R(K,IS) R(K,IS) = ACUM / A(K, K) ACUM = -ACUM / (A(K,USAVE) * A(K, K)) DO 25 I = KP1, N 25 R(I,IS) = R(I,IS) - A(I, K)*ACUM C 30 IF (IPRINT .LT. 0) GO TO 99 WRITE (OUTPUT, 209) IF (L .EQ. 0) GO TO 50 WRITE(OUTPUT,210) DO 40 I=1,L 40 WRITE(OUTPUT,212) (U(I,J), J=1,S) 50 IF (NL .GT. 0) WRITE (OUTPUT, 211) (ALF(K), K = 1, NL) WRITE(OUTPUT,214) RNORM WRITE (OUTPUT, 209) 99 RETURN C 209 FORMAT (1H0, 50(1H')) 210 FORMAT (20H0 LINEAR PARAMETERS /) 211 FORMAT (/ 23H0 NONLINEAR PARAMETERS // (7E15.7)) 212 FORMAT (7E15.7) 214 FORMAT (/ 21H0 NORM OF RESIDUAL =, E15.7) END C ============================================================== SUBROUTINE VARERR (IPRINT, IERR, K) C ============================================================== C C PRINT ERROR MESSAGES C INTEGER ERRNO, OUTPUT DATA OUTPUT /6/ C IF (IPRINT .LT. 0) GO TO 99 ERRNO = IABS(IERR) GO TO (1, 2, 99, 4, 5, 6, 7, 8), ERRNO C 1 WRITE (OUTPUT, 101) GO TO 99 2 WRITE (OUTPUT, 102) GO TO 99 4 WRITE (OUTPUT, 104) GO TO 99 5 WRITE (OUTPUT, 105) GO TO 99 6 WRITE (OUTPUT, 106) K GO TO 99 7 WRITE (OUTPUT, 107) K GO TO 99 8 WRITE (OUTPUT, 108) K C 99 RETURN 101 FORMAT (46H0 PROBLEM TERMINATED FOR EXCESSIVE ITERATIONS //) 102 FORMAT (49H0 PROBLEM TERMINATED BECAUSE OF ILL-CONDITIONING //) 104 FORMAT (/ 51H INPUT ERROR IN PARAMETER L, NL, N, LPPS1, S, LMAX,, X15H NDIM, OR NMAX. /) 105 FORMAT (69H0 ERROR -- INC MATRIX IMPROPERLY SPECIFIED, OR DISAGRE XES WITH LPPS1. /) 106 FORMAT (19H0 ERROR -- WEIGHT(, I4, 14H) IS NEGATIVE. /) 107 FORMAT (28H0 ERROR -- CONSTANT COLUMN , I3, 37H MUST BE COMPUTED XONLY WHEN ISEL = 1. /) 108 FORMAT (33H0 CATASTROPHIC FAILURE -- COLUMN , I4, 28H IS ZERO, SE XE DOCUMENTATION. /) END C ============================================================== DOUBLE PRECISION FUNCTION XNORM(N, X) C ============================================================== C C COMPUTE THE L2 (EUCLIDEAN) NORM OF A VECTOR, MAKING SURE TO C AVOID UNNECESSARY UNDERFLOWS. NO ATTEMPT IS MADE TO SUPPRESS C OVERFLOWS. C DOUBLE PRECISION X(N), RMAX, SUM, TERM, DABS, DSQRT C C FIND LARGEST (IN ABSOLUTE VALUE) ELEMENT RMAX = 0. DO 10 I = 1, N IF (DABS(X(I)) .GT. RMAX) RMAX = DABS(X(I)) 10 CONTINUE C SUM = 0. IF (RMAX .EQ. 0.) GO TO 30 DO 20 I = 1, N TERM = 0. IF (RMAX + DABS(X(I)) .NE. RMAX) TERM = X(I)/RMAX 20 SUM = SUM + TERM*TERM C 30 XNORM = RMAX*DSQRT(SUM) 99 RETURN END ********************************************************* * DATA FOR SAMPLE PROGRAM: * ********************************************************* 30 2 1 1 1 2 0.100 2.810 1.715 0.200 2.637 1.456 0.300 2.482 1.222 0.400 2.341 1.011 0.500 2.213 0.820 0.600 2.098 0.646 0.700 1.993 0.490 0.800 1.899 0.348 0.900 1.813 0.220 1.000 1.736 0.104 1.100 1.666 -0.001 1.200 1.602 -0.096 1.300 1.545 -0.182 1.400 1.493 -0.260 1.500 1.446 -0.331 1.600 1.404 -0.394 1.700 1.365 -0.452 1.800 1.331 -0.504 1.900 1.299 -0.551 2.000 1.271 -0.594 2.100 1.245 -0.633 2.200 1.222 -0.668 2.300 1.201 -0.699 2.400 1.181 -0.728 2.500 1.164 -0.754 2.600 1.149 -0.777 2.700 1.134 -0.798 2.800 1.122 -0.818 2.900 1.110 -0.835 3.000 1.100 -0.851 -0.5