module brown_m !! Contains routines for evaluating **B**._d_**W** in the SDE for Brownian Dynamics !! simulation. use iso_c_binding, only: c_loc, c_ptr, c_f_pointer use mkl_blas use mkl_lapack use constants_m use random_m use strings_m use logger_m implicit none private public :: brn_init, brn_finish, brn_calc_dw, brn_calc_bdw character(4) :: meth !!Method to use for generating Brownian terms. Possible values are !!'CHOL' or 'LANC'. integer :: sizm !!Size of the mobility matrix integer :: s !!Number of Brownian vectors to be generated; equivalent to block size in the !!block Lanczos algorithm. integer :: dimk_min = 2 !!Minimum dimension of KSP integer :: dimk_max !!Maximum dimension of KSP real(rp) :: ethres = 1.0e-3_rp !!Error threshold for KSP-based method. real(rp), dimension(:,:), allocatable :: bdw_old !!(`sizm`, `s`) array. Stores B.dw at previous iteration. real(rp), dimension(:,:), allocatable :: lncv !!(`sizm`, `s*dimk_max`) array. Stores the Lanczos vectors. real(rp), dimension(:), allocatable :: eigvlh !!(`s*dimk_max`,) array. Eigenvalues of the tridiagonal/block triadiagonal !! matrix H generated by the Lanczos algorithm. real(rp), dimension(:,:), allocatable :: eigvch !!(`s*dimk_max`, `s*dimk_max`) array. Eigenvectors of the matrix H. real(rp), dimension(:,:), allocatable :: w !!(`sizm`, `s`) array. Used for generating the Lanczos vectors. real(rp), dimension(:), allocatable :: h_d !!(`dimk_max`,) array. Diagonal of the symmetric tridiagonal matrix H !! generated in case of `s = 1`. real(rp), dimension(:), allocatable :: h_sd !!(`dimk_max`,) array. Subdiagonal of the symmetric tridiagonal matrix H !! generated in case of `s = 1`. real(rp), dimension(:), allocatable :: tau !! (`s`,) array. Used in the QR decomposition of `w`. real(rp), dimension(:,:), allocatable :: matr !! (`s`,`s`) array. Used in the QR decomposition of `w`. real(rp), dimension(:,:), allocatable :: h_bd !!(`s`, `s*dimk_max`) array. Diagonal blocks of the matrix H. !!Used for KSP with `s > 1`. real(rp), dimension(:,:), allocatable :: h_bsd !!(`s`, `s*dimk_max`) array. Subdiagonal blocks of the matrix H. !!Used for KSP with `s > 1`. real(rp), dimension(:,:), allocatable :: h_bnd !!(`2*s`, `s*dimk_max`) array. Matrix H in banded storage. !!Used for KSP with `s > 1`. real(rp), dimension(:,:), allocatable :: tmp !!(`s*dimk_max`, `2*s`) array used as temporary storage space. real(rp), dimension(:), allocatable :: work !!Workspace array for LAPACK routines. integer, dimension(:), allocatable :: iwork !!Workspace array for LAPACK routines. contains !******************************************************************************* subroutine brn_init(n, nblks, mth, f, e) !! Performs initial setup. integer, intent(in) :: n !! Size of the mobility matrix. integer, intent(in) :: nblks !! Number of **B**._d_**W** vectors to generate for the same mobility matrix. !! `1 <= nblks < n`. character(4), intent(in) :: mth !! Method for calculating the Brownian terms. `mth = 'CHOL'` for !! Cholesky decomposition; `mth = 'LANC'` for KSP-based method. integer, intent(in), optional :: f !! For KSP-based method, maximum number of iterations. !! Must be present if `mth = 'LANC'`. `f` must be less than `n`. real(rp), intent(in), optional :: e !! For KSP-based method, error threshold for convergence. !! Must be present if `mth = 'LANC'`. sizm = n; s = nblks; meth = mth if (meth == 'CHOL') then !Nothing to do here. else if (meth == 'LANC') then dimk_max = max(dimk_min, f) + 1 !Adding 1 to ensure dimk_max > dimk_min ethres = e allocate( bdw_old(sizm, s) ) allocate( lncv(sizm,s*dimk_max) ) allocate( eigvlh(s*dimk_max) ) allocate( eigvch(s*dimk_max,s*dimk_max) ) allocate( w(sizm,s)) allocate( tmp(s*dimk_max,2*s)) if (s == 1) then allocate( h_d(dimk_max) ) allocate( h_sd(dimk_max) ) else allocate( tau(s)) allocate( matr(s,s)) allocate( h_bd(s,s*dimk_max) ) allocate( h_bsd(s,s*dimk_max) ) allocate( h_bnd(2*s,s*dimk_max)) end if end if end subroutine !******************************************************************************* subroutine brn_finish() !! Performs cleanup. if (allocated(bdw_old)) deallocate(bdw_old) if (allocated(lncv)) deallocate(lncv) if (allocated(eigvlh)) deallocate(eigvlh) if (allocated(eigvch)) deallocate(eigvch) if (allocated(w)) deallocate(w) if (allocated(tmp)) deallocate(tmp) if (allocated(h_d)) deallocate(h_d) if (allocated(h_sd)) deallocate(h_sd) if (allocated(tau)) deallocate(tau) if (allocated(matr)) deallocate(matr) if (allocated(h_bd)) deallocate(h_bd) if (allocated(h_bsd)) deallocate(h_bsd) if (allocated(h_bnd)) deallocate(h_bnd) if (allocated(work)) deallocate(work) if (allocated(iwork)) deallocate(iwork) end subroutine !******************************************************************************* subroutine brn_calc_bdw(mob, bdw, ierr, lconv, f) !! Main driver routine for calculating **B**._d_**W**. real(rp), dimension(:,:), intent(in out) :: mob !! (`sizm`,`sizm`) symmetric positive definite matrix. !! On _entry_, contains the mobility matrix. On _return_, if `meth = 'CHOL'`, !! the upper triangular part is overwritten with the result of Cholesky !! decomposition. The strictly lower triangular part is left intact. real(rp), dimension(:,:), intent(in out) :: bdw !! (`sizm`,`s`) matrix. !! On _entry_, contains `s` column vectors drawn from standard normal !! distribution. On _return_, the `s` columns are overwritten with the !! result of **B**.d**W**. integer, intent(out) :: ierr !! Error flag. logical, intent(out), optional :: lconv !! If `meth = 'LANC'`, returns *true* if converged, *false* otherwise. integer, intent(out), optional :: f !! If `meth = 'LANC'`, the number of iterations performed. integer :: f_ logical :: lconv_ if (meth == 'CHOL') then call calc_bdw_cholesky(mob, bdw, ierr) else if (meth == 'LANC') then if (s == 1) then call calc_bdw_lanc(mob, bdw, lconv_, f_, ierr) else call calc_bdw_blanc(mob, bdw, lconv_, f_, ierr) end if end if if (present(lconv)) lconv = lconv_ if (present(f)) f = f_ end subroutine !******************************************************************************* subroutine calc_bdw_lanc(mob, bdw, lconv, f, ierr) !! Implements Lanczos algorithm for calculating **B**._d_**W**. real(rp), dimension(:,:), intent(in) :: mob !! (`sizm`,`sizm`) symmetric positive definite matrix. real(rp), dimension(:,:), intent(in out) :: bdw !! (`sizm`,1) array. !! On _entry_, contains a vector drawn from standard normal !! distribution. On _return_, is overwritten with the result of !! **B**._d_**W**. logical, intent(out) :: lconv !! Returns *true* if converged, *false* otherwise. integer, intent(out) :: f !! The number of iterations performed. integer, intent(out) :: ierr !! Error flag. integer :: j, info, lwork, liwork real(rp), dimension(1) :: qwork !Size 1 array for workspace query. integer, dimension(1) :: qiwork !Size 1 array for workspace query. real(rp) :: normz, err character(len=256) :: msg lconv = .false.; ierr = 0; f = 0 lncv = 0.0_rp; h_d = 0.0_rp; h_sd = 0.0_rp eigvlh = 0.0_rp; eigvch = 0.0_rp w = 0.0_rp; tmp = 0.0_rp; bdw_old = 0.0_rp !Workspace query if ((.not. allocated(work)) .or. (.not. allocated(iwork))) then call dstevd('V', dimk_max, eigvlh, tmp(1:dimk_max,1), eigvch, dimk_max, & qwork, -1, qiwork, -1, info) if (info /= 0) then write(msg,'(a,i0,a)') '<calc_bdw_lanc> dstevd err = ', info, & ' on workspace query' call logger%log_msg(msg) ierr = 1; return end if lwork = int(qwork(1)); liwork = qiwork(1) allocate( work(lwork) ); allocate( iwork(liwork) ) end if !Lanczos process normz = norm2(bdw(:,1)) lncv(:,1) = bdw(:,1)/normz do j = 1, dimk_max f = f + 1 !Iteration count call dsymv('L', sizm, 1.0_rp, mob, sizm, lncv(:,j), 1, 0.0_rp, w(:,1), 1) if (j > 1) then w(:,1) = w(:,1) - h_sd(j-1)*lncv(:,j-1) end if h_d(j) = dot_product(w(:,1), lncv(:,j)) if (j < dimk_max) then w(:,1) = w(:,1) - h_d(j)*lncv(:,j) h_sd(j) = norm2(w(:,1)) lncv(:,j+1) = w(:,1)/h_sd(j) end if !Calculate B.dW and check for convergence if (j >= dimk_min) then eigvlh(1:j) = h_d(1:j) tmp(1:j,1) = h_sd(1:j) !Eigen decomposition of H = P * d * P^T !eigvlh[1:j] <- d; eigvch[1:j,1:j] <- P call dstevd('V', j, eigvlh(1:j), tmp(1:j,1), eigvch, dimk_max, & work, size(work), iwork, size(iwork), info) if (info /= 0) then write(msg,'(a,i0,a,i0)') '<calc_bdw_lanc> dstevd err = ', info, & ' on iteration ', j call logger%log_msg(msg) ierr = 1; return end if !eigvlh(1:j) <- d^(1/2) eigvlh(1:j) = sqrt(eigvlh(1:j)) !tmp[1:j,1] <- d^(1/2) * P^T * e_1, where e_1 if the first column of !an identity matrix of size j. tmp(1:j,1) = eigvch(1,1:j) tmp(1:j,1) = eigvlh(1:j)*tmp(1:j,1) !tmp[1:j,2] <- P * tmp[1:j,1] call dgemv('N', j, j, 1.0_rp, eigvch, dimk_max, tmp(1:j,1), 1, & 0.0_rp, tmp(1:j,2), 1) !bdw <- normz * lncv[:,1:j] * tmp[1:j,2] call dgemv('N', sizm, j, normz, lncv, sizm, tmp(1:j,2), 1, & 0.0_rp, bdw, 1) !Error calculation for j > dimk_min. For j = dimk_min just save !bdw for error calculation in the next iteration. if (j > dimk_min) then err = norm2(bdw(:,1) - bdw_old(:,1))/norm2(bdw_old(:,1)) if (err <= ethres) then lconv = .true.; exit end if end if bdw_old = bdw end if end do end subroutine !******************************************************************************* subroutine calc_bdw_blanc(mob, bdw, lconv, f, ierr) !! Implements block Lanczos algorithm for calculating **B**._d_**W**. real(rp), dimension(:,:), intent(in) :: mob !! (`sizm`,`sizm`) symmetric positive definite matrix. real(rp), dimension(:,:), intent(in out) :: bdw !! (`sizm`,`s`) matrix. !! On _entry_, contains `s` column vectors drawn from standard normal !! distribution. On _return_, the `s` columns are overwritten with the !! result of **B**._d_**W**. logical, intent(out) :: lconv !! Returns *true* if converged, *false* otherwise. integer, intent(out) :: f !! The number of iterations performed. integer, intent(out) :: ierr !! Error flag. real(rp) :: err integer :: j, jb, irow, jcol integer :: jgbeg, jgend, igbeg, ig, jg integer :: info, lwork, liwork real(rp), dimension(1) :: qwork integer, dimension(1) :: qiwork character(len=256) :: msg lconv = .false.; ierr = 0; f = 0 lncv = 0.0_rp; h_bd = 0.0_rp; h_bsd = 0.0_rp eigvlh = 0.0_rp; eigvch = 0.0_rp w = 0.0_rp; tmp = 0.0_rp; bdw_old = 0.0_rp tau = 0.0_rp; matr = 0.0_rp; h_bnd = 0.0_rp !Workspace query if (.not. allocated(work)) then call dgeqp3(sizm, s, bdw, sizm, qiwork, tau, qwork, -1, info) if (info /= 0) then write(msg,'(a,i0,a)') '<calc_bdw_blanc> dgeqp3 err = ', info, & ' on workspace query' call logger%log_msg(msg) ierr = 1; return end if lwork = int( qwork(1) ) call dorgqr(sizm, s, s, bdw, sizm, tau, qwork, -1, info) if (info /= 0) then write(msg,'(a,i0,a)') '<calc_bdw_blanc> dorgqr err = ', info, & ' on workspace query' call logger%log_msg(msg) ierr = 1; return end if lwork = max( lwork, int(qwork(1)) ) call dsbevd('V', 'L', dimk_max*s, 2*s-1, h_bnd, 2*s, & eigvlh, eigvch, s*dimk_max, qwork, -1, qiwork, -1, info) if (info /= 0) then write(msg,'(a,i0,a)') '<calc_bdw_blanc> dsbevd err = ', info, & ' on workspace query' call logger%log_msg(msg) ierr = 1; return end if lwork = max( lwork, int(qwork(1)) ) liwork = qiwork(1) allocate(work(lwork)); allocate(iwork(liwork)) end if !Lanczos process !Reduced QR factorization of bdw with pivoting. The results overwrite bdw. iwork = 0; qwork = 0.0_rp call dgeqp3(sizm, s, bdw, sizm, iwork, tau, work, size(work), info) if (info /= 0) then write(msg,'(a,i0,a)') '<calc_bdw_blanc> dgeqp3 err = ', info, & ' on iteration 0' call logger%log_msg(msg) ierr = 1; return end if !Extract the (s,s) permuted upper triangular matrix R and store it in matr. do jcol = 1, s matr(1:jcol, iwork(jcol)) = bdw(1:jcol, jcol) end do !Extract the (sizm,s) matrix Q and store it in lncv[:,1:s]. call dorgqr(sizm, s, s, bdw, sizm, tau, work, size(work), info) if (info /= 0) then write(msg,'(a,i0,a)') '<calc_bdw_blanc> dorgqr err = ', info, & ' on iteration 0' call logger%log_msg(msg) ierr = 1; return end if lncv(:,1:s) = bdw do j = 1, dimk_max f = f + 1 !Iteration count !W <- mob * V_j jgbeg = (j-1)*s+1; jgend = j*s call dsymm('L', 'L', sizm, s, 1.0_rp, mob, sizm, & lncv(:,jgbeg:jgend), sizm, 0.0_rp, w, sizm ) if (j > 1) then !W -> W - V_(j-1) * H_(j-1,j) ! = W - V_(j-1) * H_(j,j-1)^T, using the subdiagonal block jgbeg = (j-2)*s+1; jgend = (j-1)*s call dgemm('N', 'T', sizm, s, s, -1.0_rp, lncv(:,jgbeg:jgend), sizm, & h_bsd(:,jgbeg:jgend), s, 1.0_rp, w, sizm) end if !H_(j,j) <- V_j^T * W jgbeg = (j-1)*s+1; jgend = j*s call dgemm('T', 'N', s, s, sizm, 1.0_rp, lncv(:,jgbeg:jgend), sizm, & w, sizm, 0.0_rp, h_bd(:,jgbeg:jgend), s) if (j < dimk_max) then !W -> W - V_(j) * H_(j,j) !Note: H_(j,j) is symmetric. jgbeg = (j-1)*s+1; jgend = j*s call dsymm('R', 'L', sizm, s, -1.0_rp, h_bd(:,jgbeg:jgend), & s, lncv(:,jgbeg:jgend), sizm, 1.0_rp, w, sizm) !Reduced QR factorization of w. iwork = 0 call dgeqp3(sizm, s, w, sizm, iwork, tau, work, size(work), info) if (info /= 0) then write(msg,'(a,i0,a,i0)') '<calc_bdw_blanc> dgeqp3 err = ', info, & ' on iteration ', j call logger%log_msg(msg) ierr = 1; return end if !Extract the (s,s) permuted upper triangular matrix R and store !it in H_(j+1,j). do jcol = 1, s h_bsd(1:jcol, (j-1)*s+iwork(jcol)) = w(1:jcol, jcol) end do !Extract the (sizm,s) matrix Q and store it in V_(j+1). call dorgqr(sizm, s, s, w, sizm, tau, work, size(work), info) if (info /= 0) then write(msg,'(a,i0,a,i0)') '<calc_bdw_blanc> dorgqr err = ', info, & ' on iteration ', j call logger%log_msg(msg) ierr = 1; return end if jgbeg = j*s+1; jgend = (j+1)*s lncv(:,jgbeg:jgend) = w end if !Calculate B.dW and check for convergence if (j >= dimk_min) then jgbeg = (j-1)*s+1; jgend = j*s !Create matrix H in band storage format, in preparation for eigen !decomposition with symmetric band matrix routine. h_bnd = 0.0_rp do jb = 1, j jgbeg = (jb-1)*s; igbeg = jgbeg do jcol = 1,s jg = jgbeg + jcol do irow = jcol,s ig = igbeg + irow h_bnd(ig-jg+1,jg) = h_bd(irow, jg) end do end do end do do jb = 1, (j-1) jgbeg = (jb-1)*s; igbeg = jb*s do jcol = 1,s jg = jgbeg + jcol do irow = 1,s ig = igbeg + irow h_bnd(ig-jg+1,jg) = h_bsd(irow, jg) end do end do end do !Eigen decomposition of H = P * d * P^T !eigvlh <- d; eigvch <- P call dsbevd('V', 'L', j*s, 2*s-1, h_bnd(:,1:j*s), 2*s, eigvlh, & eigvch, s*dimk_max, work, size(work), iwork, & size(iwork), info) if (info /= 0) then write(msg,'(a,i0,a,i0)') '<calc_bdw_blanc> dsbevd err = ', info, & ' on iteration ', j call logger%log_msg(msg) ierr = 1; return end if !Square root of the eigen values eigvlh(1:j*s) = sqrt(eigvlh(1:j*s)) !tmp[1:j*s,1:s] <- eigvch[1:s,1:j*s]^T * matr call dgemm('T', 'N', j*s, s, s, 1.0_rp, eigvch, s*dimk_max, & matr, s, 0.0_rp, tmp(:,1:s), s*dimk_max) !tmp[1:j*s,1:s] <- eigvlh[1:j*s] * tmp[1:j*s,1:s] do jcol = 1, s tmp(1:j*s,jcol) = eigvlh(1:j*s)*tmp(1:j*s,jcol) end do !tmp[1:j*s,s+1:2*s] <- eigvch[1:j*s,1:j*s] * tmp[1:j*s,1:s] call dgemm('N', 'N', j*s, s, j*s, 1.0_rp, eigvch, s*dimk_max, & tmp(:,1:s), s*dimk_max, 0.0_rp, tmp(:,s+1:2*s), s*dimk_max) !bdw <- lncv[:,1:j*s] * tmp[1:j*s,s+1:2*s] call dgemm('N', 'N', sizm, s, j*s, 1.0_rp, lncv, sizm, & tmp(:,s+1:2*s), s*dimk_max, 0.0_rp, bdw, sizm) !Error calculation for j > dimk_min. For j == dimk_min just save !bdw for error calculation in the next iteration. Using 1-norm as !the matrix norm here. if (j > dimk_min) then err = maxval(sum(abs(bdw-bdw_old), 1)) & / maxval(sum(abs(bdw_old), 1)) if (err <= ethres) then lconv = .true.; f = j exit end if end if bdw_old = bdw end if end do end subroutine !****************************************************************************** subroutine calc_bdw_cholesky(mob, bdw, ierr) !! Calculates **B**._d_**W** using Cholesky decomposition. real(rp), dimension(:,:), intent(in out) :: mob !! (`sizm`,`sizm`) symmetric positive definite matrix. !! On _entry_, contains the mobility matrix. On _return_, !! the upper triangular part is overwritten with the result of Cholesky !! decomposition. The strictly lower triangular part is left intact. real(rp), dimension(:,:), intent(in out) :: bdw !! (`sizm`,`s`) matrix. !! On _entry_, contains `s` column vectors drawn from standard normal !! distribution. On _return_, the `s` columns are overwritten with the !! result of **B**._d_**W**. integer, intent(out) :: ierr !! Error flag. integer :: info character(len=256) :: msg ierr = 0 !Cholesky decomposition: mob = U^T * U. The upper triangular part of mob is !overwritten with U. The strictly lower triangular part of mob is not !accessed. call dpotrf('U', sizm, mob, sizm, info) if (info /= 0) then write(msg,'(a,i0)') '<calc_bdw_cholesky> dpotrf err = ', info call logger%log_msg(msg) ierr = 1; return end if !Product of U (stored in the upper triangular part of mob) with each of the !random vectors: bdw = U*bdw if (s == 1) then call dtrmv('U', 'N', 'N', sizm, mob, sizm, bdw(:,1), 1) else call dtrmm('L', 'U', 'N', 'N', sizm, s, 1.0_rp, mob, sizm, bdw, sizm) end if end subroutine !******************************************************************************* subroutine brn_calc_dw(dw) !! Calculates standard normally distributed random vectors. real(rp), dimension(:,:), target, intent(out) :: dw !! (`n`,`m`) array. real(rp), dimension(:), pointer :: vec_dw type(c_ptr) :: loc_dw integer :: n, m n = size(dw,1); m = size(dw,2) dw = 0.0_rp loc_dw = c_loc(dw) call c_f_pointer(loc_dw, vec_dw, [n*m]) call get_rv_gaussian(0.0_rp, 1.0_rp, vec_dw, 1000000) end subroutine !******************************************************************************* subroutine calc_bdw_lanc_sor(mob, bdw, lconv, f, ierr) real(rp), dimension(:,:), intent(in) :: mob !! (`sizm`,`sizm`) symmetric positive definite matrix. real(rp), dimension(:,:), intent(in out) :: bdw !! (`sizm`,1) array. !! On _entry_, contains a vector drawn from standard normal !! distribution. On _return_, is overwritten with the result of !! **B**.d**W**. logical, intent(out) :: lconv !! Returns *true* if converged, *false* otherwise. integer, intent(out) :: f !! The number of iterations performed. integer, intent(out) :: ierr !! Error flag. integer :: j, info, lwork, liwork, k real(rp), dimension(1) :: qwork !Size 1 array for workspace query. integer, dimension(1) :: qiwork !Size 1 array for workspace query. real(rp), dimension(dimk_max) :: errbnd real(rp), dimension(sizm) :: ritzv real(rp) :: normz, err, normm character(len=256) :: msg lconv = .false.; ierr = 0; f = 0 tmp = 0.0_rp; eigvlh = 0.0_rp; eigvch = 0.0_rp h_d = 0.0_rp; h_sd = 0.0_rp bdw_old = 0.0_rp; lncv = 0.0_rp; w = 0.0_rp !Workspace query if ((.not. allocated(work)) .or. (.not. allocated(iwork))) then call dstevd('V', dimk_max, eigvlh, tmp(1:dimk_max,1), eigvch, dimk_max, & qwork, -1, qiwork, -1, info) if (info /= 0) then write(msg,'(a,i0,a)') '<calc_bdw_lanc_sor> dstevd err = ', info, & ' on workspace query' call logger%log_msg(msg) ierr = 1; return end if lwork = int(qwork(1)); liwork = qiwork(1) allocate( work(lwork) ); allocate( iwork(liwork) ) end if !Lanczos process normz = norm2(bdw(:,1)) normm = maxval(sum(abs(mob),1)) lncv(:,1) = bdw(:,1)/normz do j = 1, dimk_max f = f + 1 !Iteration count call dsymv('L', sizm, 1.0_rp, mob, sizm, lncv(:,j), 1, 0.0_rp, w(:,1), 1) if (j > 1) then w(:,1) = w(:,1) - h_sd(j-1)*lncv(:,j-1) end if h_d(j) = dot_product(w(:,1), lncv(:,j)) if (j < dimk_max) then w(:,1) = w(:,1) - h_d(j)*lncv(:,j) do k = 1, j if (errbnd(k) < sqrt(epsilon(1.0_rp)*normm)) then call dgemv('N', sizm, j, 1.0_rp, lncv, sizm, eigvch(:,k), 1, & 0.0_rp, ritzv, 1) w(:,1) = w(:,1) - dot_product(ritzv,w(:,1))*ritzv end if end do h_sd(j) = norm2(w(:,1)) lncv(:,j+1) = w(:,1)/h_sd(j) end if !Compute eigenvalues, eigenvectors and error bounds for H if (j > 1) then eigvlh(1:j) = h_d(1:j) tmp(1:j,1) = h_sd(1:j) !Eigen decomposition of H = P * d * P^T !eigvlh[1:j] <- d; eigvch[1:j,1:j] <- P call dstevd('V', j, eigvlh(1:j), tmp(1:j,1), eigvch, dimk_max, & work, size(work), iwork, size(iwork), info) if (info /= 0) then write(msg,'(a,i0,a,i0)') '<calc_bdw_lanc_sor> dstevd err = ', info, & ' on iteration ', j call logger%log_msg(msg) ierr = 1; return end if errbnd(1:j) = abs(h_sd(j))*eigvch(j,1:j) end if !Calculate B.dW and check for convergence if (j >= dimk_min) then !eigvlh(1:j) <- d^(1/2) eigvlh(1:j) = sqrt(eigvlh(1:j)) !tmp[1:j,1] <- d^(1/2) * P^T * e_1, where e_1 if the first column of !an identity matrix of size j. tmp(1:j,1) = eigvch(1,1:j) tmp(1:j,1) = eigvlh(1:j)*tmp(1:j,1) !tmp[1:j,2] <- P * tmp[1:j,1] call dgemv('N', j, j, 1.0_rp, eigvch, dimk_max, tmp(1:j,1), 1, & 0.0_rp, tmp(1:j,2), 1) !bdw <- normz * lncv[:,1:j] * tmp[1:j,2] call dgemv('N', sizm, j, normz, lncv, sizm, tmp(1:j,2), 1, & 0.0_rp, bdw, 1) !Error calculation for j > dimk_min. For j = dimk_min just save !bdw for error calculation in the next iteration. if (j > dimk_min) then err = norm2(bdw(:,1) - bdw_old(:,1))/norm2(bdw_old(:,1)) !write(*,*) j, 'err = ', err if (err <= ethres) then lconv = .true. exit end if end if bdw_old = bdw end if end do end subroutine !******************************************************************************* end module brown_m