summaryrefslogtreecommitdiff
path: root/modules/randlib/src/fortran/genmn.f
blob: a73982f0c2885087214508ef62b1686550be5e7f (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
      SUBROUTINE genmn(parm,x,work)
C**********************************************************************
C
C     SUBROUTINE GENMN(PARM,X,WORK)
C              GENerate Multivariate Normal random deviate
C
C
C                              Arguments
C
C
C     PARM --> Parameters needed to generate multivariate normal
C               deviates (MEANV and Cholesky decomposition of
C               COVM). Set by a previous call to SETGMN.
C               1 : 1                - size of deviate, P
C               2 : P + 1            - mean vector
C               P+2 : P*(P+3)/2 + 1  - upper half of cholesky
C                                       decomposition of cov matrix
C                                             DOUBLE PRECISION PARM(*)
C
C     X    <-- Vector deviate generated.
C                                             DOUBLE PRECISION X(P)
C
C     WORK <--> Scratch array
C                                             DOUBLE PRECISION WORK(P)
C
C
C                              Method
C
C
C     1) Generate P independent standard normal deviates - Ei ~ N(0,1)
C
C     2) Using Cholesky decomposition find A s.t. trans(A)*A = COVM
C
C     3) trans(A)E + MEANV ~ N(MEANV,COVM)
C
C**********************************************************************
C     .. Array Arguments ..
      DOUBLE PRECISION parm(*),work(*),x(*)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ae
      INTEGER i,icount,j,p
C     ..
C     .. External Functions ..
      DOUBLE PRECISION snorm
      EXTERNAL snorm
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC int
C     ..
C     .. Executable Statements ..
      p = int(parm(1))
C
C     Generate P independent normal deviates - WORK ~ N(0,1)
C
      DO 10,i = 1,p
          work(i) = snorm()
   10 CONTINUE
      DO 30,i = 1,p
C
C     PARM (P+2 : P*(P+3)/2 + 1) contains A, the Cholesky
C      decomposition of the desired covariance matrix.
C          trans(A)(1,1) = PARM(P+2)
C          trans(A)(2,1) = PARM(P+3)
C          trans(A)(2,2) = PARM(P+2+P)
C          trans(A)(3,1) = PARM(P+4)
C          trans(A)(3,2) = PARM(P+3+P)
C          trans(A)(3,3) = PARM(P+2-1+2P)  ...
C
C     trans(A)*WORK + MEANV ~ N(MEANV,COVM)
C
          icount = 0
          ae = 0.0
          DO 20,j = 1,i
              icount = icount + j - 1
              ae = ae + parm(i+ (j-1)*p-icount+p+1)*work(j)
   20     CONTINUE
          x(i) = ae + parm(i+1)
   30 CONTINUE
      RETURN
C
      END