Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- MODULE iteration
- implicit none
- contains
- SUBROUTINE iterate(alpha, beta, e, es, rank, omega, smearing, prec, max_step)
- REAL(kind=8), INTENT(in) :: omega, smearing, prec
- INTEGER :: max_step, step, rank, cnt
- COMPLEX(kind=16) :: alpha(rank,rank), beta(rank,rank), omega_mat(rank, rank), green(rank, rank)
- COMPLEX(kind=16), INTENT(inout) :: e(rank,rank), es(rank,rank)
- step = 0
- omega_mat = 0
- DO cnt=1, rank
- omega_mat(cnt, cnt) = 1.0
- ENDDO
- omega_mat = omega_mat * (omega + (0.0, 1.0) * smearing)
- DO WHILE (minval(abs(alpha)) .gt. prec .or. minval(abs(beta)) .gt. prec .and. step .lt. max_step)
- green = zInverse(rank, omega_mat - e)
- e = e + matmul(alpha, matmul(green, beta)) + matmul(beta, matmul(green, alpha))
- es = es + matmul(alpha, matmul(green, beta))
- alpha = matmul(alpha, matmul(green, alpha))
- beta = matmul(beta, matmul(green, beta))
- step = step + 1
- ENDDO
- END SUBROUTINE iterate
- FUNCTION zInverse(n, a) result(ra)
- INTEGER :: n,lda,ipiv(n),info,lwork
- COMPLEX(kind=16)::a(n,n),ra(n,n),work(n)
- ra=a
- lwork=n
- lda=n
- CALL zgetrf(n, n, ra, lda, ipiv, info)
- IF(info/=0) WRITE(0,*) 'Error occured in zgetrf!'
- CALL zgetri(n, ra, lda, ipiv, work, lwork, info)
- IF(info/=0) WRITE(0,*) 'Error occured in zgetri!'
- END FUNCTION zInverse
- END MODULE iteration
- import numpy as np
- import pyiteration
- alpha = np.array([[1,0],[0,1]], dtype='complex')
- beta = np.array([[1,0],[0,1]], dtype='complex')
- e = np.array([[1,0],[0,1]], dtype='complex')
- es = np.array([[1,0],[0,1]], dtype='complex')
- # f2py is automatically generating rank for me
- pyiteration.iteration.iterate(alpha,beta, e, es, 1.0, 0.001, 0.001, 100)
- alpha = np.array([[1,0],[0,1]], dtype='complex', order='F')
- beta = np.array([[1,0],[0,1]], dtype='complex', order='F')
- e = np.array([[1,0],[0,1]], dtype='complex', order='F')
- es = np.array([[1,0],[0,1]], dtype='complex', order='F')
- real(kind=8) dimension(:,:),intent(in) :: kernel
- integer dimension(:,:),intent(in) :: kernelmask
Add Comment
Please, Sign In to add comment