FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_estimate_condition.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5
7
8 public :: hecmw_estimate_cond_num_cg
9 public :: hecmw_estimate_cond_num_gmres
10
11contains
12
13 subroutine hecmw_estimate_condition_cg(ITER, D, E)
14 use hecmw_util
15 implicit none
16 integer(kind=kint), intent(in) :: ITER
17 real(kind=kreal), intent(in) :: d(:), e(:)
18
19#ifdef HECMW_WITH_LAPACK
20 character(len=1) :: JOBZ, range
21 ! character(len=1) :: COMPZ
22 real(kind=kreal) :: vl, vu, abstol, z(1,1)
23 integer(kind=kint) :: N, IL, IU, M, LDZ=1, isuppz(1)
24 integer(kind=kint) :: LWORK, LIWORK, INFO
25 real(kind=kreal), allocatable :: w(:), work(:)
26 integer(kind=kint), allocatable :: IWORK(:)
27 real(kind=kreal), allocatable :: d1(:), e1(:)
28 integer(kind=kint) :: i
29
30 if (iter <= 1) return
31
32 ! copy D, E
33 allocate(d1(iter),e1(iter))
34 do i=1,iter-1
35 d1(i) = d(i)
36 e1(i) = e(i)
37 enddo
38 d1(iter) = d(iter)
39
40
41 !!
42 !! dstegr version (faster than dsteqr)
43 !!
44
45 ! prepare arguments for calling dstegr
46 jobz='N'
47 range='A'
48 n=iter
49 allocate(w(iter))
50 ! estimate optimal LWORK and LIWORK
51 lwork=-1
52 liwork=-1
53 allocate(work(1),iwork(1))
54 call dstegr(jobz,range,n,d1,e1,vl,vu,il,iu,abstol, &
55 m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info)
56 if (info /= 0) then
57 write(*,*) 'ERROR: dstegr returned with INFO=',info
58 return
59 endif
60 ! calculate eigenvalues
61 lwork=work(1)
62 liwork=iwork(1)
63 deallocate(work,iwork)
64 allocate(work(lwork),iwork(liwork))
65 call dstegr(jobz,range,n,d1,e1,vl,vu,il,iu,abstol, &
66 m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info)
67 if (info /= 0) then
68 write(*,*) 'ERROR: dstegr returned with INFO=',info
69 return
70 endif
71 write(*,'("emin=",1pe13.6,", emax=",1pe13.6,", emax/emin=",1pe13.6)') &
72 w(1),w(n),w(n)/w(1)
73 deallocate(work,iwork)
74 deallocate(w)
75
76
77 ! !!
78 ! !! dsteqr version
79 ! !!
80
81 ! ! prepare arguments for calling dsteqr
82 ! COMPZ='N'
83 ! N=ITER
84 ! allocate(WORK(1))
85 ! ! calculate eigenvalues
86 ! call dsteqr(COMPZ,N,D1,E1,Z,LDZ,WORK,INFO)
87 ! if (INFO /= 0) then
88 ! write(*,*) 'ERROR: dsteqr returned with INFO=',INFO
89 ! return
90 ! endif
91 ! write(*,'("emin=",1pe13.6,", emax=",1pe13.6,", emax/emin=",1pe13.6)') &
92 ! D1(1),D1(N),D1(N)/D1(1)
93 ! deallocate(WORK)
94
95
96 deallocate(d1,e1)
97#endif
98 end subroutine hecmw_estimate_condition_cg
99
100
102 use hecmw_util
103 implicit none
104 integer(kind=kint), intent(in) :: I
105 real(kind=kreal), intent(in) :: h(:,:)
106
107#ifdef HECMW_WITH_LAPACK
108 ! character(len=1) :: JOBU, JOBVT
109 character(len=1) :: JOBZ
110 integer(kind=kint) :: N, LDH, LDZ=1, lwork, info
111 real(kind=kreal), allocatable :: wr(:), work(:), h1(:,:)
112 integer(kind=kint), allocatable :: IWORK(:)
113 real(kind=kreal) :: z(1,1)
114 integer(kind=kint) :: j, k
115
116 if (i == 0) return
117
118 ! copy H
119 n=i
120 allocate(h1(n+1,n))
121 do j = 1, n
122 do k = 1, j+1
123 h1(k,j) = h(k,j)
124 enddo
125 do k = j+2, n
126 h1(k,j) = 0.d0
127 enddo
128 enddo
129 ldh=n+1
130 allocate(wr(n))
131
132
133 ! !!
134 ! !! dgesvd version
135 ! !!
136
137 ! ! arguments for calling dgesvd
138 ! JOBU='N'
139 ! JOBVT='N'
140 ! ! estimate optimal LWORK
141 ! allocate(WORK(1))
142 ! LWORK=-1
143 ! call dgesvd(JOBU,JOBVT,N,N,H1,LDH,WR,Z,LDZ,Z,LDZ,WORK,LWORK,INFO)
144 ! if (INFO /= 0) then
145 ! write(*,*) 'ERROR: dgesvd returned with INFO=',INFO
146 ! return
147 ! endif
148 ! ! calculate singular values
149 ! LWORK=WORK(1)
150 ! deallocate(WORK)
151 ! allocate(WORK(LWORK))
152 ! call dgesvd(JOBU,JOBVT,N,N,H1,LDH,WR,Z,LDZ,Z,LDZ,WORK,LWORK,INFO)
153 ! if (INFO /= 0) then
154 ! write(*,*) 'ERROR: dgesvd returned with INFO=',INFO
155 ! return
156 ! endif
157 ! deallocate(WORK)
158
159
160 !!
161 !! dgesdd version (faster but need more workspace than dgesvd)
162 !!
163
164 ! arguments for calling dgesdd
165 jobz='N'
166 allocate(iwork(8*n))
167 ! estimate optimal LWORK
168 allocate(work(1))
169 lwork=-1
170 call dgesdd(jobz,n,n,h1,ldh,wr,z,ldz,z,ldz,work,lwork,iwork,info)
171 if (info /= 0) then
172 write(*,*) 'ERROR: dgesdd returned with INFO=',info
173 return
174 endif
175 ! calculate singular values
176 lwork=work(1)
177 deallocate(work)
178 allocate(work(lwork))
179 call dgesdd(jobz,n,n,h1,ldh,wr,z,ldz,z,ldz,work,lwork,iwork,info)
180 if (info /= 0) then
181 write(*,*) 'ERROR: dgesdd returned with INFO=',info
182 return
183 endif
184 deallocate(work)
185 deallocate(iwork)
186
187
188 write(*,'("emin=",1pe13.6,", emax=",1pe13.6,", emax/emin=",1pe13.6)') &
189 wr(n), wr(1), wr(1)/wr(n)
190
191 deallocate(wr)
192 deallocate(h1)
193#endif
194 end subroutine hecmw_estimate_condition_gmres
195
subroutine hecmw_estimate_condition_cg(iter, d, e)
subroutine hecmw_estimate_condition_gmres(i, h)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal