Advertisement
Guest User

Untitled

a guest
Sep 22nd, 2014
191
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 3.10 KB | None | 0 0
  1. complex(R_P), dimension(:,:), intent(in) :: mat
  2. complex(R_P), dimension(:), intent(inout) :: gs_vec
  3. real(R_P), intent(out) :: gs_val
  4.  
  5. integer, parameter :: nev = 1, ncv = 10
  6. integer :: iparam(11), ipntr(14)
  7. logical :: select(ncv)
  8. complex(R_P) :: d(ncv), v(size(mat,1),ncv), workd(3*size(mat,1)), workev(3*ncv), resid(size(mat,1))
  9. complex(R_P) :: workl(3*ncv*ncv+5*ncv)
  10. real(R_P) :: rwork(ncv)
  11.  
  12. character(len=1) :: bmat = 'I'
  13. character(len=2) :: which = 'SR'
  14. integer, parameter :: ishfts = 1, maxitr = 300, mode = 1
  15. integer :: ido = 0, lworkl = 3 * ncv**2+5*ncv, info = 1, n, ldv, ierr = 0
  16. complex(R_P) :: sigma = (0.0,0.0_R_P)
  17. real(R_P) :: tol = 0.0_R_P
  18. logical :: rvec = .true.
  19.  
  20. iparam = 0
  21. ipntr = 0
  22. select = .true.
  23. d = (0.0,0.0_R_P)
  24. v = (0.0,0.0_R_P)
  25. workd = (0.0,0.0_R_P)
  26. workev = (0.0,0.0_R_P)
  27. workl = (0.0,0.0_R_P)
  28. rwork = 0.0_R_P
  29. resid = gs_vec
  30. n = size(mat,1)
  31. ldv = n
  32. iparam(1) = ishfts
  33. iparam(3) = maxitr
  34. iparam(7) = mode
  35.  
  36. do while (ido == 0 .or. ido == 1 .or. ido == -1)
  37.  
  38. call znaupd( ido, bmat, n, which, nev, tol, resid, ncv,&
  39. v, ldv, iparam, ipntr, workd, workl, lworkl,&
  40. rwork,info )
  41.  
  42. if (ido .eq. -1 .or. ido .eq. 1) then
  43.  
  44. call av(workd(ipntr(1):(ipntr(1)+n-1)), workd(ipntr(2):(ipntr(2)+n-1)),mat)
  45.  
  46. end if
  47. end do
  48.  
  49. if ( info .lt. 0 ) then
  50.  
  51. print *, ' '
  52. print *, ' Error with _naupd, info = ', info
  53. print *, ' Check the documentation of _naupd'
  54. print *, ' '
  55.  
  56. else
  57.  
  58. rvec = .true.
  59.  
  60. call zneupd (rvec, 'A', select, d, v, ldv, sigma,&
  61. workev, bmat, n, which, nev, tol, resid, ncv,&
  62. v, ldv, iparam, ipntr, workd, workl, lworkl, &
  63. rwork, ierr)
  64.  
  65. if ( ierr .ne. 0) then
  66.  
  67. print *, ' '
  68. print *, ' Error with _neupd, info = ', ierr
  69. print *, ' Check the documentation of _neupd. '
  70. print *, ' '
  71.  
  72. end if
  73.  
  74. gs_vec = v(:,1)
  75. gs_val = d(1)
  76.  
  77. if ( info .eq. 1) then
  78. print *, ' '
  79. print *, ' Maximum number of iterations reached.'
  80. print *, ' '
  81. else if ( info .eq. 3) then
  82. print *, ' '
  83. print *, ' No shifts could be applied during implicit Arnoldi update, try increasing NCV.'
  84. print *, ' '
  85. end if
  86.  
  87. print *, ' '
  88. print *, '_NDRV1'
  89. print *, '====== '
  90. print *, ' '
  91. print *, ' Size of the matrix is ', n
  92. print *, ' The number of Ritz values requested is ', nev
  93. print *, ' The number of Arnoldi vectors generated (NCV) is ', ncv
  94. print *, ' What portion of the spectrum: ', which
  95. print *, ' The number of converged Ritz values is ', iparam(5)
  96. print *, ' The number of Implicit Arnoldi update iterations taken is ', iparam(3)
  97. print *, ' The number of OP*x is ', iparam(9)
  98. print *, ' The convergence criterion is ', tol
  99. print *, ' '
  100.  
  101. end if
  102.  
  103.  
  104.  
  105. contains
  106.  
  107. subroutine av ( v, w, mat)
  108.  
  109. complex(R_P), dimension(:), intent(in) :: v
  110. complex(R_P), dimension(:), intent(out) :: w
  111. complex(R_P), dimension(:,:), intent(in) :: mat
  112.  
  113. w = matmul(mat,v)
  114.  
  115. end subroutine av
  116.  
  117. end subroutine zndrv1
  118.  
  119. complex(R_P) :: sigma = (0.0,0.0_R_P)
  120.  
  121. complex(R_P) :: sigma
  122.  
  123. sigma = (0.0, 0.0_R_P)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement