hecmw_estimate_condition.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module hecmw_estimate_condition
 
  private
  public :: hecmw_estimate_condition_cg
  public :: hecmw_estimate_condition_gmres
 
contains
 
  subroutine hecmw_estimate_condition_cg(ITER, D, E)
    use hecmw_util
    implicit none
    integer(kind=kint), intent(in) :: iter
    real(kind=kreal), intent(in) :: d(:), e(:)
 
#ifdef HECMW_WITH_LAPACK
    character(len=1) :: jobz, range
    ! character(len=1) :: COMPZ
    real(kind=kreal) :: vl, vu, abstol, z(1,1)
    integer(kind=kint) :: n, il, iu, m, ldz=1, isuppz(1)
    integer(kind=kint) :: lwork, liwork, info
    real(kind=kreal), allocatable :: w(:), work(:)
    integer(kind=kint), allocatable :: iwork(:)
    real(kind=kreal), allocatable :: d1(:), e1(:)
    integer(kind=kint) :: i
 
    if (iter <= 1) return
 
    ! copy D, E
    allocate(d1(iter),e1(iter))
    do i=1,iter-1
      d1(i) = d(i)
      e1(i) = e(i)
    enddo
    d1(iter) = d(iter)
 
 
    !!
    !! dstegr version (faster than dsteqr)
    !!
 
    ! prepare arguments for calling dstegr
    jobz='N'
    range='A'
    n=iter
    allocate(w(iter))
    ! estimate optimal LWORK and LIWORK
    lwork=-1
    liwork=-1
    allocate(work(1),iwork(1))
    call dstegr(jobz,range,n,d1,e1,vl,vu,il,iu,abstol, &
      m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info)
    if (info /= 0) then
      write(*,*) 'ERROR: dstegr returned with INFO=',info
      return
    endif
    ! calculate eigenvalues
    lwork=work(1)
    liwork=iwork(1)
    deallocate(work,iwork)
    allocate(work(lwork),iwork(liwork))
    call dstegr(jobz,range,n,d1,e1,vl,vu,il,iu,abstol, &
      m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info)
    if (info /= 0) then
      write(*,*) 'ERROR: dstegr returned with INFO=',info
      return
    endif
    write(*,'("emin=",1pe13.6,", emax=",1pe13.6,", emax/emin=",1pe13.6)') &
      w(1),w(n),w(n)/w(1)
    deallocate(work,iwork)
    deallocate(w)
 
 
    ! !!
    ! !! dsteqr version
    ! !!
 
    ! ! prepare arguments for calling dsteqr
    ! COMPZ='N'
    ! N=ITER
    ! allocate(WORK(1))
    ! ! calculate eigenvalues
    ! call dsteqr(COMPZ,N,D1,E1,Z,LDZ,WORK,INFO)
    ! if (INFO /= 0) then
    !   write(*,*) 'ERROR: dsteqr returned with INFO=',INFO
    !   return
    ! endif
    ! write(*,'("emin=",1pe13.6,", emax=",1pe13.6,", emax/emin=",1pe13.6)') &
      !      D1(1),D1(N),D1(N)/D1(1)
    ! deallocate(WORK)
 
 
    deallocate(d1,e1)
#endif
  end subroutine hecmw_estimate_condition_cg
 
 
  subroutine hecmw_estimate_condition_gmres(I, H)
    use hecmw_util
    implicit none
    integer(kind=kint), intent(in) :: i
    real(kind=kreal), intent(in) :: h(:,:)
 
#ifdef HECMW_WITH_LAPACK
    ! character(len=1) :: JOBU, JOBVT
    character(len=1) :: jobz
    integer(kind=kint) :: n, ldh, ldz=1, lwork, info
    real(kind=kreal), allocatable :: wr(:), work(:), h1(:,:)
    integer(kind=kint), allocatable :: iwork(:)
    real(kind=kreal) :: z(1,1)
    integer(kind=kint) :: j, k
 
    if (i == 0) return
 
    ! copy H
    n=i
    allocate(h1(n+1,n))
    do j = 1, n
      do k = 1, j+1
        h1(k,j) = h(k,j)
      enddo
      do k = j+2, n
        h1(k,j) = 0.d0
      enddo
    enddo
    ldh=n+1
    allocate(wr(n))
 
 
    ! !!
    ! !! dgesvd version
    ! !!
 
    ! ! arguments for calling dgesvd
    ! JOBU='N'
    ! JOBVT='N'
    ! ! estimate optimal LWORK
    ! allocate(WORK(1))
    ! LWORK=-1
    ! call dgesvd(JOBU,JOBVT,N,N,H1,LDH,WR,Z,LDZ,Z,LDZ,WORK,LWORK,INFO)
    ! if (INFO /= 0) then
    !   write(*,*) 'ERROR: dgesvd returned with INFO=',INFO
    !   return
    ! endif
    ! ! calculate singular values
    ! LWORK=WORK(1)
    ! deallocate(WORK)
    ! allocate(WORK(LWORK))
    ! call dgesvd(JOBU,JOBVT,N,N,H1,LDH,WR,Z,LDZ,Z,LDZ,WORK,LWORK,INFO)
    ! if (INFO /= 0) then
    !   write(*,*) 'ERROR: dgesvd returned with INFO=',INFO
    !   return
    ! endif
    ! deallocate(WORK)
 
 
    !!
    !! dgesdd version (faster but need more workspace than dgesvd)
    !!
 
    ! arguments for calling dgesdd
    jobz='N'
    allocate(iwork(8*n))
    ! estimate optimal LWORK
    allocate(work(1))
    lwork=-1
    call dgesdd(jobz,n,n,h1,ldh,wr,z,ldz,z,ldz,work,lwork,iwork,info)
    if (info /= 0) then
      write(*,*) 'ERROR: dgesdd returned with INFO=',info
      return
    endif
    ! calculate singular values
    lwork=work(1)
    deallocate(work)
    allocate(work(lwork))
    call dgesdd(jobz,n,n,h1,ldh,wr,z,ldz,z,ldz,work,lwork,iwork,info)
    if (info /= 0) then
      write(*,*) 'ERROR: dgesdd returned with INFO=',info
      return
    endif
    deallocate(work)
    deallocate(iwork)
 
 
    write(*,'("emin=",1pe13.6,", emax=",1pe13.6,", emax/emin=",1pe13.6)') &
      wr(n), wr(1), wr(1)/wr(n)
 
    deallocate(wr)
    deallocate(h1)
#endif
  end subroutine hecmw_estimate_condition_gmres
 
end module hecmw_estimate_condition