1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
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
|
