/
dsygv_example.f90
executable file
·153 lines (123 loc) · 4.79 KB
/
dsygv_example.f90
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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
Program dsygv_example
! DSYGV Example Program Text
! Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
! For licence see
! https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md
! .. Use Statements ..
Use lapack_example_aux, Only: nagf_blas_damax_val, &
nagf_file_print_matrix_real_gen
Use lapack_interfaces, Only: ddisna, dlansy, dsygv, dtrcon
Use lapack_precision, Only: dp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Real (Kind=dp), Parameter :: zero = 0.0_dp
Integer, Parameter :: nb = 64, nin = 5, nout = 6
! .. Local Scalars ..
Real (Kind=dp) :: anorm, bnorm, eps, r, rcond, rcondb, t1, t2, t3
Integer :: i, ifail, info, k, lda, ldb, lwork, n
! .. Local Arrays ..
Real (Kind=dp), Allocatable :: a(:, :), b(:, :), eerbnd(:), rcondz(:), &
w(:), work(:), zerbnd(:)
Real (Kind=dp) :: dummy(1)
Integer, Allocatable :: iwork(:)
! .. Intrinsic Procedures ..
Intrinsic :: abs, epsilon, max, nint
! .. Executable Statements ..
Write (nout, *) 'DSYGV Example Program Results'
Write (nout, *)
! Skip heading in data file
Read (nin, *)
Read (nin, *) n
lda = n
ldb = n
Allocate (a(lda,n), b(ldb,n), eerbnd(n), rcondz(n), w(n), zerbnd(n), &
iwork(n))
! Use routine workspace query to get optimal workspace.
lwork = -1
Call dsygv(1, 'Vectors', 'Upper', n, a, lda, b, ldb, w, dummy, lwork, &
info)
! Make sure that there is enough workspace for block size nb.
lwork = max((nb+2)*n, nint(dummy(1)))
Allocate (work(lwork))
! Read the upper triangular parts of the matrices A and B
Read (nin, *)(a(i,i:n), i=1, n)
Read (nin, *)(b(i,i:n), i=1, n)
! Compute the one-norms of the symmetric matrices A and B
anorm = dlansy('One norm', 'Upper', n, a, lda, work)
bnorm = dlansy('One norm', 'Upper', n, b, ldb, work)
! Solve the generalized symmetric eigenvalue problem
! A*x = lambda*B*x (ITYPE = 1)
Call dsygv(1, 'Vectors', 'Upper', n, a, lda, b, ldb, w, work, lwork, &
info)
If (info==0) Then
! Print solution
Write (nout, *) 'Eigenvalues'
Write (nout, 100) w(1:n)
Write (nout, *)
Flush (nout)
! Normalize the eigenvectors, largest positive
Do i = 1, n
Call nagf_blas_damax_val(n, a(1,i), 1, k, r)
If (a(k,i)<zero) Then
a(1:n, i) = -a(1:n, i)
End If
End Do
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call nagf_file_print_matrix_real_gen('General', ' ', n, n, a, lda, &
'Eigenvectors', ifail)
! Call DTRCON to estimate the reciprocal condition
! number of the Cholesky factor of B. Note that:
! cond(B) = 1/rcond**2
Call dtrcon('One norm', 'Upper', 'Non-unit', n, b, ldb, rcond, work, &
iwork, info)
! Print the reciprocal condition number of B
rcondb = rcond**2
Write (nout, *)
Write (nout, *) 'Estimate of reciprocal condition number for B'
Write (nout, 110) rcondb
Flush (nout)
! Get the machine precision, eps, and if rcondb is not less
! than eps**2, compute error estimates for the eigenvalues and
! eigenvectors
eps = epsilon(1.0E0_dp)
If (rcond>=eps) Then
! Call DDISNA to estimate reciprocal condition
! numbers for the eigenvectors of (A - lambda*B)
Call ddisna('Eigenvectors', n, n, w, rcondz, info)
! Compute the error estimates for the eigenvalues and
! eigenvectors
t1 = eps/rcondb
t2 = anorm/bnorm
t3 = t2/rcond
Do i = 1, n
eerbnd(i) = t1*(t2+abs(w(i)))
zerbnd(i) = t1*(t3+abs(w(i)))/rcondz(i)
End Do
! Print the approximate error bounds for the eigenvalues
! and vectors
Write (nout, *)
Write (nout, *) 'Error estimates for the eigenvalues'
Write (nout, 110) eerbnd(1:n)
Write (nout, *)
Write (nout, *) 'Error estimates for the eigenvectors'
Write (nout, 110) zerbnd(1:n)
Else
Write (nout, *)
Write (nout, *) 'B is very ill-conditioned, error ', &
'estimates have not been computed'
End If
Else If (info>n .And. info<=2*n) Then
i = info - n
Write (nout, 120) 'The leading minor of order ', i, &
' of B is not positive definite'
Else
Write (nout, 130) 'Failure in DSYGV. INFO =', info
End If
100 Format (3X, (6F11.4))
110 Format (4X, 1P, 6E11.1)
120 Format (1X, A, I4, A)
130 Format (1X, A, I4)
End Program