parameter (n=100) real A(n,n), S(n,n),w(8*n), tiny,aln(2*n),bln(2*n) integer iw(5*n), nn * This tiny is now the relative error for ||ln(A)||. tiny = 1.0E-5 * What is the matrix, choose nn for different examples. nn = 3 call example(n,A,nn) call PADELN(n,A,S,w,iw,tiny,aln,bln) end ********************* subroutine example(n,A,nn) implicit none integer n, nn, i,j, k, p, seed real A(n,n), xp, B(500,500), rand if ( nn .EQ. 1) then do j=1, n do i=j, n A(i,j) = 1.0/(2.0 + (i-j)**2) enddo enddo elseif ( nn .EQ. 2) then xp = sqrt(real(n)) p = xp if ( p*p .NE. n) then print*, 'Is n the square of p? n=', n, ' p=', p print*, 'You should stop and try a different n' endif do j=1, n A(j,j)=4.0 do i=j+1, n A(i,j)=0.0 enddo enddo do k=1, p do i=1, p-1 j = (k-1)*p + i A(j+1,j)=-1.0 enddo enddo do j=1, n-p A(j+p,j)=-1.0 enddo elseif (nn .EQ. 3) then c c The random matrix c if ( n .GT. 500) print*, 'Wrong in Example.f' print*, 'Seed = ' read*, seed B(1,1) = rand(seed) do j=1, n do i =1 , n B(i,j) = rand(0)-0.5 enddo enddo do j=1, n do i=j,n A(i,j) = 0.0 do k=1, n A(i,j) = A(i,j) + B(i,k)*B(j,k) enddo c A(j,i) = A(i,j) enddo enddo endif return end ********************************************************* * This is the program for calculating ln(A) * where A is symmetric positive definite. * A is the input matrix, * S is the output. * w, iw, aln, bln are working space vectors. * tiny is an error control parameter. * * See the following paper for details: * * Ya Yan Lu, Computing the Logarithm of a Symmetric Positive * Definite Matrix, Applied Numerical Mathematics, * 26(4), 483-496, 1998. * * This program requires LAPACK to work. ******************************************************** subroutine PADELN(n,A,S,w,iw,tiny,aln,bln) implicit none integer n, iw(5*n) real A(n,n), S(n,n), w(8*n),tiny,aln(2*n), bln(2*n) call PDLNW(n,A,S,w,w(n+1),w(2*n+1),w(3*n+1),w(4*n+1),iw,tiny, $ aln, bln) return end ***************************************************************** subroutine PDLNW(n,A,S,d,e,tau,w,work,iw,tiny,aln,bln) implicit none integer n real A(n,n),S(n,n),d(n),e(n),tau(n),w(n),work(4*n),tiny real aln(2*n), bln(2*n) integer iw(5*n) real lams, laml, mu integer lwork, m, info, j, mfound, nsplit lwork = 4*n * Reduce the matrix A to tridiagonal form by LAPACK * routine xSYTRD. call SSYTRD( 'L', n, A, n, d, e, tau, work, lwork, info ) * Find the largest and smallest eigenvalues by LAPACK * routine SSTEBZ. call SSTEBZ( 'I', 'B', n, 0.0, 0.0, 1, 1, 0.0, d, e, $ mfound, nsplit, w, iw, iw(n+1), work, iw(2*n+1), $ info ) lams = w(1) call SSTEBZ( 'I', 'B', n, 0.0, 0.0, n, n, 0.0, d, e, $ mfound, nsplit, w, iw, iw(n+1), work, iw(2*n+1), $ info ) laml = w(1) print*, 'Eigenvalues and Condition Number:' print*, lams, laml, laml/lams * Now, we determine the number of terms needed in the * Pade approximation -- m, The scaling parameter mu * is simply sqrt(lambda(1)*lambda(n)). print*, 'The relative tiny is: ', tiny call XINGLN(laml, lams, tiny, mu, m) if ( m .GT. 2*n ) print*, 'Stop!!! m is too large!' print*, 'The absolute tiny is: ', tiny print*, 'm and mu= ', m, mu * It takes: O(m n^2) call TRILN(n,d,e,m,aln,bln,S,work,work(n+1),work(2*n+1),mu) * Given symmetric matrix S, we calculate QSQ^t * Q is given as in the SSYTRD, stored in elementary reflectors. call APPLIQ(n,A,tau,S, work, work(n+1)) do j=1, n S(j,j) = S(j,j) + log(mu) enddo return end ********************************************************* * This is the program that calculate the minimum m * such that the continued fraction of 2m stairs * will produce a satisfactory error. subroutine XINGLN(laml, lams, tiny,mu, m) implicit none real laml, lams, tiny, mu integer m real x, c, a0, a1, b0, b1, err ***************************** * This a few lines are added on March 12, 1997, to change the * input relative error tolerance to the absolute error tolerance * || ln(A)|| = max ( | ln (lambda_1)|, |lambda_n|) * The value of tiny is destroyed! ********** if ( abs(log(laml)) .GT. abs(log(lams)) ) then tiny = tiny*abs(log(laml)) else tiny = tiny*abs(log(lams)) endif ************************** mu = sqrt(laml*lams) x = sqrt(laml/lams)-1.0 c = 1.0 a0 = 0.0 a1 = c*x b0 = 1.0 b1 = 1.0 err = 1.0 m = 0 1 if ( err .GT. tiny) then m = m + 1 c = m/(2.0*(2.0*m-1.0)) b0 = b1 + c*x*b0 a0 = a1 + c*x*a0 err = log(1.0+x) - a0/b0 c print*, 'm=', m, ' err=', err c = m/(2.0*(2.0*m+1.0)) b1 = b0 + c*x*b1 a1 = a0 + c*x*a1 c c re-scale, I do not trust the results here. c b1 = b1/b0 a0 = a0/b0 a1 = a1/b0 b0 = 1.0 goto 1 endif return end ******************************************************* * This is the program to calculate * S = SUM a_j ( I + b_j X)^(-1) X * where a_j, b_j are the coefficients for ln(1+x) * The matrix T (T= mu (I+X) )is given by diagonal d, off-diagonal e. ******************************************************* subroutine TRILN(n,d,e,m,a,b,S,dx,ex,x,mu) implicit none integer n, m real d(n),e(n-1),S(n,n),dx(n),ex(n-1),x(n), mu real a(m), b(m), left, right, coef integer i,j,k,info left = 0 right = 1 call gauleg(left,right,b,a,m) do j=1, n do i=j, n S(i,j) = 0.0 enddo enddo do k=1, m * Calculate I + b(k) X based on matrix T now! coef = b(k)/mu do i=1, n dx(i) = b(m-k+1) + coef*d(i) enddo do i=1, n-1 ex(i) = coef*e(i) enddo * SPTTRF computes the factorization of a real symmetric positive * definite tridiagonal matrix A. call SPTTRF(n,dx,ex,info) do j=1, n * We now generate the right hand side vector from T directly * The subtraction is a source of round off error? do i=1, n x(i) = 0.0 enddo coef = a(k)/mu x(j) = d(j)*coef - a(k) c x(j) = a(k)*(d(j)/mu - 1.0D0) if ( j .GT. 1) x(j-1) = coef*e(j-1) if ( j .LT. n) x(j+1) = coef*e(j) call TRISPE(n, dx, ex, x, j) do i=j, n S(i,j) = S(i,j) + x(i) enddo 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 A(n,n), tau(n-1), S(n,n), v(n), w(n) integer i, j, k real 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.0 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.0 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.0 do i=k+1, n vtw = vtw + v(i)*w(i) enddo vtw = 0.5*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 *********************************************************** subroutine TRISPE(n, d, e, x, j) implicit none integer n, j real d(n), e(n-1), x(n) integer i if ( j .GT. 1) x(j) = x(j)-e(j-1)*x(j-1) if ( j .LT. n) x(j+1) = x(j+1) - e(j)*x(j) do i = j+2, n x(i) = -e(i-1)*x(i-1) enddo x(n) = x(n)/d(n) do i=n-1, j, -1 x(i) = x(i)/d(i) - e(i)*x(i+1) enddo return end SUBROUTINE gauleg(x1,x2,x,w,n) INTEGER n REAL x1,x2,x(n),w(n) DOUBLE PRECISION EPS PARAMETER (EPS=3.d-14) INTEGER i,j,m DOUBLE PRECISION p1,p2,p3,pp,xl,xm,z,z1 m=(n+1)/2 xm=0.5d0*(x2+x1) xl=0.5d0*(x2-x1) do 12 i=1,m z=cos(3.141592654d0*(i-.25d0)/(n+.5d0)) 1 continue p1=1.d0 p2=0.d0 do 11 j=1,n p3=p2 p2=p1 p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j 11 continue pp=n*(z*p1-p2)/(z*z-1.d0) z1=z z=z1-p1/pp if(abs(z-z1).gt.EPS)goto 1 x(i)=xm-xl*z x(n+1-i)=xm+xl*z w(i)=2.d0*xl/((1.d0-z*z)*pp*pp) w(n+1-i)=w(i) 12 continue return END C (C) Copr. 1986-92 Numerical Recipes Software 9!9=,5.