mvnrnd.f90

SUBROUTINE MVNRND(MU,SIGMA,R,CASES,M1,N1,M,N,SEED)

! ----------------------------------------------------------------------
!
! GENERATE RANDOM NUMBERS FROM A MULTIVARIATE RANDOM NORMAL DISTRIBUTION WITH
! SPECIFIED ERROR STATISTICS. IF ANY OFF-DIAGONAL ELEMENTS OF THE COVARIANCE
! MATRIX ARE NON-ZERO, A CHOLESKY SUBROUTINE IS CALLED TO FACTORIZE THE 
! COVARIANCE.  IF ALL OFF-DIAGONAL ELEMENTS ARE ZERO, THE FACTOR IS COMPUTED
! BY TAKING THE SQUARE ROOT OF THE DIAGONAL ELEMENTS. 

! MU - VECTOR OF MEANS OF RANDOM NUMBERS OF WHICH VARIATES WILL BE GENERATED
! SIGMA - COVARIANCES OF RANDOM NUMBERS; FOR THIS VERSION OF THE CODE, THIS 
!           MUST BE A SQUARE MATRIX WITH NO OFF-DIAGONAL NON-ZERO ELEMENTS.
! R - OUTPUT ARRAY OF VARIATES WITH MEAN MU AND VARIANCE SIGMA
! CASES - NUMBER OF VARIATES OF VECTOR OF RANDOM NUMBERS TO BE GENERATED
! M1 - NUMBER OF ROWS OF MU
! N1 - NUMBER OF COLS OF MU (SHOULD BE 1)
! M - NUMBER OF ROWS OF SIGMA
! N - NUMBER OF COLS OF SIGMA
! SEED - NUMBER TO SEED THE NORMAL RANDOM FUNCTION
!
! CALLS: RAND_NORMAL, CHOLESKY
! CALLED FROM: ENKFGEN, KUPDATE, WRITEFORCINGFILE
!
! VERSION HISTORY: ADAPTED FROM MATLAB - MD - 05/2005
!                  MODIFIED TO TAKE SEED - MD - 07/2005


IMPLICIT NONE

INTEGER,INTENT(IN)::M1,N1,M,N,CASES,SEED
REAL,INTENT(IN)::MU(M1,N1),SIGMA(M,N)
INTEGER C,I,J
REAL,DIMENSION(:),ALLOCATABLE::Z
REAL,DIMENSION(:,:),ALLOCATABLE::MU2,MU3,T,ZR,TEMP
REAL,INTENT(OUT)::R(CASES,MAX(M1,N1))
LOGICAL :: USE_CHOL,C_IDEN

!MU IS M1XN1 AND SIGMA IS MXN

C=MAX(M1,N1)

IF(M1*N1.NE.C)THEN
  PRINT *,'ERROR IN MVNRND!  MU MUST BE A VECTOR'
  STOP   
END IF

IF(M.NE.N)THEN
  PRINT *,'ERROR IN MVNRND!  SIGMA MUST BE SQUARE'
  STOP
END IF

IF(M.NE.C)THEN
  PRINT *,'ERROR IN MVNRND!  THE LENGTH OF MU MUST EQUAL THE NUMBER',&
          'OF ROWS IN SIGMA'  
  STOP
END IF

ALLOCATE(T(M,N))

IF(M==1)THEN
  T=SIGMA**0.5
ELSE
  !DETERMINE WHETHER OR NOT THE COVARIANCE MATRIX HAS ANY OFF-DIAGONAL 
  !  ELEMENTS THAT ARE NON-ZERO.  IF SO, WE MUST USE THE CHOLESKY ROUTINE
  !  TO FACTORIZE THE COVARIANCE MATRIX.  IF ALL OFF-DIAGONAL ELEMENTS ARE 
  !  ZERO, OR IF THE COVARIANCE MATRIX CONSISTS OF IDENTICAL ENTRIES, TRIVIAL
  !  SOLUTIONS OF 'T' EXIST

  !CHECK FOR IDENTICAL COVARIANCE MATRIX; IF ALL ELEMENTS ARE EQUAL TO FIRST
  !  ELEMENT, THEN I'LL ASSUME IT'S IDENTICAL. 
  C_IDEN=.TRUE.
  DO I=1,M
    DO J=1,N
      IF(SIGMA(I,J).NE.SIGMA(1,1)) THEN
        C_IDEN=.FALSE.
      END IF
    END DO
  END DO

  IF(C_IDEN.EQ..TRUE.)THEN
    !IF THE MATRIX CONSISTS OF IDENTICAL ENTRIES, DEFINITELY DON'T USE CHOLESKY
    USE_CHOL=.FALSE. 
  ELSE
    !IF THE MATRIX HAS NON-IDENTICAL ENTRIES, LOOK FOR NON-DIAGONAL ELEMENTS 
    !  THAT ARE NON-ZERO
    USE_CHOL=.FALSE. 
    DO I=1,M
      DO J=1,N
        IF (SIGMA(I,J).GT.0.AND.I.NE.J)THEN
          !SINCE AT LEAST ONE NON-DIAGONAL ELEMENT IS NOT EQUAL TO ZERO, WE MUST
          !  USE THE CHOLESKY ROUTINE
          USE_CHOL=.TRUE.
        END IF
      END DO
    END DO
  END IF
  
  IF(USE_CHOL.EQ..FALSE.)THEN
    !DON'T USE CHOLESKY...
    IF(C_IDEN.EQ..TRUE.)THEN
      !SINCE ALL ELEMENTS OF COVARIANCE MATRIX ARE IDENTICAL, 'T' IS TRIVIAL:
      !  ALL ELEMENTS OF THE FIRST ROW ARE THE SQUARE ROOT OF ANY OF THE
      !  TERMS. 
      DO I=1,M
        DO J=1,N
          IF (I.EQ.1) THEN
            T(I,J)=SIGMA(1,1)**0.5
          ELSE
            T(I,J)=0.
          END IF
        END DO
      END DO
    ELSE
      !SINCE ALL OFF-DIAGONAL TERMS ARE ZERO, 'T' IS TRIVIAL: WE CAN COMPUTE 
      !  THE SQUARE ROOT OF OF THE DIAGONAL ONLY
      DO I=1,M
        DO J=1,N
          IF (I.EQ.J) THEN
            T(I,J)=SIGMA(I,I)**0.5
          ELSE
            T(I,J)=0.
          END IF
        END DO
      END DO
    END IF
  ELSE 
    !WE MUST USE CHOLESKY
    IF(M.NE.N)THEN
      PRINT *, 'MVNRND: CHOLESKY ONLY WORKS FOR SQUARE MATRICES! ABORTING' 
    ELSE
      CALL CHOLESKY(SIGMA,T,M)
    END IF
  END IF
END IF
 
IF(M1==C.AND.C>1) THEN !IE IF MU IS A COLUMN VECTOR, TRANSPOSE IT
   ALLOCATE(MU2(N1,M1))
   MU2=TRANSPOSE(MU)
ELSE
  ALLOCATE(MU2(M1,N1))
  MU2=MU
END IF

ALLOCATE(MU3(CASES,C))
DO I=1,C
  DO J=1,CASES
    MU3(J,I)=MU2(1,I)
  END DO
END DO

ALLOCATE(Z(CASES*M))
CALL RAND_NORMAL(CASES*M,Z,SEED)

ALLOCATE(ZR(CASES,M))
IF(M>1)THEN
  ZR=RESHAPE(Z,(/CASES,M/))
ELSE
  ZR(:,1)=Z(:)
END IF

R=MATMUL(ZR,T)+MU3

DEALLOCATE(T,MU2,MU3,Z,ZR)


END SUBROUTINE MVNRND