parameter (n=5) real*8 A(n,n), S(n,n), w(14*n) integer iw(5*n), i, j do j=1, n do i=j, n A(i,j) = 1.0D0/dfloat(i+j) enddo enddo call DSYEXP(n,A,S,w,iw) do i=1, n print*, (S(i,j),j=1,i) enddo end ******************************************************************** * This is our routine for calculating exp(-A) for a symmetric * matrix A. Maybe we will change it to exp(A) later. * The theory is in the following paper: * * Ya Yan Lu, Exponentials of Symmetric Matrices through * Tridiagonal Reductions, * Linear Algebra and its Applications, Vol 279(1-3), * pp. 317-324, 1998. * * The program requires LAPACK to work. * * New prime fraction coefficients calculated by Jenny Ho * are now incorperated in this program. ****************************************************************** subroutine DSYEXP(n,A,S,w,iw) integer n integer iw(5*n) real*8 A(n,n), S(n,n), w(14*n) * call DSYEXW(n,A,S,w(1), w(n+1), w(2*n+1),w(3*n+1),iw) return end * ******************************************************************** subroutine DSYEXW(n,A,S,d,e,tau,w,iw) implicit none integer n integer iw(5*n) real*8 A(n,n),S(n,n),d(n),e(n-1),tau(n), w(11*n) real*8 lams,vl,vu,abstol integer lwork,info,mfound,nsplit,i,j,il,iu lwork = 11*n * Reduce the matrix A to tridiagonal form by LAPACK * routine xSYTRD. call DSYTRD( 'L', n, A, n, d, e, tau, w, lwork, info ) * print*, 'Best lwork = ', w(1) * print*, 'DSYTRD info = ', info * Find the largest and smallest eigenvalues by LAPACK * routine xSTEBZ. vl = 0.0D0 vu = 0.0D0 il = 1 iu = 1 abstol = 0.0D0 call DSTEBZ( 'I', 'B', n, vl,vu,il,iu,abstol, d, e, $ mfound, nsplit, w, iw, iw(n+1), w(n+1), iw(2*n+1), $ info ) * print*, 'DSTEBZ infp = ', info lams = w(1) * print*, 'The smallest eigenvalue is ', lams * Now, we calculate the exp(-(T-lambda_1)) do i=1, n d(i) = d(i) - lams enddo call DSTEXP(n,d,e,S,w,iw) * print*, 'Exp of Trid done!' * Given symmetric matrix S (lower triangular), we calculate QSQ^t * Q is given as in the xSYTRD, stored in elementary reflectors. call APPLIQ(n,A,tau,S, w, w(n+1)) * print*, 'AppliQ done!' * Scale the answer by exp(-lambda_1) lams = exp(-lams) do j=1, n do i=j,n S(i,j) = lams*S(i,j) enddo enddo return end ************************************************************ * This is our program for calculating the negative * exponential of a real symmetric positive definite * matrix T. The diagonal of T is the vector d, * the upper /lower diagonal of T is the vector e * * T = diag(d(1:n)) + diag(e(1:n-1),1)+diag(e(1:n-1),-1) * * Then, we calculate exp(-T), assuming T is symmetric * positive semi-definite. * * The approximation is given as * * B = alpha0 I + SUM alpha(i) ( T - theta(i) I)^(-1) * * for i=1, 2, ..., 7. ************************************************************* subroutine DSTEXP(n,d,e,B,w,iw) implicit none integer n, i, j, k, iw(n) real*8 d(n), e(n-1), B(n,n), alpha0, w(10*n) complex*16 alpha(7), theta(7) alpha0 = 1.832174378247001445178D-14 alpha(1)=( 5.575032388013989775304D+01, 2.0429467998037060112D+02) alpha(2)=(-9.386654897738416781398D+01,-9.1287299537244701102D+01) alpha(3)=( 4.699646418199852101538D+01, 1.1616718259322739191D+01) alpha(4)=(-9.614224197620382854756D+00, 2.6419587674922023894D+00) alpha(5)=( 7.527200775595744942100D-01,-6.7036694058767631296D-01) alpha(6)=(-1.887805062128241188701D-02, 3.4369583916841158280D-02) alpha(7)=( 1.430857612700760251239D-04,-2.8722086698954584504D-04) theta(1)=(-5.623142572742583221495D+00,-1.1940690463436614062D+00) theta(2)=(-5.089345060577111615918D+00,-3.5888240290260982584D+00) theta(3)=(-3.993369710574816069647D+00,-6.0048316422335547350D+00) theta(4)=(-2.269783829227002748129D+00,-8.4617379730382146007D+00) theta(5)=( 2.087586382547148287361D-01,-1.0991260561898794116D+01) theta(6)=( 3.703275049428656060120D+00,-1.3656371871480423129D+01) theta(7)=( 8.897773186474949120586D+00,-1.6630982619899044992D+01) do j=1, n B(j,j) = alpha0 enddo do j=1,n-1 do i=j+1,n B(i,j) = 0.0D0 enddo enddo do k=1, 7 call BUPDAT(n,d,e,alpha(k),theta(k),B,w,w(2*n+1), # w(4*n+1),w(6*n+1), w(8*n+1), iw) enddo return end * * Here we calculate B := B + alpha(k) * ( T - theta(k) I)^{-1} * subroutine BUPDAT(n,d,e,alp,tht,B,zd,zl,zu,du2,bb,ipiv) implicit none integer n real*8 d(n), e(n-1), B(n,n) complex*16 alp, tht, zd(n), zl(n-1), zu(n-1), du2(n-2), bb(n) integer ipiv(n), info, nrhs, ldb, i, j do i=1, n zd(i) = d(i) - tht enddo do i=1, n-1 zl(i) = e(i) zu(i) = e(i) enddo call ZGTTRF( n, zl, zd, zu, du2, ipiv, info ) nrhs = 1 ldb = n do j=1, n do i=1, n bb(i)=0.0D0 enddo bb(j) = alp call ZGTTRS( 'N', n, nrhs, zl, zd, zu, du2, ipiv, bb, ldb, $ info ) do i=j,n B(i,j) = B(i,j) + dreal(bb(i)) enddo enddo return end ****************************************************************** * This program requires 2n^3 flops. * * SSYTRD reduces a real symmetric matrix A to real symmetric * tridiagonal form T by an orthogonal similarity transformation: * Q**T * A * Q = T, where Q is is represented as a product of * elementary reflectors. * * It takes: 4/3 n^3 operations. * * If UPLO = 'L', the matrix Q is represented as a product of elementary * reflectors * * Q = H(1) H(2) . . . H(n-1). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a real scalar, and v is a real vector with * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), * and tau in TAU(i). * * The contents of A on exit are illustrated by the following examples * with n = 5: * * if UPLO = 'U': if UPLO = 'L': * * ( d e v2 v3 v4 ) ( d ) * ( d e v3 v4 ) ( e d ) * ( d e v4 ) ( v1 e d ) * ( d e ) ( v1 v2 e d ) * ( d ) ( v1 v2 v3 e d ) * * where d and e denote diagonal and off-diagonal elements of T, and vi * denotes an element of the vector defining H(i). subroutine APPLIQ(n,A,tau,S,v,w) implicit none integer n real*8 A(n,n), tau(n-1), S(n,n), v(n), w(n) integer i, j, k real*8 vtw do k= n-2, 1, -1 * print*, k * Redefine the vector v, as sqrt(tau)*v, H = I - v*v' v(k+1) = sqrt(tau(k)) do i=k+2, n v(i) = A(i,k)*v(k+1) enddo * Calculate the lower left corner, dimension is: (n-k) x k * Total operations: 4 (n-k) k do j=1, k w(j) = 0.0D0 do i=k+1, n w(j) = w(j) + v(i)*S(i,j) enddo enddo do j=1, k do i=k+1, n S(i,j) = S(i,j) - v(i)*w(j) enddo enddo * Calculate the lower right trailing block, only the lower * triangular part, domension (n-k) x (n-k) * Total operations: 4 (n-k)^2 * Operations: 2 (n-k)^2 do i=k+1, n w(i) = 0.0D0 do j=k+1, i w(i) = w(i) + S(i,j)*v(j) enddo do j=i+1, n w(i) = w(i) + S(j,i)*v(j) enddo enddo vtw = 0.0D0 do i=k+1, n vtw = vtw + v(i)*w(i) enddo vtw = 0.5D0*vtw do i=k+1,n w(i) = w(i) - vtw*v(i) enddo * Operations: 4/2 (n-k)^2 = 2(n-k)^2 do j= k+1, n do i=j, n S(i,j) = S(i,j) - w(i)*v(j)-v(i)*w(j) enddo enddo enddo return end