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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
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
|
