summaryrefslogtreecommitdiff
path: root/Utilib/src/AL1EIG.f
diff options
context:
space:
mode:
authorstainer_t <thomas.stainer@oecd-nea.org>2025-09-08 13:48:49 +0200
committerstainer_t <thomas.stainer@oecd-nea.org>2025-09-08 13:48:49 +0200
commit7dfcc480ba1e19bd3232349fc733caef94034292 (patch)
tree03ee104eb8846d5cc1a981d267687a729185d3f3 /Utilib/src/AL1EIG.f
Initial commit from Polytechnique Montreal
Diffstat (limited to 'Utilib/src/AL1EIG.f')
-rw-r--r--Utilib/src/AL1EIG.f67
1 files changed, 67 insertions, 0 deletions
diff --git a/Utilib/src/AL1EIG.f b/Utilib/src/AL1EIG.f
new file mode 100644
index 0000000..be699d8
--- /dev/null
+++ b/Utilib/src/AL1EIG.f
@@ -0,0 +1,67 @@
+*DECK AL1EIG
+ SUBROUTINE AL1EIG(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
+*
+*-----------------------------------------------------------------------
+*
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER N,MAXOUT,ITER,IPRINT
+ REAL A(N,N),EPSOUT,EVECT(N),EVAL
+*----
+* LOCAL VARIABLES
+*----
+ REAL, ALLOCATABLE, DIMENSION(:) :: GAR
+*----
+* POWER METHOD
+*----
+ EVECT(:N)=1.0
+ ITER=0;
+ ALLOCATE(GAR(N))
+ DO
+ ITER=ITER+1
+ IF (ITER > MAXOUT) CALL XABORT('AL1EIG: 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(14H AL1EIG: ITER=,I6,6H EVAL=,1P,E12.5,7H ERROR=,E11.4)
+ END