ksampPw.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	ksamppw (ini, coef, power, node, n, x, fx)
subroutine	ksamppw2 (coef, power, xp, cum, intgv, np, x, fx)
subroutine	ksamppwx (c, p, x0, n, xp, c2)
subroutine	ksamppwmrg (c1, p1, x1, n1, c2, p2, x2, n2, c, p, x, n)

Function/Subroutine Documentation

ksamppw()

subroutine ksamppw	(	integer	ini,
		real*8, dimension(n)	coef,
		real*8, dimension(n)	power,
		real*8, dimension(n+1)	node,
		integer	n,
		real*8	x,
		real*8	fx
	)

Definition at line 95 of file ksampPw.f.

References d0, parameter(), and rndc().

Referenced by cmbreme1(), cmbreme2(), and main().

      implicit none
!
      integer ini  ! input.  if coef, power, node, n are
                   !     different from the previous call, give 1.
                   !     If all of them are the same as the
                   !     previous call, give non 1
      integer n    ! input.  number of segments
      real*8 coef(n)  ! input.  c of cX**(-p)dX in each segment.
      real*8 power(n) ! input.  p of cX**(-p)dX in each segment.
      real*8 node(n+1)  ! input.  node(1) and node(2) are the left and
                      !         right boundaries of segment 1.,
                      !         node(2) and  node(3) are for segment 2,etc
      real*8 x        ! output. A sampled x. If the input condition 
                      !        is ng, x = 0 will result.
      real*8 fx       ! output.  function value  at  x, i.e., cx**-p
! 
      integer nmx     !  max number of nodes
      parameter( nmx = 20 )

      real*8 inte(nmx)  ! integral value of each segment region
      real*8 ci(nmx)    ! normalized cummulative integral values of inte. 
      integer nsave     ! for saving n

      save inte,  ci, nsave
      real*8 u, sum, temp
      integer i, j

      real*8  cum(*), intgv(*)
      integer np

!
      if(ini .eq. 1) then
         if(nmx .lt. n) then
            x = 0.
         else
!             get integral value of each region and sum of them
            sum = 0.
!             c/(1-p) x**(1-p) from x1 to x2
            do i = 1, n
               if(power(i) .ne. 1.d0) then
                 inte(i) = coef(i)/(1.0-power(i))*
     *           (node(i+1)**(1.-power(i)) - node(i)**(1.-power(i)))
               else
                  inte(i) = coef(i) * log(node(i+1)/node(i))
               endif
               sum = sum + inte(i)
            enddo
!             normalize and make cummulative data
            do i = 1, n
               ci(i) = inte(i)/sum
            enddo
            do i = 2, n
               ci(i) = ci(i) + ci(i-1)
            enddo
            ci(n) = 1.0   ! for safety
         endif
         nsave = n
      endif
!     --------------------------------------------------
!           assume inte and node are ready
!          choose which segment, j
      call rndc(u)
      do i = 1, n
         if(u .le. ci(i)) then
            j = i
            goto 100
         endif
      enddo
!       never come here
!
 100  continue
      

      call rndc(u)
      if(power(j) .ne. 1.d0) then
         temp = node(j+1)**(1.-power(j))  - 
     *       inte(j)*u/( coef(j)/(1. - power(j)) )       ! = x **(1-power(j))

         x = temp**(1./(1.-power(j)))
      else
         x = (node(j+1)/node(j))**u  * node(j)
      endif
      fx = coef(j) * x **(-power(j))
 1000 continue
      return
!  
!
!     ********************************
      entry ksamppwq(cum, intgv, np)
!     *******************************
!      inquire the current inner variables for the 
!      present spectrum.  
!

      np = nsave

      do i = 1, np
         cum(i) = ci(i)
         intgv(i) = inte(i) 
      enddo

Here is the call graph for this function:

Here is the caller graph for this function:

ksamppw2()

subroutine ksamppw2	(	real*8, dimension(np)	coef,
		real*8, dimension(np)	power,
		real*8, dimension(np+1)	xp,
		real*8, dimension(np)	cum,
		real*8, dimension(np)	intgv,
		integer	np,
		real*8	x,
		real*8	fx
	)

Definition at line 199 of file ksampPw.f.

References rndc().

!
!      This interface is used to skip the initial calculations
!      by using output from ksampPwQ
!
      implicit none
!

      integer np      ! input.  number of segments
      real*8 coef(np)  ! input. c of X**(-p)dX in each segment.
      real*8 power(np) ! input.  p of X**(-p)dX in each segment.
      real*8 xp(np+1)  ! input.  output from ksampPwQ
      real*8 cum(np)   ! input.  //
      real*8 intgv(np) ! //
    
      real*8 x        ! output. sampled x. 
      real*8 fx       ! output. funcion value at x. 

      real*8  u, temp
      integer i, j
!          choose which segment, j
      call rndc(u)
      do i = 1, np
         if(u .le. cum(i)) then
            j = i
            goto 100
         endif
      enddo
!
 100  continue

      call rndc(u)
      temp = xp(j+1)**(1.-power(j))  - 
     *       intgv(j)*u/( coef(j)/(1. - power(j)) )       ! = x **(1-power(j))

      x = temp**(1./(1.-power(j)))
      fx = coef(j) * x**(-power(j))

Here is the call graph for this function:

ksamppwmrg()

subroutine ksamppwmrg	(	real*8, dimension(n1)	c1,
		real*8, dimension(n1)	p1,
		real*8, dimension(n1+1)	x1,
		integer	n1,
		real*8, dimension(n2)	c2,
		real*8, dimension(n2)	p2,
		real*8, dimension(n2+1)	x2,
		integer	n2,
		real*8, dimension(*)	c,
		real*8, dimension(*)	p,
		real*8, dimension(*)	x,
		integer	n
	)

Definition at line 273 of file ksampPw.f.

Referenced by cmbreme1(), and cmbreme2().

!
      implicit none
      integer n1   ! input.         number of segments of 1 st function
                   !                c1* x**(-p1) dx 
      real*8  c1(n1)  ! input.     
      real*8  p1(n1)  ! input.
      real*8  x1(n1+1)  ! input    ! nodal point of segments of 1 st func.
      integer n2   ! input.        ! the same as the 2nd func.
      real*8  c2(n2)  ! input.   
      real*8  p2(n2)  ! input.
      real*8  x2(n2+1)  ! input
      integer n   !  output        !  those for f1 x f2
      real*8  c(*)  ! output.
      real*8  p(*)  ! output
      real*8  x(*)  ! output   (dim >= nx+1 )

!

      integer i, j
!
      i = 1
      j= 1
      n = 0

      x(1) = max(x1(1), x2(1))
      do while (i .le. n1 .and. j .le. n2)
         if(x1(i+1) .lt. x2(j+1) ) then
            n = n + 1
            x(n+1) = x1(i+1)
            c(n) = c1(i) * c2(j)
            p(n)  = p1(i) + p2(j)
            i = i + 1

         elseif(x1(i+1) .eq. x2(j+1)) then
            n = n + 1
            x(n+1) = x1(i+1)
            c(n) = c1(i) * c2(j)
            p(n)  = p1(i) + p2(j)
            i = i + 1
            j = j + 1
         else
            n = n + 1
            x(n+1) = x2(j+1)
            c(n) = c1(i) * c2(j)
            p(n)  = p1(i) + p2(j)
            j = j + 1
         endif
      enddo

Here is the caller graph for this function:

ksamppwx()

subroutine ksamppwx	(	real*8, dimension(n)	c,
		real*8, dimension(n)	p,
		real*8, dimension(n)	x0,
		integer	n,
		real*8, dimension(n-1)	xp,
		real*8, dimension(n)	c2
	)

Definition at line 240 of file ksampPw.f.

Referenced by cmbreme1(), cmbreme2(), and main().

      implicit none
!
      integer n    ! input.  number of segments
      real*8 c(n)  ! input.  c of c(X/x0)**(-p) in each segment.
      real*8 p(n) ! input.    see above
      real*8 x0(n)  ! input.  see above
      real*8 xp(n-1)  ! output.  crossing point of adjacent segments
      real*8 c2(n)  ! output. c2 = c/x0**p. That is, c2*x**(-p) becomes
                    !           the spectrum
! 
      integer i

!      get coeff.      
      do i = 1, n
         c2(i) =  c(i)* x0(i)**p(i)
      enddo
!     get crossing point of i and i+1th segmment
!              c_i X**(-p_i) = c_i+1 X**(-p_i+1)

      do i = 1, n-1
          xp(i) =
     *    (c2(i)/c2(i+1) )**(1./(p(i) - p(i+1)))
      enddo
      

Here is the caller graph for this function: