SUBROUTINE ZHECO (H,LDH,N,IPVT,RCOND,Z) INTEGER LDH,N,IPVT(N) DOUBLE PRECISION RCOND DOUBLE COMPLEX H(LDH,N),Z(N) C C *** PURPOSE: C C ZHECO FACTORS A COMPLEX UPPER HESSENBERG MATRIX BY GAUSSIAN C ELIMINATION AND ESTIMATES THE CONDITION OF THE MATRIX. THE C LINPACK ANALOGUE OF ZHECO IS ZGECO. C C IF RCOND IS NOT NEEDED, ZHEFA MAY BE SUBSTANTIALLY FASTER. TO C SOLVE THE LINEAR SYSTEM H*X = B (I.E., TO COMPUTE C (INVERSE(H)*B)), FOLLOW ZHECO BY ZHESL. C C ON ENTRY: C C H DOUBLE COMPLEX(LDH,N) C THE MATRIX TO BE FACTORED. C C LDH INTEGER C THE LEADING OR ROW DIMENSION OF THE ARRAY H AS DECLARED C IN THE MAIN CALLING PROGRAM. C C N INTEGER C THE ORDER OF THE MATRIX H. C C ON RETURN: C C H AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS WHICH WERE C USED TO OBTAIN IT. THE FACTORIZATION CAN BE WRITTEN C H = L*U WHERE L IS A PRODUCT OF PERMUTATION AND UNIT C LOWER BIDIAGONAL MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL OF THE CONDITION OF H. C FOR THE SYSTEM H*X = B , RELATIVE PERTURBATIONS C IN H AND B OF SIZE EPSILON MAY CAUSE RELATIVE PERTURBA- C TIONS IN X OF SIZE EPSILON/RCOND. IF RCOND IS SO SMALL C THAT THE LOGICAL EXPRESSION (1.0+RCOND) .EQ. 1.0 C IS TRUE, THEN H MAY BE SINGULAR TO WORKING PRECISION. C IN PARTICULAR, RCOND IS ZERO IF EXACT SINGULARITY IS C DETECTED OR THE ESTIMATE UNDERFLOWS. C C Z DOUBLE COMPLEX(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF H IS CLOSE TO A SINGULAR MATRIX, THEN Z IS AN C APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(H*Z) = RCOND*NORM(H)*NORM(Z) C C THIS VERSION DATED MAY 1980. C ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA. C C SUBROUTINES AND FUNCTIONS CALLED: C C DCL1,ZHEFA,ZL1NRM C DOUBLE PRECISION DCL1 C C INTERNAL VARIABLES: C INTEGER I,INFO,J,K,KB,KM1,KP1,L DOUBLE PRECISION HNORM,S,SM,YNORM,ZNORM DOUBLE COMPLEX EK,T,WK,WKM C C FORTRAN FUNCTIONS CALLED: C DOUBLE COMPLEX DCMPLX,DCONJG C C INTERNAL FUNCTION: C DOUBLE COMPLEX Z1,Z2,CSIGN CSIGN(Z1,Z2) = DCL1(Z1)*(Z2/DCL1(Z2)) C IF (N.GT.1) GO TO 10 C C SPECIAL CASE OF SCALAR (N = 1) SYSTEM C HNORM = DCL1(H(1,1)) YNORM = HNORM GO TO 250 10 CONTINUE C C COMPUTE COMPLEX 1-NORM OF H C HNORM = 0.0D0 DO 30 J = 1,N S = 0.0D0 DO 20 I = 1,N S = S+DCL1(H(I,J)) 20 CONTINUE IF (S .GT. HNORM) HNORM = S 30 CONTINUE C C FACTOR C CALL ZHEFA (H,LDH,N,IPVT,INFO) C C RCOND = 1/(NORM(H)*(ESTIMATE OF NORM(INVERSE(H))). C ESTIMATE = NORM(Z)/NORM(Y) WHERE H*Z =Y C AND TRANSPOSE(H)*Y = E. THE COMPONENTS OF E ARE C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS C OF W WHERE TRANSPOSE(U)*W = E. THE VECTORS ARE C FREQUENTLY RESCALED TO AVOID OVERFLOW. C C SOLVE TRANSPOSE(U)*W = E C EK = (1.0D0,0.0D0) DO 40 J = 1,N Z(J) = (0.0D0,0.0D0) 40 CONTINUE DO 130 K = 1,N IF (DCL1(Z(K)) .NE. 0.0D0) EK = CSIGN(EK,-Z(K)) IF (DCL1(EK-Z(K)) .LE. DCL1(H(K,K))) GO TO 60 S = DCL1(H(K,K))/DCL1(EK-Z(K)) DO 50 I = 1,N Z(I) = DCMPLX(S,0.0D0)*Z(I) 50 CONTINUE EK = DCMPLX(S,0.0D0)*EK 60 CONTINUE WK = EK-Z(K) WKM = -EK-Z(K) S = DCL1(WK) SM = DCL1(WKM) IF (DCL1(H(K,K)) .EQ. 0.0D0) GO TO 70 WK = WK/DCONJG(H(K,K)) WKM = WKM/DCONJG(H(K,K)) GO TO 80 70 CONTINUE WK = (1.0D0,0.0D0) WKM = (1.0D0,0.0D0) 80 CONTINUE KP1 = K+1 IF (KP1 .GT. N) GO TO 120 DO 90 J = KP1,N SM = SM+DCL1(Z(J)+WKM*DCONJG(H(K,J))) Z(J) = Z(J)+WK*DCONJG(H(K,J)) S = S+DCL1(Z(J)) 90 CONTINUE IF (S .GE. SM) GO TO 110 T = WKM-WK WK = WKM DO 100 J = KP1,N Z(J) = Z(J)+T*DCONJG(H(K,J)) 100 CONTINUE 110 CONTINUE 120 CONTINUE Z(K) = WK 130 CONTINUE CALL ZL1NRM (N,Z,ZNORM) C C SOLVE TRANSPOSE(L)*Y = W C DO 160 KB = 1,N K = N+1-KB IF (K .LT. N) Z(K) = Z(K)-DCONJG(H(K+1,K))*Z(K+1) IF (DCL1(Z(K)) .LE. 1.0D0) GO TO 150 S = 1.0D0/DCL1(Z(K)) DO 140 I = 1,N Z(I) = DCMPLX(S,0.0D0)*Z(I) 140 CONTINUE 150 CONTINUE L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T 160 CONTINUE CALL ZL1NRM (N,Z,ZNORM) C YNORM = 1.0D0 C C SOLVE L*V = Y C DO 190 K = 1,N L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T IF (K .LT. N) Z(K+1) = Z(K+1)-T*H(K+1,K) IF (DCL1(Z(K)) .LE. 1.0D0) GO TO 180 S = 1.0D0/DCL1(Z(K)) DO 170 I = 1,N Z(I) = DCMPLX(S,0.0D0)*Z(I) 170 CONTINUE YNORM = S*YNORM 180 CONTINUE 190 CONTINUE CALL ZL1NRM (N,Z,ZNORM) YNORM = YNORM/ZNORM C C SOLVE U*Z = V C DO 240 KB = 1,N K = N+1-KB IF (DCL1(Z(K)) .LE. DCL1(H(K,K))) GO TO 210 S = DCL1(H(K,K))/DCL1(Z(K)) DO 200 I = 1,N Z(I) = DCMPLX(S,0.0D0)*Z(I) 200 CONTINUE YNORM = S*YNORM 210 CONTINUE IF (DCL1(H(K,K)) .NE. 0.0D0) Z(K) = Z(K)/H(K,K) IF (DCL1(H(K,K)) .EQ. 0.0D0) Z(K) = (1.0D0,0.0D0) T = -Z(K) IF (K .LE. 1) GO TO 230 KM1 = K-1 DO 220 I = 1,KM1 Z(I) = Z(I)+T*H(I,K) 220 CONTINUE 230 CONTINUE 240 CONTINUE CALL ZL1NRM (N,Z,ZNORM) YNORM = YNORM/ZNORM 250 CONTINUE IF (HNORM .NE. 0.0D0) RCOND = YNORM/HNORM IF (HNORM .EQ. 0.0D0) RCOND = 0.0D0 RETURN END