c ****************************************************************** c subroutine mod_broydenr_VL2 c c Please feel free to use and/or modify this software as you see c fit. I simply ask that you reference the paper in which the c technique coded here is described: c Daniel Knapp et al. (2004), arXiv:cond-mat/0409701. c No guarantees are given that this program will work as advertised, c and I assume no responsibility for any problems associated with c possible (mal)function. c c This program makes extensive use of BLAS subroutines. c c written by Daniel Knapp, January 2004 c www.danielk.ca c ****************************************************************** subroutine mod_broydenr(N, nit_max, nit, alpha, psi_in, & psi_out, psi_in_old, F_old, weight, & A_mat, d_F, d_psi_in, U_mat) implicit none c ------------------------------------------------------------------ c imported variables c -> N: length of vector c -> nit: the iteration counter integer N, nit_max, nit c -> alpha: the initial inverse Jacobian is set to alpha*I double precision alpha c -> psi_in, psi_out: input and output vectors for the just completed c iteration double precision psi_in(N), psi_out(N) c -> psi_in_old: input vector for the prior iteration (nit-1) c -> F_old: the difference between the output and input vectors for c the prior iteration (nit-1) double precision psi_in_old(N), F_old(N) c -> weight: the weights associated with each iteration in the error c function that this routine seeks to minimize double precision weight(0:nit) c -> A_mat: matrix used in generating a guess for psi_in for the next c iteration double precision A_mat(nit_max, nit_max) c -> d_F: F - F_old, calculated and stored at every iteration double precision d_F(nit_max, N) c -> d_psi_in: psi_in - psi_in_old, calculated and stored at every c iteration double precision d_psi_in(nit_max, N) c -> U_mat: array of vectors used in generating a guess for psi_in for the c next iteration double precision U_mat(nit_max, N) c ------------------------------------------------------------------ c new variables integer i, j c -> B_mat: matrix, is inverse of (weight(0)^2/alpha * I - A_mat) double precision B_mat(nit, nit) c -> norm: magnitude of psi_in - psi_in_old double precision norm, temp c -> F: psi_out - psi_in, change in vector during the iteration c -> d_F_nit: F - F_old c -> d_psi_nit: psi_in - psi_in_old double precision F(N), d_F_nit(N), d_psi_in_nit(N) double precision tmp_vec1(N), tmp_vec2(N) double precision c_vec(nit), d_vec(nit) c ------------------------------------------------------------------ c functions and subroutines called double precision Ddot, Dnrm2 external Ddot, Dnrm2 c ->subroutines: all are from BLAS and handle vector and matrix c calculations in this routine external DCopy, Daxpy, Dscal, Dgemv c ================================================================== c ------------------------------------------------------------------ c initialize dummy variables do i = 1, N tmp_vec1(i) = 0.d0 tmp_vec2(i) = 0.d0 end do c ------------------------------------------------------------------ c define some quantities for intermediate steps, c -> |psi_out> = -1.d0 * |psi_in> + |psi_out> call Daxpy(N, -1.d0, psi_in, 1, psi_out, 1) c ->|F> = |psi_out> call DCopy(N, psi_out, 1, F, 1) if (nit.gt.0) then c -> |d_psi_in_nit> = |psi_in> call DCopy(N, psi_in, 1, d_psi_in_nit, 1) c -> |d_F_nit> = |F> call DCopy(N, F, 1, d_F_nit, 1) c ----------------------------------------------------------------- c iterations where the system has greatly changed are given less c weight in constructing the new Jacobian temp = Dnrm2(N, F, 1) / & ( & dble(N) & * dsqrt( ddot(N, psi_out, 1, psi_out, 1) & + ddot(N, psi_in, 1, psi_in, 1) ) & ) if(temp.lt.1) then weight(nit) = -log(temp) else weight(nit) = 1/dsqrt(temp) end if c ------------------------------------------------------------------ c |d_psi_in_nit> = ( |psi_in> - |psi_in_old> ) / norm c -> |d_psi_in_nit> = -1.d0 * |psi_in_old> + |d_psi_in_nit> call Daxpy(N, -1.d0, psi_in_old, 1, d_psi_in_nit, 1) c -> norm = DSqrt( ) norm = Dnrm2(N, d_psi_in_nit, 1) c -> |d_psi_in_nit> = weight(nit) * |d_psi_in_nit> / norm call DScal(N, weight(nit)/norm, d_psi_in_nit, 1) do i = 1, N d_psi_in(nit, i) = d_psi_in_nit(i) end do c ------------------------------------------------------------------ c |d_F_nit> = ( |F> - |F_old> ) / norm c -> |d_F_nit> = -1.d0 * |F_old> + |d_F_nit> call Daxpy(N, -1.d0, F_old, 1, d_F_nit, 1) c -> |d_F_nit> = weight(nit) * |d_F_nit> / norm call DScal(N, weight(nit)/norm, d_F_nit, 1) do i = 1, N d_F(nit, i) = d_F_nit(i) end do c ------------------------------------------------------------------ c add new elements to matrix A, c a(k,l) = do i = 1, nit - 1 do j = 1, N tmp_vec1(j) = d_F(i, j) tmp_vec2(j) = d_psi_in(i, j) end do A_mat(i, nit) = ddot(N, tmp_vec2, 1, d_F_nit, 1) A_mat(nit, i) = ddot(N, d_psi_in_nit, 1, tmp_vec1, 1) end do A_mat(nit, nit) = ddot(N, d_psi_in_nit, 1, d_F_nit, 1) c ------------------------------------------------------------------ c initialize matrix to be inverted B = [(w(0)^2 / alpha) * I + A] do j = 1, nit do i = 1, nit B_mat(i,j) = -A_mat(i,j) end do B_mat(j, j) = B_mat(j, j) + weight(0)*weight(0) / alpha end do c ------------------------------------------------------------------ c invert matrix of size nit x nit call InvertMatrix(B_mat, nit) c ------------------------------------------------------------------ c calculate c_i^m = c equivalent to a vector times a matrix: c c[nit] = (d_psi_in[nit x N]) |F(nit)[N]> c -> |c_vec[nit]> = 1.d0 * (d_psi_in[nit x N]) |F[N]> c + 0.d0 * |c_vec[nit]> call DGEMV('N', nit, N, 1.d0, d_psi_in, nit_max, & F, 1, 0.d0, c_vec, 1) c ------------------------------------------------------------------ c now calculate |u(nit)> = alpha*|d_F_nit> + |d_psi_in_nit> c -> |d_psi_in_nit> = alpha * |d_F_nit> + |d_psi_in_nit> call Daxpy(N, alpha, d_F_nit, 1, d_psi_in_nit, 1) do i = 1, N U_mat(nit, i) = d_psi_in_nit(i) end do c ------------------------------------------------------------------ c calculate d_l(m) = sum_{k=1}^{nit} c_k(m) b_{kl}, or c |d(m)> = B_mat |c(m)> c -> |d_vec> = 1.d0 * (B_mat) * |c_vec> + 0.d0 * |d_vec> call DGEMV('N', nit, nit, 1.d0, B_mat, nit, c_vec, 1, & 0.d0, d_vec, 1) c ------------------------------------------------------------------ c put all current var-s into old var-s for the next step c -> |psi_in_old> = |psi_in> call Dcopy(N, psi_in, 1, psi_in_old, 1) c -> |F_old> = |F> call Dcopy(N, F, 1, F_old, 1) c ------------------------------------------------------------------ c now update the value of next guess c blas matrix multiplication: new psi_in = psi_in + alpha*F c - sum_{i=1}^{nit} d_vec_i^m * |u_i(n)> c or |psi_in[N] = |psi_in[N]> + alpha*|F[N]> c + transpose(U_mat[nit x N])|d_vec[nit]> c -> |psi_in[N]> = |psi_in[N]> + transpose(U_mat[nit x N]) |d_vec[nit]> call DGEMV('T', nit, N, 1.d0, U_mat, nit_max, & d_vec, 1, 1.d0, psi_in, 1) c -> |psi_in> = |psi_in> + alpha |F> call DAXPY(N, alpha, F, 1, psi_in, 1) else c ------------------------------------------------------------------ c put all current var-s into old var-s for the next step c -> |psi_in_old> = |psi_in> call Dcopy(N, psi_in, 1, psi_in_old, 1) c -> |F_old> = |F> call Dcopy(N, F, 1, F_old, 1) c ------------------------------------------------------------------ c update first guess c -> |psi_in> = alpha |F> + |psi_in> call DAXPY(N, alpha, F, 1, psi_in, 1) end if end c ******************************************************************* c End of subroutine broydenr c ******************************************************************* c ******************************************************************* subroutine InvertMatrix(A, N) c ******************************************************************* c Uses LAPACK routines to invert a matrix implicit none integer N double precision A(N,N) double precision work(N) integer ipiv(N) integer info call dgetrf(n,n,A,n,ipiv,info) ! get LU of a if(info.ne.0) stop 'DGETRF: error in LU factorization' call dgetri(n,A,n,ipiv,work,n,info) ! now get the inverse if(info /= 0) stop 'DGETRI: error in InvertMatrixing a' end