### Blame view

m_cholesky.F90 1.55 KB
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62``` ``````module m_cholesky contains subroutine cholesky(Aa,nn) integer, intent(in) :: nn real, intent(inout), dimension(1:nn,1:nn) :: Aa integer :: i, j, k real :: sum=0 ! print*, '<><><><><><><> Cholesky is working for you<><><><><><><>' do i=1, nn do j=1, i-1 sum=0 do k=1, j-1 sum=sum + Aa(i,k)*Aa(j,k) end do Aa(i,j)=(1/Aa(j,j))*(Aa(i,j)-sum) end do sum=0 do k=1, i-1 sum=sum+(Aa(i,k))**2 end do Aa(i,i)=sqrt(Aa(i,i)-sum) do j=i+1, nn ! Lower triangle sqrt matrix Aa(i,j)=0 end do end do ! print*, '<><><><><><><> Cholesky made good job!! <><><><><><><>' return end subroutine cholesky ! Exponential covariance function function exponential(h) real, intent(in) :: h real exponential if(h<0) stop 'm_cholesky exponential covariance: negative distance' exponential=exp(-h) end function exponential ! Gaussian covariance function function gaussian(h) real, intent(in) :: h real gaussian if(h<0) stop 'm_cholesky gaussian covariance: negative distance' gaussian=0.999*exp(-h*h) if (h.eq.0.0) gaussian=1.0 end function gaussian ! spherical covariance function function spherical(h) real, intent(in) :: h real spherical if(h<0) stop 'm_cholesky spherical covariance: negative distance' if(h>1.) then spherical=0. else spherical=1.-1.5*h+0.5*h*h*h endif end function spherical end module m_cholesky ``````