00001 SUBROUTINE setgmn(meanv,covm,ldcovm,p,parm) 00002 C SUBROUTINE setgmn(meanv,covm,p,parm) 00003 C JJV changed this routine to take leading dimension of COVM 00004 C JJV argument and pass it to SPOTRF, making it easier to use 00005 C JJV if the COVM which is used is contained in a larger matrix 00006 C JJV and to make the routine more consistent with LAPACK. 00007 C JJV Changes are in comments, declarations, and the call to SPOTRF. 00008 C********************************************************************** 00009 C 00010 C SUBROUTINE SETGMN( MEANV, COVM, LDCOVM, P, PARM) 00011 C SET Generate Multivariate Normal random deviate 00012 C 00013 C 00014 C Function 00015 C 00016 C 00017 C Places P, MEANV, and the Cholesky factoriztion of COVM 00018 C in PARM for GENMN. 00019 C 00020 C 00021 C Arguments 00022 C 00023 C 00024 C MEANV --> Mean vector of multivariate normal distribution. 00025 C REAL MEANV(P) 00026 C 00027 C COVM <--> (Input) Covariance matrix of the multivariate 00028 C normal distribution. This routine uses only the 00029 C (1:P,1:P) slice of COVM, but needs to know LDCOVM. 00030 C 00031 C (Output) Destroyed on output 00032 C REAL COVM(LDCOVM,P) 00033 C 00034 C LDCOVM --> Leading actual dimension of COVM. 00035 C INTEGER LDCOVM 00036 C 00037 C P --> Dimension of the normal, or length of MEANV. 00038 C INTEGER P 00039 C 00040 C PARM <-- Array of parameters needed to generate multivariate 00041 C normal deviates (P, MEANV and Cholesky decomposition 00042 C of COVM). 00043 C 1 : 1 - P 00044 C 2 : P + 1 - MEANV 00045 C P+2 : P*(P+3)/2 + 1 - Cholesky decomposition of COVM 00046 C REAL PARM(P*(P+3)/2 + 1) 00047 C 00048 C********************************************************************** 00049 C .. Scalar Arguments .. 00050 C INTEGER p 00051 INTEGER p, ldcovm 00052 C .. 00053 C .. Array Arguments .. 00054 C REAL covm(p,p),meanv(p),parm(p* (p+3)/2+1) 00055 REAL covm(ldcovm,p),meanv(p),parm(p* (p+3)/2+1) 00056 C .. 00057 C .. Local Scalars .. 00058 INTEGER i,icount,info,j 00059 C .. 00060 C .. External Subroutines .. 00061 EXTERNAL spotrf 00062 C .. 00063 C .. Executable Statements .. 00064 C 00065 C 00066 C TEST THE INPUT 00067 C 00068 IF (.NOT. (p.LE.0)) GO TO 10 00069 WRITE (*,*) 'P nonpositive in SETGMN' 00070 WRITE (*,*) 'Value of P: ',p 00071 CALL XSTOPX ('P nonpositive in SETGMN') 00072 00073 10 parm(1) = p 00074 C 00075 C PUT P AND MEANV INTO PARM 00076 C 00077 DO 20,i = 2,p + 1 00078 parm(i) = meanv(i-1) 00079 20 CONTINUE 00080 C 00081 C Cholesky decomposition to find A s.t. trans(A)*(A) = COVM 00082 C 00083 C CALL spofa(covm,p,p,info) 00084 C CALL spofa(covm,ldcovm,p,info) 00085 CALL spotrf ( 'Upper', p, covm, ldcovm, info) 00086 IF (.NOT. (info.NE.0)) GO TO 30 00087 WRITE (*,*) ' COVM not positive definite in SETGMN' 00088 CALL XSTOPX (' COVM not positive definite in SETGMN') 00089 00090 30 icount = p + 1 00091 C 00092 C PUT UPPER HALF OF A, WHICH IS NOW THE CHOLESKY FACTOR, INTO PARM 00093 C COVM(1,1) = PARM(P+2) 00094 C COVM(1,2) = PARM(P+3) 00095 C : 00096 C COVM(1,P) = PARM(2P+1) 00097 C COVM(2,2) = PARM(2P+2) ... 00098 C 00099 DO 50,i = 1,p 00100 DO 40,j = i,p 00101 icount = icount + 1 00102 parm(icount) = covm(i,j) 00103 40 CONTINUE 00104 50 CONTINUE 00105 RETURN 00106 C 00107 END
1.7.1
