home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / ranlib / setgmn.f < prev    next >
Text File  |  1996-09-28  |  3KB  |  93 lines
  1.       SUBROUTINE setgmn(meanv,covm,p,parm)
  2. C**********************************************************************
  3. C
  4. C     SUBROUTINE SETGMN( MEANV, COVM, P, PARM)
  5. C            SET Generate Multivariate Normal random deviate
  6. C
  7. C
  8. C                              Function
  9. C
  10. C
  11. C      Places P, MEANV, and the Cholesky factoriztion of COVM
  12. C      in GENMN.
  13. C
  14. C
  15. C                              Arguments
  16. C
  17. C
  18. C     MEANV --> Mean vector of multivariate normal distribution.
  19. C                                        REAL MEANV(P)
  20. C
  21. C     COVM   <--> (Input) Covariance   matrix    of  the  multivariate
  22. C                 normal distribution
  23. C                 (Output) Destroyed on output
  24. C                                        REAL COVM(P,P)
  25. C
  26. C     P     --> Dimension of the normal, or length of MEANV.
  27. C                                        INTEGER P
  28. C
  29. C     PARM <-- Array of parameters needed to generate multivariate norma
  30. C                deviates (P, MEANV and Cholesky decomposition of
  31. C                COVM).
  32. C                1 : 1                - P
  33. C                2 : P + 1            - MEANV
  34. C                P+2 : P*(P+3)/2 + 1  - Cholesky decomposition of COVM
  35. C                                             REAL PARM(P*(P+3)/2 + 1)
  36. C
  37. C**********************************************************************
  38. C     .. Scalar Arguments ..
  39.       INTEGER p
  40. C     ..
  41. C     .. Array Arguments ..
  42.       REAL covm(p,p),meanv(p),parm(p* (p+3)/2+1)
  43. C     ..
  44. C     .. Local Scalars ..
  45.       INTEGER i,icount,info,j
  46. C     ..
  47. C     .. External Subroutines ..
  48.       EXTERNAL spofa
  49. C     ..
  50. C     .. Executable Statements ..
  51. C
  52. C
  53. C     TEST THE INPUT
  54. C
  55.       IF (.NOT. (p.LE.0)) GO TO 10
  56.       WRITE (*,*) 'P nonpositive in SETGMN'
  57.       WRITE (*,*) 'Value of P: ',p
  58.       CALL XSTOPX ('P nonpositive in SETGMN')
  59.  
  60.    10 parm(1) = p
  61. C
  62. C     PUT P AND MEANV INTO PARM
  63. C
  64.       DO 20,i = 2,p + 1
  65.           parm(i) = meanv(i-1)
  66.    20 CONTINUE
  67. C
  68. C      Cholesky decomposition to find A s.t. trans(A)*(A) = COVM
  69. C
  70.       CALL spofa(covm,p,p,info)
  71.       IF (.NOT. (info.NE.0)) GO TO 30
  72.       WRITE (*,*) ' COVM not positive definite in SETGMN'
  73.       CALL XSTOPX (' COVM not positive definite in SETGMN')
  74.  
  75.    30 icount = p + 1
  76. C
  77. C     PUT UPPER HALF OF A, WHICH IS NOW THE CHOLESKY FACTOR, INTO PARM
  78. C          COVM(1,1) = PARM(P+2)
  79. C          COVM(1,2) = PARM(P+3)
  80. C                    :
  81. C          COVM(1,P) = PARM(2P+1)
  82. C          COVM(2,2) = PARM(2P+2)  ...
  83. C
  84.       DO 50,i = 1,p
  85.           DO 40,j = i,p
  86.               icount = icount + 1
  87.               parm(icount) = covm(i,j)
  88.    40     CONTINUE
  89.    50 CONTINUE
  90.       RETURN
  91. C
  92.       END
  93.