GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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