C certain help routines for arithmetic operations with C expressions of the form dexp(e)*c C sven ivansson 1992 SUBROUTINE SCALEEC ( alim1,alim2, E,C ) IMPLICIT NONE REAL*8 alim1,alim2 REAL*8 E COMPLEX*16 C REAL*8 SLASK SLASK = DMAX1 ( DABS(DREAL(C)), DABS(DIMAG(C)) ) IF ( SLASK.LT.alim1 .OR. alim2.LT.SLASK ) THEN IF (SLASK.GT.1.0D0) THEN SLASK = 2/SLASK E = E - DLOG(SLASK) C = SLASK*C ELSE IF (SLASK.GT.0.0D0) THEN SLASK = 0.5D0/SLASK E = E - DLOG(SLASK) C = SLASK*C ELSE E = 0.0D0 C = (0.0D0,0.0D0) END IF END IF RETURN END SUBROUTINE SCALE2EC ( alim1,alim2, E,C1,C2 ) IMPLICIT NONE REAL*8 alim1,alim2 REAL*8 E COMPLEX*16 C1,C2 REAL*8 SLASK SLASK = DMAX1 ( DABS(DREAL(C1)), DABS(DIMAG(C1)), $ DABS(DREAL(C2)), DABS(DIMAG(C2)) ) IF ( SLASK.LT.alim1 .OR. alim2.LT.SLASK ) THEN IF (SLASK.GT.1.0D0) THEN SLASK = 2/SLASK E = E - DLOG(SLASK) C1 = SLASK*C1 C2 = SLASK*C2 ELSE IF (SLASK.GT.0.0D0) THEN SLASK = 0.5D0/SLASK E = E - DLOG(SLASK) C1 = SLASK*C1 C2 = SLASK*C2 ELSE E = 0.0D0 C1 = (0.0D0,0.0D0) C2 = (0.0D0,0.0D0) END IF END IF RETURN END SUBROUTINE SCALE3EC ( alim1,alim2, E,C1,C2,C3 ) IMPLICIT NONE REAL*8 alim1,alim2 REAL*8 E COMPLEX*16 C1,C2,C3 REAL*8 SLASK SLASK = DMAX1 ( DABS(DREAL(C1)), DABS(DIMAG(C1)), $ DABS(DREAL(C2)), DABS(DIMAG(C2)), DABS(DREAL(C3)), DABS(DIMAG(C3)) ) IF ( SLASK.LT.alim1 .OR. alim2.LT.SLASK ) THEN IF (SLASK.GT.1.0D0) THEN SLASK = 2/SLASK E = E - DLOG(SLASK) C1 = SLASK*C1 C2 = SLASK*C2 C3 = SLASK*C3 ELSE IF (SLASK.GT.0.0D0) THEN SLASK = 0.5D0/SLASK E = E - DLOG(SLASK) C1 = SLASK*C1 C2 = SLASK*C2 C3 = SLASK*C3 ELSE E = 0.0D0 C1 = (0.0D0,0.0D0) C2 = (0.0D0,0.0D0) C3 = (0.0D0,0.0D0) END IF END IF RETURN END SUBROUTINE SCALE4EC ( alim1,alim2, E,C1,C2,C3,C4 ) IMPLICIT NONE REAL*8 alim1,alim2 REAL*8 E COMPLEX*16 C1,C2,C3,C4 REAL*8 SLASK SLASK = DMAX1 ( DABS(DREAL(C1)), DABS(DIMAG(C1)), $ DABS(DREAL(C2)), DABS(DIMAG(C2)), DABS(DREAL(C3)), DABS(DIMAG(C3)), $ DABS(DREAL(C4)), DABS(DIMAG(C4)) ) IF ( SLASK.LT.alim1 .OR. alim2.LT.SLASK ) THEN IF (SLASK.GT.1.0D0) THEN SLASK = 2/SLASK E = E - DLOG(SLASK) C1 = SLASK*C1 C2 = SLASK*C2 C3 = SLASK*C3 C4 = SLASK*C4 ELSE IF (SLASK.GT.0.0D0) THEN SLASK = 0.5D0/SLASK E = E - DLOG(SLASK) C1 = SLASK*C1 C2 = SLASK*C2 C3 = SLASK*C3 C4 = SLASK*C4 ELSE E = 0.0D0 C1 = (0.0D0,0.0D0) C2 = (0.0D0,0.0D0) C3 = (0.0D0,0.0D0) C4 = (0.0D0,0.0D0) END IF END IF RETURN END SUBROUTINE SCALE5EC ( alim1,alim2, E,C1,C2,C3,C4,C6 ) IMPLICIT NONE REAL*8 alim1,alim2 REAL*8 E COMPLEX*16 C1,C2,C3,C4,C6 REAL*8 SLASK SLASK = DMAX1 ( DABS(DREAL(C1)), DABS(DIMAG(C1)), $ DABS(DREAL(C2)), DABS(DIMAG(C2)), DABS(DREAL(C3)), DABS(DIMAG(C3)), $ DABS(DREAL(C4)), DABS(DIMAG(C4)), DABS(DREAL(C6)), DABS(DIMAG(C6)) ) IF ( SLASK.LT.alim1 .OR. alim2.LT.SLASK ) THEN IF (SLASK.GT.1.0D0) THEN SLASK = 2/SLASK E = E - DLOG(SLASK) C1 = SLASK*C1 C2 = SLASK*C2 C3 = SLASK*C3 C4 = SLASK*C4 C6 = SLASK*C6 ELSE IF (SLASK.GT.0.0D0) THEN SLASK = 0.5D0/SLASK E = E - DLOG(SLASK) C1 = SLASK*C1 C2 = SLASK*C2 C3 = SLASK*C3 C4 = SLASK*C4 C6 = SLASK*C6 ELSE E = 0.0D0 C1 = (0.0D0,0.0D0) C2 = (0.0D0,0.0D0) C3 = (0.0D0,0.0D0) C4 = (0.0D0,0.0D0) C6 = (0.0D0,0.0D0) END IF END IF RETURN END SUBROUTINE ADDEC ( EA,CA, EB,CB, E,C ) C dexp(e)*c will be dexp(ea)*ca + dexp(eb)*cb C "overwrite" possible IMPLICIT NONE REAL*8 EA,EB COMPLEX*16 CA,CB REAL*8 E COMPLEX*16 C IF (EA.EQ.EB) THEN E = EA C = CA + CB ELSE IF (CA.EQ.(0.0D0,0.0D0)) THEN E = EB C = CB ELSE IF (CB.EQ.(0.0D0,0.0D0)) THEN E = EA C = CA ELSE IF (EB.GE.EA) THEN C = DEXP(EA-EB)*CA + CB E = EB ELSE C = CA + DEXP(EB-EA)*CB E = EA END IF RETURN END SUBROUTINE ADD2EC ( EA,C1A,C2A, EB,C1B,C2B, E,C1,C2 ) C dexp(e)*c1 will be dexp(ea)*c1a + dexp(eb)*c1b C dexp(e)*c2 will be dexp(ea)*c2a + dexp(eb)*c2b C "overwrite" possible IMPLICIT NONE REAL*8 EA,EB COMPLEX*16 C1A,C2A, C1B,C2B REAL*8 E COMPLEX*16 C1,C2 REAL*8 SLASK IF (EA.EQ.EB) THEN E = EA C1 = C1A + C1B C2 = C2A + C2B RETURN END IF IF (C1A.EQ.(0.0D0,0.0D0)) THEN IF (C2A.EQ.(0.0D0,0.0D0)) THEN E = EB C1 = C1B C2 = C2B RETURN END IF END IF IF (C1B.EQ.(0.0D0,0.0D0)) THEN IF (C2B.EQ.(0.0D0,0.0D0)) THEN E = EA C1 = C1A C2 = C2A RETURN END IF END IF IF (EB.GE.EA) THEN SLASK = DEXP(EA-EB) C1 = SLASK*C1A + C1B C2 = SLASK*C2A + C2B E = EB ELSE SLASK = DEXP(EB-EA) C1 = C1A + SLASK*C1B C2 = C2A + SLASK*C2B E = EA END IF RETURN END SUBROUTINE ADD4EC ( EA,C1A,C2A,C3A,C4A, EB,C1B,C2B, $ C3B,C4B, E,C1,C2,C3,C4 ) C dexp(e)*c1 will be dexp(ea)*c1a + dexp(eb)*c1b C dexp(e)*c2 will be dexp(ea)*c2a + dexp(eb)*c2b C dexp(e)*c3 will be dexp(ea)*c3a + dexp(eb)*c3b C dexp(e)*c4 will be dexp(ea)*c4a + dexp(eb)*c4b C "overwrite" possible IMPLICIT NONE REAL*8 EA,EB COMPLEX*16 C1A,C2A,C3A,C4A, C1B,C2B,C3B,C4B REAL*8 E COMPLEX*16 C1,C2,C3,C4 REAL*8 SLASK IF (EA.EQ.EB) THEN E = EA C1 = C1A + C1B C2 = C2A + C2B C3 = C3A + C3B C4 = C4A + C4B RETURN END IF IF (C1A.EQ.(0.0D0,0.0D0)) THEN IF (C2A.EQ.(0.0D0,0.0D0)) THEN IF (C3A.EQ.(0.0D0,0.0D0)) THEN IF (C4A.EQ.(0.0D0,0.0D0)) THEN E = EB C1 = C1B C2 = C2B C3 = C3B C4 = C4B RETURN END IF END IF END IF END IF IF (C1B.EQ.(0.0D0,0.0D0)) THEN IF (C2B.EQ.(0.0D0,0.0D0)) THEN IF (C3B.EQ.(0.0D0,0.0D0)) THEN IF (C4B.EQ.(0.0D0,0.0D0)) THEN E = EA C1 = C1A C2 = C2A C3 = C3A C4 = C4A RETURN END IF END IF END IF END IF IF (EB.GE.EA) THEN SLASK = DEXP(EA-EB) C1 = SLASK*C1A + C1B C2 = SLASK*C2A + C2B C3 = SLASK*C3A + C3B C4 = SLASK*C4A + C4B E = EB ELSE SLASK = DEXP(EB-EA) C1 = C1A + SLASK*C1B C2 = C2A + SLASK*C2B C3 = C3A + SLASK*C3B C4 = C4A + SLASK*C4B E = EA END IF RETURN END SUBROUTINE EQ2EC ( E1,C1, E2,C2 ) C equalize the exponents e1,e2 IMPLICIT NONE REAL*8 E1,E2 COMPLEX*16 C1,C2 IF (E1.EQ.E2) RETURN IF (C1.EQ.(0.0D0,0.0D0)) THEN E1 = E2 ELSE IF (C2.EQ.(0.0D0,0.0D0)) THEN E2 = E1 ELSE IF (E2.GE.E1) THEN C1 = DEXP(E1-E2)*C1 E1 = E2 ELSE C2 = DEXP(E2-E1)*C2 E2 = E1 END IF RETURN END SUBROUTINE EQ22EC ( E1,C11,C12, E2,C21,C22 ) C equalize the exponents e1 (for c11,c12) and e2 (for c21,c22) IMPLICIT NONE REAL*8 E1, E2 COMPLEX*16 C11,C12, C21,C22 REAL*8 SLASK IF (E1.EQ.E2) RETURN IF (C11.EQ.(0.0D0,0.0D0)) THEN IF (C12.EQ.(0.0D0,0.0D0)) THEN E1 = E2 RETURN END IF END IF IF (C21.EQ.(0.0D0,0.0D0)) THEN IF (C22.EQ.(0.0D0,0.0D0)) THEN E2 = E1 RETURN END IF END IF IF (E2.GE.E1) THEN SLASK = DEXP(E1-E2) C11 = SLASK*C11 C12 = SLASK*C12 E1 = E2 ELSE SLASK = DEXP(E2-E1) C21 = SLASK*C21 C22 = SLASK*C22 E2 = E1 END IF RETURN END SUBROUTINE EQ3EC ( E1,C1, E2,C2, E3,C3 ) C equalize the exponents e1,e2,e3 IMPLICIT NONE REAL*8 E1,E2,E3 COMPLEX*16 C1,C2,C3 IF (E1.EQ.E2) THEN IF (E2.EQ.E3) RETURN END IF IF (C1.EQ.(0.0D0,0.0D0)) THEN CALL EQ2EC ( E2,C2, E3,C3 ) E1 = E2 ! .eq.e3 ELSE IF (C2.EQ.(0.0D0,0.0D0)) THEN CALL EQ2EC ( E1,C1, E3,C3 ) E2 = E1 ! .eq.e3 ELSE IF (C3.EQ.(0.0D0,0.0D0)) THEN CALL EQ2EC ( E1,C1, E2,C2 ) E3 = E1 ! .eq.e2 ELSE IF (E1.GE.E2) THEN IF (E1.GE.E3) THEN GOTO 1 ELSE GOTO 3 END IF ELSE IF (E2.GE.E3) THEN GOTO 2 ELSE GOTO 3 END IF END IF 1 C2 = DEXP(E2-E1)*C2 E2 = E1 C3 = DEXP(E3-E1)*C3 E3 = E1 RETURN 2 C1 = DEXP(E1-E2)*C1 E1 = E2 C3 = DEXP(E3-E2)*C3 E3 = E2 RETURN 3 C1 = DEXP(E1-E3)*C1 E1 = E3 C2 = DEXP(E2-E3)*C2 E2 = E3 RETURN END SUBROUTINE EQ4EC ( E1,C1, E2,C2, E3,C3, E4,C4 ) C equalize the exponents e1,e2,e3,e4 IMPLICIT NONE REAL*8 E1,E2,E3,E4 COMPLEX*16 C1,C2,C3,C4 IF (E1.EQ.E2) THEN IF (E2.EQ.E3) THEN IF (E3.EQ.E4) RETURN END IF END IF IF (C1.EQ.(0.0D0,0.0D0)) THEN CALL EQ3EC ( E2,C2, E3,C3, E4,C4 ) E1 = E2 ! .eq.e3.eq.e4 ELSE IF (C2.EQ.(0.0D0,0.0D0)) THEN CALL EQ3EC ( E1,C1, E3,C3, E4,C4 ) E2 = E1 ! .eq.e3.eq.e4 ELSE IF (C3.EQ.(0.0D0,0.0D0)) THEN CALL EQ3EC ( E1,C1, E2,C2, E4,C4 ) E3 = E1 ! .eq.e2.eq.e4 ELSE IF (C4.EQ.(0.0D0,0.0D0)) THEN CALL EQ3EC ( E1,C1, E2,C2, E3,C3 ) E4 = E1 ! .eq.e2.eq.e3 ELSE IF (E1.GE.E2) THEN IF (E1.GE.E3) THEN IF (E1.GE.E4) THEN GOTO 1 ELSE GOTO 4 END IF ELSE IF (E3.GE.E4) THEN GOTO 3 ELSE GOTO 4 END IF END IF ELSE IF (E2.GE.E3) THEN IF (E2.GE.E4) THEN GOTO 2 ELSE GOTO 4 END IF ELSE IF (E3.GE.E4) THEN GOTO 3 ELSE GOTO 4 END IF END IF END IF 1 C2 = DEXP(E2-E1)*C2 E2 = E1 C3 = DEXP(E3-E1)*C3 E3 = E1 C4 = DEXP(E4-E1)*C4 E4 = E1 RETURN 2 C1 = DEXP(E1-E2)*C1 E1 = E2 C3 = DEXP(E3-E2)*C3 E3 = E2 C4 = DEXP(E4-E2)*C4 E4 = E2 RETURN 3 C1 = DEXP(E1-E3)*C1 E1 = E3 C2 = DEXP(E2-E3)*C2 E2 = E3 C4 = DEXP(E4-E3)*C4 E4 = E3 RETURN 4 C1 = DEXP(E1-E4)*C1 E1 = E4 C2 = DEXP(E2-E4)*C2 E2 = E4 C3 = DEXP(E3-E4)*C3 E3 = E4 RETURN END SUBROUTINE EQ24EC ( E1,C11,C12,C13,C14, $ E2,C21,C22,C23,C24 ) C equalize the exponents e1 (for c11,c12,c13,c14) and c e2 (for c21,c22,c23,c24) IMPLICIT NONE REAL*8 E1, E2 COMPLEX*16 C11,C12,C13,C14, C21,C22,C23,C24 REAL*8 SLASK IF (E1.EQ.E2) RETURN IF (C11.EQ.(0.0D0,0.0D0)) THEN IF (C12.EQ.(0.0D0,0.0D0)) THEN IF (C13.EQ.(0.0D0,0.0D0)) THEN IF (C14.EQ.(0.0D0,0.0D0)) THEN E1 = E2 RETURN END IF END IF END IF END IF IF (C21.EQ.(0.0D0,0.0D0)) THEN IF (C22.EQ.(0.0D0,0.0D0)) THEN IF (C23.EQ.(0.0D0,0.0D0)) THEN IF (C24.EQ.(0.0D0,0.0D0)) THEN E2 = E1 RETURN END IF END IF END IF END IF IF (E2.GE.E1) THEN SLASK = DEXP(E1-E2) C11 = SLASK*C11 C12 = SLASK*C12 C13 = SLASK*C13 C14 = SLASK*C14 E1 = E2 ELSE SLASK = DEXP(E2-E1) C21 = SLASK*C21 C22 = SLASK*C22 C23 = SLASK*C23 C24 = SLASK*C24 E2 = E1 END IF RETURN END