C ALGORITHM 635 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.11, NO. 3, C SEPT., 1985, P. 242-249. PROGRAM CPROG C CPROG IS A MAIN ROUTINE WHICH COMPUTES C EXAMPLES FOR ACM ALGORITHM CAPROX. C AUGUST 16, 1983 C REVISED: JULY 30, 1984, MAY 10, 1985 C INTEGER P1 PARAMETER (NDIM=8,NDIM21=2*NDIM+1,NDIM24=2*NDIM+4) PARAMETER (MDIM=285) PARAMETER (L=2,LDIM=2) PARAMETER (LOGP=10,P1=2**LOGP+1) C INTEGER IJSWT(3,NDIM21),IOEXIT(LOGP),ITLOG(LOGP),ICOUNT(NDIM) REAL C(NDIM),D(LDIM),CHEB(LOGP),ESTORE(NDIM) DOUBLE PRECISION BASINV(NDIM21,NDIM24),COSSIN(2,P1), + PI,DELTA,DUM COMPLEX A(NDIM,MDIM),Z(NDIM,LOGP),F(MDIM),B(NDIM,LDIM), + H(NDIM),G(LDIM),DATA(MDIM) C STANDARD FILE NUMBER OF THE OUTPUT STREAM ... FOR VAX 11/780 DATA NOUT/6/ C PI=4.0D0*DATAN(1.0D0) C C DEFINE THE APPROXIMATION DOMAIN C ... GENERATE MCIRC EQUISPACED DATA POINTS ON A C ... CIRCLE OF RADIUS 1 AND CENTER AT CMPLX(0,2). C ... THE FIRST POINT IS CMPLX(0,3). MCIRC=125 DELTA=2.D0*PI/FLOAT(MCIRC) DATA(1)=CMPLX(0.,3.) DO 10 K=2,MCIRC DUM=(PI/2.D0)+FLOAT(K-1)*DELTA SX=DCOS(DUM) SY=2.D0+DSIN(DUM) DATA(K)=CMPLX(SX,SY) 10 CONTINUE C ... GENERATE MSQUAR EQUISPACED DATA POINTS ON A SQUARE WITH C ... CORNERS AT CMPLX(1,-1), CMPLX(-1,-1), CMPLX(-1,-3), AND C ... CMPLX(1,-3). ASSUMES MSQUAR IS DIVISIBLE BY 4. MSIDE=40 MSQUAR=4*MSIDE DELTA=2.D0/FLOAT(MSIDE) DO 20 K=1,MSIDE DUM=FLOAT(K-1)*DELTA K1=K+MCIRC DELTAK=1.D0-DUM DATA(K1)=CMPLX(DELTAK,-1.) K1=K1+MSIDE DELTAK=-1.D0-DUM DATA(K1)=CMPLX(-1.,DELTAK) K1=K1+MSIDE DELTAK=-1.D0+DUM DATA(K1)=CMPLX(DELTAK,-3.) K1=K1+MSIDE DELTAK=-3.D0+DUM DATA(K1)=CMPLX(1.,DELTAK) 20 CONTINUE C ... SET THE TOTAL NUMBER OF POINTS M=MCIRC+MSQUAR C C DEFINE THE PROBLEM DATA FOR THE LARGEST ORDER C ... THE MATRIX A AND THE VECTOR F DO 40 I=1,M A(1,I)=DATA(I) DO 30 J=2,NDIM A(J,I)=DATA(I)**J 30 CONTINUE F(I)=CMPLX(1.0,0.0) 40 CONTINUE C ... THE MATRIX B AND THE VECTORS G AND D DO 50 I=1,NDIM B(I,1)=CMPLX(FLOAT(I),0.0) B(I,2)=CMPLX(FLOAT(I*(I-1)),0.0) 50 CONTINUE G(1)=CMPLX(0.0,0.0) G(2)=CMPLX(0.0,0.0) D(1)=.75 D(2)=.75 C ... SIMPLE BOUNDS WITH VECTORS H AND C DO 60 I=1,NDIM H(I)=CMPLX(0.0,0.0) C(I)=.2 60 CONTINUE C C SOLVE ALL PROBLEMS FROM ORDER 1 THROUGH NDIM DO 70 K=1,NDIM C ... GET CURRENT TIMING AND PAGE FAULT DATA: VAX 11/780 C CALL GETJPI(KTIME1,KPAGE1) C ... SET THE PROBLEM ORDER N=K C ... SET MAXIMUM ITERATION COUNT ITLOG(1)=100*(2*N+1) C ... SET INTERNAL PRINT-OUT OPTION IOEXIT(1)=-1 C ... SOLVE THE PROBLEM OF ORDER N CALL CAPROX(A,Z,F,B,N,NDIM,M,L,LOGP,H,C,G, + D,ITLOG,NDIM21,IJSWT,BASINV,COSSIN,CHEB,IOEXIT) C ... GET CURRENT TIMING AND PAGE FAULT DATA: VAX 11/780 C CALL GETJPI(KTIME2,KPAGE2) C SX=FLOAT(KTIME2-KTIME1)/100. C KP=KPAGE2-KPAGE1 C WRITE(NOUT,80)N,SX,KP C ... STORE SUMMARY RESULTS ESTORE(N)=CHEB(LOGP) ICOUNT(N)=ITLOG(LOGP) C ... PRINT SOLUTION VECTOR Z WRITE(NOUT,90) WRITE(NOUT,100)(J,Z(J,LOGP),J=1,N) WRITE(NOUT,130) 70 CONTINUE C C PRINT SUMMARY TABLE WRITE(NOUT,110) WRITE(NOUT,120)(J,ESTORE(J),ICOUNT(J),J=1,NDIM) C 80 FORMAT(9H ORDER = ,I3,16H, TIME IN SEC = ,F8.2, C + 16H, PAGE FAULTS = ,I7) 90 FORMAT(40H COMPONENT REAL PART IMAG PART) 100 FORMAT(I10,2F15.8) 110 FORMAT(//40H ORDER CHEBYSHEV ERROR ITERATIONS) 120 FORMAT(I10,F15.8,I10) 130 FORMAT(//////) END C SUBROUTINE CAPROX SOLVES SYSTEMS OF COMPLEX LINEAR EQUATIONS C IN THE CHEBYSHEV NORM WITH CONSTRAINTS ON THE UNKNOWNS. C C =============================================================== C PROBLEM: MIN (OVER Z) THE MAX (OVER K=1,...,M) OF C C ABS( Z * A(K) - F(K) ) C C SUBJECT TO: C C ABS( Z(K) - H(K) ) .LE. C(K) FOR K=1,...,N C C ABS( Z * B(K) - G(K) ) .LE. D(K) FOR K=1,...,L C C WHERE A(K) AND B(K) DENOTE COLUMN K OF MATRICES A AND B, RESP., C AND WHERE Z(K) DENOTES COMPONENT K OF THE ROW VECTOR Z. C =============================================================== C C WRITTEN BY: ROY L. STREIT C NAVAL UNDERWATER SYSTEMS CENTER C NEW LONDON LABORATORY C NEW LONDON, CONNECTICUT 06320 C JUNE 11, 1982 C REVISED: AUGUST 16, 1983, JULY 31, 1984, MAY 10, 1985 C C C REFERENCES: C 1. R.L. STREIT AND A.H. NUTTALL, "LINEAR CHEBYSHEV COMPLEX C FUNCTION APPROXIMATION," NUSC TECHNICAL REPORT 6403, C NAVAL UNDERWATER SYSTEMS CENTER, NEW LONDON LABORATORY, C NEW LONDON, CONNECTICUT, 26 FEBRUARY 1981. C 2. R.L. STREIT AND A.H. NUTTALL, "A GENERAL CHEBYSHEV C COMPLEX FUNCTION APPROXIMATION PROCEDURE AND AN APPLICATION C TO BEAMFORMING," JOURNAL OF THE ACOUSTICAL SOCIETY OF C AMERICA, VOL. 72(1), JULY, 1982, PP.181-190. C 3. R.L. STREIT AND A.H. NUTTALL, "A NOTE ON THE SEMI-INFINITE C PROGRAMMING APPROACH TO COMPLEX APPROXIMATION," MATHEMATICS C OF COMPUTATION, 40(1983), 599-605. C 4. R.L. STREIT, "SOLUTION OF SYSTEMS OF COMPLEX LINEAR C EQUATIONS IN THE L-SUB-INFINITY NORM WITH CONSTRAINTS C ON THE UNKNOWNS," SIAM JOURNAL ON SCIENTIFIC AND STATISTICAL C COMPUTATION, TO APPEAR (VOL 7, NUMBER 1, JANUARY 1986). C C C ALGORITHM SKETCH: C 1. SET P=4. THIS IS THE START OF DISCRETIZATION STAGE R:=1. C 2. APPROXIMATE THE ABSOLUTE VALUES OF ALL COMPLEX NUMBERS C WITH P PROJECTIONS; I.E., DISCRETIZE ABSOLUTE VALUES. C 3. SOLVE THE DUAL OF THE RESULTING PRIMAL LINEAR PROGRAM C USING THE REVISED SIMPLEX METHOD WITH EXPLICIT INVERSE. C (A) THE STORAGE REQUIRED IS REDUCED BY A FACTOR OF P C BY CONSTRUCTING THE INCOMING COLUMNS AS NEEDED. C (B) THE NUMBER OF INNER PRODUCTS REQUIRED IN EACH C SIMPLEX STEP TO COMPUTE THE REDUCED COSTS IS C DECREASED BY A FACTOR OF P BY UTILIZING THE C INTRINSIC COMPLEX ARITHMETIC PROBLEM STRUCTURE. C 4. DOUBLE THE NUMBER P OF PROJECTIONS AND LOOP BACK TO C STEP 2 TO START DISCRETIZATION STAGE R:=R+1. C 5. NOTE: THE TOTAL NUMBER OF STAGES NORMALLY EQUALS LOGP-1 C C C ARGUMENT LIST: C COMPLEX A(NDIM,*) = COEFFICIENT MATRIX IS CONTAINED IN THE C FIRST N ROWS. MUST HAVE .GE. M COLUMNS C COMPLEX Z(NDIM,*) = ON OUTPUT, CONTAINS THE SOLUTION VECTORS. C THE R-TH COLUMN IS SOLUTION VECTOR OF THE C R-TH STAGE. MUST HAVE .GE. LOGP COLUMNS. C ON OUTPUT, THE BEST COMPUTED SOLUTION IS C STORED IN THE LAST (I.E., LOGP) COLUMN. C COMPLEX F(*) = CONSTANT VECTOR. MUST BE OF LENGTH .GE. M C COMPLEX B(NDIM,*) = GENERAL BOUND MATRIX COEFFICIENTS ARE C CONTAINED IN THE FIRST N ROWS. C MUST HAVE AT LEAST ONE COLUMN WHEN L=0. C INTEGER N = NUMBER OF UNKNOWNS. REQUIRES N.GE.1 C INTEGER NDIM = NUMBER OF ROWS OF MATRICES A, Z, AND B IN C THE DIMENSION STATEMENT IN THE CALLING C PROGRAM. REQUIRES NDIM .GE. N C INTEGER M = NUMBER OF EQUATIONS. REQUIRES M.GE.1 C INTEGER L = NUMBER OF GENERAL BOUND CONSTRAINTS; L.GE.0 C INTEGER LOGP = THE BASE 2 LOG OF THE MAXIMUM NUMBER OF C PROJECTIONS: P=2**LOGP. REQUIRES LOGP.GE.1 C WARNING: LOGP=1 WORKS ONLY IF ALL THE C PROBLEM DATA ARE REAL-VALUED. C COMPLEX H(*) = CENTERS OF CIRCLES BOUNDING Z COMPONENTS. C MUST BE OF LENGTH .GE. N C REAL C(*) = RADII OF CIRCLES BOUNDING Z COMPONENTS. C REQUIRES C(K) .GT. 0 FOR ALL K. C MUST BE OF LENGTH .GE. N C COMPLEX G(*) = CONSTANT TERM OF GENERAL BOUND CONSTRAINTS. C MUST BE OF LENGTH .GE. MAX(1,L) C REAL D(*) = RADII OF GENERAL BOUND CONSTRAINTS. C REQUIRES D(K) .GT. 0 FOR ALL K. C MUST BE OF LENGTH .GE. MAX(1,L) C INTEGER ITLOG(*) = INPUT/OUTPUT ARRAY OF LENGTH .GE. LOGP C INPUT: ITLOG(1) = MAXIMUM ITERATION COUNT C OUTPUT: ITLOG(R) = ITERATION COUNT TO C FINISH STAGE R. SPECIAL CASE: C ITLOG(LOGP)=FINAL ITERATION COUNT. C INTEGER NDIM21 = 2*NDIM + 1 C INTEGER IJSWT(3,*)= CONTAINS THE BASIC VARIABLE NAMES. C MUST HAVE .GE. (2*N+1) COLUMNS. C THE J-TH BASIC VARIABLE IS SPECIFIED BY: C IZSWT(1,J)=1, 2, OR 3 IF AN S, W, OR T C IJSWT(2,J)=EQUATION/CONSTRAINT NUMBER C IJSWT(3,J)=PROJECTION NUMBER C DOUBLE BASINV(NDIM21,*) = MUST HAVE .GE. (2*N+4) COLUMNS. C PRECISION IDENTIFIED BY COLUMNS AS FOLLOWS: C COL. 1 TO 2*N+1=CURRENT BASIS INVERSE C COL. 2*N+1 + 1 =CURRENT B.F.S. C COL. 2*N+1 + 2 =INCOMING COL OF NEW B.F.S. C COL. 2*N+1 + 3 =COSTS OF CURRENT B.F.S. C DOUBLE COSSIN(2,*)= COSINES AND SINES OF THE P PROJECTION C PRECISION ANGLES. DEFINED IN ROUTINE PABS. C MUST HAVE .GE. (P+1) COLUMNS C REAL CHEB(*) = AT END OF STAGE R, CHEBYSHEV MULTIPLIER C IS CHEB(R). MUST BE OF LENGTH .GE. LOGP C SPECIAL CASE: CHEB(LOGP)=CHEBYSHEV MULTI- C PLIER OF BEST SOLUTION. C INTEGER IOEXIT(*) = INPUT/OUTPUT ARRAY OF LENGTH .GE. LOGP C ..........INPUT: SETS PRINT OPTION ACCORDING TO C IOEXIT(1) = -1 IF ONLY BEST SOLUTION IS PRINTED C = 0 IF NOTHING IS TO BE PRINTED C = 1 IF ONLY B.F.S. AND SIMPLEX MULTIPLIERS C AT END OF EACH STAGE ARE PRINTED C = 2 IF INCOMING AND OUTGOING BASIC VARIABLES C AT EACH ITERATION ARE ALSO PRINTED C ..........OUTPUT: AT THE END OF THE R-TH STAGE, C IOEXIT(R) = 4 IF MAXIMUM ITERATION COUNT REACHED C = 3 IF CHEBYSHEV ERROR ESTIMATE DECREASED C THREE TIMES IN A ROW C = 2 IF INCOMING AND OUTGOING BASIC VARIABLES C ARE IDENTICAL C = 1 IF ALL REDUCED COSTS ARE .GE. -TOLP C = 0 IF ALL REDUCED COSTS ARE .GE. ZERO C = -1 IF PRIMAL PROBLEM IS INFEASIBLE C = -2 IF IMPOSSIBLE BASIS FOUND TWICE IN A ROW C = -3 IF SINGULAR BASIS IS DETECTED C = -10000-J IF J INPUT ERRORS WERE FOUND. C (IN THIS CASE, R=1 AND PROGRAM EXITS.) C = -999 IF THE R-TH STAGE WAS NOT ATTEMPTED C SPECIAL CASE: IOEXIT(LOGP) = MAXIMUM (OVER ALL C STAGES R) OF THE NUMBERS IOEXIT(R). C C THE CONTENTS OF THE FOLLOWING ARRAYS ARE NOT ALTERED BY THIS C SUBROUTINE... A, F, B, H, C, G, D. C SUBROUTINE CAPROX(A,Z,F,B,N,NDIM,M,L,LOGP,H,C, + G,D,ITLOG,NDIM21,IJSWT, + BASINV,COSSIN,CHEB,IOEXIT) C C TYPE AND DIMENSION ALL SUBROUTINE ARGUMENTS INTEGER N,NDIM,M,L,LOGP,ITLOG(1),NDIM21,IJSWT(3,*),IOEXIT(1) COMPLEX A(NDIM,*),Z(NDIM,*),F(1),B(NDIM,*),H(1),G(1) REAL C(1),D(1),CHEB(1) DOUBLE PRECISION COSSIN(2,*),BASINV(NDIM21,*) C C TYPE ALL INTERNAL VARIABLE NAMES COMPLEX ZDUM REAL ZERO DOUBLE PRECISION X,Y,DZERO,DONE,DRANGE,DRELPR,CHEBBO, + CHEBER,RCC,TOLP,ZR,ZI,PNORM INTEGER I,J,K,PARG,KPT,JPT,IWHICH,ITER,ITMAX,LOOP,R, + NOUT,ISTAR,KSTAR,JSTAR,IMPOSS,IPRT,SLACK,ICHEB,CYCLE, + N2,NP1,N2M1,N2P1,N2P2,N2P3,N2P4,LOGM1,T(3),TSAVE(3) C DATA DZERO,DONE/+0.0D+0,+1.0D+0/ DATA ZERO/+0.0E+0/ C C MACHINE DEPENDENT CONSTANTS --- THESE ARE FOR THE DEC VAX 11/780 C ... LARGEST DOUBLE PRECISION NUMBER DATA DRANGE/+1.0D+38/ C ... RELATIVE PRECISION OF SINGLE PRECISION ARITHMETIC DATA DRELPR/+0.6D-7/ C ... STANDARD FILE NUMBER OF THE OUTPUT STREAM DATA NOUT/6/ C C CHECK INPUT ARGUMENTS J=0 IF(LOGP.LT.1)J=J+1 IF(M.LT.1)J=J+1 IF(N.LT.1)J=J+1 IF(L.LT.0)J=J+1 IF(NDIM21.NE.(NDIM+NDIM+1))J=J+1 IF(ITLOG(1).LT.1)J=J+1 IF(IOEXIT(1).LT.(-1).OR.IOEXIT(1).GT.2)J=J+1 DO 10 I=1,N IF(C(I).LE.ZERO)J=J+1 10 CONTINUE IF(L.LE.0)GO TO 30 DO 20 I=1,L IF(D(I).LE.ZERO)J=J+1 20 CONTINUE 30 IF(J.GT.0)GO TO 710 C C INITIALIZE IPRT=IOEXIT(1) ITMAX=ITLOG(1) ITER=-1 NP1=N+1 N2=N+N N2M1=N2-1 N2P1=N2+1 N2P2=N2P1+1 N2P3=N2P2+1 N2P4=N2P3+1 LOGM1=LOGP-1 IF(LOGM1.LT.1)LOGM1=1 CHEBBO=DZERO SLACK=1 CYCLE=1 DO 50 I=1,LOGP DO 40 K=1,N Z(K,I)=CMPLX(ZERO,ZERO) 40 CONTINUE ITLOG(I)=0 CHEB(I)=ZERO IOEXIT(I)=-999 50 CONTINUE T(1)=0 T(2)=0 T(3)=0 C C DEFINE THE FIRST BASIS INVERSE AND B.F.S. DO 80 I=1,N2P3 DO 70 J=1,N2P1 IF(I.NE.J)GO TO 60 BASINV(J,J)=DONE GO TO 70 60 BASINV(J,I)=DZERO 70 CONTINUE 80 CONTINUE BASINV(N2P1,N2P2)=DONE C C IDENTIFY THE B.F.S. C ... THE W VARIABLES DO 90 I=1,N2M1,2 K=(I+1)/2 IJSWT(1,I)=2 IJSWT(2,I)=K IJSWT(3,I)=1 IJSWT(1,I+1)=2 IJSWT(2,I+1)=K IJSWT(3,I+1)=2 90 CONTINUE C ... THE SLACK VARIABLE IJSWT(1,N2P1)=0 IJSWT(2,N2P1)=0 IJSWT(3,N2P1)=0 C C DEFINE THE COST COEFFICIENTS C ... OF THE W VARIABLES DO 100 I=1,N2M1,2 K=(I+1)/2 BASINV(I ,N2P4)=DBLE(C(K))+DBLE( REAL(H(K))) BASINV(I+1,N2P4)=DBLE(C(K))+DBLE(AIMAG(H(K))) 100 CONTINUE C ... OF THE SLACK VARIABLE BASINV(N2P1,N2P4)=DZERO C C ========================= C BEGIN DISCRETIZATION LOOP C ========================= C DO 690 R=1,LOGM1 LOOP=R+1 IMPOSS=0 ICHEB=0 IF(LOOP.EQ.2)GO TO 110 IF(IOEXIT(R-1).GT.2)GO TO 150 IF(IOEXIT(R-1).LT.3)GO TO 180 C C --------------------------------------------------------- C BEGIN THE REVISED SIMPLEX ALGORITHM WITH EXPLICIT INVERSE C --------------------------------------------------------- C C COMPUTE THE CHEBYSHEV SIMPLEX MULTIPLIER 110 CHEBER=DZERO DO 120 J=1,N2P1 CHEBER=CHEBER-BASINV(J,N2P4)*BASINV(J,N2P1) 120 CONTINUE C C TESTS FOR MINIMUM INCREMENT TERMINATION IF(SLACK.GT.0)GO TO 150 IF(CHEBER.LT.CHEBBO)GO TO 130 GO TO 140 130 ICHEB=ICHEB+1 IF(ICHEB.LT.4)GO TO 150 GO TO 630 140 ICHEB=0 150 RCC=CHEBER-CHEBBO CHEBBO=CHEBER C C COMPUTE THE COMPLEX SIMPLEX MULTIPLIERS DO 170 I=1,N2M1,2 K=(I+1)/2 X=DZERO Y=DZERO DO 160 J=1,N2P1 X=X+BASINV(J,N2P4)*BASINV(J,I) Y=Y+BASINV(J,N2P4)*BASINV(J,I+1) 160 CONTINUE Z(K,R)=CMPLX(SNGL(X),SNGL(Y)) 170 CONTINUE C C TEST FOR MAXIMUM ITERATION COUNT ITER=ITER+1 IF(ITER.GT.ITMAX)GO TO 580 C C PRINT CHEBYSHEV MULTIPLIER THIS ITERATION, IF REQUESTED IF(IPRT.LT.2)GO TO 180 WRITE(NOUT,730)LOOP,ITER,CHEBER,RCC C C DETERMINE THE ENTERING BASIC VARIABLE := (IWHICH,KPT,JPT). C TIES ARE BROKEN BY CHOOSING THE MOST NEGATIVE REDUCED C COST COEFFICIENT WITH THE SMALLEST INDEX. EXCEPTIONS: C BRING IN AN S VARIABLE IF POSSIBLE WHEN SLACK VARIABLE C IS STILL IN THE BASIS. THIS SHOULD AVOID CYCLING. 180 IWHICH=0 RCC=DRANGE C ... CHECK THE S VARIABLES DO 210 K=1,M ZDUM=CMPLX(ZERO,ZERO) DO 190 J=1,N ZDUM=ZDUM+Z(J,R)*A(J,K) 190 CONTINUE ZDUM=ZDUM-F(K) CALL PABS(ZDUM,LOOP,PNORM,PARG,COSSIN) X=-PNORM+CHEBER IF(X.GE.RCC)GO TO 210 IF(IMPOSS.NE.1)GO TO 200 IF(K.NE.KSTAR)GO TO 200 IF(ISTAR.EQ.1)GO TO 210 200 JPT=PARG KPT=K RCC=X IWHICH=1 210 CONTINUE IF(CYCLE.EQ.0)GO TO 220 IF(RCC.GT.(-DRELPR))CYCLE=0 IF(CYCLE.EQ.1)GO TO 290 C ... CHECK THE W VARIABLES 220 DO 240 K=1,N ZDUM=Z(K,R)-H(K) CALL PABS(ZDUM,LOOP,PNORM,PARG,COSSIN) X=-PNORM+C(K) IF(X.GE.RCC)GO TO 240 IF(IMPOSS.NE.1)GO TO 230 IF(K.NE.KSTAR)GO TO 230 IF(ISTAR.EQ.2)GO TO 240 230 JPT=PARG KPT=K RCC=X IWHICH=2 240 CONTINUE IF(L.EQ.0)GO TO 280 C ... CHECK THE T VARIABLES DO 270 K=1,L ZDUM=CMPLX(ZERO,ZERO) DO 250 J=1,N ZDUM=ZDUM+Z(J,R)*B(J,K) 250 CONTINUE ZDUM=ZDUM-G(K) CALL PABS(ZDUM,LOOP,PNORM,PARG,COSSIN) X=-PNORM+D(K) IF(X.GE.RCC)GO TO 270 IF(IMPOSS.NE.1)GO TO 260 IF(K.NE.KSTAR)GO TO 260 IF(ISTAR.EQ.3)GO TO 270 260 JPT=PARG KPT=K RCC=X IWHICH=3 270 CONTINUE C C TEST FOR TERMINATION 280 IF(IWHICH.GT.0)GO TO 290 IF(IMPOSS.EQ.0)GO TO 610 GO TO 560 290 IF(RCC.GE.DZERO.AND.CYCLE.EQ.0)GO TO 610 C C GENERATE INCOMING S, W, OR T COLUMN X=COSSIN(1,JPT) Y=COSSIN(2,JPT) IF(IWHICH.NE.2)GO TO 310 C ... FINAL FORM OF ANY W COLUMN K=2*KPT-1 DO 300 I=1,N2P1 BASINV(I,N2P3)=X*BASINV(I,K)+Y*BASINV(I,K+1) 300 CONTINUE GO TO 380 C BASIS INVERSE TIMES ORIGINAL S OR T COLUMN 310 DO 320 I=1,N2P1 BASINV(I,N2P3)=DZERO 320 CONTINUE DO 360 J=1,N IF(IWHICH.GT.2)GO TO 330 ZR=DBLE( REAL(A(J,KPT)))*X+DBLE(AIMAG(A(J,KPT)))*Y ZI=DBLE(AIMAG(A(J,KPT)))*X-DBLE( REAL(A(J,KPT)))*Y GO TO 340 330 ZR=DBLE( REAL(B(J,KPT)))*X+DBLE(AIMAG(B(J,KPT)))*Y ZI=DBLE(AIMAG(B(J,KPT)))*X-DBLE( REAL(B(J,KPT)))*Y 340 K=2*J-1 DO 350 I=1,N2P1 BASINV(I,N2P3)=BASINV(I,N2P3)+ZR*BASINV(I,K) + -ZI*BASINV(I,K+1) 350 CONTINUE 360 CONTINUE C ... FINAL FORM OF ANY T COLUMN IF(IWHICH.GT.2)GO TO 380 C ... FINAL FORM OF ANY S COLUMN DO 370 I=1,N2P1 BASINV(I,N2P3)=BASINV(I,N2P3)+BASINV(I,N2P1) 370 CONTINUE C C DETERMINE NUMERICAL TOLERANCE FOR PIVOT AND REDUCED COST 380 X=DABS(BASINV(1,N2P3)) DO 390 I=2,N2P1 X=X+DABS(BASINV(I,N2P3)) 390 CONTINUE TOLP=X*DRELPR C C REDUCED COST TOLERANCE TESTS FOR TERMINATION IF(TOLP.EQ.DZERO)GO TO 570 IF(RCC.LE.(-TOLP))GO TO 400 IF(CYCLE.EQ.0)GO TO 590 CYCLE=0 GO TO 220 C C DETERMINE THE OUTGOING BASIC VARIABLE := K-TH IN THE LIST. C TIES IN THE MINIMUM RATIO TEST ARE BROKEN BY SELECTING THE C MAXIMUM MAGNITUDE PIVOT WITH THE SMALLEST INDEX. C EXCEPTIONS: PIVOT OUT A W VARIABLE IF POSSIBLE WHEN C SLACK VARIABLE STILL IN BASIS. THIS SHOULD AVOID CYCLING. 400 X=DRANGE K=0 DO 440 I=1,N2P1 IF(BASINV(I,N2P3).LT.TOLP)GO TO 440 PNORM=BASINV(I,N2P2)/BASINV(I,N2P3) IF(PNORM.LT.DZERO)PNORM=DZERO IF(PNORM.GT.X)GO TO 440 IF(PNORM.LT.X)GO TO 420 IF(K.LT.1)GO TO 430 IF(SLACK.LT.1)GO TO 410 IF(IJSWT(1,K).NE.2.AND.IJSWT(1,I).EQ.2)GO TO 430 IF(IJSWT(1,K).NE.2.AND.IJSWT(1,I).NE.2)GO TO 410 IF(IJSWT(1,I).NE.2)GO TO 440 410 IF(BASINV(I,N2P3).LE.Y)GO TO 440 GO TO 430 420 X=PNORM 430 Y=BASINV(I,N2P3) K=I 440 CONTINUE C C TERMINATE IF SOLUTION IS UNBOUNDED BELOW IF(K.EQ.0)GO TO 600 C C UPDATE THE BASIC VARIABLE LIST IF(K.EQ.N2P1)SLACK=0 IF(K.EQ.N2P1)CYCLE=0 ISTAR=IJSWT(1,K) KSTAR=IJSWT(2,K) JSTAR=IJSWT(3,K) IJSWT(1,K)=IWHICH IJSWT(2,K)=KPT IJSWT(3,K)=JPT C ... TEST FOR IDENTICAL INCOMING AND OUTGOING BASIS VARIABLES IF(KSTAR.NE. KPT)GO TO 450 IF(JSTAR.NE. JPT)GO TO 450 IF(ISTAR.NE.IWHICH)GO TO 450 IF(CYCLE.EQ.0)GO TO 620 CYCLE=0 GO TO 220 C ... TERMINATE IF IMPOSSIBLE BASIS TWICE IN A ROW 450 J=0 DO 460 I=1,N2P1 IF(IJSWT(2,I).NE.KPT)GO TO 460 IF(IJSWT(1,I).NE.IWHICH)GO TO 460 J=J+1 460 CONTINUE IF(J.LE.2)GO TO 470 IF(IMPOSS.EQ.1)GO TO 560 IMPOSS=1 CYCLE=0 IJSWT(1,K)=ISTAR IJSWT(2,K)=KSTAR IJSWT(3,K)=JSTAR ISTAR=IWHICH KSTAR=KPT GO TO 180 470 IMPOSS=0 C C PRINT INCOMING AND OUTGOING BASIC VARIABLE, IF REQUESTED IF(IPRT.LT.2)GO TO 480 WRITE(NOUT,740)LOOP,ITER,IJSWT(1,K),IJSWT(2,K),IJSWT(3,K), + ISTAR,KSTAR,JSTAR C C UPDATE COST COEFFICIENT 480 IF(IWHICH.GT.1)GO TO 490 BASINV(K,N2P4)=DBLE( REAL(F(KPT)))*COSSIN(1,JPT) + +DBLE(AIMAG(F(KPT)))*COSSIN(2,JPT) GO TO 510 490 IF(IWHICH.GT.2)GO TO 500 BASINV(K,N2P4)=DBLE(C(KPT))+DBLE(REAL(H(KPT)))* + COSSIN(1,JPT)+DBLE(AIMAG(H(KPT)))*COSSIN(2,JPT) GO TO 510 500 BASINV(K,N2P4)=DBLE(D(KPT))+DBLE(REAL(G(KPT)))* + COSSIN(1,JPT)+DBLE(AIMAG(G(KPT)))*COSSIN(2,JPT) C C PIVOT TO UPDATE B.F.S. AND BASIS INVERSE 510 X=BASINV(K,N2P3) KPT=K-1 KSTAR=K+1 DO 550 J=1,N2P2 Y=BASINV(K,J)/X IF(K.EQ.1)GO TO 530 DO 520 I=1,KPT BASINV(I,J)=BASINV(I,J)-Y*BASINV(I,N2P3) 520 CONTINUE 530 BASINV(K,J)=Y IF(K.EQ.N2P1)GO TO 550 DO 540 I=KSTAR,N2P1 BASINV(I,J)=BASINV(I,J)-Y*BASINV(I,N2P3) 540 CONTINUE 550 CONTINUE C GO TO 110 C C ------------------------------------- C END OF REVISED SIMPLEX ALGORITHM LOOP C ------------------------------------- C C DEFINE ALL EXIT CONDITIONS FROM SIMPLEX ALGORITHM C ... IMPOSSIBLE BASIS; EITHER TWICE IN A ROW... FROM C NEAR 460, OR ONLY POSSIBLE PIVOT... FROM NEAR 290 560 IOEXIT(R)=-2 GO TO 640 C ... SINGULAR BASIS INVERSE, IMPOSSIBLE... FROM NEAR 390 570 IOEXIT(R)=-3 GO TO 640 C ... MAXIMUM ITERATION COUNT ACHIEVED... FROM NEAR 170 580 IOEXIT(R)=4 GO TO 640 C ... REDUCED COSTS EXCEED TOLERANCE... FROM NEAR 400 590 IOEXIT(R)=1 GO TO 640 C ... SOLUTION UNBOUNDED BELOW; PRIMAL INFEASIBLE... NEAR 440 600 IOEXIT(R)=-1 GO TO 640 C ... REDUCED COSTS ARE NON-NEGATIVE... FROM NEAR 280,290 610 IOEXIT(R)=0 GO TO 640 C ... IN- AND OUT-GOING BASIC VARIABLES THE SAME... NEAR 450 620 IOEXIT(R)=2 GO TO 640 C ... CONSECUTIVE DECREASES IN CHEBYSHEV MULT ... FROM NEAR 130 630 IOEXIT(R)=3 CHEB(R)=CHEBBO ITLOG(R)=ITER GO TO 650 640 CHEB(R)=CHEBER ITLOG(R)=ITER-1 C C PRINT B.F.S. AND SIMPLEX MULTIPLIERS, IF REQUESTED 650 K=IOEXIT(R) J=2**LOOP IF(LOGP.EQ.1)J=2 TSAVE(1)=T(1) TSAVE(2)=T(2) TSAVE(3)=T(3) CALL CPAIRS(IJSWT,BASINV,NDIM21,N2P1,N2P2,K,J,IPRT,T) IF(IPRT.LT.1)GO TO 660 WRITE(NOUT,760)(K,Z(K,R),K=1,N) WRITE(NOUT,770)NP1,CHEB(R) WRITE(NOUT,780)J,LOOP,IOEXIT(R),ITLOG(R) C C TERMINATE IF AND ONLY IF SIMPLEX ALGORITHM ENDED BADLY 660 IF(IOEXIT(R).LT.0.OR.IOEXIT(R).GT.3)GO TO 700 C C INITIALIZE NEXT STAGE, UNLESS FINISHED IF(R.EQ.LOGM1)GO TO 700 C ... BASIS NAMES DO 670 K=1,N2P1 IJSWT(3,K)=2*IJSWT(3,K)-1 670 CONTINUE C ... COMPLEX SIMPLEX MULTIPLIERS DO 680 K=1,N Z(K,LOOP)=Z(K,R) 680 CONTINUE 690 CONTINUE C C ========================== C END OF DISCRETIZATION LOOP C ========================== C C C STORE THE BEST COMPUTED SOLUTION AND RELEVANT DATA. C (ALTERNATIVELY, AT THIS POINT ONE COULD CALL A C SUBROUTINE DESIGNED TO IMPROVE THE BEST SOLUTION C PROVIDED BY THE DISCRETIZATION PROCESS.) 700 CALL CEND(Z,CHEB,LOGP,LOGM1,N,NDIM,IOEXIT,ITLOG,K) IF(IPRT.EQ.0)GO TO 720 IF(K.EQ.0)WRITE(NOUT,800) IF(K.EQ.0)GO TO 720 WRITE(NOUT,790) WRITE(NOUT,760)(I,Z(I,LOGP),I=1,N) WRITE(NOUT,770)NP1,CHEB(LOGP) IF(K.EQ.R)WRITE(NOUT,820)T(1),T(2),T(3) IF(K.EQ.(R-1))WRITE(NOUT,820)TSAVE(1),TSAVE(2),TSAVE(3) K=K+1 IF(LOGP.EQ.1)K=1 J=2**K WRITE(NOUT,780)J,K,IOEXIT(LOGP),ITLOG(LOGP) GO TO 720 C C INPUT ERRORS... FROM 30 710 IOEXIT(1)=-10000-J WRITE(NOUT,810)J C C UNIQUE PROGRAM EXIT 720 RETURN C 730 FORMAT(I4,I7,23H CHEBYSHEV MULTIPLIER =,D15.8,10H CHANGE =,D15.8) 740 FORMAT(I4,I7,6H NEW=(,I2,1H,,I6,1H,,I6,1H),9H OLD=(,I2, + 1H,,I6,1H,,I6,1H)) 750 FORMAT(17H BASIS. NAME = (,I2,1H,,I6,1H,,I6,10H) VALUE =,D15.8) 760 FORMAT(23H COMPLEX MULTIPLIER (,I6,4H) = ,E15.8,5H + I ,E15.8) 770 FORMAT(23H CHEBYSHEV MULTIPLIER (,I6,4H) = ,E15.8) 780 FORMAT(4H P =,I8,9H, LOG P =,I3,14H, IOEXIT(R) = ,I2, + 21H, ITERATION NUMBER = ,I7//) 790 FORMAT(//24H BEST COMPUTED SOLUTION.) 800 FORMAT(//29H NO FEASIBLE SOLUTION EXISTS.) 810 FORMAT(//4H ***,I6,38H INPUT ERRORS TO SUBROUTINE CAPROX ***) 820 FORMAT(22H ACTIVE CONSTRAINTS. ,I6,22H OF APPROXIMATION TYPE/ + 22H ACTIVE CONSTRAINTS. ,I6,22H OF SIMPLE BOUND TYPE/ + 22H ACTIVE CONSTRAINTS. ,I6,22H OF GENERAL BOUND TYPE) END C SUBROUTINE CPAIRS PRINTS THE BASIS NAMES SO THAT C NATURAL PAIRINGS IN OPTIMAL BASES ARE MADE C IMMEDIATELY APPARENT. NONE OF THE SUBROUTINE C ARGUMENTS ARE ALTERED IN ANY WAY UPON RETURN. C ASTERISKS ARE PRINTED AFTER SUSPICIOUS DATA. C ALSO, EACH ACTIVE CONSTRAINT TYPE IS COUNTED. C C WRITTEN BY: ROY L. STREIT C NAVAL UNDERWATER SYSTEMS CENTER C NEW LONDON LABORATORY C NEW LONDON, CONNECTICUT 06320 C OCTOBER 29, 1982 C REVISED: AUGUST 16, 1983, MAY 10, 1985 C C ARGUMENT LIST: C IJSWT = SAME AS CAPROX ARGUMENT LIST C BASINV = SAME AS CAPROX ARGUMENT LIST C NDIM21 = SAME AS CAPROX ARGUMENT LIST C N2P1 = 2*N+1 C N2P2 = 2*N+2 C KJUMP = IOEXIT(R), WHERE R IS CURRENT DISCRETIZATION C STAGE. (SEE IOEXIT IN CAPROX ARGUMENT LIST.) C P = CURRENT NUMBER OF PROJECTION ANGLES C IPRT = SAME AS CAPROX INPUT ARGUMENT IOEXIT(1) C T = COUNT OF EACH ACTIVE CONSTRAINT TYPE. C LINEAR EQUALITIES ARE COUNTED AS ONE WHEN C THEY HAVE THE SAME FIRST TWO NAMES. C SUBROUTINE CPAIRS(IJSWT,BASINV,NDIM21,N2P1,N2P2,KJUMP,P,IPRT,T) C INTEGER IJSWT(3,*),N2P1,N2P2,J,K,J1,LOOK1,LOOK2,KJUMP,P,L,P1, + NOUT,NDIM21,BTYPE,ICOUNT,IPRT,T(3) DOUBLE PRECISION BASINV(NDIM21,1),DZERO DATA DZERO/+0.0D+0/ C STANDARD FILE NUMBER OF THE OUTPUT STREAM ... FOR VAX 11/780 DATA NOUT/6/ C P1=P-1 ICOUNT=0 DO 80 BTYPE=1,3 T(BTYPE)=0 DO 70 J=1,N2P1 LOOK1=IJSWT(1,J) IF(IABS(LOOK1).NE.BTYPE)GO TO 70 IF(LOOK1.GT.0)GO TO 10 IJSWT(1,J)=-IJSWT(1,J) GO TO 70 10 LOOK2=IJSWT(2,J) T(BTYPE)=T(BTYPE)+1 IF(IPRT.LT.1)GO TO 30 ICOUNT=ICOUNT+1 IF(BASINV(J,N2P2).GE.DZERO)GO TO 20 WRITE(NOUT,100)ICOUNT,LOOK1,LOOK2,IJSWT(3,J),BASINV(J,N2P2) GO TO 30 20 WRITE(NOUT,90)ICOUNT,LOOK1,LOOK2,IJSWT(3,J),BASINV(J,N2P2) 30 J1=J+1 IF(J1.GT.N2P1)GO TO 70 DO 60 K=J1,N2P1 IF(IJSWT(1,K).LT.0)GO TO 60 IF(LOOK1.NE.IJSWT(1,K).OR.LOOK2.NE.IJSWT(2,K))GO TO 60 IF(IPRT.LT.1)GO TO 50 ICOUNT=ICOUNT+1 L=IABS(IJSWT(3,J)-IJSWT(3,K)) IF((L.EQ.1.OR.L.EQ.P1).AND.(BASINV(K,N2P2).GE.DZERO)) + GO TO 40 WRITE(NOUT,100)ICOUNT,IJSWT(1,K),IJSWT(2,K),IJSWT(3,K), + BASINV(K,N2P2) GO TO 50 40 WRITE(NOUT,90)ICOUNT,IJSWT(1,K),IJSWT(2,K),IJSWT(3,K), + BASINV(K,N2P2) 50 IJSWT(1,K)=-IJSWT(1,K) IF(KJUMP.GT.(-1))GO TO 70 60 CONTINUE 70 CONTINUE 80 CONTINUE RETURN 90 FORMAT(7H BASIS.,I6,11H NAME = (,I2,1H,,I6,1H,,I6,10H) VALUE =, + D15.8) 100 FORMAT(7H BASIS.,I6,11H NAME = (,I2,1H,,I6,1H,,I6,10H) VALUE =, + D15.8,3H *) END C SUBROUTINE CEND STORES THE BEST COMPUTED SOLUTION C IN THE LAST COLUMN OF Z, AND SETS RELEVANT DATA C IN IOEXIT(LOGP), ITLOG(LOGP), AND K. NONE OF THE C OTHER SUBROUTINE ARGUMENTS ARE ALTERED. C C WRITTEN BY: ROY L. STREIT C NAVAL UNDERWATER SYSTEMS CENTER C NEW LONDON LABORATORY C NEW LONDON, CT 06320 C MARCH 9, 1983 C REVISED: JUNE 7, 1983, MAY 10, 1985 C C ARGUMENT LIST: C Z, CHEB, LOGP, N, NDIM, IOEXIT, ITLOG C = SAME AS SUBROUTINE CAPROX ARGUMENT LIST C LOGM1 = LOGP-1 C K = ON OUTPUT, THE COLUMN OF Z WHICH CONTAINS THE C BEST COMPUTED SOLUTION. (THE K-TH COLUMN IS C IDENTICAL TO THE LAST COLUMN OF Z, ON RETURN.) C EXCEPTION: K=0 MEANS NO GOOD SOLUTION COMPUTED. C SUBROUTINE CEND(Z,CHEB,LOGP,LOGM1,N,NDIM,IOEXIT,ITLOG,K) C INTEGER I,K,LOGP,LOGM1,N,NDIM,IOEXIT(1),ITLOG(1) COMPLEX Z(NDIM,*) REAL CHEB(1) C K=0 DO 10 I=1,LOGM1 IF(IOEXIT(I).LT.0)GO TO 10 IF(IOEXIT(I).EQ.4)IOEXIT(LOGP)=4 IF(IOEXIT(I).EQ.4)GO TO 20 K=I IOEXIT(LOGP)=MAX0(IOEXIT(LOGP),IOEXIT(I)) 10 CONTINUE 20 IF(K.EQ.0)GO TO 40 ITLOG(LOGP)=ITLOG(K) CHEB(LOGP)=CHEB(K) DO 30 I=1,N Z(I,LOGP)=Z(I,K) 30 CONTINUE 40 RETURN END C SUBROUTINE PABS COMPUTES THE DISCRETIZED ABSOLUTE VALUE C OF A COMPLEX NUMBER AND ITS MINIMAL CLOCKWISE INDEX. C IN ADDITION, PABS INITIALIZES THE COSINES AND SINES C REQUIRED FOR ALL THE PROJECTIONS. C C WRITTEN BY: ROY L. STREIT C NAVAL UNDERWATER SYSTEMS CENTER C NEW LONDON LABORATORY C NEW LONDON, CONNECTICUT 06320 C MAY 28, 1982 C REVISED: APRIL 27, 1983, MAY 10, 1985 C C ARGUMENT LIST: C COMPLEX Z = THE GIVEN COMPLEX NUMBER C INTEGER LOGP = THE BASE 2 LOG OF THE NUMBER OF C PROJECTIONS, P=2**LOGP; LOGP.GE.1 C DOUBLE PNORM = THE MAXIMUM (OVER J=1,...,P) OF: C PRECISION REAL(Z)*COSSIN(1,J)+AIMAG(Z)*COSSIN(2,J) C INTEGER PARG = THE MINIMAL INDEX J FOR WHICH PNORM IS C ATTAINED WHEN TIES ARE RESOLVED IN THE C CLOCKWISE DIRECTION. IN PARTICULAR, C PARG = 1 FOR Z = ZERO. C DOUBLE COSSIN(2,P+1) = THE COSINES AND SINES OF THE ANGLES C PRECISION ASSOCIATED WITH THE PROJECTIONS. THEY C ARE INITIALIZED IN THIS ROUTINE AND MUST C NOT BE ALTERED IN THE CALLING ROUTINE. C COSSIN(1,J) = COS(J-TH PROJECTION ANGLE) C COSSIN(2,J) = SIN(J-TH PROJECTION ANGLE) C NOTE: THE FIRST AND (P+1)-ST PROJECTIONS C ARE IDENTICAL; I.E., BOTH ANGLES = 0. C C COMMENTS: C 1. THE PROJECTION ANGLES ARE TAKEN TO BE EQUISPACED AROUND C THE CIRCLE BEGINNING ON THE POSITIVE X-AXIS. C 2. FOR ALL LOGP .GT. 2, THIS SUBROUTINE REQUIRES THE C CALCULATION OF AT MOST 2*LOGP - 5 MULTIPLICATIONS. C NO MULTIPLICATIONS ARE REQUIRED FOR LOGP = 1,2. C SUBROUTINE PABS(Z,LOGP,PNORM,PARG,COSSIN) COMPLEX Z DOUBLE PRECISION PNORM,COSSIN(2,*) INTEGER LOGP,PARG DOUBLE PRECISION X,Y,PTEST1,PTEST2,DZERO,DONE,PI,COS45 INTEGER P,P2,P4,P8,IDEL,I1,I2,P34,STORPN,P18,P28, + P38,P48,P58,P68,P78,P88 DATA DZERO,DONE/+0.0D+0,+1.0D+0/ DATA STORPN/-1/ C C MACHINE DEPENDENT CONSTANT C DATA PI/+3.141592653589793238D+0/ IF(STORPN.LT.0)PI=4.0D0*DATAN(DONE) C IF(LOGP.GT.1)GO TO 5 C C THE SPECIAL CASE P=2 PNORM=DBLE(ABS(REAL(Z))) PARG=1 IF(DBLE(REAL(Z)).LT.DZERO)PARG=2 IF(DBLE(REAL(Z)).EQ.DZERO.AND.DBLE(AIMAG(Z)).LT.DZERO)PARG=2 IF(STORPN.EQ.1)GO TO 330 STORPN=1 COSSIN(1,1)= DONE COSSIN(2,1)= DZERO COSSIN(1,2)=-DONE COSSIN(2,2)= DZERO COSSIN(1,3)= DONE COSSIN(2,3)= DZERO GO TO 330 C C INITIALIZE IF NECESSARY 5 IF(LOGP.EQ.STORPN)GO TO 40 STORPN=LOGP IF(LOGP.EQ.2)GO TO 20 P=2**LOGP P2=P/2 P4=P2/2 P8=P4/2 P34=P2+P4 P18= 1 +P8 P28=P18+P8 P38=P28+P8 P48=P38+P8 P58=P48+P8 P68=P58+P8 P78=P68+P8 P88=P78+P8 C COMPUTE COSINES AND SINES OF ALL PROJECTIONS X=DBLE(FLOAT(P2)) PTEST1=PI/X DO 10 I1=2,P18 I2=I1-1 X=DBLE(FLOAT(I2))*PTEST1 Y=DSIN(X) X=DCOS(X) C ...FIRST QUADRANT COSSIN(1,I1 )= X COSSIN(2,I1 )= Y IDEL=P28-I2 COSSIN(1,IDEL)= Y COSSIN(2,IDEL)= X C ... SECOND QUADRANT IDEL=I1+P4 COSSIN(1,IDEL)=-Y COSSIN(2,IDEL)= X IDEL=P48-I2 COSSIN(1,IDEL)=-X COSSIN(2,IDEL)= Y C ... THIRD QUADRANT IDEL=I1+P2 COSSIN(1,IDEL)=-X COSSIN(2,IDEL)=-Y IDEL=P68-I2 COSSIN(1,IDEL)=-Y COSSIN(2,IDEL)=-X C ... FOURTH QUADRANT IDEL=I1+P34 COSSIN(1,IDEL)= Y COSSIN(2,IDEL)=-X IDEL=P88-I2 COSSIN(1,IDEL)= X COSSIN(2,IDEL)=-Y 10 CONTINUE COS45=COSSIN(1,P18) GO TO 30 20 P=4 P28=2 P48=3 P68=4 P88=5 30 COSSIN(1,1 )= DONE COSSIN(2,1 )= DZERO COSSIN(1,P28)= DZERO COSSIN(2,P28)= DONE COSSIN(1,P48)=-DONE COSSIN(2,P48)= DZERO COSSIN(1,P68)= DZERO COSSIN(2,P68)=-DONE COSSIN(1,P88)= DONE COSSIN(2,P88)= DZERO C C LOCATE THE COMPLEX NUMBER IN THE PLANE 40 X=DBLE( REAL(Z)) Y=DBLE(AIMAG(Z)) IF(Y.GT.DZERO)GO TO 50 IF(Y.LT.DZERO)GO TO 60 IF(X.GT.DZERO)GO TO 70 IF(X.LT.DZERO)GO TO 80 C ... AT THE ORIGIN PNORM=DZERO PARG=1 GO TO 330 50 IF(X.GT.DZERO)GO TO 90 IF(X.LT.DZERO)GO TO 110 C ... ON THE POSITIVE Y-AXIS PNORM=Y PARG=P28 GO TO 330 60 IF(X.LT.DZERO)GO TO 130 IF(X.GT.DZERO)GO TO 150 C ... ON THE NEGATIVE Y-AXIS PNORM=-Y PARG=P68 GO TO 330 C ... ON THE POSITIVE X-AXIS 70 PNORM=X PARG=1 GO TO 330 C ... ON THE NEGATIVE X-AXIS 80 PNORM=-X PARG=P48 GO TO 330 C C SET SEARCH LIMITS AND FIRST TWO PROJECTIONS IF P.NE.4 C ... IN THE FIRST QUADRANT 90 IF(Y.GT.X)GO TO 100 IF(P.EQ.4)GO TO 170 I1=1 I2=P18 PTEST1=X PTEST2=(X+Y)*COS45 GO TO 250 100 IF(P.EQ.4)GO TO 180 I1=P18 I2=P28 PTEST1=(X+Y)*COS45 PTEST2=Y GO TO 250 C ...IN THE SECOND QUADRANT 110 IF((-X).GT.Y)GO TO 120 IF(P.EQ.4)GO TO 190 I1=P28 I2=P38 PTEST1=Y PTEST2=(Y-X)*COS45 GO TO 250 120 IF(P.EQ.4)GO TO 200 I1=P38 I2=P48 PTEST1=(Y-X)*COS45 PTEST2=-X GO TO 250 C ... IN THE THIRD QUADRANT 130 IF((-Y).GT.(-X))GO TO 140 IF(P.EQ.4)GO TO 210 I1=P48 I2=P58 PTEST1=-X PTEST2=-(X+Y)*COS45 GO TO 250 140 IF(P.EQ.4)GO TO 220 I1=P58 I2=P68 PTEST1=-(X+Y)*COS45 PTEST2=-Y GO TO 250 C ... IN THE FOURTH QUADRANT 150 IF(X.GT.(-Y))GO TO 160 IF(P.EQ.4)GO TO 230 I1=P68 I2=P78 PTEST1=-Y PTEST2=(X-Y)*COS45 GO TO 250 160 IF(P.EQ.4)GO TO 240 I1=P78 I2=P88 PTEST1=(X-Y)*COS45 PTEST2=X GO TO 250 C C FINISH THE SPECIAL CASE P = 4 170 PNORM=X PARG=1 GO TO 330 180 PNORM=Y PARG=2 GO TO 330 190 PNORM=Y PARG=2 GO TO 330 200 PNORM=-X PARG=3 GO TO 330 210 PNORM=-X PARG=3 GO TO 330 220 PNORM=-Y PARG=4 GO TO 330 230 PNORM=-Y PARG=4 GO TO 330 240 PNORM=X PARG=1 GO TO 330 C C USE BISECTION TO GET LARGEST PROJECTION AND MINIMAL INDEX 250 IDEL=P8 260 IF(PTEST1.GT.PTEST2)GO TO 270 IF(PTEST1.LT.PTEST2)GO TO 290 GO TO 310 270 IF(IDEL.GT.1)GO TO 280 PNORM=PTEST1 PARG=I1 GO TO 330 280 IDEL=IDEL/2 I2=I1+IDEL PTEST2=X*COSSIN(1,I2)+Y*COSSIN(2,I2) GO TO 260 290 IF(IDEL.GT.1)GO TO 300 PNORM=PTEST2 IF(I2.EQ.P88)I2=1 PARG=I2 GO TO 330 300 IDEL=IDEL/2 I1=I1+IDEL PTEST1=X*COSSIN(1,I1)+Y*COSSIN(2,I1) GO TO 260 310 IF(IDEL.GT.1)GO TO 320 PNORM=PTEST1 PARG=I1 GO TO 330 320 PARG=(I1+I2)/2 PNORM=X*COSSIN(1,PARG)+Y*COSSIN(2,PARG) C C UNIQUE PROGRAM EXIT 330 RETURN END PROGRAM KPROG C KPROG IS A MAIN ROUTINE WHICH COMPUTES C EXAMPLES FOR ACM ALGORITHM KAPROX. C AUGUST 16, 1983 C REVISED: JULY 30, 1984, MAY 10, 1985 C INTEGER P1 PARAMETER (NDIM=8,NDIMP1=NDIM+1,NDIMP4=NDIM+4) PARAMETER (MDIM=285) PARAMETER (L=2,LDIM=2) PARAMETER (LOGP=10,P1=2**LOGP+1) C INTEGER IJSWT(3,NDIMP1),IOEXIT(LOGP),ITLOG(LOGP),ICOUNT(NDIM) REAL C(NDIM),D(LDIM),CHEB(LOGP),ESTORE(NDIM), + Z(NDIM,LOGP),H(NDIM) DOUBLE PRECISION BASINV(NDIMP1,NDIMP4),COSSIN(2,P1), + PI,DELTA,DUM COMPLEX A(NDIM,MDIM),F(MDIM),B(NDIM,LDIM), + G(LDIM),DATA(MDIM) C C STANDARD FILE NUMBER OF THE OUTPUT STREAM ... FOR VAX 11/780 DATA NOUT/6/ C PI=4.0D0*DATAN(1.0D0) C C DEFINE THE APPROXIMATION DOMAIN C ... GENERATE MCIRC EQUISPACED DATA POINTS ON A C ... CIRCLE OF RADIUS 1 AND CENTER AT CMPLX(0,2). C ... THE FIRST POINT IS CMPLX(0,3). MCIRC=125 DELTA=2.D0*PI/FLOAT(MCIRC) DATA(1)=CMPLX(0.,3.) DO 10 K=2,MCIRC DUM=(PI/2.D0)+FLOAT(K-1)*DELTA SX=DCOS(DUM) SY=2.D0+DSIN(DUM) DATA(K)=CMPLX(SX,SY) 10 CONTINUE C ... GENERATE MSQUAR EQUISPACED DATA POINTS ON A SQUARE WITH C ... CORNERS AT CMPLX(1,-1), CMPLX(-1,-1), CMPLX(-1,-3), AND C ... CMPLX(1,-3). ASSUMES MSQUAR IS DIVISIBLE BY 4. MSIDE=40 MSQUAR=4*MSIDE DELTA=2.D0/FLOAT(MSIDE) DO 20 K=1,MSIDE DUM=FLOAT(K-1)*DELTA K1=K+MCIRC DELTAK=1.D0-DUM DATA(K1)=CMPLX(DELTAK,-1.) K1=K1+MSIDE DELTAK=-1.D0-DUM DATA(K1)=CMPLX(-1.,DELTAK) K1=K1+MSIDE DELTAK=-1.D0+DUM DATA(K1)=CMPLX(DELTAK,-3.) K1=K1+MSIDE DELTAK=-3.D0+DUM DATA(K1)=CMPLX(1.,DELTAK) 20 CONTINUE C ... SET THE TOTAL NUMBER OF POINTS M=MCIRC+MSQUAR C C DEFINE THE PROBLEM DATA FOR THE LARGEST ORDER C ... THE MATRIX A AND THE VECTOR F DO 40 I=1,M A(1,I)=DATA(I) DO 30 J=2,NDIM A(J,I)=DATA(I)**J 30 CONTINUE F(I)=CMPLX(1.0,0.0) 40 CONTINUE C ... THE MATRIX B AND THE VECTORS G AND D DO 50 I=1,NDIM B(I,1)=CMPLX(FLOAT(I),0.0) B(I,2)=CMPLX(FLOAT(I*(I-1)),0.0) 50 CONTINUE G(1)=CMPLX(0.0,0.0) G(2)=CMPLX(0.0,0.0) D(1)=.75 D(2)=.75 C ... SIMPLE BOUNDS WITH VECTORS H AND C DO 60 I=1,NDIM H(I)=0.0 C(I)=.2 60 CONTINUE C C SOLVE ALL PROBLEMS FROM ORDER 1 THROUGH NDIM DO 70 K=1,NDIM C ... GET CURRENT TIMING AND PAGE FAULT DATA: VAX 11/780 C CALL GETJPI(KTIME1,KPAGE1) C ... SET THE PROBLEM ORDER N=K C ... SET MAXIMUM ITERATION COUNT ITLOG(1)=100*(N+1) C ... SET INTERNAL PRINT-OUT OPTION IOEXIT(1)=-1 C ... SOLVE THE PROBLEM OF ORDER N CALL KAPROX(A,Z,F,B,N,NDIM,M,L,LOGP,H,C,G, + D,ITLOG,NDIMP1,IJSWT,BASINV,COSSIN,CHEB,IOEXIT) C ... GET CURRENT TIMING AND PAGE FAULT DATA: VAX 11/780 C CALL GETJPI(KTIME2,KPAGE2) C SX=FLOAT(KTIME2-KTIME1)/100. C KP=KPAGE2-KPAGE1 C WRITE(NOUT,80)N,SX,KP C ... STORE SUMMARY RESULTS ESTORE(N)=CHEB(LOGP) ICOUNT(N)=ITLOG(LOGP) C ... PRINT SOLUTION VECTOR Z WRITE(NOUT,90) WRITE(NOUT,100)(J,Z(J,LOGP),J=1,N) WRITE(NOUT,130) 70 CONTINUE C C PRINT SUMMARY TABLE WRITE(NOUT,110) WRITE(NOUT,120)(J,ESTORE(J),ICOUNT(J),J=1,NDIM) C 80 FORMAT(9H ORDER = ,I3,16H, TIME IN SEC = ,F8.2, C + 16H, PAGE FAULTS = ,I7) 90 FORMAT(25H COMPONENT REAL COEF) 100 FORMAT(I10,F15.8) 110 FORMAT(//40H ORDER CHEBYSHEV ERROR ITERATIONS) 120 FORMAT(I10,F15.8,I10) 130 FORMAT(//////) END C SUBROUTINE KAPROX SOLVES SYSTEMS OF COMPLEX LINEAR EQUATIONS C IN THE CHEBYSHEV NORM WITH CONSTRAINTS ON THE UNKNOWNS, WHEN C THE SOLUTION VECTOR IS REQUIRED TO BE REAL-VALUED. C C =============================================================== C PROBLEM: MIN (OVER Z) THE MAX (OVER K=1,...,M) OF C C ABS( Z * A(K) - F(K) ) C C SUBJECT TO: C C ABS( Z(K) - H(K) ) .LE. C(K) FOR K=1,...,N C C ABS( Z * B(K) - G(K) ) .LE. D(K) FOR K=1,...,L C C WHERE A(K) AND B(K) DENOTE COLUMN K OF MATRICES A AND B, RESP., C AND WHERE Z(K) DENOTES COMPONENT K OF THE ROW VECTOR Z. C =============================================================== C C WRITTEN BY: ROY L. STREIT C NAVAL UNDERWATER SYSTEMS CENTER C NEW LONDON LABORATORY C NEW LONDON, CONNECTICUT 06320 C NOVEMBER 10, 1982 C REVISED: AUGUST 16, 1983, JULY 31, 1984, MAY 10, 1985 C C C REFERENCES: C 1. R.L. STREIT AND A.H. NUTTALL, "LINEAR CHEBYSHEV COMPLEX C FUNCTION APPROXIMATION," NUSC TECHNICAL REPORT 6403, C NAVAL UNDERWATER SYSTEMS CENTER, NEW LONDON LABORATORY, C NEW LONDON, CONNECTICUT, 26 FEBRUARY 1981. C 2. R.L. STREIT AND A.H. NUTTALL, "A GENERAL CHEBYSHEV C COMPLEX FUNCTION APPROXIMATION PROCEDURE AND AN APPLICATION C TO BEAMFORMING," JOURNAL OF THE ACOUSTICAL SOCIETY OF C AMERICA, VOL. 72(1), JULY, 1982, PP.181-190. C 3. R.L. STREIT AND A.H. NUTTALL, "A NOTE ON THE SEMI-INFINITE C PROGRAMMING APPROACH TO COMPLEX APPROXIMATION," MATHEMATICS C OF COMPUTATION, 40(1983), 599-605. C 4. R.L. STREIT, "SOLUTION OF SYSTEMS OF COMPLEX LINEAR C EQUATIONS IN THE L-SUB-INFINITY NORM WITH CONSTRAINTS C ON THE UNKNOWNS," SIAM JOURNAL ON SCIENTIFIC AND STATISTICAL C COMPUTATION, TO APPEAR (VOL 7, NUMER 1, JANUARY 1986). C C C ALGORITHM SKETCH: C 1. SET P=4. THIS IS THE START OF DISCRETIZATION STAGE R:=1. C 2. APPROXIMATE THE ABSOLUTE VALUES OF ALL COMPLEX NUMBERS C WITH P PROJECTIONS; I.E., DISCRETIZE ABSOLUTE VALUES. C 3. SOLVE THE DUAL OF THE RESULTING PRIMAL LINEAR PROGRAM C USING THE REVISED SIMPLEX METHOD WITH EXPLICIT INVERSE. C (A) THE STORAGE REQUIRED IS REDUCED BY A FACTOR OF P C BY CONSTRUCTING THE INCOMING COLUMNS AS NEEDED. C (B) THE NUMBER OF INNER PRODUCTS REQUIRED IN EACH C SIMPLEX STEP TO COMPUTE THE REDUCED COSTS IS C DECREASED BY A FACTOR OF P BY UTILIZING THE C INTRINSIC COMPLEX ARITHMETIC PROBLEM STRUCTURE. C 4. DOUBLE THE NUMBER P OF PROJECTIONS AND LOOP BACK TO C STEP 2 TO START DISCRETIZATION STAGE R:=R+1. C 5. NOTE: THE TOTAL NUMBER OF STAGES NORMALLY EQUALS LOGP-1 C C C ARGUMENT LIST: C COMPLEX A(NDIM,*) = COEFFICIENT MATRIX IS CONTAINED IN THE C FIRST N ROWS. MUST HAVE .GE. M COLUMNS C REAL Z(NDIM,*) = ON OUTPUT, CONTAINS THE SOLUTION VECTORS. C THE R-TH COLUMN IS SOLUTION VECTOR OF THE C R-TH STAGE. MUST HAVE .GE. LOGP COLUMNS. C ON OUTPUT, THE BEST COMPUTED SOLUTION IS C STORED IN THE LAST (I.E., LOGP) COLUMN. C COMPLEX F(*) = CONSTANT VECTOR. MUST BE OF LENGTH .GE. M C COMPLEX B(NDIM,*) = GENERAL BOUND MATRIX COEFFICIENTS ARE C CONTAINED IN THE FIRST N ROWS. C MUST HAVE AT LEAST ONE COLUMN WHEN L=0. C INTEGER N = NUMBER OF UNKNOWNS. REQUIRES N.GE.1 C INTEGER NDIM = NUMBER OF ROWS OF MATRICES A, Z, AND B IN C THE DIMENSION STATEMENT IN THE CALLING C PROGRAM. REQUIRES NDIM .GE. N C INTEGER M = NUMBER OF EQUATIONS. REQUIRES M.GE.1 C INTEGER L = NUMBER OF GENERAL BOUND CONSTRAINTS; L.GE.0 C INTEGER LOGP = THE BASE 2 LOG OF THE MAXIMUM NUMBER OF C PROJECTIONS: P=2**LOGP. REQUIRES LOGP.GE.1 C WARNING: LOGP=1 WORKS ONLY IF ALL THE C PROBLEM DATA ARE REAL-VALUED. C REAL H(*) = CENTERS OF INTERVALS BOUNDING Z COMPONENTS C MUST BE OF LENGTH .GE. N C REAL C(*) = RADII OF INTERVALS BOUNDING Z COMPONENTS C REQUIRES C(K) .GT. 0 FOR ALL K. C MUST BE OF LENGTH .GE. N C COMPLEX G(*) = CONSTANT TERM OF GENERAL BOUND CONSTRAINTS. C MUST BE OF LENGTH .GE. MAX(1,L) C REAL D(*) = RADII OF GENERAL BOUND CONSTRAINTS. C REQUIRES D(K) .GT. 0 FOR ALL K. C MUST BE OF LENGTH .GE. MAX(1,L) C INTEGER ITLOG(*) = INPUT/OUTPUT ARRAY OF LENGTH .GE. LOGP C INPUT: ITLOG(1) = MAXIMUM ITERATION COUNT C OUTPUT: ITLOG(R) = ITERATION COUNT TO C FINISH STAGE R. SPECIAL CASE: C ITLOG(LOGP)=FINAL ITERATION COUNT. C INTEGER NDIMP1 = NDIM + 1 C INTEGER IJSWT(3,*)= CONTAINS THE BASIC VARIABLE NAMES. C MUST HAVE .GE. (N+1) COLUMNS. C THE J-TH BASIC VARIABLE IS SPECIFIED BY: C IZSWT(1,J)=1, 2, OR 3 IF AN S, W, OR T C IJSWT(2,J)=EQUATION/CONSTRAINT NUMBER C IJSWT(3,J)=PROJECTION NUMBER C DOUBLE BASINV(NDIMP1,*) = MUST HAVE .GE. (N+4) COLUMNS. C PRECISION IDENTIFIED BY COLUMNS AS FOLLOWS: C COL. 1 TO N+1=CURRENT BASIS INVERSE C COL. N+1 + 1 =CURRENT B.F.S. C COL. N+1 + 2 =INCOMING COL OF NEW B.F.S. C COL. N+1 + 3 =COSTS OF CURRENT B.F.S. C DOUBLE COSSIN(2,*)= COSINES AND SINES OF THE P PROJECTION C PRECISION ANGLES. DEFINED IN ROUTINE PABS. C MUST HAVE .GE. (P+1) COLUMNS C REAL CHEB(*) = AT END OF STAGE R, CHEBYSHEV MULTIPLIER C IS CHEB(R). MUST BE OF LENGTH .GE. LOGP C SPECIAL CASE: CHEB(LOGP)=CHEBYSHEV MULTI- C PLIER OF BEST SOLUTION. C INTEGER IOEXIT(*) = INPUT/OUTPUT ARRAY OF LENGTH .GE. LOGP C ..........INPUT: SETS PRINT OPTION ACCORDING TO C IOEXIT(1) = -1 IF ONLY BEST SOLUTION IS PRINTED C = 0 IF NOTHING IS TO BE PRINTED C = 1 IF ONLY B.F.S. AND SIMPLEX MULTIPLIERS C AT END OF EACH STAGE ARE PRINTED C = 2 IF INCOMING AND OUTGOING BASIC VARIABLES C ARE ALSO PRINTED C ..........OUTPUT: AT THE END OF THE R-TH STAGE, C IOEXIT(R) = 4 IF MAXIMUM ITERATION COUNT REACHED C = 3 IF CHEBYSHEV ERROR ESTIMATE DECREASED C THREE TIMES IN A ROW C = 2 IF INCOMING AND OUTGOING BASIC VARIABLES C ARE IDENTICAL C = 1 IF ALL REDUCED COSTS ARE .GE. -TOLP C = 0 IF ALL REDUCED COSTS ARE .GE. ZERO C = -1 IF PRIMAL PROBLEM IS INFEASIBLE C = -2 IF IMPOSSIBLE BASIS FOUND TWICE IN A ROW C = -3 IF SINGULAR BASIS IS DETECTED C = -10000-J IF J INPUT ERRORS WERE FOUND. C (IN THIS CASE, R=1 AND PROGRAM EXITS.) C = -999 IF THE R-TH STAGE WAS NOT ATTEMPTED C SPECIAL CASE: IOEXIT(LOGP) = MAXIMUM (OVER ALL C STAGES R) OF THE NUMBERS IOEXIT(R). C C THE CONTENTS OF THE FOLLOWING ARRAYS ARE NOT ALTERED BY THIS C SUBROUTINE... A, F, B, H, C, G, D. C SUBROUTINE KAPROX(A,Z,F,B,N,NDIM,M,L,LOGP,H,C, + G,D,ITLOG,NDIMP1,IJSWT, + BASINV,COSSIN,CHEB,IOEXIT) C C TYPE AND DIMENSION ALL SUBROUTINE ARGUMENTS INTEGER N,NDIM,M,L,LOGP,ITLOG(1),NDIMP1,IJSWT(3,*),IOEXIT(1) COMPLEX A(NDIM,*),F(1),B(NDIM,*),G(1) REAL Z(NDIM,*),H(1),C(1),D(1),CHEB(1) DOUBLE PRECISION COSSIN(2,*),BASINV(NDIMP1,*) C C TYPE ALL INTERNAL VARIABLE NAMES COMPLEX ZDUM REAL ZERO DOUBLE PRECISION X,Y,DZERO,DONE,DRANGE,DRELPR,CHEBBO, + CHEBER,RCC,TOLP,ZR,ZI,PNORM INTEGER I,J,K,PARG,KPT,JPT,IWHICH,ITER,ITMAX,LOOP,R, + NOUT,ISTAR,KSTAR,JSTAR,IMPOSS,IPRT,SLACK,ICHEB, + CYCLE,NP1,NP2,NP3,NP4,LOGM1,T(3),TSAVE(3) C DATA DZERO,DONE/+0.0D+0,+1.0D+0/ DATA ZERO/+0.0E+0/ C C MACHINE DEPENDENT CONSTANTS --- THESE ARE FOR THE DEC VAX 11/780 C ... LARGEST DOUBLE PRECISION NUMBER DATA DRANGE/+1.0D+38/ C ... RELATIVE PRECISION OF SINGLE PRECISION ARITHMETIC DATA DRELPR/+0.6D-7/ C ... STANDARD FILE NUMBER OF THE OUTPUT STREAM DATA NOUT/6/ C C CHECK INPUT ARGUMENTS J=0 IF(LOGP.LT.1)J=J+1 IF(M.LT.1)J=J+1 IF(N.LT.1)J=J+1 IF(L.LT.0)J=J+1 IF(NDIMP1.NE.(NDIM+1))J=J+1 IF(ITLOG(1).LT.1)J=J+1 IF(IOEXIT(1).LT.(-1).OR.IOEXIT(1).GT.2)J=J+1 DO 10 I=1,N IF(C(I).LE.ZERO)J=J+1 10 CONTINUE IF(L.LE.0)GO TO 30 DO 20 I=1,L IF(D(I).LE.ZERO)J=J+1 20 CONTINUE 30 IF(J.GT.0)GO TO 710 C C INITIALIZE IPRT=IOEXIT(1) ITMAX=ITLOG(1) ITER=-1 NP1=N+1 NP2=NP1+1 NP3=NP2+1 NP4=NP3+1 LOGM1=LOGP-1 IF(LOGM1.LT.1)LOGM1=1 CHEBBO=DZERO SLACK=1 CYCLE=1 DO 50 I=1,LOGP DO 40 K=1,N Z(K,I)=ZERO 40 CONTINUE ITLOG(I)=0 CHEB(I)=ZERO IOEXIT(I)=-999 50 CONTINUE T(1)=0 T(2)=0 T(3)=0 C C DEFINE THE FIRST BASIS INVERSE AND B.F.S. DO 80 I=1,NP3 DO 70 J=1,NP1 IF(I.NE.J)GO TO 60 BASINV(J,J)=DONE GO TO 70 60 BASINV(J,I)=DZERO 70 CONTINUE 80 CONTINUE BASINV(NP1,NP2)=DONE C C IDENTIFY THE B.F.S. C ... THE W VARIABLES DO 90 K=1,N IJSWT(1,K)=2 IJSWT(2,K)=K IJSWT(3,K)=1 90 CONTINUE C ... THE SLACK VARIABLE IJSWT(1,NP1)=0 IJSWT(2,NP1)=0 IJSWT(3,NP1)=0 C C DEFINE THE COST COEFFICIENTS C ... OF THE W VARIABLES DO 100 K=1,N BASINV(K,NP4)=DBLE(C(K))+DBLE(H(K)) 100 CONTINUE C ... OF THE SLACK VARIABLE BASINV(NP1,NP4)=DZERO C C ========================= C BEGIN DISCRETIZATION LOOP C ========================= C DO 690 R=1,LOGM1 LOOP=R+1 IMPOSS=0 ICHEB=0 IF(LOOP.EQ.2)GO TO 110 IF(IOEXIT(R-1).GT.2)GO TO 150 IF(IOEXIT(R-1).LT.3)GO TO 180 C C --------------------------------------------------------- C BEGIN THE REVISED SIMPLEX ALGORITHM WITH EXPLICIT INVERSE C --------------------------------------------------------- C C COMPUTE THE CHEBYSHEV SIMPLEX MULTIPLIER 110 CHEBER=DZERO DO 120 J=1,NP1 CHEBER=CHEBER-BASINV(J,NP4)*BASINV(J,NP1) 120 CONTINUE C C TESTS FOR MINIMUM INCREMENT TERMINATION IF(SLACK.GT.0)GO TO 150 IF(CHEBER.LT.CHEBBO)GO TO 130 GO TO 140 130 ICHEB=ICHEB+1 IF(ICHEB.LT.4)GO TO 150 GO TO 630 140 ICHEB=0 150 RCC=CHEBER-CHEBBO CHEBBO=CHEBER C C COMPUTE THE REAL SIMPLEX MULTIPLIERS DO 170 I=1,N X=DZERO DO 160 J=1,NP1 X=X+BASINV(J,NP4)*BASINV(J,I) 160 CONTINUE Z(I,R)=SNGL(X) 170 CONTINUE C C TEST FOR MAXIMUM ITERATION COUNT ITER=ITER+1 IF(ITER.GT.ITMAX)GO TO 580 C C PRINT CHEBYSHEV MULTIPLIER THIS ITERATION, IF REQUESTED IF(IPRT.LT.2)GO TO 180 WRITE(NOUT,730)LOOP,ITER,CHEBER,RCC C C DETERMINE THE ENTERING BASIC VARIABLE := (IWHICH,KPT,JPT). C TIES ARE BROKEN BY CHOOSING THE MOST NEGATIVE REDUCED C COST COEFFICIENT WITH THE SMALLEST INDEX. EXCEPTIONS: C BRING IN AN S VARIABLE IF POSSIBLE WHEN SLACK VARIABLE C IS STILL IN THE BASIS. THIS SHOULD AVOID CYCLING. 180 IWHICH=0 RCC=DRANGE C ... CHECK THE S VARIABLES DO 210 K=1,M ZDUM=CMPLX(ZERO,ZERO) DO 190 J=1,N ZDUM=ZDUM+CMPLX(Z(J,R)*REAL(A(J,K)), + Z(J,R)*AIMAG(A(J,K))) 190 CONTINUE ZDUM=ZDUM-F(K) CALL PABS(ZDUM,LOOP,PNORM,PARG,COSSIN) X=-PNORM+CHEBER IF(X.GE.RCC)GO TO 210 IF(IMPOSS.NE.1)GO TO 200 IF(K.NE.KSTAR)GO TO 200 IF(ISTAR.EQ.1)GO TO 210 200 JPT=PARG KPT=K RCC=X IWHICH=1 210 CONTINUE IF(CYCLE.EQ.0)GO TO 220 IF(RCC.GT.(-DRELPR))CYCLE=0 IF(CYCLE.EQ.1)GO TO 290 C ... CHECK THE W VARIABLES 220 DO 240 K=1,N ZDUM=CMPLX(Z(K,R)-H(K),ZERO) CALL PABS(ZDUM,LOOP,PNORM,PARG,COSSIN) X=-PNORM+C(K) IF(X.GE.RCC)GO TO 240 IF(IMPOSS.NE.1)GO TO 230 IF(K.NE.KSTAR)GO TO 230 IF(ISTAR.EQ.2)GO TO 240 230 JPT=PARG KPT=K RCC=X IWHICH=2 240 CONTINUE IF(L.EQ.0)GO TO 280 C ... CHECK THE T VARIABLES DO 270 K=1,L ZDUM=CMPLX(ZERO,ZERO) DO 250 J=1,N ZDUM=ZDUM+CMPLX(Z(J,R)*REAL(B(J,K)), + Z(J,R)*AIMAG(B(J,K))) 250 CONTINUE ZDUM=ZDUM-G(K) CALL PABS(ZDUM,LOOP,PNORM,PARG,COSSIN) X=-PNORM+D(K) IF(X.GE.RCC)GO TO 270 IF(IMPOSS.NE.1)GO TO 260 IF(K.NE.KSTAR)GO TO 260 IF(ISTAR.EQ.3)GO TO 270 260 JPT=PARG KPT=K RCC=X IWHICH=3 270 CONTINUE C C TEST FOR TERMINATION 280 IF(IWHICH.GT.0)GO TO 290 IF(IMPOSS.EQ.0)GO TO 610 GO TO 560 290 IF(RCC.GE.DZERO.AND.CYCLE.EQ.0)GO TO 610 C C GENERATE INCOMING S, W, OR T COLUMN X=COSSIN(1,JPT) Y=COSSIN(2,JPT) IF(IWHICH.NE.2)GO TO 310 C ... FINAL FORM OF ANY W COLUMN DO 300 I=1,NP1 BASINV(I,NP3)=X*BASINV(I,KPT) 300 CONTINUE GO TO 380 C BASIS INVERSE TIMES ORIGINAL S OR T COLUMN 310 DO 320 I=1,NP1 BASINV(I,NP3)=DZERO 320 CONTINUE DO 360 J=1,N IF(IWHICH.GT.2)GO TO 330 ZR=DBLE(REAL(A(J,KPT)))*X+DBLE(AIMAG(A(J,KPT)))*Y GO TO 340 330 ZR=DBLE(REAL(B(J,KPT)))*X+DBLE(AIMAG(B(J,KPT)))*Y 340 DO 350 I=1,NP1 BASINV(I,NP3)=BASINV(I,NP3)+ZR*BASINV(I,J) 350 CONTINUE 360 CONTINUE C ... FINAL FORM OF ANY T COLUMN IF(IWHICH.GT.2)GO TO 380 C ... FINAL FORM OF ANY S COLUMN DO 370 I=1,NP1 BASINV(I,NP3)=BASINV(I,NP3)+BASINV(I,NP1) 370 CONTINUE C C DETERMINE NUMERICAL TOLERANCE FOR PIVOT AND REDUCED COST 380 X=DABS(BASINV(1,NP3)) DO 390 I=2,NP1 X=X+DABS(BASINV(I,NP3)) 390 CONTINUE TOLP=X*DRELPR C C REDUCED COST TOLERANCE TESTS FOR TERMINATION IF(TOLP.EQ.DZERO)GO TO 570 IF(RCC.LE.(-TOLP))GO TO 400 IF(CYCLE.EQ.0)GO TO 590 CYCLE=0 GO TO 220 C C DETERMINE THE OUTGOING BASIC VARIABLE := K-TH IN THE LIST. C TIES IN THE MINIMUM RATIO TEST ARE BROKEN BY SELECTING THE C MAXIMUM MAGNITUDE PIVOT WITH THE SMALLEST INDEX. C EXCEPTIONS: PIVOT OUT A W VARIABLE IF POSSIBLE WHEN C SLACK VARIABLE STILL IN BASIS. THIS SHOULD AVOID CYCLING. 400 X=DRANGE K=0 DO 440 I=1,NP1 IF(BASINV(I,NP3).LT.TOLP)GO TO 440 PNORM=BASINV(I,NP2)/BASINV(I,NP3) IF(PNORM.LT.DZERO)PNORM=DZERO IF(PNORM.GT.X)GO TO 440 IF(PNORM.LT.X)GO TO 420 IF(K.LT.1)GO TO 430 IF(SLACK.LT.1)GO TO 410 IF(IJSWT(1,K).NE.2.AND.IJSWT(1,I).EQ.2)GO TO 430 IF(IJSWT(1,K).NE.2.AND.IJSWT(1,I).NE.2)GO TO 410 IF(IJSWT(1,I).NE.2)GO TO 440 410 IF(BASINV(I,NP3).LE.Y)GO TO 440 GO TO 430 420 X=PNORM 430 Y=BASINV(I,NP3) K=I 440 CONTINUE C C TERMINATE IF SOLUTION IS UNBOUNDED BELOW IF(K.EQ.0)GO TO 600 C C UPDATE THE BASIC VARIABLE LIST IF(K.EQ.NP1)SLACK=0 IF(K.EQ.NP1)CYCLE=0 ISTAR=IJSWT(1,K) KSTAR=IJSWT(2,K) JSTAR=IJSWT(3,K) IJSWT(1,K)=IWHICH IJSWT(2,K)=KPT IJSWT(3,K)=JPT C ... TEST FOR IDENTICAL INCOMING AND OUTGOING BASIS VARIABLES IF(KSTAR.NE. KPT)GO TO 450 IF(JSTAR.NE. JPT)GO TO 450 IF(ISTAR.NE.IWHICH)GO TO 450 IF(CYCLE.EQ.0)GO TO 620 CYCLE=0 GO TO 220 C ... TERMINATE IF IMPOSSIBLE BASIS TWICE IN A ROW 450 J=0 DO 460 I=1,NP1 IF(IJSWT(2,I).NE.KPT)GO TO 460 IF(IJSWT(1,I).NE.IWHICH)GO TO 460 J=J+1 460 CONTINUE IF(J.LE.2)GO TO 470 IF(IMPOSS.EQ.1)GO TO 560 IMPOSS=1 CYCLE=0 IJSWT(1,K)=ISTAR IJSWT(2,K)=KSTAR IJSWT(3,K)=JSTAR ISTAR=IWHICH KSTAR=KPT GO TO 180 470 IMPOSS=0 C C PRINT INCOMING AND OUTGOING BASIC VARIABLE, IF REQUESTED IF(IPRT.LT.2)GO TO 480 WRITE(NOUT,740)LOOP,ITER,IJSWT(1,K),IJSWT(2,K),IJSWT(3,K), + ISTAR,KSTAR,JSTAR C C UPDATE COST COEFFICIENT 480 IF(IWHICH.GT.1)GO TO 490 BASINV(K,NP4)=DBLE( REAL(F(KPT)))*COSSIN(1,JPT) + +DBLE(AIMAG(F(KPT)))*COSSIN(2,JPT) GO TO 510 490 IF(IWHICH.GT.2)GO TO 500 BASINV(K,NP4)=DBLE(C(KPT))+DBLE(H(KPT))* + COSSIN(1,JPT)+DBLE(H(KPT))*COSSIN(2,JPT) GO TO 510 500 BASINV(K,NP4)=DBLE(D(KPT))+DBLE(REAL(G(KPT)))* + COSSIN(1,JPT)+DBLE(AIMAG(G(KPT)))*COSSIN(2,JPT) C C PIVOT TO UPDATE B.F.S. AND BASIS INVERSE 510 X=BASINV(K,NP3) KPT=K-1 KSTAR=K+1 DO 550 J=1,NP2 Y=BASINV(K,J)/X IF(K.EQ.1)GO TO 530 DO 520 I=1,KPT BASINV(I,J)=BASINV(I,J)-Y*BASINV(I,NP3) 520 CONTINUE 530 BASINV(K,J)=Y IF(K.EQ.NP1)GO TO 550 DO 540 I=KSTAR,NP1 BASINV(I,J)=BASINV(I,J)-Y*BASINV(I,NP3) 540 CONTINUE 550 CONTINUE C GO TO 110 C C ------------------------------------- C END OF REVISED SIMPLEX ALGORITHM LOOP C ------------------------------------- C C DEFINE ALL EXIT CONDITIONS FROM SIMPLEX ALGORITHM C ... IMPOSSIBLE BASIS; EITHER TWICE IN A ROW... FROM C NEAR 460, OR ONLY POSSIBLE PIVOT... FROM NEAR 290 560 IOEXIT(R)=-2 GO TO 640 C ... SINGULAR BASIS INVERSE, IMPOSSIBLE... FROM NEAR 390 570 IOEXIT(R)=-3 GO TO 640 C ... MAXIMUM ITERATION COUNT ACHIEVED... FROM NEAR 170 580 IOEXIT(R)=4 GO TO 640 C ... REDUCED COSTS EXCEED TOLERANCE... FROM NEAR 400 590 IOEXIT(R)=1 GO TO 640 C ... SOLUTION UNBOUNDED BELOW; PRIMAL INFEASIBLE...NEAR 440 600 IOEXIT(R)=-1 GO TO 640 C ... REDUCED COSTS ARE NON-NEGATIVE... FROM NEAR 280,290 610 IOEXIT(R)=0 GO TO 640 C ... IN- AND OUT-GOING BASIC VARIABLES THE SAME... NEAR 450 620 IOEXIT(R)=2 GO TO 640 C ... CONSECUTIVE DECREASES IN CHEBYSHEV MULT ... FROM NEAR 130 630 IOEXIT(R)=3 CHEB(R)=CHEBBO ITLOG(R)=ITER GO TO 650 640 CHEB(R)=CHEBER ITLOG(R)=ITER-1 C C PRINT B.F.S. AND SIMPLEX MULTIPLIERS, IF REQUESTED 650 K=IOEXIT(R) J=2**LOOP IF(LOGP.EQ.1)J=2 TSAVE(1)=T(1) TSAVE(2)=T(2) TSAVE(3)=T(3) CALL KPAIRS(IJSWT,BASINV,NDIMP1,NP1,NP2,K,J,IPRT,T) IF(IPRT.LT.1)GO TO 660 WRITE(NOUT,760)(K,Z(K,R),K=1,N) WRITE(NOUT,770)NP1,CHEB(R) WRITE(NOUT,780)J,LOOP,IOEXIT(R),ITLOG(R) C C TERMINATE IF AND ONLY IF THE SIMPLEX ALGORITHM ENDED BADLY 660 IF(IOEXIT(R).LT.0.OR.IOEXIT(R).GT.3)GO TO 700 C C INITIALIZE NEXT STAGE, UNLESS FINISHED IF(R.EQ.LOGM1)GO TO 700 C ... BASIS NAMES DO 670 K=1,NP1 IJSWT(3,K)=2*IJSWT(3,K)-1 670 CONTINUE C ... REAL SIMPLEX MULTIPLIERS DO 680 K=1,N Z(K,LOOP)=Z(K,R) 680 CONTINUE 690 CONTINUE C C ========================== C END OF DISCRETIZATION LOOP C ========================== C C STORE THE BEST COMPUTED SOLUTION AND RELEVANT DATA. C (ALTERNATIVELY, AT THIS POINT ONE COULD CALL A C SUBROUTINE DESIGNED TO IMPROVE THE BEST SOLUTION C PROVIDED BY THE DISCRETIZATION PROCESS.) 700 CALL KEND(Z,CHEB,LOGP,LOGM1,N,NDIM,IOEXIT,ITLOG,K) IF(IPRT.EQ.0)GO TO 720 IF(K.EQ.0)WRITE(NOUT,800) IF(K.EQ.0)GO TO 720 WRITE(NOUT,790) WRITE(NOUT,760)(I,Z(I,LOGP),I=1,N) WRITE(NOUT,770)NP1,CHEB(LOGP) IF(K.EQ.R)WRITE(NOUT,820)T(1),T(2),T(3) IF(K.EQ.(R-1))WRITE(NOUT,820)TSAVE(1),TSAVE(2),TSAVE(3) K=K+1 IF(LOGP.EQ.1)K=1 J=2**K WRITE(NOUT,780)J,K,IOEXIT(LOGP),ITLOG(LOGP) GO TO 720 C C INPUT ERRORS... FROM 30 710 IOEXIT(1)=-10000-J WRITE(NOUT,810)J C C UNIQUE PROGRAM EXIT 720 RETURN C 730 FORMAT(I4,I7,23H CHEBYSHEV MULTIPLIER =,D15.8,10H CHANGE =,D15.8) 740 FORMAT(I4,I7,6H NEW=(,I2,1H,,I6,1H,,I6,1H),9H OLD=(,I2, + 1H,,I6,1H,,I6,1H)) 750 FORMAT(17H BASIS. NAME = (,I2,1H,,I6,1H,,I6,10H) VALUE =,D15.8) 760 FORMAT(23H REAL MULTIPLIER (,I6,4H) = ,E15.8) 770 FORMAT(23H CHEBYSHEV MULTIPLIER (,I6,4H) = ,E15.8) 780 FORMAT(4H P =,I8,9H, LOG P =,I3,14H, IOEXIT(R) = ,I2, + 21H, ITERATION NUMBER = ,I7//) 790 FORMAT(//24H BEST COMPUTED SOLUTION.) 800 FORMAT(//29H NO FEASIBLE SOLUTION EXISTS.) 810 FORMAT(//4H ***,I6,38H INPUT ERRORS TO SUBROUTINE KAPROX ***) 820 FORMAT(22H ACTIVE CONSTRAINTS. ,I6,22H OF APPROXIMATION TYPE/ + 22H ACTIVE CONSTRAINTS. ,I6,22H OF SIMPLE BOUND TYPE/ + 22H ACTIVE CONSTRAINTS. ,I6,22H OF GENERAL BOUND TYPE) END C SUBROUTINE KPAIRS PRINTS THE BASIS NAMES SO THAT C NATURAL PAIRINGS IN OPTIMAL BASES ARE MADE C IMMEDIATELY APPARENT. NONE OF THE SUBROUTINE C ARGUMENTS ARE ALTERED IN ANY WAY UPON RETURN. C ASTERISKS ARE PRINTED AFTER SUSPICIOUS DATA. C ALSO, EACH ACTIVE CONSTRAINT TYPE IS COUNTED. C C WRITTEN BY: ROY L. STREIT C NAVAL UNDERWATER SYSTEMS CENTER C NEW LONDON LABORATORY C NEW LONDON, CONNECTICUT 06320 C NOVEMBER 10, 1982 C REVISED: AUGUST 16, 1983, MAY 10, 1985 C C ARGUMENT LIST: C IJSWT = SAME AS KAPROX ARGUMENT LIST C BASINV = SAME AS KAPROX ARGUMENT LIST C NDIMP1 = SAME AS KAPROX ARGUMENT LIST C NP1 = N+1 C NP2 = N+2 C KJUMP = IOEXIT(R), WHERE R IS CURRENT DISCRETIZATION C STAGE. (SEE IOEXIT IN KAPROX ARGUMENT LIST.) C P = CURRENT NUMBER OF PROJECTION ANGLES C IPRT = SAME AS KAPROX INPUT ARGUMENT IOEXIT(1) C T = COUNT OF EACH ACTIVE CONSTRAINT TYPE. C LINEAR EQUALITIES ARE COUNTED AS ONE WHEN C THEY HAVE THE SAME FIRST TWO NAMES. C SUBROUTINE KPAIRS(IJSWT,BASINV,NDIMP1,NP1,NP2,KJUMP,P,IPRT,T) C INTEGER IJSWT(3,*),NP1,NP2,J,K,J1,LOOK1,LOOK2,KJUMP,P,L,P1, + NOUT,NDIMP1,BTYPE,ICOUNT,IPRT,T(3) DOUBLE PRECISION BASINV(NDIMP1,1),DZERO C DATA DZERO/+0.0D+0/ C STANDARD FILE NUMBER OF THE OUTPUT STREAM ... FOR VAX 11/780 DATA NOUT/6/ C P1=P-1 ICOUNT=0 DO 80 BTYPE=1,3 T(BTYPE)=0 DO 70 J=1,NP1 LOOK1=IJSWT(1,J) IF(IABS(LOOK1).NE.BTYPE)GO TO 70 IF(LOOK1.GT.0)GO TO 10 IJSWT(1,J)=-IJSWT(1,J) GO TO 70 10 LOOK2=IJSWT(2,J) T(BTYPE)=T(BTYPE)+1 IF(IPRT.LT.1)GO TO 30 ICOUNT=ICOUNT+1 IF(BASINV(J,NP2).GE.DZERO)GO TO 20 WRITE(NOUT,100)ICOUNT,LOOK1,LOOK2,IJSWT(3,J),BASINV(J,NP2) GO TO 30 20 WRITE(NOUT,90)ICOUNT,LOOK1,LOOK2,IJSWT(3,J),BASINV(J,NP2) 30 J1=J+1 IF(J1.GT.NP1)GO TO 70 DO 60 K=J1,NP1 IF(IJSWT(1,K).LT.0)GO TO 60 IF(LOOK1.NE.IJSWT(1,K).OR.LOOK2.NE.IJSWT(2,K))GO TO 60 IF(IPRT.LT.1)GO TO 50 ICOUNT=ICOUNT+1 L=IABS(IJSWT(3,J)-IJSWT(3,K)) IF((L.EQ.1.OR.L.EQ.P1).AND.(BASINV(K,NP2).GE.DZERO)) + GO TO 40 WRITE(NOUT,100)ICOUNT,IJSWT(1,K),IJSWT(2,K),IJSWT(3,K), + BASINV(K,NP2) GO TO 50 40 WRITE(NOUT,90)ICOUNT,IJSWT(1,K),IJSWT(2,K),IJSWT(3,K), + BASINV(K,NP2) 50 IJSWT(1,K)=-IJSWT(1,K) IF(KJUMP.GT.(-1))GO TO 70 60 CONTINUE 70 CONTINUE 80 CONTINUE RETURN 90 FORMAT(7H BASIS.,I6,11H NAME = (,I2,1H,,I6,1H,,I6,10H) VALUE =, + D15.8) 100 FORMAT(7H BASIS.,I6,11H NAME = (,I2,1H,,I6,1H,,I6,10H) VALUE =, + D15.8,3H *) END C SUBROUTINE KEND STORES THE BEST COMPUTED SOLUTION C IN THE LAST COLUMN OF Z, AND SETS RELEVANT DATA C IN IOEXIT(LOGP), ITLOG(LOGP), AND K. NONE OF THE C OTHER SUBROUTINE ARGUMENTS ARE ALTERED. C C WRITTEN BY: ROY L. STREIT C NAVAL UNDERWATER SYSTEMS CENTER C NEW LONDON LABORATORY C NEW LONDON, CT 06320 C MARCH 9, 1983 C REVISED: JUNE 7, 1983, MAY 10, 1985 C C ARGUMENT LIST: C Z, CHEB, LOGP, N, NDIM, IOEXIT, ITLOG C = SAME AS SUBROUTINE KAPROX ARGUMENT LIST C LOGM1 = LOGP-1 C K = ON OUTPUT, THE COLUMN OF Z WHICH CONTAINS THE C BEST COMPUTED SOLUTION. (THE K-TH COLUMN IS C IDENTICAL TO THE LAST COLUMN OF Z, ON RETURN.) C EXCEPTION: K=0 MEANS NO GOOD SOLUTION COMPUTED. C SUBROUTINE KEND(Z,CHEB,LOGP,LOGM1,N,NDIM,IOEXIT,ITLOG,K) C INTEGER I,K,LOGP,LOGM1,N,NDIM,IOEXIT(1),ITLOG(1) REAL Z(NDIM,*),CHEB(1) C K=0 DO 10 I=1,LOGM1 IF(IOEXIT(I).LT.0)GO TO 10 IF(IOEXIT(I).EQ.4)IOEXIT(LOGP)=4 IF(IOEXIT(I).EQ.4)GO TO 20 K=I IOEXIT(LOGP)=MAX0(IOEXIT(LOGP),IOEXIT(I)) 10 CONTINUE 20 IF(K.EQ.0)GO TO 40 ITLOG(LOGP)=ITLOG(K) CHEB(LOGP)=CHEB(K) DO 30 I=1,N Z(I,LOGP)=Z(I,K) 30 CONTINUE 40 RETURN END C SUBROUTINE PABS COMPUTES THE DISCRETIZED ABSOLUTE VALUE C OF A COMPLEX NUMBER AND ITS MINIMAL CLOCKWISE INDEX. C IN ADDITION, PABS INITIALIZES THE COSINES AND SINES C REQUIRED FOR ALL THE PROJECTIONS. C C WRITTEN BY: ROY L. STREIT C NAVAL UNDERWATER SYSTEMS CENTER C NEW LONDON LABORATORY C NEW LONDON, CONNECTICUT 06320 C MAY 28, 1982 C REVISED: APRIL 27, 1983, MAY 10, 1985 C C ARGUMENT LIST: C COMPLEX Z = THE GIVEN COMPLEX NUMBER C INTEGER LOGP = THE BASE 2 LOG OF THE NUMBER OF C PROJECTIONS, P=2**LOGP; LOGP.GE.1 C DOUBLE PNORM = THE MAXIMUM (OVER J=1,...,P) OF: C PRECISION REAL(Z)*COSSIN(1,J)+AIMAG(Z)*COSSIN(2,J) C INTEGER PARG = THE MINIMAL INDEX J FOR WHICH PNORM IS C ATTAINED WHEN TIES ARE RESOLVED IN THE C CLOCKWISE DIRECTION. IN PARTICULAR, C PARG = 1 FOR Z = ZERO. C DOUBLE COSSIN(2,P+1) = THE COSINES AND SINES OF THE ANGLES C PRECISION ASSOCIATED WITH THE PROJECTIONS. THEY C ARE INITIALIZED IN THIS ROUTINE AND MUST C NOT BE ALTERED IN THE CALLING ROUTINE. C COSSIN(1,J) = COS(J-TH PROJECTION ANGLE) C COSSIN(2,J) = SIN(J-TH PROJECTION ANGLE) C NOTE: THE FIRST AND (P+1)-ST PROJECTIONS C ARE IDENTICAL; I.E., BOTH ANGLES = 0. C C COMMENTS: C 1. THE PROJECTION ANGLES ARE TAKEN TO BE EQUISPACED AROUND C THE CIRCLE BEGINNING ON THE POSITIVE X-AXIS. C 2. FOR ALL LOGP .GT. 2, THIS SUBROUTINE REQUIRES THE C CALCULATION OF AT MOST 2*LOGP - 5 MULTIPLICATIONS. C NO MULTIPLICATIONS ARE REQUIRED FOR LOGP = 1,2. C SUBROUTINE PABS(Z,LOGP,PNORM,PARG,COSSIN) COMPLEX Z DOUBLE PRECISION PNORM,COSSIN(2,*) INTEGER LOGP,PARG DOUBLE PRECISION X,Y,PTEST1,PTEST2,DZERO,DONE,PI,COS45 INTEGER P,P2,P4,P8,IDEL,I1,I2,P34,STORPN,P18,P28, + P38,P48,P58,P68,P78,P88 DATA DZERO,DONE/+0.0D+0,+1.0D+0/ DATA STORPN/-1/ C C MACHINE DEPENDENT CONSTANT C DATA PI/+3.141592653589793238D+0/ IF(STORPN.LT.0)PI=4.0D0*DATAN(DONE) C IF(LOGP.GT.1)GO TO 5 C C THE SPECIAL CASE P=2 PNORM=DBLE(ABS(REAL(Z))) PARG=1 IF(DBLE(REAL(Z)).LT.DZERO)PARG=2 IF(DBLE(REAL(Z)).EQ.DZERO.AND.DBLE(AIMAG(Z)).LT.DZERO)PARG=2 IF(STORPN.EQ.1)GO TO 330 STORPN=1 COSSIN(1,1)= DONE COSSIN(2,1)= DZERO COSSIN(1,2)=-DONE COSSIN(2,2)= DZERO COSSIN(1,3)= DONE COSSIN(2,3)= DZERO GO TO 330 C C INITIALIZE IF NECESSARY 5 IF(LOGP.EQ.STORPN)GO TO 40 STORPN=LOGP IF(LOGP.EQ.2)GO TO 20 P=2**LOGP P2=P/2 P4=P2/2 P8=P4/2 P34=P2+P4 P18= 1 +P8 P28=P18+P8 P38=P28+P8 P48=P38+P8 P58=P48+P8 P68=P58+P8 P78=P68+P8 P88=P78+P8 C COMPUTE COSINES AND SINES OF ALL PROJECTIONS X=DBLE(FLOAT(P2)) PTEST1=PI/X DO 10 I1=2,P18 I2=I1-1 X=DBLE(FLOAT(I2))*PTEST1 Y=DSIN(X) X=DCOS(X) C ...FIRST QUADRANT COSSIN(1,I1 )= X COSSIN(2,I1 )= Y IDEL=P28-I2 COSSIN(1,IDEL)= Y COSSIN(2,IDEL)= X C ... SECOND QUADRANT IDEL=I1+P4 COSSIN(1,IDEL)=-Y COSSIN(2,IDEL)= X IDEL=P48-I2 COSSIN(1,IDEL)=-X COSSIN(2,IDEL)= Y C ... THIRD QUADRANT IDEL=I1+P2 COSSIN(1,IDEL)=-X COSSIN(2,IDEL)=-Y IDEL=P68-I2 COSSIN(1,IDEL)=-Y COSSIN(2,IDEL)=-X C ... FOURTH QUADRANT IDEL=I1+P34 COSSIN(1,IDEL)= Y COSSIN(2,IDEL)=-X IDEL=P88-I2 COSSIN(1,IDEL)= X COSSIN(2,IDEL)=-Y 10 CONTINUE COS45=COSSIN(1,P18) GO TO 30 20 P=4 P28=2 P48=3 P68=4 P88=5 30 COSSIN(1,1 )= DONE COSSIN(2,1 )= DZERO COSSIN(1,P28)= DZERO COSSIN(2,P28)= DONE COSSIN(1,P48)=-DONE COSSIN(2,P48)= DZERO COSSIN(1,P68)= DZERO COSSIN(2,P68)=-DONE COSSIN(1,P88)= DONE COSSIN(2,P88)= DZERO C C LOCATE THE COMPLEX NUMBER IN THE PLANE 40 X=DBLE( REAL(Z)) Y=DBLE(AIMAG(Z)) IF(Y.GT.DZERO)GO TO 50 IF(Y.LT.DZERO)GO TO 60 IF(X.GT.DZERO)GO TO 70 IF(X.LT.DZERO)GO TO 80 C ... AT THE ORIGIN PNORM=DZERO PARG=1 GO TO 330 50 IF(X.GT.DZERO)GO TO 90 IF(X.LT.DZERO)GO TO 110 C ... ON THE POSITIVE Y-AXIS PNORM=Y PARG=P28 GO TO 330 60 IF(X.LT.DZERO)GO TO 130 IF(X.GT.DZERO)GO TO 150 C ... ON THE NEGATIVE Y-AXIS PNORM=-Y PARG=P68 GO TO 330 C ... ON THE POSITIVE X-AXIS 70 PNORM=X PARG=1 GO TO 330 C ... ON THE NEGATIVE X-AXIS 80 PNORM=-X PARG=P48 GO TO 330 C C SET SEARCH LIMITS AND FIRST TWO PROJECTIONS IF P.NE.4 C ... IN THE FIRST QUADRANT 90 IF(Y.GT.X)GO TO 100 IF(P.EQ.4)GO TO 170 I1=1 I2=P18 PTEST1=X PTEST2=(X+Y)*COS45 GO TO 250 100 IF(P.EQ.4)GO TO 180 I1=P18 I2=P28 PTEST1=(X+Y)*COS45 PTEST2=Y GO TO 250 C ...IN THE SECOND QUADRANT 110 IF((-X).GT.Y)GO TO 120 IF(P.EQ.4)GO TO 190 I1=P28 I2=P38 PTEST1=Y PTEST2=(Y-X)*COS45 GO TO 250 120 IF(P.EQ.4)GO TO 200 I1=P38 I2=P48 PTEST1=(Y-X)*COS45 PTEST2=-X GO TO 250 C ... IN THE THIRD QUADRANT 130 IF((-Y).GT.(-X))GO TO 140 IF(P.EQ.4)GO TO 210 I1=P48 I2=P58 PTEST1=-X PTEST2=-(X+Y)*COS45 GO TO 250 140 IF(P.EQ.4)GO TO 220 I1=P58 I2=P68 PTEST1=-(X+Y)*COS45 PTEST2=-Y GO TO 250 C ... IN THE FOURTH QUADRANT 150 IF(X.GT.(-Y))GO TO 160 IF(P.EQ.4)GO TO 230 I1=P68 I2=P78 PTEST1=-Y PTEST2=(X-Y)*COS45 GO TO 250 160 IF(P.EQ.4)GO TO 240 I1=P78 I2=P88 PTEST1=(X-Y)*COS45 PTEST2=X GO TO 250 C C FINISH THE SPECIAL CASE P = 4 170 PNORM=X PARG=1 GO TO 330 180 PNORM=Y PARG=2 GO TO 330 190 PNORM=Y PARG=2 GO TO 330 200 PNORM=-X PARG=3 GO TO 330 210 PNORM=-X PARG=3 GO TO 330 220 PNORM=-Y PARG=4 GO TO 330 230 PNORM=-Y PARG=4 GO TO 330 240 PNORM=X PARG=1 GO TO 330 C C USE BISECTION TO GET LARGEST PROJECTION AND MINIMAL INDEX 250 IDEL=P8 260 IF(PTEST1.GT.PTEST2)GO TO 270 IF(PTEST1.LT.PTEST2)GO TO 290 GO TO 310 270 IF(IDEL.GT.1)GO TO 280 PNORM=PTEST1 PARG=I1 GO TO 330 280 IDEL=IDEL/2 I2=I1+IDEL PTEST2=X*COSSIN(1,I2)+Y*COSSIN(2,I2) GO TO 260 290 IF(IDEL.GT.1)GO TO 300 PNORM=PTEST2 IF(I2.EQ.P88)I2=1 PARG=I2 GO TO 330 300 IDEL=IDEL/2 I1=I1+IDEL PTEST1=X*COSSIN(1,I1)+Y*COSSIN(2,I1) GO TO 260 310 IF(IDEL.GT.1)GO TO 320 PNORM=PTEST1 PARG=I1 GO TO 330 320 PARG=(I1+I2)/2 PNORM=X*COSSIN(1,PARG)+Y*COSSIN(2,PARG) C C UNIQUE PROGRAM EXIT 330 RETURN END