Guest User

Untitled

a guest
Apr 25th, 2018
74
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.03 KB | None | 0 0
  1. MODULE iteration
  2. implicit none
  3. contains
  4.  
  5. SUBROUTINE iterate(alpha, beta, e, es, rank, omega, smearing, prec, max_step)
  6. REAL(kind=8), INTENT(in) :: omega, smearing, prec
  7. INTEGER :: max_step, step, rank, cnt
  8. COMPLEX(kind=16) :: alpha(rank,rank), beta(rank,rank), omega_mat(rank, rank), green(rank, rank)
  9. COMPLEX(kind=16), INTENT(inout) :: e(rank,rank), es(rank,rank)
  10. step = 0
  11. omega_mat = 0
  12. DO cnt=1, rank
  13. omega_mat(cnt, cnt) = 1.0
  14. ENDDO
  15. omega_mat = omega_mat * (omega + (0.0, 1.0) * smearing)
  16. DO WHILE (minval(abs(alpha)) .gt. prec .or. minval(abs(beta)) .gt. prec .and. step .lt. max_step)
  17. green = zInverse(rank, omega_mat - e)
  18. e = e + matmul(alpha, matmul(green, beta)) + matmul(beta, matmul(green, alpha))
  19. es = es + matmul(alpha, matmul(green, beta))
  20. alpha = matmul(alpha, matmul(green, alpha))
  21. beta = matmul(beta, matmul(green, beta))
  22. step = step + 1
  23. ENDDO
  24. END SUBROUTINE iterate
  25.  
  26. FUNCTION zInverse(n, a) result(ra)
  27. INTEGER :: n,lda,ipiv(n),info,lwork
  28. COMPLEX(kind=16)::a(n,n),ra(n,n),work(n)
  29. ra=a
  30. lwork=n
  31. lda=n
  32. CALL zgetrf(n, n, ra, lda, ipiv, info)
  33. IF(info/=0) WRITE(0,*) 'Error occured in zgetrf!'
  34. CALL zgetri(n, ra, lda, ipiv, work, lwork, info)
  35. IF(info/=0) WRITE(0,*) 'Error occured in zgetri!'
  36. END FUNCTION zInverse
  37. END MODULE iteration
  38.  
  39. import numpy as np
  40. import pyiteration
  41. alpha = np.array([[1,0],[0,1]], dtype='complex')
  42. beta = np.array([[1,0],[0,1]], dtype='complex')
  43. e = np.array([[1,0],[0,1]], dtype='complex')
  44. es = np.array([[1,0],[0,1]], dtype='complex')
  45. # f2py is automatically generating rank for me
  46. pyiteration.iteration.iterate(alpha,beta, e, es, 1.0, 0.001, 0.001, 100)
  47.  
  48. alpha = np.array([[1,0],[0,1]], dtype='complex', order='F')
  49. beta = np.array([[1,0],[0,1]], dtype='complex', order='F')
  50. e = np.array([[1,0],[0,1]], dtype='complex', order='F')
  51. es = np.array([[1,0],[0,1]], dtype='complex', order='F')
  52.  
  53. real(kind=8) dimension(:,:),intent(in) :: kernel
  54. integer dimension(:,:),intent(in) :: kernelmask
Add Comment
Please, Sign In to add comment