var/da/da_tools/da_1d_eigendecomposition.inc

   1 subroutine da_1d_eigendecomposition( bx, e, l )
   2
   3    !------------------------------------------------------------------------------
   4    !  Purpose: Compute eigenvectors E and eigenvalues L of vertical covariance matrix
   5    !           B_{x} defined by equation:  E^{T} B_{x} E = L, given input 3D field of
   6    !           errors (sum over all horizontal locations).
   7    !------------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    real, intent(in)         :: bx(:,:)          ! Global vert. background error.
  12
  13    real, intent(out)        :: e(:,:)           ! Eigenvectors of Bx.
  14    real, intent(out)        :: l(:)             ! Global eigenvalues of Bx.
  15
  16    integer                  :: kz               ! Size of 3rd dimension.
  17    integer                  :: m                ! Loop counters
  18    integer                  :: work             ! Size of work array.
  19    integer                  :: info             ! Info code.
  20
  21    real*8, allocatable        :: ecopy(:,:)
  22    real*8, allocatable        :: lcopy(:)
  23    real*8, allocatable        :: work_array(:)
  24
  25    if (trace_use_dull) call da_trace_entry("da_1d_eigendecomposition")
  26
  27    !-------------------------------------------------------------------------
  28    ! [1.0]: Initialise:
  29    !-------------------------------------------------------------------------
  30
  31    kz = size(bx, dim=1)
  32
  33    !-------------------------------------------------------------------------
  34    ! [5.0]: Perform global eigenvalue decomposition using LAPACK software:
  35    !-------------------------------------------------------------------------
  36
  37    work = 3 * kz - 1
  38    allocate( work_array(1:work) )
  39
  40    allocate( ecopy(1:kz,1:kz) )
  41    allocate( lcopy(1:kz) )
  42
  43    ecopy(:,:) = bx(:,:)
  44
  45    lcopy = 0.0
  46
  47    call dsyev( 'V', 'U', kz, ecopy, kz, lcopy, &
  48               work_array, work, info )
  49
  50    if ( info /= 0 ) then
  51       write(unit=message(1),fmt='(A,I4,A)') &
  52          ' da_1d_eigendecomposition: info = ', &
  53          info,' - error in decomposition.'
  54       call da_error(__FILE__,__LINE__,message(1:1))
  55    end if
  56
  57    !--Swap order of eigenvalues, vectors so 1st is one with most
  58    !  variance, etc:
  59
  60    do m=1,kz
  61       l(m) = lcopy(kz+1-m)
  62       e(1:kz,m) = ecopy(1:kz,kz+1-m)
  63    end do
  64
  65    deallocate (work_array)
  66    deallocate (ecopy)
  67    deallocate (lcopy)
  68
  69    if (trace_use_dull) call da_trace_exit("da_1d_eigendecomposition")
  70
  71 end subroutine da_1d_eigendecomposition
  72
  73