Initial commit from Polytechnique Montreal

1 files changed, 132 insertions, 0 deletions

diff --git a/Utilib/src/ALROOT.f b/Utilib/src/ALROOT.f
new file mode 100644
index 0000000..0d65793
--- /dev/null
+++ b/Utilib/src/ALROOT.f
@@ -0,0 +1,132 @@
+*DECK ALROOT
+      SUBROUTINE ALROOT(A,M,ROOTS,LFAIL)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* find the roots of a polynomial.
+*
+*Copyright:
+* Copyright (C) 1993 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): A. Hebert
+*
+*Parameters: input
+* A       polynomial coefficients.
+* M       polynomial order.
+*
+*Parameters: output
+* ROOTS   complex roots.
+* LFAIL   flag set to .true. in case of failure.
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER M
+      DOUBLE PRECISION A(M+1)
+      COMPLEX*16 ROOTS(M)
+      LOGICAL LFAIL
+*----
+*  LOCAL VARIABLES
+*----
+      COMPLEX*16 AAA,BBB,SQ1,TEST,SQRTM3
+      PARAMETER (SQRTM3=(0.0,1.73205080756888))
+      DOUBLE PRECISION EPS
+      PARAMETER (EPS=1.D-6,MAXM=101)
+      COMPLEX*16 AD(MAXM),X,B,C
+*
+      LFAIL=.FALSE.
+      IF(M+1.GT.MAXM) CALL XABORT('ALROOT: INSUFFICIENT STORAGE.')
+      IF(A(M+1).EQ.0.0D0) CALL XABORT('ALROOT: INVALID COEFFICIENT.')
+      IF(M.EQ.1) THEN
+         ROOTS(1)=-A(1)/A(2)
+      ELSE IF(M.EQ.2) THEN
+         CQ=A(2)/A(3)
+         DQ=A(1)/A(3)
+         AAA=CQ*CQ-4.0D0*DQ
+         AAA=SQRT(AAA)
+         ROOTS(1)=-0.5D0*(CQ+AAA)
+         ROOTS(2)=-0.5D0*(CQ-AAA)
+      ELSE IF(M.EQ.3) THEN
+         IF(A(1).EQ.0.0) THEN
+            CQ=A(3)/A(4)
+            DQ=A(2)/A(4)
+            AAA=CQ*CQ-4.0D0*DQ
+            AAA=SQRT(AAA)
+            ROOTS(1)=0.0
+            ROOTS(2)=-0.5D0*(CQ+AAA)
+            ROOTS(3)=-0.5D0*(CQ-AAA)
+         ELSE
+            BQ=A(3)/A(4)
+            CQ=A(2)/A(4)
+            DQ=A(1)/A(4)
+            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
+         ENDIF
+      ELSE IF(M.EQ.4) THEN
+         CALL ALQUAR(A,ROOTS)
+      ELSE
+         DO 10 J=1,M+1
+         AD(J)=CMPLX(A(J),0.D0,KIND=KIND(AD))
+   10    CONTINUE
+         DO 25 J=M,1,-1
+         X=CMPLX(0.D0,0.D0,KIND=KIND(X))
+         CALL ALGUER(AD,J,X,ITS,LFAIL)
+         IF(LFAIL) RETURN
+         IF(ABS(DIMAG(X)).LE.2.D0*EPS**2*ABS(DBLE(X)))
+     1   X=CMPLX(DBLE(X),0.D0,KIND=KIND(X))
+         ROOTS(J)=X
+         B=AD(J+1)
+         DO 20 JJ=J,1,-1
+         C=AD(JJ)
+         AD(JJ)=B
+         B=X*B+C
+   20    CONTINUE
+   25    CONTINUE
+         DO 30 J=1,M+1
+         AD(J)=CMPLX(A(J),0.D0,KIND=KIND(AD))
+   30    CONTINUE
+         DO 40 J=1,M
+         CALL ALGUER(AD,M,ROOTS(J),ITS,LFAIL)
+         IF(LFAIL) RETURN
+   40    CONTINUE
+      ENDIF
+*
+      DO 70 J=2,M
+      X=ROOTS(J)
+      DO 50 I=J-1,1,-1
+      IF(DBLE(ROOTS(I)).LE.DBLE(X)) GOTO 60
+      ROOTS(I+1)=ROOTS(I)
+   50 CONTINUE
+      I=0
+   60 ROOTS(I+1)=X
+   70 CONTINUE
+      RETURN
+      END