Fortran 90 SUBROUTINE eof_scalar
!>============================================================================<
! FILE: /fs/cgd/home0/davestep/f90_routines/SUB_eof_scalar.f90
!
! AUTHOR: David Stepaniak, NCAR/CGD/CAS
! DATE INITIATED: Mon Jan 26 13:15:21 MST 1998
! LAST MODIFIED: Fri Mar 13 08:31:36 MST 1998
!
! SYNOPSIS: Computes EOFs, eigenvalues, and coefficients of EOFs for
! a scalar data field.
!
! DESCRIPTION: To fix ideas, suppose that one has constructed a data
! matrix consisting of M rows, each row corresponding to
! a spatial grid point, and N columns, each column
! corresponding to a temporal observation. Each row is
! simply a time series of observations for a particular
! grid point, and each column is simply a set of grid
! point observations at a particular time. We may call
! this data matrix F.
!
! The time series of F may consist of raw data, or anomalies
! computed using a climatology, or some other user-defined
! form.
!
! In what follows it will be assumed that, as a final step,
! the overall mean of each time series of F has been subtracted
! from its corresponding time series so as to form deviations
! from the mean. (This allows efficient matrix based
! calculation of the variance-covariance matrix and its
! relatives.) The matrix of deviations will be denoted by Fd.
!
! The M Empirical Orthogonal Functions (EOFs) of this data
! Fd may be obtained by computing the eigenvectors of the MxM
! variance-covariance matix R = Fd * Fd_T / N, where Fd_T is
! the transpose of Fd. Since it is most often the case that
! M >> N with M numbering in the thousands, the task of
! computing the eigenvectors (and eigenvalues) of R is
! computationally prohibitive.
!
! An alternative and efficient method put forth by von
! Storch and Hannoschock (1984) computes the N eigenvalues
! and eigenvectors of the NxN matrix formed by Fd_T * Fd / N,
! a much less computationally intensive exercise. From this
! NxN eigensystem it is straight forward to obtain the first
! N of M eigenvalues and eigenvectors of the original MxM
! system where the corresponding pairs of eigenvalues and
! eigenvectors are arranged in descending order according to
! eigenvalue. This subroutine implements the von Storch and
! Hannoschock algorithm.
!
! INPUT: ucM. Number of spatial grid points.
! ucN. Number of temporal observations.
! ucFd. Data matrix in which rows are time series of
! deviations from their corresponding overall time
! mean.
!
! OUTPUT: ucE. EOFs arranged column by column in order of their
! corresponding set of eigenvalues.
! ucL_N1. Eigenvalues arranged in descending order.
! ucC. Coefficent matrix of EOFs.
!
! NOTES: Input and output arguments are more fully described in the
! argument list below.
!
! Time series containing missing values are not allowed. It
! is left to the user to make provisions for time series
! containing missing values. This usually involves employing
! a grid point mask to assist in the packing and bookkeeping
! of the input data array and the unpacking and bookkeeping
! of the output EOF array.
!
! It is assumed that M > N.
!
! Utilizes SSYEV routine from LAPACK to compute eigenvalues
! and eigenvectors of a real symmetric matrix. On winterpark
! compile with -r4 and -lcomplib.sgimath flags where r4
! indicates 32 bit (4 byte) storage of real variables and
! and complib.sgimath is the Silicon Graphics Scientific
! Mathematical Library which contains LAPACK routines among
! others.
!
! REFERENCE: von Storch, H., and G. Hannoschock (1984) Comments on
! "Emipirical Orthogonal Function Analysis of Wind Vectors
! Over the Tropical Pacific Ocean". Bulletin of the
! Meteorological Society of America, 65, 162.
!
! [Appeared as a letter to the editor concerning: Legler,
! D. M. (1983) Empirical Orthogonal Function Analysis of Wind
! Vectors over the Tropical Pacific Region. Bulletin of the
! Meteorological Society of America, 64, 234-241.]
!
!>============================================================================<
SUBROUTINE eof_scalar( ucM, ucN, ucFd, ucE, ucL_N1, ucC )
IMPLICIT NONE
!>----------------------------------------------------------------------------<
! Argument list in order of appearance:
INTEGER, INTENT(IN) :: ucM
! Number of spatial grid points to be
! represented in MxN data matrix. (ucM
! refers to upper case M.)
INTEGER, INTENT(IN) :: ucN
! Number of temporal observations of a
! time series. (ucN refers to uppercase
! N.) It is assumed that M > N.
REAL, DIMENSION( 1:ucM, 1:ucN ), INTENT(IN) :: ucFd
! MxN data array in which there are M
! rows (M grid points) each of which
! contains N columns (N temporal
! observations). Thus each row is
! simply a time series of observations
! at a particular grid point.
!---------------- NOTE! ---------------+
! Irregardless of how the data (in |
! particular the time series) are |
! prepared, AS A FINAL STEP the overall|
! mean of each time series should be |
! subtracted from its corresponding |
! time series so as to form deviations |
! from the mean. This allows for |
! efficient matrix-based calculation of|
! the variance-covariance matrix and |
! its relatives. |
!--------------------------------------+
REAL, DIMENSION( 1:ucM, 1:ucN ), INTENT(OUT) :: ucE
! MxN array containing N EOFs
! (eigenvectors.)
! Each EOF has M components which
! correspond to the M spatial grid
! points. The EOFs are stored in
! columns in the same order as their
! corresponding eigenvalues arranged
! from greatest to least.
REAL, DIMENSION( 1:ucN ), INTENT(OUT) :: ucL_N1
! Nx1 one-dimensional array containing
! eigenvalues arranged from greatest
! to least.
REAL, DIMENSION( 1:ucN, 1:ucN ), INTENT(OUT) :: ucC
! NxN array containing raw coefficients
! of EOFs. The raw coefficients
! form time series of the temporal
! behavoir of each EOF, and are stored
! row by row with the first row
! corresponding to the first EOF and so
! on.
!>----------------------------------------------------------------------------<
! Declare useful transposes:
REAL, DIMENSION(:,:), ALLOCATABLE :: ucFd_T
! NxM transpose array of ucFd.
REAL, DIMENSION(:,:), ALLOCATABLE :: ucE_T
! NxM transpose array of ucE.
!>----------------------------------------------------------------------------<
! Declare arrays and scalars to implement von Storch and Hannoschock (1984)
! algorithm:
REAL, DIMENSION(:,:), ALLOCATABLE :: ucQ
! NxN array formed by the multiplication
! ucFd_T * ucFd / ucN => ucQ. Thus
! [NxM] * [MxN] => [NxN].
! Note that ucQ is a real symmetric
! matrix.
REAL, DIMENSION(:,:), ALLOCATABLE :: ucZ
! NxN array containing N eigenvectors
! of ucQ. The eigenvectors each have N
! components, and are stored in
! columns in the same order as their
! corresponding eigenvalues arranged
! from greatest to least.
INTEGER :: col
! Simply a DO loop control variable for
! incrementing over columns.
REAL, DIMENSION(:), ALLOCATABLE :: ucFdz
! Mx1 column vector (array) formed by the
! multiplication of ucFd and a single
! column (eigenvector) of ucZ. Thus
! ucFd * ucZ(:,col) => ucFdz,
! [MxN]* [Nx1] => [Mx1]
REAL :: euc_norm_ucFdz
! Scalar. Euclidean norm of ucFdz. That
! is, the square root of the sum of the
! element by element multiplication of
! ucFdz and ucFdz.
!>----------------------------------------------------------------------------<
! Arguments for LAPACK SSYEV routine. SSYEV computes the eigenvectors and
! eigenvalues of a real symmetric matrix.
CHARACTER(1) :: arg_JOBZ
! Flag for computing only eigenvalues,
! or eigenvalues and eigenvectors.
CHARACTER(1) :: arg_UPLO
! Flag indicating storage of upper or
! lower triangle of arg_A.
INTEGER :: arg_N
! Number of rows (or columns) of arg_A.
REAL, DIMENSION(:,:), ALLOCATABLE :: arg_A
! On entry, real symmetric matrix
! arg_A. On exit, and when computing
! both eigenvalues and eigenvectors,
! contains the orthonormal
! eigenvectors of arg_A with the ith
! column of arg_A holding the
! eigenvector associated with
! arg_W(i).
INTEGER :: arg_LDA
! A parameter (typically N).
REAL, DIMENSION(:), ALLOCATABLE :: arg_W
! Contains the eigenvalues of arg_A in
! ascending order
REAL, DIMENSION(:), ALLOCATABLE :: arg_WORK
! One-dimensional work array.
INTEGER :: arg_LWORK
! A parameter (typically 3*N-1).
INTEGER :: arg_INFO
! Error trapping parameter diagnosed in
! code as below.
! +------------------------+
!>------------------------| BEGIN EXECUTABLE CODE: |--------------------------<
! +------------------------+
!
! Form the NxN array Q = ucFd_T * ucFd / ucN:
ALLOCATE( ucFd_T( 1:ucN, 1:ucM ) )
ucFd_T = TRANSPOSE( ucFd )
ALLOCATE( ucQ( 1:ucN, 1:ucN ) )
ucQ = MATMUL( ucFd_T, ucFd ) / REAL( ucN )
DEALLOCATE( ucFd_T )
!>----------------------------------------------------------------------------<
! Assign arguments for LAPACK SSYEV routine:
arg_JOBZ = "V" ! Compute both eigenvalues and
! eigenvectors.
arg_UPLO = "L" ! Store lower triangle of arg_A.
arg_N = ucN ! arg_N is simply ucN.
ALLOCATE( arg_A( 1:ucN, 1:ucN ) )
arg_A = ucQ ! Assign arg_A as ucQ.
DEALLOCATE( ucQ )
arg_LDA = ucN ! arg_LDA is simply ucN.
ALLOCATE( arg_W( 1:ucN ))
arg_LWORK = 3*ucN - 1
! arg_LWORK is simply 3*ucN - 1.
ALLOCATE( arg_WORK( 1:arg_LWORK ) )
arg_INFO = 0 ! arg_info initialized as zero.
!>----------------------------------------------------------------------------<
!> Compute eigenvalues and eigenvectors of ucQ:
CALL SSYEV( arg_JOBZ, &
arg_UPLO, &
arg_N, &
arg_A, &
arg_LDA, &
arg_W, &
arg_WORK, &
arg_LWORK, &
arg_INFO &
)
IF ( arg_INFO < 0 ) THEN
WRITE (*,* ) "FATAL ERROR CONDITION: "
WRITE (*,100) -arg_INFO,"th argument of SSYEV "
WRITE (*,*) "had an illegal value. "
WRITE (*,*) "Execution halted in subroutine eof_scalar."
STOP
ELSE IF ( arg_INFO > 0 ) THEN
WRITE (*,*) "FATAL ERROR CONDITION: "
WRITE (*,*) "LAPACK SSYEV algorithm failed to converge "
WRITE (*,200) "after finding only ",arg_INFO-1, &
" eigenvalues. "
WRITE (*,*) "Execution halted in subroutine eof_scalar."
STOP
END IF
DEALLOCATE( arg_WORK )
!>----------------------------------------------------------------------------<
! Extract ucE (EOFs) from ucZ:
ALLOCATE( ucZ( 1:ucN, 1:ucN ) )
ucZ(:,1:ucN) = arg_A(:,ucN:1:-1)
! Arrange eigenvectors of ucZ so as to
! have same order as corresponding set
! of descending eigenvalues.
DEALLOCATE( arg_A )
ALLOCATE( ucFdz( 1:ucM ) )
DO col = 1, ucN
! Loop over all eigenvectors of ucZ.
ucFdz(:) = MATMUL( ucFd(:,:), ucZ(:,col) )
! Compute ucFd * ucZ for a single
! eigenvector of ucZ.
euc_norm_ucFdz = SQRT( &
DOT_PRODUCT( ucFdz(:),ucFdz(:) ) &
)
! Compute Euclidean norm of ucFdz.
ucE(:,col) = ucFdz(:) / euc_norm_ucFdz
! Normalize ucFdz by its Euclidean
! norm to get corresponding Mx1
! eigenvector (EOF) of ucE.
END DO
DEALLOCATE( ucZ )
DEALLOCATE( ucFdz )
ALLOCATE( ucE_T( 1:ucN,1:ucM ) )
ucE_T = TRANSPOSE( ucE )
ucC = MATMUL( ucE_T, ucFd )
! Compute coefficient matrix.
DEALLOCATE( ucE_T )
ucL_N1(:) = arg_W(ucN:1:-1)
! Assign eigenvalues in descending order.
DEALLOCATE( arg_W )
!>----------------------------------------------------------------------------<
! Format statements:
100 FORMAT( I1,A23 )
200 FORMAT( A19,I4,A37 )
!>----------------------------------------------------------------------------<
END SUBROUTINE eof_scalar
!>----------------------------------------------------------------------------<