Initial commit from Polytechnique Montreal

1 files changed, 237 insertions, 0 deletions

diff --git a/Utilib/src/ALQUAR.f b/Utilib/src/ALQUAR.f
new file mode 100644
index 0000000..0dcbf41
--- /dev/null
+++ b/Utilib/src/ALQUAR.f
@@ -0,0 +1,237 @@
+*DECK ALQUAR
+      SUBROUTINE ALQUAR(A,ROOTS)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* compute the roots of the real quartic polynomial defined as
+* A(1)+A(2)*Z + ... + A(5)*Z**4.
+* NOTE: It is assumed that A(5) is non-zero. No test is made here.
+*
+*Author(s): A. H. Morris, W. L. Davis, A. Miller, and R. L. Carmichael
+*
+*Parameters: input
+* A      polynomial coefficients
+*
+*Parameters: output
+* ROOTS  complex roots
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      DOUBLE PRECISION, INTENT(IN) :: A(5)
+      COMPLEX*16, INTENT(OUT) :: ROOTS(4)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER :: I, J
+      INTEGER,DIMENSION(1) :: K
+      DOUBLE PRECISION :: B, B2, C, D, E, H, P, Q, R, T, WORK, U, V,
+     > V1, V2, X, XX(3), Y, BQ, CQ, DQ, AA, BB
+      COMPLEX*16 :: W, AAA, BBB, SQ1, TEST, SQRTM3
+      DOUBLE PRECISION,DIMENSION(4) :: TEMP
+      PARAMETER (SQRTM3=(0.0,1.73205080756888))
+*
+      IF(A(1)==0.0) THEN
+        IF(A(2).EQ.0.0) THEN
+           CQ=A(4)/A(5)
+           DQ=A(3)/A(5)
+           AAA=CQ*CQ-4.0D0*DQ
+           AAA=SQRT(AAA)
+           ROOTS(1)=0.0
+           ROOTS(2)=0.0
+           ROOTS(3)=-0.5D0*(CQ+AAA)
+           ROOTS(4)=-0.5D0*(CQ-AAA)
+        ELSE
+           BQ=A(4)/A(5)
+           CQ=A(3)/A(5)
+           DQ=A(2)/A(5)
+           AA=(3.0D0*CQ-BQ**2)/3.0D0
+           BB=(2.0D0*BQ**3-9.0D0*BQ*CQ+27.0D0*DQ)/27.0D0
+           SQ1=BB**2/4.0D0+AA**3/27.0D0
+           TEST=BB/2.0D0-SQRT(SQ1)
+           IF(DBLE(TEST).EQ.0.0) THEN
+              AAA=0.0D0
+           ELSE IF(DBLE(TEST).GT.0.0) THEN
+             AAA=-(TEST)**(1.0D0/3.0D0)
+           ELSE
+             AAA=(-TEST)**(1.0D0/3.0D0)
+           ENDIF
+           TEST=BB/2.0D0+SQRT(SQ1)
+           IF(DBLE(TEST).EQ.0.0) THEN
+              BBB=0.0D0
+           ELSE IF(DBLE(TEST).GT.0.0) THEN
+              BBB=-(TEST)**(1.0D0/3.0D0)
+           ELSE
+              BBB=(-TEST)**(1.0D0/3.0D0)
+           ENDIF
+           ROOTS(1)=0.0
+           ROOTS(2)=AAA+BBB-BQ/3.0D0
+           ROOTS(3)=-(AAA+BBB)/2.0D0+(AAA-BBB)*SQRTM3/2.0D0-BQ/3.0D0
+           ROOTS(4)=-(AAA+BBB)/2.0D0-(AAA-BBB)*SQRTM3/2.0D0-BQ/3.0D0
+        ENDIF
+        RETURN
+      ENDIF
+*----
+*  Solve a quartic equation
+*----
+      B = A(4)/(4.0D0*A(5))
+      C = A(3)/A(5)
+      D = A(2)/A(5)
+      E = A(1)/A(5)
+      B2 = B*B
+
+      P = 0.5D0*(C - 6.0D0*B2)
+      Q = D - 2.0D0*B*(C - 4.0D0*B2)
+      R = B2*(C - 3.0D0*B2) - B*D + E
+*----
+*  Solve the resolvent cubic equation. the cubic has at least one
+*  nonnegative real root.  if W1, W2, W3 are the roots of the cubic
+*  then the roots of the original equation are
+*     ROOTS = -B + CSQRT(W1) + CSQRT(W2) + CSQRT(W3)
+*  where the signs of the square roots are chosen so
+*  that CSQRT(W1) * CSQRT(W2) * CSQRT(W3) = -Q/8.
+*----
+      TEMP(1) = -Q*Q/64.0D0
+      TEMP(2) = 0.25D0*(P*P - R)
+      TEMP(3) =  P
+      TEMP(4) = 1.0D0
+      BQ=TEMP(3)
+      CQ=TEMP(2)
+      DQ=TEMP(1)
+      AA=(3.0D0*CQ-BQ**2)/3.0D0
+      BB=(2.0D0*BQ**3-9.0D0*BQ*CQ+27.0D0*DQ)/27.0D0
+      SQ1=BB**2/4.0D0+AA**3/27.0D0
+      TEST=BB/2.0D0-SQRT(SQ1)
+      IF(DBLE(TEST).EQ.0.0) THEN
+         AAA=0.0D0
+      ELSE IF(DBLE(TEST).GT.0.0) THEN
+         AAA=-(TEST)**(1.0D0/3.0D0)
+      ELSE
+         AAA=(-TEST)**(1.0D0/3.0D0)
+      ENDIF
+      TEST=BB/2.0D0+SQRT(SQ1)
+      IF(DBLE(TEST).EQ.0.0) THEN
+         BBB=0.0D0
+      ELSE IF(DBLE(TEST).GT.0.0) THEN
+         BBB=-(TEST)**(1.0D0/3.0D0)
+      ELSE
+         BBB=(-TEST)**(1.0D0/3.0D0)
+      ENDIF
+      ROOTS(1)=AAA+BBB-BQ/3.0D0
+      ROOTS(2)=-(AAA+BBB)/2.0D0+(AAA-BBB)*SQRTM3/2.0D0-BQ/3.0D0
+      ROOTS(3)=-(AAA+BBB)/2.0D0-(AAA-BBB)*SQRTM3/2.0D0-BQ/3.0D0
+      IF(AIMAG(ROOTS(2)).NE.0.0D0) GO TO 60
+*----
+*  The resolvent cubic has only real roots.
+*  Reorder the roots in increasing order.
+*----
+      XX(1) = DBLE(ROOTS(1))
+      XX(2) = DBLE(ROOTS(2))
+      XX(3) = DBLE(ROOTS(3))
+      DO 25 J=2,3
+      X=XX(J)
+      DO 10 I=J-1,1,-1
+      IF(XX(I).LE.X) GOTO 20
+      XX(I+1)=XX(I)
+   10 CONTINUE
+      I=0
+   20 XX(I+1)=X
+   25 CONTINUE
+
+      U = 0.0D0
+      IF(XX(3).GT.0.0D0) U = SQRT(XX(3))
+      IF(XX(2).LE.0.0D0) GO TO 41
+      IF(XX(1).GE.0.0D0) GO TO 30
+      IF(ABS(XX(1)).GT.XX(2)) GO TO 40
+      XX(1) = 0.0D0
+
+   30 XX(1) = SQRT(XX(1))
+      XX(2) = SQRT(XX(2))
+      IF(Q.GT.0.0D0) XX(1) = -XX(1)
+      TEMP(1) = (( XX(1) + XX(2)) + U) - B
+      TEMP(2) = ((-XX(1) - XX(2)) + U) - B
+      TEMP(3) = (( XX(1) - XX(2)) - U) - B
+      TEMP(4) = ((-XX(1) + XX(2)) - U) - B
+
+      DO J=1,3
+        K=MINLOC(TEMP(J:))
+        IF(J.NE.K(1)) THEN
+          WORK = TEMP(J)
+          TEMP(J) = TEMP(K(1))
+          TEMP(K(1)) = WORK
+        ENDIF
+      ENDDO
+
+      IF(ABS(TEMP(1)).GE.0.1D0*ABS(TEMP(4))) GO TO 31
+      T = TEMP(2)*TEMP(3)*TEMP(4)
+      IF(T.NE.0.0D0) TEMP(1) = E/T
+   31 ROOTS(1) = CMPLX(TEMP(1), 0.0D0, KIND=KIND(ROOTS))
+      ROOTS(2) = CMPLX(TEMP(2), 0.0D0, KIND=KIND(ROOTS))
+      ROOTS(3) = CMPLX(TEMP(3), 0.0D0, KIND=KIND(ROOTS))
+      ROOTS(4) = CMPLX(TEMP(4), 0.0D0, KIND=KIND(ROOTS))
+      RETURN
+
+   40 V1 = SQRT(ABS(XX(1)))
+      V2 = 0.0D0
+      GO TO 50
+   41 V1 = SQRT(ABS(XX(1)))
+      V2 = SQRT(ABS(XX(2)))
+      IF(Q < 0.0D0) U = -U
+
+   50 X = -U - B
+      Y = V1 - V2
+      ROOTS(1) = CMPLX(X, Y, KIND=KIND(ROOTS))
+      ROOTS(2) = CMPLX(X,-Y, KIND=KIND(ROOTS))
+      X =  U - B
+      Y = V1 + V2
+      ROOTS(3) = CMPLX(X, Y, KIND=KIND(ROOTS))
+      ROOTS(4) = CMPLX(X,-Y, KIND=KIND(ROOTS))
+      RETURN
+*----
+*  The resolvent cubic has complex roots.
+*----
+   60 T = DBLE(ROOTS(1))
+      X = 0.0D0
+      IF(T < 0.0D0) THEN
+        GO TO 61
+      ELSE IF(T.EQ.0.0D0) THEN
+        GO TO 70
+      ELSE
+        GO TO 62
+      ENDIF
+   61 H = ABS(DBLE(ROOTS(2))) + ABS(AIMAG(ROOTS(2)))
+      IF(ABS(T).LE.H) GO TO 70
+      GO TO 80
+   62 X = SQRT(T)
+      IF(Q.GT.0.0D0) X = -X
+
+   70 W = SQRT(ROOTS(2))
+      U = 2.0D0*DBLE(W)
+      V = 2.0D0*ABS(AIMAG(W))
+      T =  X - B
+      XX(1) = T + U
+      XX(2) = T - U
+      IF(ABS(XX(1)).LE.ABS(XX(2))) GO TO 71
+      T = XX(1)
+      XX(1) = XX(2)
+      XX(2) = T
+   71 U = -X - B
+      H = U*U + V*V
+      IF(XX(1)*XX(1) < 0.01D0*MIN(XX(2)*XX(2),H)) XX(1) = E/(XX(2)*H)
+      ROOTS(1) = CMPLX(XX(1), 0.0D0, KIND=KIND(ROOTS))
+      ROOTS(2) = CMPLX(XX(2), 0.0D0, KIND=KIND(ROOTS))
+      ROOTS(3) = CMPLX(U, V, KIND=KIND(ROOTS))
+      ROOTS(4) = CMPLX(U,-V, KIND=KIND(ROOTS))
+      RETURN
+
+   80 V = SQRT(ABS(T))
+      ROOTS(1) = CMPLX(-B, V, KIND=KIND(ROOTS))
+      ROOTS(2) = CMPLX(-B,-V, KIND=KIND(ROOTS))
+      ROOTS(3) = ROOTS(1)
+      ROOTS(4) = ROOTS(2)
+      RETURN
+      END