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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
SUBROUTINE setgmn(meanv,covm,ldcovm,p,parm,ierr)
C SUBROUTINE setgmn(meanv,covm,p,parm)
C JJV changed this routine to take leading dimension of COVM
C JJV argument and pass it to SPOFA, making it easier to use
C JJV if the COVM which is used is contained in a larger matrix
C JJV and to make the routine more consistent with LINPACK.
C JJV Changes are in comments, declarations, and the call to SPOFA.
C**********************************************************************
C
C SUBROUTINE SETGMN( MEANV, COVM, LDCOVM, P, PARM)
C SET Generate Multivariate Normal random deviate
C
C
C Function
C
C
C Places P, MEANV, and the Cholesky factoriztion of COVM
C in PARM for GENMN.
C
C
C Arguments
C
C
C MEANV --> Mean vector of multivariate normal distribution.
C DOUBLE PRECISION MEANV(P)
C
C COVM <--> (Input) Covariance matrix of the multivariate
C normal distribution. This routine uses only the
C (1:P,1:P) slice of COVM, but needs to know LDCOVM.
C
C (Output) Destroyed on output
C DOUBLE PRECISION COVM(LDCOVM,P)
C
C LDCOVM --> Leading actual dimension of COVM.
C INTEGER LDCOVM
C
C P --> Dimension of the normal, or length of MEANV.
C INTEGER P
C
C PARM <-- Array of parameters needed to generate multivariate
C normal deviates (P, MEANV and Cholesky decomposition
C of COVM).
C 1 : 1 - P
C 2 : P + 1 - MEANV
C P+2 : P*(P+3)/2 + 1 - Cholesky decomposition of COVM
C DOUBLE PRECISION PARM(P*(P+3)/2 + 1)
C
C**********************************************************************
C .. Scalar Arguments ..
C INTEGER p
include 'stack.h'
INTEGER p, ldcovm
C ..
C .. Array Arguments ..
C DOUBLE PRECISION covm(p,p),meanv(p),parm(p* (p+3)/2+1)
DOUBLE PRECISION covm(ldcovm,p),meanv(p),parm(p* (p+3)/2+1)
C ..
C .. Local Scalars ..
INTEGER i,icount,info,j
C ..
C .. External Subroutines ..
EXTERNAL dpofa
C ..
C .. Executable Statements ..
C
C
C TEST THE INPUT
C
10 parm(1) = p
C
C PUT P AND MEANV INTO PARM
C
DO 20,i = 2,p + 1
parm(i) = meanv(i-1)
20 CONTINUE
C
C Cholesky decomposition to find A s.t. trans(A)*(A) = COVM
C
C CALL dpofa(covm,p,p,info)
CALL dpofa(covm,ldcovm,p,info)
ierr=0
IF (.NOT. (info.NE.0)) GO TO 30
call basout(io,wte,"Rand: COV not positive definite")
ierr=1
return
30 icount = p + 1
C
C PUT UPPER HALF OF A, WHICH IS NOW THE CHOLESKY FACTOR, INTO PARM
C COVM(1,1) = PARM(P+2)
C COVM(1,2) = PARM(P+3)
C :
C COVM(1,P) = PARM(2P+1)
C COVM(2,2) = PARM(2P+2) ...
C
DO 50,i = 1,p
DO 40,j = i,p
icount = icount + 1
parm(icount) = covm(i,j)
40 CONTINUE
50 CONTINUE
RETURN
C
END
|
