Initial commit from Polytechnique Montreal

1 files changed, 68 insertions, 0 deletions

diff --git a/Utilib/src/AL1EIGD.f b/Utilib/src/AL1EIGD.f
new file mode 100644
index 0000000..4e3a411
--- /dev/null
+++ b/Utilib/src/AL1EIGD.f
@@ -0,0 +1,68 @@
+*DECK AL1EIGD
+      SUBROUTINE AL1EIGD(N,A,EPSOUT,MAXOUT,ITER,EVECT,EVAL,IPRINT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Find the fundamental eigenvalue and corresponding eigenvector of
+* equation (A-EVAL)*EVECT=0 using the power method.
+*
+*Copyright:
+* Copyright (C) 2021 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
+* N       number of unknowns
+* A       coefficient matrix
+* EPSOUT  convergence epsilon for the power method
+* MAXOUT  maximum number of iterations for the power method
+* EVECT   initial estimate
+* IPRINT  print parameter
+*
+*Parameters: output
+* ITER    number of iterations
+* EVECT   corresponding eigenvector
+* EVAL    fondamental eigenvalue
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT REAL*8(A-H,O-Z)
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER N,MAXOUT,ITER,IPRINT
+      REAL(KIND=8) A(N,N),EPSOUT,EVECT(N),EVAL
+*----
+*  LOCAL VARIABLES
+*----
+      REAL(KIND=8), ALLOCATABLE, DIMENSION(:) :: GAR
+*----
+*  POWER METHOD
+*----
+      EVECT(:N)=1.D0
+      ITER=0;
+      ALLOCATE(GAR(N))
+      DO
+        ITER=ITER+1
+        IF (ITER > MAXOUT) CALL XABORT('AL1EIGD: UNABLE TO CONVERGE.')
+        GAR(:)=EVECT(:)
+        EVECT(:)=MATMUL(A(:,:),EVECT(:))
+        EVAL=SQRT(DOT_PRODUCT(EVECT(:),EVECT(:)))
+        EVECT(:)=EVECT(:)/EVAL
+        ERR1=MAXVAL(ABS(EVECT))
+        ERR2=MAXVAL(ABS(GAR(:)-EVECT(:)))
+        IF(IPRINT.GT.1) THEN
+          IF(MOD(ITER,5) == 1) WRITE(6,10) ITER,EVAL,ERR2
+        ENDIF
+        IF(ERR2 <= ERR1*EPSOUT) EXIT
+      ENDDO
+      IF(IPRINT.GT.1) WRITE(6,10) ITER,EVAL,ERR2
+      DEALLOCATE(GAR)
+      RETURN
+   10 FORMAT(15H AL1EIGD: ITER=,I6,6H EVAL=,1P,E13.6,7H ERROR=,E11.4)
+      END