c================================================================ c routines from zlinpack that are used in envelope estimation c================================================================ SUBROUTINE ZQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB) INTEGER LDX,N,P,JOB INTEGER JPVT(1) COMPLEX*16 X(LDX,1),QRAUX(1),WORK(1) INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU DOUBLE PRECISION MAXNRM,DZNRM2,TT COMPLEX*16 ZDOTC,NRMXL,T LOGICAL NEGJ,SWAPJ COMPLEX*16 CSIGN,ZDUM,ZDUM1,ZDUM2 DOUBLE PRECISION CABS1 DOUBLE PRECISION DREAL,DIMAG COMPLEX*16 ZDUMR,ZDUMI DREAL(ZDUMR) = ZDUMR DIMAG(ZDUMI) = (0.0D0,-1.0D0)*ZDUMI CSIGN(ZDUM1,ZDUM2) = CDABS(ZDUM1)*(ZDUM2/CDABS(ZDUM2)) CABS1(ZDUM) = DABS(DREAL(ZDUM)) + DABS(DIMAG(ZDUM)) PL = 1 PU = 0 IF (JOB .EQ. 0) GO TO 60 DO 20 J = 1, P SWAPJ = JPVT(J) .GT. 0 NEGJ = JPVT(J) .LT. 0 JPVT(J) = J IF (NEGJ) JPVT(J) = -J IF (.NOT.SWAPJ) GO TO 10 IF (J .NE. PL) CALL ZSWAP(N,X(1,PL),1,X(1,J),1) JPVT(J) = JPVT(PL) JPVT(PL) = J PL = PL + 1 10 CONTINUE 20 CONTINUE PU = P DO 50 JJ = 1, P J = P - JJ + 1 IF (JPVT(J) .GE. 0) GO TO 40 JPVT(J) = -JPVT(J) IF (J .EQ. PU) GO TO 30 CALL ZSWAP(N,X(1,PU),1,X(1,J),1) JP = JPVT(PU) JPVT(PU) = JPVT(J) JPVT(J) = JP 30 CONTINUE PU = PU - 1 40 CONTINUE 50 CONTINUE 60 CONTINUE IF (PU .LT. PL) GO TO 80 DO 70 J = PL, PU QRAUX(J) = DCMPLX(DZNRM2(N,X(1,J),1),0.0D0) WORK(J) = QRAUX(J) 70 CONTINUE 80 CONTINUE LUP = MIN0(N,P) DO 200 L = 1, LUP IF (L .LT. PL .OR. L .GE. PU) GO TO 120 MAXNRM = 0.0D0 MAXJ = L DO 100 J = L, PU IF (DREAL(QRAUX(J)) .LE. MAXNRM) GO TO 90 MAXNRM = DREAL(QRAUX(J)) MAXJ = J 90 CONTINUE 100 CONTINUE IF (MAXJ .EQ. L) GO TO 110 CALL ZSWAP(N,X(1,L),1,X(1,MAXJ),1) QRAUX(MAXJ) = QRAUX(L) WORK(MAXJ) = WORK(L) JP = JPVT(MAXJ) JPVT(MAXJ) = JPVT(L) JPVT(L) = JP 110 CONTINUE 120 CONTINUE QRAUX(L) = (0.0D0,0.0D0) IF (L .EQ. N) GO TO 190 NRMXL = DCMPLX(DZNRM2(N-L+1,X(L,L),1),0.0D0) IF (CABS1(NRMXL) .EQ. 0.0D0) GO TO 180 IF (CABS1(X(L,L)) .NE. 0.0D0) * NRMXL = CSIGN(NRMXL,X(L,L)) CALL ZSCAL(N-L+1,(1.0D0,0.0D0)/NRMXL,X(L,L),1) X(L,L) = (1.0D0,0.0D0) + X(L,L) LP1 = L + 1 IF (P .LT. LP1) GO TO 170 DO 160 J = LP1, P T = -ZDOTC(N-L+1,X(L,L),1,X(L,J),1)/X(L,L) CALL ZAXPY(N-L+1,T,X(L,L),1,X(L,J),1) IF (J .LT. PL .OR. J .GT. PU) GO TO 150 IF (CABS1(QRAUX(J)) .EQ. 0.0D0) GO TO 150 TT = 1.0D0 - (CDABS(X(L,J))/DREAL(QRAUX(J)))**2 TT = DMAX1(TT,0.0D0) T = DCMPLX(TT,0.0D0) TT = 1.0D0 * + 0.05D0*TT * *(DREAL(QRAUX(J))/DREAL(WORK(J)))**2 IF (TT .EQ. 1.0D0) GO TO 130 QRAUX(J) = QRAUX(J)*CDSQRT(T) GO TO 140 130 CONTINUE QRAUX(J) = DCMPLX(DZNRM2(N-L,X(L+1,J),1),0.0D0) WORK(J) = QRAUX(J) 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE QRAUX(L) = X(L,L) X(L,L) = -NRMXL 180 CONTINUE 190 CONTINUE 200 CONTINUE RETURN END SUBROUTINE ZQRSL(X,LDX,N,K,QRAUX,Y,QY,QTY,B,RSD,XB,JOB,INFO) INTEGER LDX,N,K,JOB,INFO COMPLEX*16 X(LDX,1),QRAUX(1),Y(1),QY(1),QTY(1),B(1),RSD(1),XB(1) INTEGER I,J,JJ,JU,KP1 COMPLEX*16 ZDOTC,T,TEMP LOGICAL CB,CQY,CQTY,CR,CXB COMPLEX*16 ZDUM DOUBLE PRECISION CABS1 DOUBLE PRECISION DREAL,DIMAG COMPLEX*16 ZDUMR,ZDUMI DREAL(ZDUMR) = ZDUMR DIMAG(ZDUMI) = (0.0D0,-1.0D0)*ZDUMI CABS1(ZDUM) = DABS(DREAL(ZDUM)) + DABS(DIMAG(ZDUM)) INFO = 0 CQY = JOB/10000 .NE. 0 CQTY = MOD(JOB,10000) .NE. 0 CB = MOD(JOB,1000)/100 .NE. 0 CR = MOD(JOB,100)/10 .NE. 0 CXB = MOD(JOB,10) .NE. 0 JU = MIN0(K,N-1) IF (JU .NE. 0) GO TO 40 IF (CQY) QY(1) = Y(1) IF (CQTY) QTY(1) = Y(1) IF (CXB) XB(1) = Y(1) IF (.NOT.CB) GO TO 30 IF (CABS1(X(1,1)) .NE. 0.0D0) GO TO 10 INFO = 1 GO TO 20 10 CONTINUE B(1) = Y(1)/X(1,1) 20 CONTINUE 30 CONTINUE IF (CR) RSD(1) = (0.0D0,0.0D0) GO TO 250 40 CONTINUE IF (CQY) CALL ZCOPY(N,Y,1,QY,1) IF (CQTY) CALL ZCOPY(N,Y,1,QTY,1) IF (.NOT.CQY) GO TO 70 DO 60 JJ = 1, JU J = JU - JJ + 1 IF (CABS1(QRAUX(J)) .EQ. 0.0D0) GO TO 50 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -ZDOTC(N-J+1,X(J,J),1,QY(J),1)/X(J,J) CALL ZAXPY(N-J+1,T,X(J,J),1,QY(J),1) X(J,J) = TEMP 50 CONTINUE 60 CONTINUE 70 CONTINUE IF (.NOT.CQTY) GO TO 100 DO 90 J = 1, JU IF (CABS1(QRAUX(J)) .EQ. 0.0D0) GO TO 80 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -ZDOTC(N-J+1,X(J,J),1,QTY(J),1)/X(J,J) CALL ZAXPY(N-J+1,T,X(J,J),1,QTY(J),1) X(J,J) = TEMP 80 CONTINUE 90 CONTINUE 100 CONTINUE IF (CB) CALL ZCOPY(K,QTY,1,B,1) KP1 = K + 1 IF (CXB) CALL ZCOPY(K,QTY,1,XB,1) IF (CR .AND. K .LT. N) CALL ZCOPY(N-K,QTY(KP1),1,RSD(KP1),1) IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 120 DO 110 I = KP1, N XB(I) = (0.0D0,0.0D0) 110 CONTINUE 120 CONTINUE IF (.NOT.CR) GO TO 140 DO 130 I = 1, K RSD(I) = (0.0D0,0.0D0) 130 CONTINUE 140 CONTINUE IF (.NOT.CB) GO TO 190 DO 170 JJ = 1, K J = K - JJ + 1 IF (CABS1(X(J,J)) .NE. 0.0D0) GO TO 150 INFO = J GO TO 180 150 CONTINUE B(J) = B(J)/X(J,J) IF (J .EQ. 1) GO TO 160 T = -B(J) CALL ZAXPY(J-1,T,X(1,J),1,B,1) 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE IF (.NOT.CR .AND. .NOT.CXB) GO TO 240 DO 230 JJ = 1, JU J = JU - JJ + 1 IF (CABS1(QRAUX(J)) .EQ. 0.0D0) GO TO 220 TEMP = X(J,J) X(J,J) = QRAUX(J) IF (.NOT.CR) GO TO 200 T = -ZDOTC(N-J+1,X(J,J),1,RSD(J),1)/X(J,J) CALL ZAXPY(N-J+1,T,X(J,J),1,RSD(J),1) 200 CONTINUE IF (.NOT.CXB) GO TO 210 T = -ZDOTC(N-J+1,X(J,J),1,XB(J),1)/X(J,J) CALL ZAXPY(N-J+1,T,X(J,J),1,XB(J),1) 210 CONTINUE X(J,J) = TEMP 220 CONTINUE 230 CONTINUE 240 CONTINUE 250 CONTINUE RETURN END SUBROUTINE ZAXPY(N,ZA,ZX,INCX,ZY,INCY) COMPLEX*16 ZX(1),ZY(1),ZA DOUBLE PRECISION DCABS1 IF(N.LE.0)RETURN IF (DCABS1(ZA) .EQ. 0.0D0) RETURN IF (INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N ZY(IY) = ZY(IY) + ZA*ZX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN 20 DO 30 I = 1,N ZY(I) = ZY(I) + ZA*ZX(I) 30 CONTINUE RETURN END SUBROUTINE ZCOPY(N,ZX,INCX,ZY,INCY) COMPLEX*16 ZX(1),ZY(1) INTEGER I,INCX,INCY,IX,IY,N IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N ZY(IY) = ZX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN 20 DO 30 I = 1,N ZY(I) = ZX(I) 30 CONTINUE RETURN END SUBROUTINE ZSCAL(N,ZA,ZX,INCX) COMPLEX*16 ZA,ZX(1) IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 IX = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 DO 10 I = 1,N ZX(IX) = ZA*ZX(IX) IX = IX + INCX 10 CONTINUE RETURN 20 DO 30 I = 1,N ZX(I) = ZA*ZX(I) 30 CONTINUE RETURN END SUBROUTINE ZSWAP (N,ZX,INCX,ZY,INCY) COMPLEX*16 ZX(1),ZY(1),ZTEMP IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N ZTEMP = ZX(IX) ZX(IX) = ZY(IY) ZY(IY) = ZTEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN 20 DO 30 I = 1,N ZTEMP = ZX(I) ZX(I) = ZY(I) ZY(I) = ZTEMP 30 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DZNRM2( N, ZX, INCX) LOGICAL IMAG, SCALE INTEGER NEXT DOUBLE PRECISION CUTLO, CUTHI, HITEST, SUM, XMAX, ABSX, ZERO, ONE COMPLEX*16 ZX(1) DOUBLE PRECISION DREAL,DIMAG COMPLEX*16 ZDUMR,ZDUMI DREAL(ZDUMR) = ZDUMR DIMAG(ZDUMI) = (0.0D0,-1.0D0)*ZDUMI DATA ZERO, ONE /0.0D0, 1.0D0/ DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / IF(N .GT. 0) GO TO 10 DZNRM2 = ZERO GO TO 300 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX DO 210 I=1,NN,INCX ABSX = DABS(DREAL(ZX(I))) IMAG = .FALSE. GO TO NEXT,(30, 50, 70, 90, 110) 30 IF( ABSX .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT SCALE = .FALSE. 50 IF( ABSX .EQ. ZERO) GO TO 200 IF( ABSX .GT. CUTLO) GO TO 85 ASSIGN 70 TO NEXT GO TO 105 100 ASSIGN 110 TO NEXT SUM = (SUM / ABSX) / ABSX 105 SCALE = .TRUE. XMAX = ABSX GO TO 115 70 IF( ABSX .GT. CUTLO ) GO TO 75 110 IF( ABSX .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / ABSX)**2 XMAX = ABSX GO TO 200 115 SUM = SUM + (ABSX/XMAX)**2 GO TO 200 75 SUM = (SUM * XMAX) * XMAX 85 ASSIGN 90 TO NEXT SCALE = .FALSE. HITEST = CUTHI/FLOAT( N ) 90 IF(ABSX .GE. HITEST) GO TO 100 SUM = SUM + ABSX**2 200 CONTINUE IF(IMAG) GO TO 210 ABSX = DABS(DIMAG(ZX(I))) IMAG = .TRUE. GO TO NEXT,( 50, 70, 90, 110 ) 210 CONTINUE DZNRM2 = DSQRT(SUM) IF(SCALE) DZNRM2 = DZNRM2 * XMAX 300 CONTINUE RETURN END COMPLEX FUNCTION ZDOTC*16(N,ZX,INCX,ZY,INCY) COMPLEX*16 ZX(1),ZY(1),ZTEMP ZTEMP = (0.0D0,0.0D0) ZDOTC = (0.0D0,0.0D0) IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N ZTEMP = ZTEMP + DCONJG(ZX(IX))*ZY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE ZDOTC = ZTEMP RETURN 20 DO 30 I = 1,N ZTEMP = ZTEMP + DCONJG(ZX(I))*ZY(I) 30 CONTINUE ZDOTC = ZTEMP RETURN END DOUBLE PRECISION FUNCTION DCABS1(Z) COMPLEX*16 Z,ZZ DOUBLE PRECISION T(2) EQUIVALENCE (ZZ,T(1)) ZZ = Z DCABS1 = DABS(T(1)) + DABS(T(2)) RETURN END