FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
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 
11 contains
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 
196 end module hecmw_estimate_condition
subroutine hecmw_estimate_condition_gmres(I, H)
subroutine hecmw_estimate_condition_cg(ITER, D, E)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal