*DECK DLSEIT SUBROUTINE DLSEIT (LUN, KPRINT, IPASS) C***BEGIN PROLOGUE DLSEIT C***PURPOSE Quick check for DLSEI. C***LIBRARY SLATEC C***TYPE DOUBLE PRECISION (LSEIQX-S, DLSEIT-D) C***KEYWORDS QUICK CHECK C***AUTHOR Hanson, R. J., (SNLA) C Haskell, Karen, (SNLA) C***DESCRIPTION C C The sample problem solved is from a paper by J. Stoer, in C SIAM Journal of Numerical Analysis, June 1971. C C***ROUTINES CALLED D1MACH, DAXPY, DCOPY, DDOT, DLSEI, DNRM2, DVOUT C***REVISION HISTORY (YYMMDD) C 790216 DATE WRITTEN C 890618 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 901010 Restructured using IF-THEN-ELSE-ENDIF, modified tolerances C to use D1MACH(4) rather than D1MACH(3) and cleaned up C FORMATs. (RWC) C 920722 Initialized IP(1) and IP(2) for CALL to DLSEI. (BKS, WRB) C 930214 Declarations sections added, code revised to test error C returns for all values of KPRINT and code polished. (WRB) C***END PROLOGUE DDLSEIT C .. Scalar Arguments .. INTEGER IPASS, KPRINT, LUN C .. Local Scalars .. DOUBLE PRECISION CNORM, RELERR, RELNRM, RESNRM, RNORME, RNORML, * TNORM INTEGER I, IDIGIT, JDIGIT, KONTRL, MA, MDD, ME, MEAP1, MEP1, MG, * MODE, N, NERR, NP1 LOGICAL FATAL C .. Local Arrays .. DOUBLE PRECISION A(6,5), D(11,6), ERR(5), F(6), G(5,5), H(5), * PRGOPT(4), SOL(5), WORK(105), X(5) INTEGER IP(17) C .. External Functions .. DOUBLE PRECISION D1MACH, DDOT, DNRM2 INTEGER NUMXER EXTERNAL NUMXER, D1MACH, DDOT, DNRM2 C .. External Subroutines .. EXTERNAL DAXPY, DCOPY, DLSEI, DVOUT, XGETF, XSETF C .. Intrinsic Functions .. INTRINSIC SQRT C .. Data statements .. C C Define the data arrays for the example. The array A contains C the least squares equations. (There are no equality constraints C in this example). C DATA A(1,1),A(1,2),A(1,3),A(1,4),A(1,5) * /-74.,80.,18.,-11.,-4./ DATA A(2,1),A(2,2),A(2,3),A(2,4),A(2,5) * /14.,-69.,21.,28.,0./ DATA A(3,1),A(3,2),A(3,3),A(3,4),A(3,5) * /66.,-72.,-5.,7.,1./ DATA A(4,1),A(4,2),A(4,3),A(4,4),A(4,5) * /-12.,66.,-30.,-23.,3./ DATA A(5,1),A(5,2),A(5,3),A(5,4),A(5,5) * /3.,8.,-7.,-4.,1./ DATA A(6,1),A(6,2),A(6,3),A(6,4),A(6,5) * /4.,-12.,4.,4.,0./ C C The array G contains the inequality constraint equations, C written in the sense C (row vector)*(solution vector) .GE. (given value). C DATA G(1,1),G(1,2),G(1,3),G(1,4),G(1,5) * /-1.,-1.,-1.,-1.,-1./ DATA G(2,1),G(2,2),G(2,3),G(2,4),G(2,5) * /10.,10.,-3.,5.,4./ DATA G(3,1),G(3,2),G(3,3),G(3,4),G(3,5) * /-8.,1.,-2.,-5.,3./ DATA G(4,1),G(4,2),G(4,3),G(4,4),G(4,5) * /8.,-1.,2.,5.,-3./ DATA G(5,1),G(5,2),G(5,3),G(5,4),G(5,5) * /-4.,-2.,3.,-5.,1./ C C Define the least squares right-side vector. C DATA F(1),F(2),F(3),F(4),F(5),F(6) * /-5.,-9.,708.,4165.,-13266.,8409./ C C Define the inequality constraint right-side vector. C DATA H(1),H(2),H(3),H(4),H(5) * /-5.,20.,-40.,11.,-30./ C C Define the vector that is the known solution. C DATA SOL(1),SOL(2),SOL(3),SOL(4),SOL(5) * /1.,2.,-1.,3.,-4./ C***FIRST EXECUTABLE STATEMENT DDLSEIT IF (KPRINT .GE. 2) WRITE (LUN, 9000) C C Define the matrix dimensions, number of least squares equations, C number of equality constraints, total number of equations, and C number of variables. Set ME=0 to indicate there are no equality C constraints. C MDD = 11 MA = 6 MG = 5 N = 5 ME = 0 C IP(1) = 105 IP(2) = 17 C NP1 = N + 1 MEP1 = ME + 1 MEAP1 = ME + MA + 1 C C Copy the problem matrices. C DO 10 I = 1, N C C Copy the i-th column of the inequality constraint matrix into C the work array. C CALL DCOPY(MG, G(1,I), 1, D(MEAP1,I), 1) C C Copy the i-th column of the least squares matrix into the work C array. C CALL DCOPY(MA, A(1,I), 1, D(MEP1,I), 1) 10 CONTINUE C C Copy the right-side vectors into the work array in compatible C order. C CALL DCOPY(MG, H, 1, D(MEAP1,NP1), 1) CALL DCOPY(MA, F, 1, D(MEP1,NP1), 1) C C Use default program options in DLSEI, and set matrix-vector C printing accuracy parameters. C PRGOPT(1) = 1 IDIGIT = -4 JDIGIT = -11 C C Compute residual norm of known least squares solution. C (to be used to check computed residual norm = RNORML.) C DO 20 I = 1, MA WORK(I) = DDOT(N,D(I,1),MDD,SOL,1) - F(I) 20 CONTINUE RESNRM = DNRM2(MA,WORK,1) C C Call DLSEI to get solution in X(*), least squares residual in C RNORML. C CALL DLSEI(D, MDD, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML, MODE, * WORK, IP) C C Compute relative error in problem variable solution and residual C norm computation. C TNORM = DNRM2(N,SOL,1) CALL DCOPY(N, SOL, 1, ERR, 1) CALL DAXPY(N, -1.0D0, X, 1, ERR, 1) CNORM = DNRM2(N, ERR, 1) RELERR = CNORM/TNORM RELNRM = (RESNRM-RNORML)/RESNRM C IF (RELERR .LE. 70.0D0*SQRT(D1MACH(4)) .AND. * RELNRM .LE. 5.0D0*D1MACH(4)) THEN IPASS = 1 IF (KPRINT .GE. 3) WRITE (LUN, 9010) ELSE IPASS = 0 IF (KPRINT .GE. 2) WRITE (LUN, 9020) RELERR, RELNRM ENDIF C C Print out known and computed solutions. C IF (KPRINT .GE. 3) THEN CALL DVOUT(N, ERR, * '('' RESIDUALS FROM KNOWN LEAST SQUARES SOLUTION'')', * IDIGIT) CALL DVOUT(N, X, '(/'' SOLUTION COMPUTED BY DLSEI'')', JDIGIT) ENDIF C IF (KPRINT .GE. 2) THEN IF (.NOT.(KPRINT.EQ.2 .AND. IPASS.NE.0)) THEN C C Print out the known and computed residual norms. C CALL DVOUT(1, RESNRM, * '(/'' RESIDUAL NORM OF KNOWN LEAST SQUARES SOLUTION'')', * JDIGIT) CALL DVOUT(1, RNORML, * '(/'' RESIDUAL NORM COMPUTED BY DLSEI'')', JDIGIT) C C Print out the computed solution relative error. C CALL DVOUT(1, RELERR, * '(/'' COMPUTED SOLUTION RELATIVE ERROR'')', IDIGIT) C C Print out the computed relative error in residual norm. C CALL DVOUT(1, RELNRM, * '(/'' COMPUTED RELATIVE ERROR IN RESIDUAL NORM'')', IDIGIT) ENDIF ENDIF C C Check calls to error processor. C CALL XGETF (KONTRL) IF (KPRINT .LE. 2) THEN CALL XSETF (0) ELSE CALL XSETF (1) ENDIF FATAL = .FALSE. CALL XERCLR C IF (KPRINT .GE. 3) WRITE (LUN, 9030) C CALL DLSEI (D, 0, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML, * MODE, WORK, IP) IF (NUMXER(NERR) .NE. 2) THEN IPASS = 0 FATAL = .TRUE. ENDIF CALL XERCLR C PRGOPT(1) = -1 CALL DLSEI (D, MDD, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML, * MODE, WORK, IP) IF (NUMXER(NERR) .NE. 2) THEN IPASS = 0 FATAL = .TRUE. ENDIF CALL XERCLR C C Restore KONTRL and check to see if the tests of error detection C passed. C CALL XSETF (KONTRL) IF (FATAL) THEN IF (KPRINT .GE. 2) THEN WRITE (LUN, 9040) ENDIF ELSE IF (KPRINT .GE. 3) THEN WRITE (LUN, 9050) ENDIF ENDIF C C Print PASS/FAIL message. C IF (IPASS.EQ.1 .AND. KPRINT.GE.2) WRITE (LUN, 9100) IF (IPASS.EQ.0 .AND. KPRINT.GE.1) WRITE (LUN, 9110) RETURN C 9000 FORMAT ('1TEST OF SUBROUTINE DLSEI') 9010 FORMAT (/' DLSEI PASSED TEST') 9020 FORMAT (/' DLSEI FAILED TEST'/' RELERR = ',1P,D20.6/' RELNRM = ', * D20.6) 9030 FORMAT (/ ' 2 ERROR MESSAGES EXPECTED') 9040 FORMAT (/ ' AT LEAST ONE INCORRECT ARGUMENT TEST FAILED') 9050 FORMAT (/ ' ALL INCORRECT ARGUMENT TESTS PASSED') 9100 FORMAT (/' ****************DLSEI PASSED ALL TESTS***************') 9110 FORMAT (/' ****************DLSEI FAILED SOME TESTS**************') END