1 files changed, 97 insertions, 0 deletions
diff --git a/Utilib/src/ALEIGD.f b/Utilib/src/ALEIGD.f
new file mode 100644
index 0000000..350387f
--- /dev/null
+++ b/Utilib/src/ALEIGD.f
@@ -0,0 +1,97 @@
+*DECK ALEIGD
+      SUBROUTINE ALEIGD (A,B,N,EVAL,EVECT,EPS,ITER)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* find the fondamental eigenvalue and corresponding eigenvector of
+* equation (A-EVAL*B)*EVECT=0 using the inverse power method.
+*
+*Copyright:
+* Copyright (C) 2002 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      first coefficient matrix
+* B      second coefficient matrix
+* N      number of unknowns
+* EVECT  initial estimate
+* EPS2   stopping criterion
+*
+*Parameters: output
+* EVAL   fondamental eigenvalue
+* EVECT  corresponding eigenvector
+* ITER   number of iterations
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER N,ITER
+      DOUBLE PRECISION A(N,N),B(N,N),EVAL,EVECT(N),EPS
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (MMAX=1000)
+      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: GAR
+*----
+*  COMPUTE THE ITERATIVE MATRIX
+*----
+      CALL ALINVD (N,A,N,IER)
+      IF(IER.EQ.1) CALL XABORT('ALEIGD: SINGULAR MATRIX.')
+      ALLOCATE(GAR(N))
+      DO 30 I=1,N
+      DO 10 J=1,N
+      GAR(J)=A(I,J)
+10    CONTINUE
+      DO 25 J=1,N
+      A(I,J)=0.0D0
+      DO 20 K=1,N
+      A(I,J)=A(I,J)+GAR(K)*B(K,J)
+20    CONTINUE
+25    CONTINUE
+30    CONTINUE
+*----
+*  PERFORM POWER ITERATIONS
+*----
+      TEST=0.0D0
+      ITER=0
+      EVAL=0.0D0
+40    ITER=ITER+1
+      IF(ITER.GT.MMAX) CALL XABORT('ALEIGD: UNABLE TO CONVERGE(1).')
+      S1=0.0D0
+      S2=0.0D0
+      DO 60 I=1,N
+      GAR(I)=0.0D0
+      DO 50 J=1,N
+      GAR(I)=GAR(I)+A(I,J)*EVECT(J)
+50    CONTINUE
+      S1=S1+GAR(I)*EVECT(I)
+      S2=S2+GAR(I)**2
+60    CONTINUE
+      IF(S2.EQ.0.0D0) CALL XABORT('ALEIGD: DIVIDE CHECK.')
+      ZZ=ABS(EVAL-S1/S2)
+      EVAL=S1/S2
+      ERR1=0.0D0
+      ERR2=0.0D0
+      DO 70 I=1,N
+      ERR1=MAX(ERR1,ABS(GAR(I)*EVAL))
+      ERR2=MAX(ERR2,ABS(GAR(I)*EVAL-EVECT(I)))
+      EVECT(I)=GAR(I)*EVAL
+70    CONTINUE
+      IF((ZZ.LE.EPS).AND.(ERR2.LE.ERR1*EPS)) THEN
+         DEALLOCATE(GAR)
+         RETURN
+      ENDIF
+      IF(ITER.EQ.1) TEST=ZZ
+      IF((ITER.GE.10).AND.(ZZ.GT.TEST)) CALL XABORT('ALEIGD: UNABLE TO'
+     1 //' CONVERGE(2).')
+      GO TO 40
+      END