m_pseudo.F90 8.01 KB
module m_pseudo
contains
  subroutine pseudo(A,n,nrfields,lowersqrt)
    use mod_dimensions
    use m_random
    implicit none
    real, intent(out)   :: A(n,nrfields) 
    integer, intent(in) :: nrfields,n
    real, intent(in)    :: lowersqrt(n,n) ! cholesky square root
    
    real nugget(n)
    integer l,i
    do l=1,nrfields
       call random(nugget,n)
       do i=1,n
          A(i,l)=dot_product(lowersqrt(:,i),nugget)  ! matmul(transpose(lowersqrt),nugget)
       end do
    end do
    
  end subroutine pseudo
end module m_pseudo
!
!   real, save :: sigg    ! sigg unchanged by subroutines if its name is repeated
!
!
!   integer l,p,j,n,m,i
!   real sigg2,c
!   real kappa2,lambda2,kappa,lambda
!   real pi2,deltak,sum
!   real aa,bb,tol
!
!   real fampl(0:nfft/2,2)
!   real phi(0:nfft/2)
!
!#ifdef SGI
!   integer sign
!   real(kind=8), save :: coeff(15+4*nfft) ! 19+2*nfft for imsl routines
!   real cpy(2*nfft) ! real *8 only for imsl df2tcb routine
!!!   real array(0:nfft+1)
!!   real array(2,nfft/2+1)
!!   complex(kind=8) arrayC(0:nfft-1)  ! nfft/2 for scomplib
!   complex arrayC(0:nfft-1)  ! nfft/2 for scomplib
!!   complex arrayC(0:nfft-1)
!   real tmp(0:nfft-1)
!#endif
!   real tt
!   integer, save :: iprepfft=0
!   real,    save :: rh_save=0.0
!   logical, save :: diag=.true.
!
!   real dx,rh,fval
!
!   integer ii
!   character(len=3) tag
!
!!!!!!!!!!!!!!!!!!
!
!   dx=1.0 ; rh=rhh/dxx
!
!   pi2=2.0*pi
!   deltak=pi2/(float(nfft)*dx)
!   kappa=pi2/(float(nfft)*dx)
!   kappa2=kappa**2
!
!   if (iprepfft == 0) then
!#ifdef SGI
!!      call dzfft1dui(nfft,coeff)
!      call fftci(nfft,coeff) !for imsl only nfft/2+1 ??
!!      print '(a,99f9.4)','ini coeff : ',coeff
!#endif
!      iprepfft=1
!   endif
!
!   if (rhh /= rh_save) then
!      rh_save=rhh
!      aa=0.1e-01
!      bb=0.1e+01
!      tol=0.1e-10
!      call zeroin(sigg,aa,bb,tol,rh,dx,fval)
!      print *,'zeroin: sigma and fval: ',sigg,fval
!      diag=.true.
!   endif
!
!   sigg2=sigg**2
!   sum=0.0
!   do l=-nfft/2+1,nfft/2
!      sum=sum+exp(-2.0*(kappa2*float(l*l))/sigg2)
!   enddo
!   c=sqrt(1.0/(deltak*sum))
!
!   if (diag) then
!      print *,'pseudo: sum ',sum
!      print *,'pseudo: rh ',rh
!      print *,'pseudo: sigg  ',sigg
!      print *,'pseudo: c=     ',c
!      print *,'pseudo: sqrt(deltak)*c ',sqrt(deltak)*c
!      diag=.false.
!   endif
!
!   do i=1,nrfields
!!     Calculating the random wave phases 
!      call random_number(phi)
!! NOW:
!      write(tag,'(i3.3)')i
!      open(23,file='random'//tag//'.dat')
!      do ii=0,nfft/2
!        write(23,*)ii,phi(ii)
!      enddo
!      close(23)
!
!      phi=pi2*phi
!
!!     Calculating the wave amplitues 
!      do l=0,nfft/2 
!         tt=kappa2*real(l*l)/sigg2
!         fampl(l,1)=exp(-tt)*cos(phi(l))*sqrt(deltak)*c
!         fampl(l,2)=exp(-tt)*sin(phi(l))*sqrt(deltak)*c
!      enddo
!      fampl(0,2)=0.0
!
!! Here are two parallel ways of computing the  inverse fft
!#ifdef FFT_CASEA
!!  1'st case
!!  real arrayC(2,nfft/2+1)
!      stop 'error : case a not changed from scomplib to imsl'
!      print *,'SJALABAIS I 1'st!'
!      array(1,:)=fampl(:,1)
!      array(2,:)=fampl(:,2)
!      print *,'ARRAY intermediate (real)'
!      print '(10e12.4)',array
!      sign=-1 
!      call zdfft1du(sign,nfft,array,1,coeff)
!      print *,'ARRAY out'
!      print '(10e12.4)',array
!      print '(10e12.4)',arrayC
!#endif
!
!!  2'st case is easier when copying the data back to A
!!   complex arrayC(0:nfft/2)
!      arrayC(0:nfft/2)=cmplx(fampl(0:nfft/2,1),fampl(0:nfft/2,2))
!      do l=1,nfft/2-1
!         arrayC(nfft-l)=conjg(arrayC(l)) ! for real values sequences
!      end do
!      arrayC(nfft/2)=cmplx(real(arrayC(nfft/2)),0.0)
!      sign=-1 
!!      print '(a,99f5.2)', 'fourierC :', arrayC
!!      call zdfft1du(sign,nfft,arrayC,1,coeff)
!      call f2tcb(nfft,arrayC,arrayC,coeff,cpy)  ! for imsl only (nfft/2+1)
!      tmp=real(arrayC) 
!!      print '(a,99f5.2)', 'arrayC :', arrayC
!!      print '(a,5f5.2)','tmp=', tmp(0:5)
!      A(1:n,i)=tmp(0:n-1)
!
!#ifdef DEBUG_FFT
!      stop 'error : debug_fft not changed to imsl'
!      array=0.0
!      call random_number(array(0:nfft-1))
!      print *,'ARRAY in'
!      print '(10e12.4)',array
!      sign=1 
!      call dzfft1dui(nfft,coeff)
!      call dzfft1du(sign,nfft,array,1,coeff)
!      print *,'ARRAY intermediate'
!      print '(10e12.4)',array
!      sign=-1 
!      call zdfft1du(sign,nfft,array,1,coeff)
!      call dscal1d(nfft,1.0/real(nfft),array,1)
!      print *,'ARRAY out'
!      print '(10e12.4)',array
!#endif
!
!   enddo
!
!!   do i=1,n
!!      write(21,'(i4,10e12.4)')i,A(i,1:min(nrfields,10))
!!   enddo
!
!end subroutine
!
!subroutine zeroin(zeropkt,ax,bx,tol,length,dx,fval)
!! Finds zero of function f.
!! A zero of the function  $func(x,length,dx)$
!! is computed in the interval $[ax,bx]$.
!! Zeroin| returns a zero $x$ in the given interval
!! to within a tolerance  $4*macheps*abs(x) + tol$, where macheps
!! is the relative machine precision.
!
!! This function subprogram is a slightly  modified  translation  of
!! the algol 60 procedure  zero  given in  richard brent, algorithms for
!! minimization without derivatives, prentice - hall, inc. (1973).
!
!   real zeropkt,length,dx
!   real ax   ! left endpoint of initial interval
!   real bx   ! right endpoint of initial interval
!   real tol  !  desired length of the interval of uncertainty of the
!   real  a,b,c,d,e,eps,fa,fb,fc,tol1,xm,p,q,r,s
!   real  abs,sign,fval
!
!!  compute eps, the relative machine precision
!
!#ifdef DEBUG
!   write(*,*)'A=',ax,bx,tol,length,dx
!#endif
!   icorr=0
!
!   eps = 1.0
! 10 eps = eps/2.0
!   tol1 = 1.0 + eps
!   if (tol1 .gt. 1.0) go to 10
!
!! initialization
! 77 a = ax
!   b = bx
!   fa = func(a,length,dx)
!   fb = func(b,length,dx)
!
!
!   if (fa*fb.gt.0.0) then
!      write(*,*)'fa=',fa
!      write(*,*)'fb=',fb
!      write(*,*)'fa*fb =',fa*fb,'is greater than zero'
!      ax=0.1*ax
!      bx=10.0*bx
!      icorr=icorr+1
!      if (icorr < 20) then
!         goto 77
!      else
!         write(*,'(2(a,g13.5))')'zeroin: No convergence, ax=',ax,' bx=',bx
!         stop
!      endif
!   endif
!
!! begin step
!
! 20 c = a
!   fc = fa
!   d = b - a
!   e = d
! 30 if (abs(fc) .ge. abs(fb)) go to 40
!   a = b
!   b = c
!   c = a
!   fa = fb
!   fb = fc
!   fc = fa
!
!! convergence test
!
! 40 tol1 = 2.0*eps*abs(b) + 0.5*tol
!   xm = .5*(c - b)
!   if (abs(xm) .le. tol1) go to 90
!   if (fb .eq. 0.0) go to 90
!
!! is bisection necessary
!
!   if (abs(e) .lt. tol1) go to 70
!   if (abs(fa) .le. abs(fb)) go to 70
!
!! is quadratic interpolation possible
!
!   if (a .ne. c) go to 50
!
!! linear interpolation
!
!   s = fb/fa
!   p = 2.0*xm*s
!   q = 1.0 - s
!   go to 60
!
!! inverse quadratic interpolation
!
! 50 q = fa/fc
!   r = fb/fc
!   s = fb/fa
!   p = s*(2.0*xm*q*(q - r) - (b - a)*(r - 1.0))
!   q = (q - 1.0)*(r - 1.0)*(s - 1.0)
!
!! adjust signs
!
! 60 if (p .gt. 0.0) q = -q
!   p = abs(p)
!
!! is interpolation acceptable
!
!   if ((2.0*p) .ge. (3.0*xm*q - abs(tol1*q))) go to 70
!   if (p .ge. abs(0.5*e*q)) go to 70
!   e = d
!   d = p/q
!   go to 80
!
!! bisection
!
! 70 d = xm
!   e = d
!
!! complete step
!
! 80 a = b
!   fa = fb
!   if (abs(d) .gt. tol1) b = b + d
!   if (abs(d) .le. tol1) b = b + sign(tol1, xm)
!   fb = func(b,length,dx)
!   if ((fb*(fc/abs(fc))) .gt. 0.0) go to 20
!   go to 30
!
!! done
!
! 90 zeropkt = b
!   fval=func(b,length,dx) 
!
!
!   contains
!
!      real function func(sigg,length,dx)
!      use mod_dimensions
!      implicit none
!      real sum1,sum2,sigg,length
!      real sigg2,pi2,kappa,kappa2,lambda,lambda2,dx
!      integer l,p
!
!      sigg2=sigg**2
!
!      pi2=2.0*pi
!      kappa=pi2/(float(2*kdim)*dx)
!      kappa2=kappa**2
!
!!      /* Calculate sum1 */
!      sum1=0.0
!      do l=-nfft/2+1,nfft/2
!         sum1=sum1+&
!         exp(-2.0*(kappa2*float(l*l))/sigg2)*cos(kappa*float(l)*length)
!      enddo
!
!!      /* Calculate sum2 */
!      sum2=0.0
!      do l=-nfft/2+1,nfft/2
!         sum2=sum2+exp(-2.0*(kappa2*float(l*l))/sigg2)
!      enddo
!
!      func = sum1/sum2 - exp(-1.0)
!
!      end function
!end subroutine
!end module