PROGRAM HEART C REVISED 890628-0930, 950228-1300 C REVISED YYMMDD-HHMM C ILLUSTRATE THE USE OF THE DQED, HANSON-KROGH NONLINEAR LEAST C SQUARES SOLVER: SOLVING A CERTAIN HEART DIPOLE MOMENT EQUATION. C REFERENCE: A NEW NONLINEAR EQUATIONS TEST PROBLEM, RICE UNIV. C DEPT. OF MATH. REPORT 83-16, 6/83, 6/85, BY J. DENNIS, C JR., DAVID GAY, AND PHUONG AHN VU. C C COMMENTS ON THE MODULES DISTRIBUTED WITH THIS PACKAGE. C HEART - A TEST PROGRAM THAT READS AN INPUT FILE, CALLS C THE NONLINEAR TENSOR-QUADRATIC MODEL SOLVER, AND C WRITES THE OUTPUT FILE. OUTPUT RESULTS WILL VARY C IN ENVIRONMENTS DIFFERENT FROM THE PC/AT WITH THE C LAHEY F77L FORTRAN 77 COMPILER, v. 3.0. C DQED, DQEDMN, DQEDGN, DQEDIP, DQEDEV - C THE NONLINEAR SOLVING SOFTWARE. C DBOCLS, DBOLS, DBOLSM - C LINEARLY CONSTRAINED LINEAR LEAST SQUARES SOLVER. C DMOUT, DVOUT, IVOUT - C MATRIX-VECTOR PRINTING SUBPROGRAMS. C XERROR, XERRWV, CHRCNT - C A SMALL AND PARTIAL IMPLEMENTATION OF THE SLATEC C LIBRARY ERROR PROCESSING SUBPROGRAMS. DISCARD THESE C IF YOUR LOCAL LIBRARY HAS THIS PACKAGE. C D1MACH, I1MACH - C THE BELL LABS./PORT LIBRARY ENVIRONMENT INQUIRY C SUBPROGRAMS. DISCARD THESE IF YOUR LOCAL LIBRARY C HAS THIS PACKAGE. THE CURRENT VERSION IS SET FOR C THE LAHEY F77L FORTRAN 77 COMPILER ON THE PC/AT. C YOU MUST MODIFY THESE SUBROUTINES FOR YOUR OWN SYSTEM C IF THEY ARE NOT AVAILABLE IN YOUR LOCAL LIBRARY. C DGECO, DGEFA, DGESL - LINPACK LINEAR EQUATION SOLVING SOFTWARE. C (NOT INCLUDED IN THIS PACKAGE.) C DCOPY, DASUM, DDOT, DROTG, DAXPY, IDAMAX, DROT, DSWAP, DNRM2, C DSCAL - C BLAS, BASIC LINEAR ALGEBRA SUBPROGRAMS. C (NOT INCLUDED IN THIS PACKAGE.) C C EDIT on 950228-1300: ***LINE BELOW IS THE START OF THE INPUT FILE. Make this *** a separate file and remove the * in column 1 for all lines. *EXP. 0121C *-.807 -.021 -2.379 -3.64 -10.541 -1.961 -51.551 21.053 *-.074 -.733 0.013 -.034 -3.632 3.632 -.289 .289 *EXP. 0121B *-.809 -.021 -2.04 -.614 -6.903 -2.934 -26.328 18.639 *-.056 -.753 .026 -.047 -2.991 2.991 -.568 .568 *EXP. 0121A *-.816 -.017 -1.826 -.754 -4.839 -3.259 -14.023 15.467 *-.041 -.775 0.03 -.047 -2.565 2.565 -.754 .754 *EXP. 791226 *-.69 -.044 -1.57 -1.31 -2.65 2.0 -12.6 9.48 *-.3 -.39 .3 -.344 -1.2 2.69 1.59 -1.5 *EXP. 791129 *.485 -.0019 -.0581 .015 .105 .0406 .167 -.399 *.299 .186 -.0273 .0254 -.474 .474 -.0892 .0892 ***LINE ABOVE IS END OF THE INPUT FILE. C ***LINE BELOW IS THE START OF THE OUTPUT FILE. * PROBLEM TITLE IS = EXP. 0121C * HEART DIPOLE MODEL, * X(1),...,X(8) = * -0.692099 -0.114901 0.710964 -0.731964 * -0.815740 3.426596 1.184719 -2.332848 * RESIDUAL AFTER THE FIT = 1.5561E-10 * OUTPUT FLAG FROM SOLVER = 2 * NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS = 19 * * * PROBLEM TITLE IS = EXP. 0121B * HEART DIPOLE MODEL, * X(1),...,X(8) = * 0.009035 -0.818035 -0.000445 -0.020555 * 2.773429 2.529477 -14.800972 0.522047 * RESIDUAL AFTER THE FIT = 7.7707E-08 * OUTPUT FLAG FROM SOLVER = 6 * NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS = 19 * * * PROBLEM TITLE IS = EXP. 0121A * HEART DIPOLE MODEL, * X(1),...,X(8) = * 0.003100 -0.819100 -0.000224 -0.016776 * 2.681514 2.250216 -20.241705 0.797098 * RESIDUAL AFTER THE FIT = 9.6877E-13 * OUTPUT FLAG FROM SOLVER = 2 * NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS = 20 * * * PROBLEM TITLE IS = EXP. 791226 * HEART DIPOLE MODEL, * X(1),...,X(8) = * -0.311627 -0.378373 0.328244 -0.372244 * -1.282227 2.494300 1.554866 -1.384638 * RESIDUAL AFTER THE FIT = 3.3973E-08 * OUTPUT FLAG FROM SOLVER = 6 * NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS = 4 * * * PROBLEM TITLE IS = EXP. 791129 * HEART DIPOLE MODEL, * X(1),...,X(8) = * 0.491321 -0.006321 0.000098 -0.001998 * -0.100315 0.122657 -0.020718 -4.023518 * RESIDUAL AFTER THE FIT = 1.0066E-13 * OUTPUT FLAG FROM SOLVER = 2 * NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS = 18 ***LINE ABOVE IS THE END OF OUTPUT FILE. C C DIMENSION FOR THE NONLINEAR SOLVER. DOUBLE PRECISION FJ(08,9),BL(10),BU(10),X(8),ROPT(001),WA(2000) DOUBLE PRECISION SIGMA(8),Y(08),DF(08),T,D1MACH,RNORM INTEGER IND(10),IOPT(28),IWA(150) INTEGER LDFJ,LWA,LIWA,NVARS,NIN,I,J,NFEV,MCON,MEQUA,NOUT INTEGER IGO,K,MODE CHARACTER*12 TITLE LOGICAL WANTPR COMMON /SIGMA/SIGMA COMMON /MODE/MODE,NFEV EXTERNAL DQEDHD DATA LDFJ,LWA,LIWA/08,2000,150/ C C TO CHECK PARTIAL DERIVATIVES NUMERICALLY SET WANTPR=.TRUE. WANTPR = .FALSE. NVARS = 8 NIN = 1 open ( unit = NIN, file = 'dqed_prb1_input.txt', status = 'old' ) 10 CONTINUE C C READ THE PROBLEM TITLE. c READ (NIN,'(A)',END=60) TITLE C C READ IN THE SUMS OF VALUES AND INITIAL PARAMETER VALUES. READ (NIN,*,END=60) SIGMA,X C C TEST THE PARTIAL DERIVATIVE COMPUTATION. IF (WANTPR) THEN T = SQRT(D1MACH(4)) DO 30 J = 1,NVARS CALL DCOPY(NVARS,X,1,Y,1) Y(J) = X(J) + T*X(J) CALL DQEDHD(Y,FJ,LDFJ,IGO,IOPT,ROPT) CALL DCOPY(NVARS,FJ(1,NVARS+1),1,DF,1) Y(J) = X(J) - T*X(J) CALL DQEDHD(Y,FJ,LDFJ,IGO,IOPT,ROPT) DO 20 I = 1,NVARS DF(I) = (DF(I)-FJ(I,NVARS+1))/ (2.*T*X(J)) 20 CONTINUE CALL DQEDHD(X,FJ,LDFJ,IGO,IOPT,ROPT) CALL DVOUT(NVARS,DF,'(///,'' NUMERICAL DERIVATIVE'')',-4) CALL IVOUT(1,J,'('' = VARIABLE NUMBER'')',-4) CALL DVOUT(NVARS,FJ(1,J),'('' ANALYTIC PARTIAL'')',-4) 30 CONTINUE END IF C C SOLVE THE PROBLEM WITH THE FIRST TWO EQUATIONS AS CONSTRAINTS. DO 50 K = 1,1 C C THE LINEAR EQUATIONS WILL NOT BE CONSTRAINTS IF MODE=0. C CONVERGENCE IS NORMALLY FASTER IF THE LINEAR EQUATIONS ARE USED C AS CONSTRAINTS, MODE=1. MODE = K NFEV = 0 C TELL HOW MUCH STORAGE THE SOLVER HAS. IWA(1) = LWA IWA(2) = LIWA C C SET UP PRINT OPTION. (NOT USED HERE). IOPT(1) = -1 IOPT(2) = 1 C C SET UP FOR LINEAR MODEL WITHOUT ANY QUADRATIC TERMS. (NOT USED). IOPT(3) = -14 IOPT(4) = 1 C C DO NOT ALLOW CONVERGENCE TO BE CLAIMED ON SMALL STEPS. IOPT(5) = 17 IOPT(6) = 1 C C ALLOW UP TO NVARS = 8 QUADRATIC MODEL TERMS. IOPT(7)=15 IOPT(8)=NVARS C C CHANGE CONDITION NUMBER FOR QUADRATIC MODEL DEGREE. IOPT(9)=10 IOPT(10)=1 ROPT(01)=1.D4 C C NO MORE OPTIONS. IOPT(11) = 99 IF (MODE.EQ.1) THEN C SOLVE WITH TWO CONSTRAINT EQUATIONS. MCON = 2 IND(NVARS+1) = 3 IND(NVARS+2) = 3 BL(NVARS+1) = SIGMA(1) BU(NVARS+1) = SIGMA(1) BL(NVARS+2) = SIGMA(2) BU(NVARS+2) = SIGMA(2) ELSE C SOLVE WITH ALL REGRESSION EQUATIONS. MCON = 0 END IF MEQUA = NVARS - MCON C C ALL VARIABLES ARE OTHERWISE FREE. DO 40 I = 1,NVARS IND(I) = 4 40 CONTINUE CALL DQED(DQEDHD,MEQUA,NVARS,MCON,IND,BL,BU,X,FJ,LDFJ,RNORM, . IGO,IOPT,ROPT,IWA,WA) NOUT = 6 WRITE (NOUT,9031) TITLE WRITE (NOUT,9001) (X(J),J=1,NVARS) WRITE (NOUT,9011) RNORM WRITE (NOUT,9021) IGO,NFEV 50 CONTINUE GO TO 10 60 continue close ( unit = nin ) STOP 9001 FORMAT (' HEART DIPOLE MODEL,',/,' X(1),...,X(8) = ',/,4F12.6,/, . 4F12.6) 9011 FORMAT (' RESIDUAL AFTER THE FIT = ',1PE12.4) 9021 FORMAT (' OUTPUT FLAG FROM SOLVER =',17X,I6,/, . ' NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS =',08X,I6,//) 9031 FORMAT (' PROBLEM TITLE IS = ',A) END SUBROUTINE DQEDHD(X,FJ,LDFJ,IGO,IOPT,ROPT) C REVISED 860113-0930 C REVISED YYMMDD-HHMM C THIS IS THE SUBPROGRAM FOR EVALUATING THE FUNCTIONS C AND DERIVATIVES FOR THE NONLINEAR SOLVER, DQED. C C THE USER PROBLEM HAS MCON CONSTRAINT FUNCTIONS, C MEQUA LEAST SQUARES EQUATIONS, AND INVOLVES NVARS C UNKNOWN VARIABLES. C C WHEN THIS SUBPROGRAM IS ENTERED, THE GENERAL (NEAR) C LINEAR CONSTRAINT PARTIAL DERIVATIVES, THE DERIVATIVES C FOR THE LEAST SQUARES EQUATIONS, AND THE ASSOCIATED C FUNCTION VALUES ARE PLACED INTO THE ARRAY FJ(*,*). C ALL PARTIALS AND FUNCTIONS ARE EVALUATED AT THE POINT C IN X(*). THEN THE SUBPROGRAM RETURNS TO THE CALLING C PROGRAM UNIT. TYPICALLY ONE COULD DO THE FOLLOWING C STEPS: C C IF(IGO.NE.0) THEN C STEP 1. PLACE THE PARTIALS OF THE I-TH CONSTRAINT C FUNCTION WITH REPECT TO VARIABLE J IN THE C ARRAY FJ(I,J), I=1,...,MCON, J=1,...,NVARS. C END IF C STEP 2. PLACE THE VALUES OF THE I-TH CONSTRAINT C EQUATION INTO FJ(I,NVARS+1). C IF(IGO.NE.0) THEN C STEP 3. PLACE THE PARTIALS OF THE I-TH LEAST SQUARES C EQUATION WITH RESPECT TO VARIABLE J IN THE C ARRAY FJ(I,J), I=1,...,MEQUA, C J=1,...,NVARS. C END IF C STEP 4. PLACE THE VALUE OF THE I-TH LEAST SQUARES C EQUATION INTO FJ(I,NVARS+1). C STEP 5. RETURN TO THE CALLING PROGRAM UNIT. INTEGER LDFJ,NVARS,I,J,NFEV,MODE DOUBLE PRECISION FJ(LDFJ,*),X(*),ROPT(*),ONE,TWO,THREE,ZERO DOUBLE PRECISION SIGMA(08) COMMON /SIGMA/SIGMA COMMON /MODE/MODE,NFEV INTEGER IGO,IOPT(*) DATA ONE,TWO,THREE/1.D0,2.D0,3.D0/,ZERO/0.D0/ C C HEART DIPOLE EQUATIONS. FIRST TWO EQUATIONS MAY BE CONSTRAINTS. C C DEFINE THE CONSTRAINTS AND THEIR DERIVATIVES. NVARS = 8 DO 20 J = 1,NVARS DO 10 I = 1,2 FJ(I,J) = ZERO 10 CONTINUE 20 CONTINUE FJ(1,1) = ONE FJ(1,2) = ONE FJ(1,NVARS+1) = X(1) + X(2) FJ(2,3) = ONE FJ(2,4) = ONE FJ(2,NVARS+1) = X(3) + X(4) C C DEFINE THE REMAINING NONLINEAR EQUATIONS AND THEIR DERIVATIVES. C (START EQUATION ONE). FJ(3,1) = X(5) FJ(3,2) = X(6) FJ(3,3) = -X(7) FJ(3,4) = -X(8) FJ(3,5) = X(1) FJ(3,6) = X(2) FJ(3,7) = -X(3) FJ(3,8) = -X(4) FJ(3,NVARS+1) = X(1)*X(5) + X(2)*X(6) - X(3)*X(7) - X(4)*X(8) - . SIGMA(3) C (END EQUATION ONE). C (START EQUATION TWO). FJ(4,1) = X(7) FJ(4,2) = X(8) FJ(4,3) = X(5) FJ(4,4) = X(6) FJ(4,5) = X(3) FJ(4,6) = X(4) FJ(4,7) = X(1) FJ(4,8) = X(2) FJ(4,NVARS+1) = X(1)*X(7) + X(2)*X(8) + X(3)*X(5) + X(4)*X(6) - . SIGMA(4) C (END EQUATION TWO). C (START EQUATION THREE). FJ(5,1) = X(5)**2 - X(7)**2 FJ(5,2) = X(6)**2 - X(8)**2 FJ(5,3) = -TWO*X(5)*X(7) FJ(5,4) = -TWO*X(6)*X(8) FJ(5,5) = TWO* (X(1)*X(5)-X(3)*X(7)) FJ(5,6) = TWO* (X(2)*X(6)-X(4)*X(8)) FJ(5,7) = -TWO* (X(1)*X(7)+X(3)*X(5)) FJ(5,8) = -TWO* (X(2)*X(8)+X(4)*X(6)) FJ(5,NVARS+1) = X(1)*FJ(5,1) + X(2)*FJ(5,2) + X(3)*FJ(5,3) + . X(4)*FJ(5,4) - SIGMA(5) C (END EQUATION THREE). C (START EQUATION FOUR). C FJ(6,1) = TWO*X(5)*X(7) C FJ(6,2) = TWO*X(6)*X(8) C FJ(6,3) = X(5)**2 - X(7)**2 C FJ(6,4) = X(6)**2 - X(8)**2 FJ(6,1) = -FJ(5,3) FJ(6,2) = -FJ(5,4) FJ(6,3) = FJ(5,1) FJ(6,4) = FJ(5,2) C FJ(6,5) = TWO* (X(3)*X(5)+X(1)*X(7)) C FJ(6,6) = TWO* (X(4)*X(6)+X(2)*X(8)) C FJ(6,7) = -TWO* (X(3)*X(7)-X(1)*X(5)) C FJ(6,8) = -TWO* (X(4)*X(8)-X(2)*X(6)) FJ(6,5) = -FJ(5,7) FJ(6,6) = -FJ(5,8) FJ(6,7) = FJ(5,5) FJ(6,8) = FJ(5,6) FJ(6,NVARS+1) = X(1)*FJ(6,1) + X(2)*FJ(6,2) + X(3)*FJ(6,3) + . X(4)*FJ(6,4) - SIGMA(6) C (END EQUATION FOUR). C (START EQUATION FIVE). FJ(7,1) = X(5)* (X(5)**2-THREE*X(7)**2) FJ(7,2) = X(6)* (X(6)**2-THREE*X(8)**2) FJ(7,3) = X(7)* (X(7)**2-THREE*X(5)**2) FJ(7,4) = X(8)* (X(8)**2-THREE*X(6)**2) FJ(7,5) = THREE* (X(1)*FJ(5,1)+X(3)*FJ(5,3)) FJ(7,6) = THREE* (X(2)*FJ(5,2)+X(4)*FJ(5,4)) FJ(7,7) = -THREE* (X(1)*FJ(6,1)+X(3)*FJ(6,3)) FJ(7,8) = -THREE* (X(2)*FJ(6,2)+X(4)*FJ(6,4)) FJ(7,NVARS+1) = X(1)*FJ(7,1) + X(2)*FJ(7,2) + X(3)*FJ(7,3) + . X(4)*FJ(7,4) - SIGMA(7) C (END EQUATION FIVE). C (START EQUATION SIX). C FJ(8,1) = -X(7)* (X(7)**2-THREE*X(5)**2) C FJ(8,2) = -X(8)* (X(8)**2-THREE*X(6)**2) C FJ(8,3) = X(5)* (X(5)**2-THREE*X(7)**2) C FJ(8,4) = X(6)* (X(6)**2-THREE*X(8)**2) FJ(8,1) = -FJ(7,3) FJ(8,2) = -FJ(7,4) FJ(8,3) = FJ(7,1) FJ(8,4) = FJ(7,2) C FJ(8,5) = THREE* (X(1)*FJ(6,1)+X(3)*FJ(6,3)) C FJ(8,6) = THREE* (X(2)*FJ(6,2)+X(4)*FJ(6,4)) C FJ(8,7) = THREE* (X(1)*FJ(5,1)+X(3)*FJ(5,3)) C FJ(8,8) = THREE* (X(2)*FJ(5,2)+X(4)*FJ(5,4)) FJ(8,5) = THREE* (X(1)*FJ(6,1)+X(3)*FJ(6,3)) FJ(8,6) = THREE* (X(2)*FJ(6,2)+X(4)*FJ(6,4)) FJ(8,7) = FJ(7,5) FJ(8,8) = FJ(7,6) FJ(8,NVARS+1) = X(1)*FJ(8,1) + X(2)*FJ(8,2) + X(3)*FJ(8,3) + . X(4)*FJ(8,4) - SIGMA(8) C (END EQUATION SIX). IF (MODE.NE.1) THEN FJ(1,NVARS+1) = FJ(1,NVARS+1) - SIGMA(1) FJ(2,NVARS+1) = FJ(2,NVARS+1) - SIGMA(2) END IF NFEV = NFEV + 1 RETURN END