GNU Octave  4.0.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
setgmn.f
Go to the documentation of this file.
1  SUBROUTINE setgmn(meanv,covm,ldcovm,p,parm)
2 C SUBROUTINE setgmn(meanv,covm,p,parm)
3 C JJV changed this routine to take leading dimension of COVM
4 C JJV argument and pass it to SPOTRF, making it easier to use
5 C JJV if the COVM which is used is contained in a larger matrix
6 C JJV and to make the routine more consistent with LAPACK.
7 C JJV Changes are in comments, declarations, and the call to SPOTRF.
8 C**********************************************************************
9 C
10 C SUBROUTINE SETGMN( MEANV, COVM, LDCOVM, P, PARM)
11 C SET Generate Multivariate Normal random deviate
12 C
13 C
14 C Function
15 C
16 C
17 C Places P, MEANV, and the Cholesky factoriztion of COVM
18 C in PARM for GENMN.
19 C
20 C
21 C Arguments
22 C
23 C
24 C MEANV --> Mean vector of multivariate normal distribution.
25 C REAL MEANV(P)
26 C
27 C COVM <--> (Input) Covariance matrix of the multivariate
28 C normal distribution. This routine uses only the
29 C (1:P,1:P) slice of COVM, but needs to know LDCOVM.
30 C
31 C (Output) Destroyed on output
32 C REAL COVM(LDCOVM,P)
33 C
34 C LDCOVM --> Leading actual dimension of COVM.
35 C INTEGER LDCOVM
36 C
37 C P --> Dimension of the normal, or length of MEANV.
38 C INTEGER P
39 C
40 C PARM <-- Array of parameters needed to generate multivariate
41 C normal deviates (P, MEANV and Cholesky decomposition
42 C of COVM).
43 C 1 : 1 - P
44 C 2 : P + 1 - MEANV
45 C P+2 : P*(P+3)/2 + 1 - Cholesky decomposition of COVM
46 C REAL PARM(P*(P+3)/2 + 1)
47 C
48 C**********************************************************************
49 C .. Scalar Arguments ..
50 C INTEGER p
51  INTEGER p, ldcovm
52 C ..
53 C .. Array Arguments ..
54 C REAL covm(p,p),meanv(p),parm(p* (p+3)/2+1)
55  REAL covm(ldcovm,p),meanv(p),parm(p* (p+3)/2+1)
56 C ..
57 C .. Local Scalars ..
58  INTEGER i,icount,info,j
59 C ..
60 C .. External Subroutines ..
61  EXTERNAL spotrf
62 C ..
63 C .. Executable Statements ..
64 C
65 C
66 C TEST THE INPUT
67 C
68  IF (.NOT. (p.LE.0)) go to 10
69  WRITE (*,*) 'P nonpositive in SETGMN'
70  WRITE (*,*) 'Value of P: ',p
71  CALL xstopx('P nonpositive in SETGMN')
72
73  10 parm(1) = p
74 C
75 C PUT P AND MEANV INTO PARM
76 C
77  DO 20,i = 2,p + 1
78  parm(i) = meanv(i-1)
79  20 CONTINUE
80 C
81 C Cholesky decomposition to find A s.t. trans(A)*(A) = COVM
82 C
83 C CALL spofa(covm,p,p,info)
84 C CALL spofa(covm,ldcovm,p,info)
85  CALL spotrf( 'Upper', p, covm, ldcovm, info)
86  IF (.NOT. (info.NE.0)) go to 30
87  WRITE (*,*) ' COVM not positive definite in SETGMN'
88  CALL xstopx(' COVM not positive definite in SETGMN')
89
90  30 icount = p + 1
91 C
92 C PUT UPPER HALF OF A, WHICH IS NOW THE CHOLESKY FACTOR, INTO PARM
93 C COVM(1,1) = PARM(P+2)
94 C COVM(1,2) = PARM(P+3)
95 C :
96 C COVM(1,P) = PARM(2P+1)
97 C COVM(2,2) = PARM(2P+2) ...
98 C
99  DO 50,i = 1,p
100  DO 40,j = i,p
101  icount = icount + 1
102  parm(icount) = covm(i,j)
103  40 CONTINUE
104  50 CONTINUE
105  RETURN
106 C
107  END
subroutine setgmn(meanv, covm, ldcovm, p, parm)
Definition: setgmn.f:1