*DECK DLPDP SUBROUTINE DLPDP (A, MDA, M, N1, N2, PRGOPT, X, WNORM, MODE, WS, + IS) C***BEGIN PROLOGUE DLPDP C***SUBSIDIARY C***PURPOSE Subsidiary to DLSEI C***LIBRARY SLATEC C***TYPE DOUBLE PRECISION (LPDP-S, DLPDP-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION C C **** Double Precision version of LPDP **** C DIMENSION A(MDA,N+1),PRGOPT(*),X(N),WS((M+2)*(N+7)),IS(M+N+1), C where N=N1+N2. This is a slight overestimate for WS(*). C C Determine an N1-vector W, and C an N2-vector Z C which minimizes the Euclidean length of W C subject to G*W+H*Z .GE. Y. C This is the least projected distance problem, LPDP. C The matrices G and H are of respective C dimensions M by N1 and M by N2. C C Called by subprogram DLSI( ). C C The matrix C (G H Y) C C occupies rows 1,...,M and cols 1,...,N1+N2+1 of A(*,*). C C The solution (W) is returned in X(*). C (Z) C C The value of MODE indicates the status of C the computation after returning to the user. C C MODE=1 The solution was successfully obtained. C C MODE=2 The inequalities are inconsistent. C C***SEE ALSO DLSEI C***ROUTINES CALLED DCOPY, DDOT, DNRM2, DSCAL, DWNNLS C***REVISION HISTORY (YYMMDD) C 790701 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 891214 Prologue converted to Version 4.0 format. (BAB) C 900328 Added TYPE section. (WRB) C 910408 Updated the AUTHOR section. (WRB) C***END PROLOGUE DLPDP C INTEGER I, IS(*), IW, IX, J, L, M, MDA, MODE, MODEW, N, N1, N2, * NP1 DOUBLE PRECISION A(MDA,*), DDOT, DNRM2, FAC, ONE, * PRGOPT(*), RNORM, SC, WNORM, WS(*), X(*), YNORM, ZERO SAVE ZERO, ONE, FAC DATA ZERO,ONE /0.0D0,1.0D0/, FAC /0.1D0/ C***FIRST EXECUTABLE STATEMENT DLPDP N = N1 + N2 MODE = 1 IF (M .GT. 0) GO TO 20 IF (N .LE. 0) GO TO 10 X(1) = ZERO CALL DCOPY(N,X,0,X,1) 10 CONTINUE WNORM = ZERO GO TO 200 20 CONTINUE C BEGIN BLOCK PERMITTING ...EXITS TO 190 NP1 = N + 1 C C SCALE NONZERO ROWS OF INEQUALITY MATRIX TO HAVE LENGTH ONE. DO 40 I = 1, M SC = DNRM2(N,A(I,1),MDA) IF (SC .EQ. ZERO) GO TO 30 SC = ONE/SC CALL DSCAL(NP1,SC,A(I,1),MDA) 30 CONTINUE 40 CONTINUE C C SCALE RT.-SIDE VECTOR TO HAVE LENGTH ONE (OR ZERO). YNORM = DNRM2(M,A(1,NP1),1) IF (YNORM .EQ. ZERO) GO TO 50 SC = ONE/YNORM CALL DSCAL(M,SC,A(1,NP1),1) 50 CONTINUE C C SCALE COLS OF MATRIX H. J = N1 + 1 60 IF (J .GT. N) GO TO 70 SC = DNRM2(M,A(1,J),1) IF (SC .NE. ZERO) SC = ONE/SC CALL DSCAL(M,SC,A(1,J),1) X(J) = SC J = J + 1 GO TO 60 70 CONTINUE IF (N1 .LE. 0) GO TO 130 C C COPY TRANSPOSE OF (H G Y) TO WORK ARRAY WS(*). IW = 0 DO 80 I = 1, M C C MOVE COL OF TRANSPOSE OF H INTO WORK ARRAY. CALL DCOPY(N2,A(I,N1+1),MDA,WS(IW+1),1) IW = IW + N2 C C MOVE COL OF TRANSPOSE OF G INTO WORK ARRAY. CALL DCOPY(N1,A(I,1),MDA,WS(IW+1),1) IW = IW + N1 C C MOVE COMPONENT OF VECTOR Y INTO WORK ARRAY. WS(IW+1) = A(I,NP1) IW = IW + 1 80 CONTINUE WS(IW+1) = ZERO CALL DCOPY(N,WS(IW+1),0,WS(IW+1),1) IW = IW + N WS(IW+1) = ONE IW = IW + 1 C C SOLVE EU=F SUBJECT TO (TRANSPOSE OF H)U=0, U.GE.0. THE C MATRIX E = TRANSPOSE OF (G Y), AND THE (N+1)-VECTOR C F = TRANSPOSE OF (0,...,0,1). IX = IW + 1 IW = IW + M C C DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF C DWNNLS( ). IS(1) = 0 IS(2) = 0 CALL DWNNLS(WS,NP1,N2,NP1-N2,M,0,PRGOPT,WS(IX),RNORM, * MODEW,IS,WS(IW+1)) C C COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY W. SC = ONE - DDOT(M,A(1,NP1),1,WS(IX),1) IF (ONE + FAC*ABS(SC) .EQ. ONE .OR. RNORM .LE. ZERO) * GO TO 110 SC = ONE/SC DO 90 J = 1, N1 X(J) = SC*DDOT(M,A(1,J),1,WS(IX),1) 90 CONTINUE C C COMPUTE THE VECTOR Q=Y-GW. OVERWRITE Y WITH THIS C VECTOR. DO 100 I = 1, M A(I,NP1) = A(I,NP1) - DDOT(N1,A(I,1),MDA,X,1) 100 CONTINUE GO TO 120 110 CONTINUE MODE = 2 C .........EXIT GO TO 190 120 CONTINUE 130 CONTINUE IF (N2 .LE. 0) GO TO 180 C C COPY TRANSPOSE OF (H Q) TO WORK ARRAY WS(*). IW = 0 DO 140 I = 1, M CALL DCOPY(N2,A(I,N1+1),MDA,WS(IW+1),1) IW = IW + N2 WS(IW+1) = A(I,NP1) IW = IW + 1 140 CONTINUE WS(IW+1) = ZERO CALL DCOPY(N2,WS(IW+1),0,WS(IW+1),1) IW = IW + N2 WS(IW+1) = ONE IW = IW + 1 IX = IW + 1 IW = IW + M C C SOLVE RV=S SUBJECT TO V.GE.0. THE MATRIX R =(TRANSPOSE C OF (H Q)), WHERE Q=Y-GW. THE (N2+1)-VECTOR S =(TRANSPOSE C OF (0,...,0,1)). C C DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF C DWNNLS( ). IS(1) = 0 IS(2) = 0 CALL DWNNLS(WS,N2+1,0,N2+1,M,0,PRGOPT,WS(IX),RNORM,MODEW, * IS,WS(IW+1)) C C COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY Z. SC = ONE - DDOT(M,A(1,NP1),1,WS(IX),1) IF (ONE + FAC*ABS(SC) .EQ. ONE .OR. RNORM .LE. ZERO) * GO TO 160 SC = ONE/SC DO 150 J = 1, N2 L = N1 + J X(L) = SC*DDOT(M,A(1,L),1,WS(IX),1)*X(L) 150 CONTINUE GO TO 170 160 CONTINUE MODE = 2 C .........EXIT GO TO 190 170 CONTINUE 180 CONTINUE C C ACCOUNT FOR SCALING OF RT.-SIDE VECTOR IN SOLUTION. CALL DSCAL(N,YNORM,X,1) WNORM = DNRM2(N1,X,1) 190 CONTINUE 200 CONTINUE RETURN END