c c Copyright (C) 2002 Vanderbilt group c This file is distributed under the terms of the GNU General Public c License as described in the file 'License' in the current directory. c * ====================================================================== * NIST Guide to Available Math Software. * Fullsource for module RSG from package EISPACK. * Retrieved from NETLIB on Fri Jun 21 11:57:40 1996. * ====================================================================== subroutine rsg(nm,n,a,b,w,matz,z,fv1,fv2,ierr) c integer n,nm,ierr,matz double precision a(nm,n),b(nm,n),w(n),z(nm,n),fv1(n),fv2(n) c c this subroutine calls the recommended sequence of c subroutines from the eigensystem subroutine package (eispack) c to find the eigenvalues and eigenvectors (if desired) c for the real symmetric generalized eigenproblem ax = (lambda)bx. c c on input c c nm must be set to the row dimension of the two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrices a and b. c c a contains a real symmetric matrix. c c b contains a positive definite real symmetric matrix. c c matz is an integer variable set equal to zero if c only eigenvalues are desired. otherwise it is set to c any non-zero integer for both eigenvalues and eigenvectors. c c on output c c w contains the eigenvalues in ascending order. c c z contains the eigenvectors if matz is not zero. c c ierr is an integer output variable set equal to an error c completion code described in the documentation for tqlrat c and tql2. the normal completion code is zero. c c fv1 and fv2 are temporary storage arrays. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c if (n .le. nm) go to 10 ierr = 10 * n go to 50 c 10 call reduc(nm,n,a,b,fv2,ierr) if (ierr .ne. 0) go to 50 if (matz .ne. 0) go to 20 c .......... find eigenvalues only .......... call tred1(nm,n,a,w,fv1,fv2) call tqlrat(n,w,fv2,ierr) go to 50 c .......... find both eigenvalues and eigenvectors .......... 20 call tred2(nm,n,a,w,fv1,z) call tql2(nm,n,w,fv1,z,ierr) if (ierr .ne. 0) go to 50 call rebak(nm,n,b,fv2,n,z) 50 return end double precision function epslon (x) double precision x c c estimate unit roundoff in quantities of size x. c double precision a,b,c,eps c c this program should function properly on all systems c satisfying the following two assumptions, c 1. the base used in representing floating point c numbers is not a power of three. c 2. the quantity a in statement 10 is represented to c the accuracy used in floating point variables c that are stored in memory. c the statement number 10 and the go to 10 are intended to c force optimizing compilers to generate code satisfying c assumption 2. c under these assumptions, it should be true that, c a is not exactly equal to four-thirds, c b has a zero for its last bit or digit, c c is not exactly equal to one, c eps measures the separation of 1.0 from c the next larger floating point number. c the developers of eispack would appreciate being informed c about any systems where these assumptions do not hold. c c this version dated 4/6/83. c a = 4.0d0/3.0d0 10 b = a - 1.0d0 c = b + b + b eps = dabs(c-1.0d0) if (eps .eq. 0.0d0) go to 10 epslon = eps*dabs(x) return end subroutine tqlrat(n,d,e2,ierr) c integer i,j,l,m,n,ii,l1,mml,ierr double precision d(n),e2(n) double precision b,c,f,g,h,p,r,s,t,epslon,pythag c c this subroutine is a translation of the algol procedure tqlrat, c algorithm 464, comm. acm 16, 689(1973) by reinsch. c c this subroutine finds the eigenvalues of a symmetric c tridiagonal matrix by the rational ql method. c c on input c c n is the order of the matrix. c c d contains the diagonal elements of the input matrix. c c e2 contains the squares of the subdiagonal elements of the c input matrix in its last n-1 positions. e2(1) is arbitrary. c c on output c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct and c ordered for indices 1,2,...ierr-1, but may not be c the smallest eigenvalues. c c e2 has been destroyed. c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 if (n .eq. 1) go to 1001 c do 100 i = 2, n 100 e2(i-1) = e2(i) c f = 0.0d0 t = 0.0d0 e2(n) = 0.0d0 c do 290 l = 1, n j = 0 h = dabs(d(l)) + dsqrt(e2(l)) if (t .gt. h) go to 105 t = h b = epslon(t) c = b * b c .......... look for small squared sub-diagonal element .......... 105 do 110 m = l, n if (e2(m) .le. c) go to 120 c .......... e2(n) is always zero, so there is no exit c through the bottom of the loop .......... 110 continue c 120 if (m .eq. l) go to 210 130 if (j .eq. 30) go to 1000 j = j + 1 c .......... form shift .......... l1 = l + 1 s = dsqrt(e2(l)) g = d(l) p = (d(l1) - g) / (2.0d0 * s) r = pythag(p,1.0d0) d(l) = s / (p + dsign(r,p)) h = g - d(l) c do 140 i = l1, n 140 d(i) = d(i) - h c f = f + h c .......... rational ql transformation .......... g = d(m) if (g .eq. 0.0d0) g = b h = g s = 0.0d0 mml = m - l c .......... for i=m-1 step -1 until l do -- .......... do 200 ii = 1, mml i = m - ii p = g * h r = p + e2(i) e2(i+1) = s * r s = e2(i) / r d(i+1) = h + s * (h + d(i)) g = d(i) - e2(i) / g if (g .eq. 0.0d0) g = b h = g * p / r 200 continue c e2(l) = s * g d(l) = h c .......... guard against underflow in convergence test .......... if (h .eq. 0.0d0) go to 210 if (dabs(e2(l)) .le. dabs(c/h)) go to 210 e2(l) = h * e2(l) if (e2(l) .ne. 0.0d0) go to 130 210 p = d(l) + f c .......... order eigenvalues .......... if (l .eq. 1) go to 250 c .......... for i=l step -1 until 2 do -- .......... do 230 ii = 2, l i = l + 2 - ii if (p .ge. d(i-1)) go to 270 d(i) = d(i-1) 230 continue c 250 i = 1 270 d(i) = p 290 continue c go to 1001 c .......... set error -- no convergence to an c eigenvalue after 30 iterations .......... 1000 ierr = l 1001 return end double precision function pythag(a,b) double precision a,b c c finds dsqrt(a**2+b**2) without overflow or destructive underflow c double precision p,r,s,t,u p = dmax1(dabs(a),dabs(b)) if (p .eq. 0.0d0) go to 20 r = (dmin1(dabs(a),dabs(b))/p)**2 10 continue t = 4.0d0 + r if (t .eq. 4.0d0) go to 20 s = r/t u = 1.0d0 + 2.0d0*s p = u*p r = (s/u)**2 * r go to 10 20 pythag = p return end subroutine rebak(nm,n,b,dl,m,z) c integer i,j,k,m,n,i1,ii,nm double precision b(nm,n),dl(n),z(nm,m) double precision x c c this subroutine is a translation of the algol procedure rebaka, c num. math. 11, 99-110(1968) by martin and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 303-314(1971). c c this subroutine forms the eigenvectors of a generalized c symmetric eigensystem by back transforming those of the c derived symmetric matrix determined by reduc. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix system. c c b contains information about the similarity transformation c (cholesky decomposition) used in the reduction by reduc c in its strict lower triangle. c c dl contains further information about the transformation. c c m is the number of eigenvectors to be back transformed. c c z contains the eigenvectors to be back transformed c in its first m columns. c c on output c c z contains the transformed eigenvectors c in its first m columns. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c if (m .eq. 0) go to 200 c do 100 j = 1, m c .......... for i=n step -1 until 1 do -- .......... do 100 ii = 1, n i = n + 1 - ii i1 = i + 1 x = z(i,j) if (i .eq. n) go to 80 c do 60 k = i1, n 60 x = x - b(k,i) * z(k,j) c 80 z(i,j) = x / dl(i) 100 continue c 200 return end subroutine reduc(nm,n,a,b,dl,ierr) c integer i,j,k,n,i1,j1,nm,nn,ierr double precision a(nm,n),b(nm,n),dl(n) double precision x,y c c this subroutine is a translation of the algol procedure reduc1, c num. math. 11, 99-110(1968) by martin and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 303-314(1971). c c this subroutine reduces the generalized symmetric eigenproblem c ax=(lambda)bx, where b is positive definite, to the standard c symmetric eigenproblem using the cholesky factorization of b. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrices a and b. if the cholesky c factor l of b is already available, n should be prefixed c with a minus sign. c c a and b contain the real symmetric input matrices. only the c full upper triangles of the matrices need be supplied. if c n is negative, the strict lower triangle of b contains, c instead, the strict lower triangle of its cholesky factor l. c c dl contains, if n is negative, the diagonal elements of l. c c on output c c a contains in its full lower triangle the full lower triangle c of the symmetric matrix derived from the reduction to the c standard form. the strict upper triangle of a is unaltered. c c b contains in its strict lower triangle the strict lower c triangle of its cholesky factor l. the full upper c triangle of b is unaltered. c c dl contains the diagonal elements of l. c c ierr is set to c zero for normal return, c 7*n+1 if b is not positive definite. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 nn = iabs(n) if (n .lt. 0) go to 100 c .......... form l in the arrays b and dl .......... do 80 i = 1, n i1 = i - 1 c do 80 j = i, n x = b(i,j) if (i .eq. 1) go to 40 c do 20 k = 1, i1 20 x = x - b(i,k) * b(j,k) c 40 if (j .ne. i) go to 60 if (x .le. 0.0d0) go to 1000 y = dsqrt(x) dl(i) = y go to 80 60 b(j,i) = x / y 80 continue c .......... form the transpose of the upper triangle of inv(l)*a c in the lower triangle of the array a .......... 100 do 200 i = 1, nn i1 = i - 1 y = dl(i) c do 200 j = i, nn x = a(i,j) if (i .eq. 1) go to 180 c do 160 k = 1, i1 160 x = x - b(i,k) * a(j,k) c 180 a(j,i) = x / y 200 continue c .......... pre-multiply by inv(l) and overwrite .......... do 300 j = 1, nn j1 = j - 1 c do 300 i = j, nn x = a(i,j) if (i .eq. j) go to 240 i1 = i - 1 c do 220 k = j, i1 220 x = x - a(k,j) * b(i,k) c 240 if (j .eq. 1) go to 280 c do 260 k = 1, j1 260 x = x - a(j,k) * b(i,k) c 280 a(i,j) = x / dl(i) 300 continue c go to 1001 c .......... set error -- b is not positive definite .......... 1000 ierr = 7 * n + 1 1001 return end subroutine tql2(nm,n,d,e,z,ierr) c integer i,j,k,l,m,n,ii,l1,l2,nm,mml,ierr double precision d(n),e(n),z(nm,n) double precision c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2,pythag c c this subroutine is a translation of the algol procedure tql2, c num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and c wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). c c this subroutine finds the eigenvalues and eigenvectors c of a symmetric tridiagonal matrix by the ql method. c the eigenvectors of a full symmetric matrix can also c be found if tred2 has been used to reduce this c full matrix to tridiagonal form. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c d contains the diagonal elements of the input matrix. c c e contains the subdiagonal elements of the input matrix c in its last n-1 positions. e(1) is arbitrary. c c z contains the transformation matrix produced in the c reduction by tred2, if performed. if the eigenvectors c of the tridiagonal matrix are desired, z must contain c the identity matrix. c c on output c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct but c unordered for indices 1,2,...,ierr-1. c c e has been destroyed. c c z contains orthonormal eigenvectors of the symmetric c tridiagonal (or full) matrix. if an error exit is made, c z contains the eigenvectors associated with the stored c eigenvalues. c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 if (n .eq. 1) go to 1001 c do 100 i = 2, n 100 e(i-1) = e(i) c f = 0.0d0 tst1 = 0.0d0 e(n) = 0.0d0 c do 240 l = 1, n j = 0 h = dabs(d(l)) + dabs(e(l)) if (tst1 .lt. h) tst1 = h c .......... look for small sub-diagonal element .......... do 110 m = l, n tst2 = tst1 + dabs(e(m)) if (tst2 .eq. tst1) go to 120 c .......... e(n) is always zero, so there is no exit c through the bottom of the loop .......... 110 continue c 120 if (m .eq. l) go to 220 130 if (j .eq. 30) go to 1000 j = j + 1 c .......... form shift .......... l1 = l + 1 l2 = l1 + 1 g = d(l) p = (d(l1) - g) / (2.0d0 * e(l)) r = pythag(p,1.0d0) d(l) = e(l) / (p + dsign(r,p)) d(l1) = e(l) * (p + dsign(r,p)) dl1 = d(l1) h = g - d(l) if (l2 .gt. n) go to 145 c do 140 i = l2, n 140 d(i) = d(i) - h c 145 f = f + h c .......... ql transformation .......... p = d(m) c = 1.0d0 c2 = c el1 = e(l1) s = 0.0d0 mml = m - l c .......... for i=m-1 step -1 until l do -- .......... do 200 ii = 1, mml c3 = c2 c2 = c s2 = s i = m - ii g = c * e(i) h = c * p r = pythag(p,e(i)) e(i+1) = s * r s = e(i) / r c = p / r p = c * d(i) - s * g d(i+1) = h + s * (c * g + s * d(i)) c .......... form vector .......... do 180 k = 1, n h = z(k,i+1) z(k,i+1) = s * z(k,i) + c * h z(k,i) = c * z(k,i) - s * h 180 continue c 200 continue c p = -s * s2 * c3 * el1 * e(l) / dl1 e(l) = s * p d(l) = c * p tst2 = tst1 + dabs(e(l)) if (tst2 .gt. tst1) go to 130 220 d(l) = d(l) + f 240 continue c .......... order eigenvalues and eigenvectors .......... do 300 ii = 2, n i = ii - 1 k = i p = d(i) c do 260 j = ii, n if (d(j) .ge. p) go to 260 k = j p = d(j) 260 continue c if (k .eq. i) go to 300 d(k) = d(i) d(i) = p c do 280 j = 1, n p = z(j,i) z(j,i) = z(j,k) z(j,k) = p 280 continue c 300 continue c go to 1001 c .......... set error -- no convergence to an c eigenvalue after 30 iterations .......... 1000 ierr = l 1001 return end subroutine tred1(nm,n,a,d,e,e2) c integer i,j,k,l,n,ii,nm,jp1 double precision a(nm,n),d(n),e(n),e2(n) double precision f,g,h,scale c c this subroutine is a translation of the algol procedure tred1, c num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). c c this subroutine reduces a real symmetric matrix c to a symmetric tridiagonal matrix using c orthogonal similarity transformations. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c a contains the real symmetric input matrix. only the c lower triangle of the matrix need be supplied. c c on output c c a contains information about the orthogonal trans- c formations used in the reduction in its strict lower c triangle. the full upper triangle of a is unaltered. c c d contains the diagonal elements of the tridiagonal matrix. c c e contains the subdiagonal elements of the tridiagonal c matrix in its last n-1 positions. e(1) is set to zero. c c e2 contains the squares of the corresponding elements of e. c e2 may coincide with e if the squares are not needed. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c do 100 i = 1, n d(i) = a(n,i) a(n,i) = a(i,i) 100 continue c .......... for i=n step -1 until 1 do -- .......... do 300 ii = 1, n i = n + 1 - ii l = i - 1 h = 0.0d0 scale = 0.0d0 if (l .lt. 1) go to 130 c .......... scale row (algol tol then not needed) .......... do 120 k = 1, l 120 scale = scale + dabs(d(k)) c if (scale .ne. 0.0d0) go to 140 c do 125 j = 1, l d(j) = a(l,j) a(l,j) = a(i,j) a(i,j) = 0.0d0 125 continue c 130 e(i) = 0.0d0 e2(i) = 0.0d0 go to 300 c 140 do 150 k = 1, l d(k) = d(k) / scale h = h + d(k) * d(k) 150 continue c e2(i) = scale * scale * h f = d(l) g = -dsign(dsqrt(h),f) e(i) = scale * g h = h - f * g d(l) = f - g if (l .eq. 1) go to 285 c .......... form a*u .......... do 170 j = 1, l 170 e(j) = 0.0d0 c do 240 j = 1, l f = d(j) g = e(j) + a(j,j) * f jp1 = j + 1 if (l .lt. jp1) go to 220 c do 200 k = jp1, l g = g + a(k,j) * d(k) e(k) = e(k) + a(k,j) * f 200 continue c 220 e(j) = g 240 continue c .......... form p .......... f = 0.0d0 c do 245 j = 1, l e(j) = e(j) / h f = f + e(j) * d(j) 245 continue c h = f / (h + h) c .......... form q .......... do 250 j = 1, l 250 e(j) = e(j) - h * d(j) c .......... form reduced a .......... do 280 j = 1, l f = d(j) g = e(j) c do 260 k = j, l 260 a(k,j) = a(k,j) - f * e(k) - g * d(k) c 280 continue c 285 do 290 j = 1, l f = d(j) d(j) = a(l,j) a(l,j) = a(i,j) a(i,j) = f * scale 290 continue c 300 continue c return end subroutine tred2(nm,n,a,d,e,z) c integer i,j,k,l,n,ii,nm,jp1 double precision a(nm,n),d(n),e(n),z(nm,n) double precision f,g,h,hh,scale c c this subroutine is a translation of the algol procedure tred2, c num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). c c this subroutine reduces a real symmetric matrix to a c symmetric tridiagonal matrix using and accumulating c orthogonal similarity transformations. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c a contains the real symmetric input matrix. only the c lower triangle of the matrix need be supplied. c c on output c c d contains the diagonal elements of the tridiagonal matrix. c c e contains the subdiagonal elements of the tridiagonal c matrix in its last n-1 positions. e(1) is set to zero. c c z contains the orthogonal transformation matrix c produced in the reduction. c c a and z may coincide. if distinct, a is unaltered. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c do 100 i = 1, n c do 80 j = i, n 80 z(j,i) = a(j,i) c d(i) = a(n,i) 100 continue c if (n .eq. 1) go to 510 c .......... for i=n step -1 until 2 do -- .......... do 300 ii = 2, n i = n + 2 - ii l = i - 1 h = 0.0d0 scale = 0.0d0 if (l .lt. 2) go to 130 c .......... scale row (algol tol then not needed) .......... do 120 k = 1, l 120 scale = scale + dabs(d(k)) c if (scale .ne. 0.0d0) go to 140 130 e(i) = d(l) c do 135 j = 1, l d(j) = z(l,j) z(i,j) = 0.0d0 z(j,i) = 0.0d0 135 continue c go to 290 c 140 do 150 k = 1, l d(k) = d(k) / scale h = h + d(k) * d(k) 150 continue c f = d(l) g = -dsign(dsqrt(h),f) e(i) = scale * g h = h - f * g d(l) = f - g c .......... form a*u .......... do 170 j = 1, l 170 e(j) = 0.0d0 c do 240 j = 1, l f = d(j) z(j,i) = f g = e(j) + z(j,j) * f jp1 = j + 1 if (l .lt. jp1) go to 220 c do 200 k = jp1, l g = g + z(k,j) * d(k) e(k) = e(k) + z(k,j) * f 200 continue c 220 e(j) = g 240 continue c .......... form p .......... f = 0.0d0 c do 245 j = 1, l e(j) = e(j) / h f = f + e(j) * d(j) 245 continue c hh = f / (h + h) c .......... form q .......... do 250 j = 1, l 250 e(j) = e(j) - hh * d(j) c .......... form reduced a .......... do 280 j = 1, l f = d(j) g = e(j) c do 260 k = j, l 260 z(k,j) = z(k,j) - f * e(k) - g * d(k) c d(j) = z(l,j) z(i,j) = 0.0d0 280 continue c 290 d(i) = h 300 continue c .......... accumulation of transformation matrices .......... do 500 i = 2, n l = i - 1 z(n,l) = z(l,l) z(l,l) = 1.0d0 h = d(i) if (h .eq. 0.0d0) go to 380 c do 330 k = 1, l 330 d(k) = z(k,i) / h c do 360 j = 1, l g = 0.0d0 c do 340 k = 1, l 340 g = g + z(k,i) * z(k,j) c do 360 k = 1, l z(k,j) = z(k,j) - g * d(k) 360 continue c 380 do 400 k = 1, l 400 z(k,i) = 0.0d0 c 500 continue c 510 do 520 i = 1, n d(i) = z(n,i) z(n,i) = 0.0d0 520 continue c z(n,n) = 1.0d0 e(1) = 0.0d0 return end