From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Utilib/src/ALROOT.f | 132 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 132 insertions(+) create mode 100644 Utilib/src/ALROOT.f (limited to 'Utilib/src/ALROOT.f') 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 -- cgit v1.2.3