! Copyright (c) 2004-2024 Lars Nerger ! ! This file is part of PDAF. ! ! PDAF is free software: you can redistribute it and/or modify ! it under the terms of the GNU Lesser General Public License ! as published by the Free Software Foundation, either version ! 3 of the License, or (at your option) any later version. ! ! PDAF is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public ! License along with PDAF. If not, see <http://www.gnu.org/licenses/>. ! !$Id$ !BOP ! ! !ROUTINE: PDAF_enkf_Omega --- Generate random matrix with special properties ! ! !INTERFACE: SUBROUTINE PDAF_enkf_Omega(seed, r, dim_ens, Omega, norm, & otype, screen) ! !DESCRIPTION: ! Generate a random matrix OMEGA for the initilization ! of ensembles for the EnKF. ! ! The following properties can be set:\\ ! 1. Simply fill the matrix with random numbers from a ! Gaussian distribution with mean zero and unit ! variance. (This corresponds to the simple random ! initialization of EnKF.)\\ ! 2. Constrain the columns of OMEGA to be of unit norm ! (This corrects error in the variance estimates caused ! by the finite number of random numbers.)\\ ! 3. Constrain the columns of OMEGA to be of norm dim\_ens$^{-1/2}$ ! (This corrects variance errors as in 2)\\ ! 4. Project columns of OMEGA to be orthogonal to the vector ! $(1,....,1)^T$ by Householder reflections. (This assures ! that the mean of the generated ensemble equals the ! prescribed state estimate.)\\ ! Property 4 can be combined with either property 2 or 3. ! ! ! This is a core routine of PDAF and ! should not be changed by the user ! ! ! !REVISION HISTORY: ! 2005-04 - Lars Nerger - Initial code ! Later revisions - see svn log ! ! !USES: ! Include definitions for real type of different precision ! (Defines BLAS/LAPACK routines and MPI_REALTYPE) #include "typedefs.h" IMPLICIT NONE ! !ARGUMENTS: INTEGER, INTENT(in) :: seed(4) ! Seed for random number generation INTEGER, INTENT(in) :: r ! Approximated rank of covar matrix INTEGER, INTENT(in) :: dim_ens ! Ensemble size REAL, INTENT(inout) :: Omega(dim_ens,r) ! Random matrix REAL, INTENT(inout) :: norm ! Norm for ensemble transformation INTEGER, INTENT(in) :: otype ! Type of Omega: ! (1) Simple Gaussian random matrix ! (2) Columns of unit norm ! (3) Columns of norm dim_ens^(-1/2) ! (4) Projection orthogonal (1,..,1)^T ! (6) Combination of 2 and 4 ! (7) Combination of 3 and 4 ! (8) Rows of sum 0 and variance 1 INTEGER, INTENT(in) :: screen ! Control verbosity ! !CALLING SEQUENCE: ! Called by: Used-provided ensemble initialization routine for EnKF ! Calls: gemmTYPE (BLAS; dgemm or sgemm dependent on precision) ! Calls: larnvTYPE (BLAS; dlarnv or slarnv dependent on precision) !EOP ! *** local variables *** INTEGER :: i, j ! Counters INTEGER, SAVE :: iseed(4) ! seed array for random number routine REAL, ALLOCATABLE :: rndvec(:) ! Vector of random numbers REAL, ALLOCATABLE :: house(:,:) ! Householder matrix REAL, ALLOCATABLE :: Omega_tmp(:,:) ! Temporary OMEGA for projection REAL :: colnorm ! Norm of matrix column REAL :: rownorm ! Norm of matrix row ! ********************** ! *** INITIALIZATION *** ! ********************** IF (seed(1) >= 0) THEN ! Use given seed iseed=seed END IF ALLOCATE(rndvec(r)) ! *** Initialize random matrix *** DO i = 1, dim_ens ! Fill row-wise to be consistent with old sampling formulation CALL larnvTYPE(3, seed, r, rndvec) Omega(i, :) = rndvec END DO ! *** Normalize columns of Omega *** normcols: IF (otype == 2 .OR. otype == 3 .OR. otype == 6 .OR. otype == 7) THEN IF ((screen > 0) .AND. (otype == 2 .OR. otype == 6)) THEN WRITE (*, '(a, 5x, a)') 'PDAF', '--- EnKF_Omega: Normalize columns of random matrix' ELSE IF (screen > 0) THEN WRITE (*, '(a, 5x,a)') & 'PDAF', '--- EnKF_Omega: Normalize columns of random matrix to dim_ens^(-1/2)' END IF DO j = 1, r ! Compute norm colnorm = 0.0 DO i = 1, dim_ens colnorm = colnorm + Omega(i, j)**2 END DO IF (otype == 3 .OR. otype == 7) THEN ! Set column norm to 1/sqrt(dim_ens) colnorm = colnorm / REAL(dim_ens) END IF colnorm = SQRT(colnorm) ! Perform normalization DO i = 1, dim_ens Omega(i, j) = Omega(i, j) / colnorm END DO END DO END IF normcols ! *** Project columns orthogonal to (1,1,...,1)^T *** doproject: IF (otype == 4 .OR. otype == 6 .OR. otype == 7) THEN IF (screen > 0) & WRITE (*, '(a, 5x, a)') & 'PDAF', '--- EnKF_Omega: Project columns orthogonal to (1,...,1)^T' ALLOCATE(Omega_tmp(dim_ens, r)) ALLOCATE(house(dim_ens, dim_ens)) ! Store Omega Omega_tmp = Omega ! Initialize Householder matrix housecolb: DO j = 1, dim_ens houserowb: DO i = 1, dim_ens house(i, j) = -1.0 / REAL(dim_ens) END DO houserowb END DO housecolb DO j = 1, dim_ens house(j, j) = house(j, j) + 1.0 END DO ! Perform reflection CALL gemmTYPE ('n', 'n', dim_ens, r, dim_ens, & 1.0, house, dim_ens, Omega_tmp, dim_ens, & 0.0, Omega, dim_ens) DEALLOCATE(Omega_tmp, house) END IF doproject rowzero: IF (otype == 8) THEN IF (screen > 0) & WRITE (*, '(a, 5x, a)') & 'PDAF', '--- EnKF_Omega: Ensure that row sums are zero' DO i = 1, dim_ens rownorm = 0.0 DO j = 1, r rownorm = rownorm + Omega(i,j) ENDDO rownorm = rownorm / REAL(r) DO j = 1, r Omega(i,j) = Omega(i,j) - rownorm ENDDO END DO IF (screen > 0) & WRITE (*, '(a, 5x, a)') & 'PDAF', '--- EnKF_Omega: Ensure that variance in rows is one' DO i = 1, dim_ens rownorm = 0.0 DO j = 1, r rownorm = rownorm + Omega(i,j)*Omega(i,j) ENDDO rownorm = rownorm / REAL(r-1) rownorm = SQRT(rownorm) DO j = 1, r Omega(i,j) = Omega(i,j) / rownorm ENDDO END DO END IF rowzero ! *** Initialize norm for ensemble transformation *** IF (otype == 1) THEN norm = 1.0 ELSEIF (otype == 2) THEN norm = SQRT(REAL(dim_ens - 1)) ELSEIF (otype == 3) THEN norm = 1.0 ELSEIF (otype == 4) THEN norm = 1.0 ELSEIF (otype == 6) THEN norm = SQRT(REAL(dim_ens - 1)) ELSEIF (otype == 7) THEN norm = 1.0 END IF ! ******************** ! *** Finishing up *** ! ******************** DEALLOCATE(rndvec) END SUBROUTINE PDAF_enkf_Omega