Mercurial > octave-nkf
changeset 3181:3f6a813eb09e
[project @ 1998-09-25 02:50:29 by jwe]
| author | jwe |
|---|---|
| date | Fri, 25 Sep 1998 02:53:39 +0000 |
| parents | c17387059fd3 |
| children | 936b9ae8f7d2 |
| files | libcruft/ChangeLog libcruft/Makefile.in libcruft/balgen/Makefile.in libcruft/balgen/balgen.f libcruft/balgen/gradeq.f libcruft/balgen/reduce.f libcruft/balgen/scaleg.f libcruft/eispack/Makefile.in libcruft/eispack/epslon.f libcruft/eispack/qzhes.f libcruft/eispack/qzit.f libcruft/eispack/qzval.f src/DLD-FUNCTIONS/balance.cc src/DLD-FUNCTIONS/qzval.cc |
| diffstat | 14 files changed, 182 insertions(+), 1849 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog Thu Sep 24 19:00:19 1998 +0000 +++ b/libcruft/ChangeLog Fri Sep 25 02:53:39 1998 +0000 @@ -1,5 +1,8 @@ Thu Sep 24 11:59:02 1998 John W. Eaton <jwe@bevo.che.wisc.edu> + * balgen, eispack: Delete directories and unnecesary files. + * Makefile.in (CRUFT_DIRS): Delete eispack and balgen from the list. + * lapack/xdlamch.f: New file. * ordered-qz: New directory.
--- a/libcruft/Makefile.in Thu Sep 24 19:00:19 1998 +0000 +++ b/libcruft/Makefile.in Fri Sep 25 02:53:39 1998 +0000 @@ -24,9 +24,8 @@ # generate a new configure script in the top-level directory (edit # configure.in and run autoconf). -CRUFT_DIRS = balgen blas dassl eispack fftpack lapack linpack \ - minpack misc odepack ordered-qz quadpack ranlib slatec-err \ - slatec-fn specfun villad +CRUFT_DIRS = blas dassl fftpack lapack linpack minpack misc odepack \ + ordered-qz quadpack ranlib slatec-err slatec-fn specfun villad SUBDIRS = $(CRUFT_DIRS)
--- a/libcruft/balgen/Makefile.in Thu Sep 24 19:00:19 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,19 +0,0 @@ -# -# Makefile for octave's libcruft/balgen directory -# -# John W. Eaton -# jwe@bevo.che.wisc.edu -# University of Wisconsin-Madison -# Department of Chemical Engineering - -TOPDIR = ../.. - -srcdir = @srcdir@ -top_srcdir = @top_srcdir@ -VPATH = @srcdir@ - -EXTERNAL_DISTFILES = $(DISTFILES) - -include $(TOPDIR)/Makeconf - -include ../Makerules
--- a/libcruft/balgen/balgen.f Thu Sep 24 19:00:19 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,102 +0,0 @@ - subroutine balgen (n,ma,a,mb,b,low,igh,cscale,cperm,wk) -c -c *****parameters: - integer igh,low,ma,mb,n - double precision a(ma,n),b(mb,n),cperm(n),cscale(n),wk(n,6) -c -c *****local variables: -c none -c -c *****functions: -c none -c -c *****subroutines called: -c reduce, scaleg, gradeq -c -c --------------------------------------------------------------- -c -c *****purpose: -c this subroutine balances the matrices a and b to improve the -c accuracy of computing the eigensystem of the generalized -c eigenproblem a*x = (lambda)*b*x. the algorithm is specifically -c designed to precede qz type algorithms, but improved performance -c is expected from most eigensystem solvers. -c ref.: ward, r. c., balancing the generalized eigenvalue -c problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981, -c 141-152. -c -c *****parameter description: -c -c on input: -c -c ma,mb integer -c row dimensions of the arrays containing matrices -c a and b respectively, as declared in the main calling -c program dimension statement; -c -c n integer -c order of the matrices a and b; -c -c a real(ma,n) -c contains the a matrix of the generalized eigenproblem -c defined above; -c -c b real(mb,n) -c contains the b matrix of the generalized eigenproblem -c defined above; -c -c wk real(n,6) -c work array that must contain at least 6*n storage -c locations. wk is altered by this subroutine. -c -c on output: -c -c a,b contain the balanced a and b matrices; -c -c low integer -c beginning -1 of the submatrices of a and b -c containing the non-isolated eigenvalues; -c -c igh integer -c ending -1 of the submatrices of a and b -c containing the non-isolated eigenvalues. if -c igh = 1 (low = 1 also), the a and b matrices have -c been permuted into upper triangular form and have -c not been balanced; -c -c cscale real(n) -c contains the exponents of the column scaling factors -c in its low through igh locations and the reducing -c column permutations in its first low-1 and its -c igh+1 through n locations; -c -c cperm real(n) -c contains the column permutations applied in grading -c the a and b submatrices in its low through igh -c locations; -c -c wk contains the exponents of the row scaling factors -c in its low through igh locations, the reducing row -c permutations in its first low-1 and its igh+1 -c through n locations, and the row permutations -c applied in grading the a and b submatrices in its -c n+low through n+igh locations. -c -c *****algorithm notes: -c none -c -c *****history: -c written by r. c. ward....... -c -c --------------------------------------------------------------- -c - call reduce (n,ma,a,mb,b,low,igh,cscale,wk) - if (low .eq. igh) go to 10 - call scaleg (n,ma,a,mb,b,low,igh,cscale,cperm,wk) - call gradeq (n,ma,a,mb,b,low,igh,cperm,wk(1,2)) - 10 continue - return -c -c last line of balgen -c - end
--- a/libcruft/balgen/gradeq.f Thu Sep 24 19:00:19 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,177 +0,0 @@ - subroutine gradeq (n,ma,a,mb,b,low,igh,cperm,wk) -c -c *****parameters: - integer igh,low,ma,mb,n - double precision a(ma,n),b(mb,n),cperm(n),wk(n,2) -c -c *****local variables: - integer i,ighm1,im,ip1,j,jm,jp1,k - double precision cmax,rmax,suma,sumb,temp -c -c *****fortran functions: - double precision dabs -c -c *****subroutines called: -c none -c -c --------------------------------------------------------------- -c -c *****purpose: -c this subroutine grades the submatrices of a and b given by -c starting -1 low and ending -1 igh in the generalized -c eigenvalue problem a*x = (lambda)*b*x by permuting rows and -c columns such that the norm of the i-th row (column) of the -c a submatrix divided by the norm of the i-th row (column) of -c the b submatrix becomes smaller as i increases. -c ref.: ward, r. c., balancing the generalized eigenvalue -c problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981, -c 141-152. -c -c *****parameter description: -c -c on input: -c -c ma,mb integer -c row dimensions of the arrays containing matrices -c a and b respectively, as declared in the main calling -c program dimension statement; -c -c n integer -c order of the matrices a and b; -c -c a real(ma,n) -c contains the a matrix of the generalized eigenproblem -c defined above; -c -c b real(mb,n) -c contains the b matrix of the generalized eigenproblem -c defined above; -c -c low integer -c specifies the beginning -1 for the rows and -c columns of a and b to be graded; -c -c igh integer -c specifies the ending -1 for the rows and columns -c of a and b to be graded; -c -c wk real(n,2) -c work array that must contain at least 2*n locations. -c only locations low through igh and n+low through -c n+igh are referenced by this subroutine. -c -c on output: -c -c a,b contain the permuted and graded a and b matrices; -c -c cperm real(n) -c contains in its low through igh locations the -c column permutations applied in grading the -c submatrices. the other locations are not referenced -c by this subroutine; -c -c wk contains in its low through igh locations the row -c permutations applied in grading the submatrices. -c -c *****algorithm notes: -c none. -c -c *****history: -c written by r. c. ward....... -c -c --------------------------------------------------------------- -c - if (low .eq. igh) go to 510 - ighm1 = igh-1 -c -c compute column norms of a / those of b -c - do 420 j = low,igh - suma = 0.0d0 - sumb = 0.0d0 - do 410 i = low,igh - suma = suma + dabs(a(i,j)) - sumb = sumb + dabs(b(i,j)) - 410 continue - if (sumb .eq. 0.0d0) go to 415 - wk(j,2) = suma / sumb - go to 420 - 415 continue - wk(j,2) = 1.0d38 - 420 continue -c -c permute columns to order them by decreasing quotients -c - do 450 j = low,ighm1 - cmax = wk(j,2) - jm = j - jp1 = j+1 - do 430 k = jp1,igh - if (cmax .ge. wk(k,2)) go to 430 - jm = k - cmax = wk(k,2) - 430 continue - cperm(j) = jm - if (jm .eq. j) go to 450 - temp = wk(j,2) - wk(j,2) = wk(jm,2) - wk(jm,2) = temp - do 440 i = 1,igh - temp = b(i,j) - b(i,j) = b(i,jm) - b(i,jm) = temp - temp = a(i,j) - a(i,j) = a(i,jm) - a(i,jm) = temp - 440 continue - 450 continue - cperm(igh) = igh -c -c compute row norms of a / those of b -c - do 470 i = low,igh - suma = 0.0d0 - sumb = 0.0d0 - do 460 j = low,igh - suma = suma + dabs(a(i,j)) - sumb = sumb + dabs(b(i,j)) - 460 continue - if (sumb .eq. 0.0d0) go to 465 - wk(i,2) = suma / sumb - go to 470 - 465 continue - wk(i,2) = 1.0d38 -c -c permute rows to order them by decreasing quotients -c - 470 continue - do 500 i = low,ighm1 - rmax = wk(i,2) - im = i - ip1 = i+1 - do 480 k = ip1,igh - if (rmax .ge. wk(k,2)) go to 480 - im = k - rmax = wk(k,2) - 480 continue - wk(i,1) = im - if (im .eq. i) go to 500 - temp = wk(i,2) - wk(i,2) = wk(im,2) - wk(im,2) = temp - do 490 j = low,n - temp = b(i,j) - b(i,j) = b(im,j) - b(im,j) = temp - temp = a(i,j) - a(i,j) = a(im,j) - a(im,j) = temp - 490 continue - 500 continue - wk(igh,1) = igh - 510 continue - return -c -c last line of gradeq -c - end
--- a/libcruft/balgen/reduce.f Thu Sep 24 19:00:19 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,181 +0,0 @@ - subroutine reduce (n,ma,a,mb,b,low,igh,cscale,wk) -c -c *****parameters: - integer igh,low,ma,mb,n - double precision a(ma,n),b(mb,n),cscale(n),wk(n) -c -c *****local variables: - integer i,iflow,ii,ip1,is,j,jp1,k,l,lm1,m - double precision f -c -c *****functions: -c none -c -c *****subroutines called: -c none -c -c --------------------------------------------------------------- -c -c *****purpose: -c this subroutine reduces, if possible, the order of the -c generalized eigenvalue problem a*x = (lambda)*b*x by permuting -c the rows and columns of a and b so that they each have the -c form -c u x y -c 0 c z -c 0 0 r -c -c where u and r are upper triangular and c, x, y, and z are -c arbitrary. thus, the isolated eigenvalues corresponding to -c the triangular matrices are obtained by a division, leaving -c only eigenvalues corresponding to the center matrices to be -c computed. -c ref.: ward, r. c., balancing the generalized eigenvalue -c problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981, -c 141-152. -c -c *****parameter description: -c -c on input: -c -c ma,mb integer -c row dimensions of the arrays containing matrices -c a and b respectively, as declared in the main calling -c program dimension statement; -c -c n integer -c order of the matrices a and b; -c -c a real(ma,n) -c contains the a matrix of the generalized eigenproblem -c defined above; -c -c b real(mb,n) -c contains the b matrix of the generalized eigenproblem -c defined above. -c -c on output: -c -c a,b contain the permuted a and b matrices; -c -c low integer -c beginning -1 of the submatrices of a and b -c containing the non-isolated eigenvalues; -c -c igh integer -c ending -1 of the submatrices of a and b -c containing the non-isolated eigenvalues. if -c igh = 1 (low = 1 also), the permuted a and b -c matrices are upper triangular; -c -c cscale real(n) -c contains the required column permutations in its -c first low-1 and its igh+1 through n locations; -c -c wk real(n) -c contains the required row permutations in its first -c low-1 and its igh+1 through n locations. -c -c *****algorithm notes: -c none -c -c *****history: -c written by r. c. ward....... -c -c --------------------------------------------------------------- -c - k = 1 - l = n - go to 20 -c -c find row with one nonzero in columns 1 through l -c - 10 continue - l = lm1 - if (l .ne. 1) go to 20 - wk(1) = 1 - cscale(1) = 1 - go to 200 - 20 continue - lm1 = l-1 - do 70 ii = 1,l - i = l+1-ii - do 30 j = 1,lm1 - jp1 = j+1 - if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 40 - 30 continue - j = l - go to 60 - 40 continue - do 50 j = jp1,l - if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 70 - 50 continue - j = jp1-1 - 60 continue - m = l - iflow = 1 - go to 150 - 70 continue - go to 90 -c -c find column with one nonzero in rows k through n -c - 80 continue - k = k+1 - 90 continue - do 140 j = k,l - do 100 i = k,lm1 - ip1 = i+1 - if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 110 - 100 continue - i = l - go to 130 - 110 continue - do 120 i = ip1,l - if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 140 - 120 continue - i = ip1-1 - 130 continue - m = k - iflow = 2 - go to 150 - 140 continue - go to 200 -c -c permute rows m and i -c - 150 continue - wk(m) = i - if (i .eq. m) go to 170 - do 160 is = k,n - f = a(i,is) - a(i,is) = a(m,is) - a(m,is) = f - f = b(i,is) - b(i,is) = b(m,is) - b(m,is) = f - 160 continue -c -c permute columns m and j -c - 170 continue - cscale(m) = j - if (j .eq. m) go to 190 - do 180 is = 1,l - f = a(is,j) - a(is,j) = a(is,m) - a(is,m) = f - f = b(is,j) - b(is,j) = b(is,m) - b(is,m) = f - 180 continue - 190 continue - go to (10,80), iflow - 200 continue - low = k - igh = l - return -c -c last line of reduce -c - end
--- a/libcruft/balgen/scaleg.f Thu Sep 24 19:00:19 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,236 +0,0 @@ - subroutine scaleg (n,ma,a,mb,b,low,igh,cscale,cperm,wk) -c -c *****parameters: - integer igh,low,ma,mb,n - double precision a(ma,n),b(mb,n),cperm(n),cscale(n),wk(n,6) -c -c *****local variables: - integer i,ir,it,j,jc,kount,nr,nrp2 - double precision alpha,basl,beta,cmax,coef,coef2,coef5,cor, - * ew,ewc,fi,fj,gamma,pgamma,sum,t,ta,tb,tc -c -c *****fortran functions: - double precision dabs, dlog10, dsign -c float -c -c *****subroutines called: -c none -c -c --------------------------------------------------------------- -c -c *****purpose: -c scales the matrices a and b in the generalized eigenvalue -c problem a*x = (lambda)*b*x such that the magnitudes of the -c elements of the submatrices of a and b (as specified by low -c and igh) are close to unity in the least squares sense. -c ref.: ward, r. c., balancing the generalized eigenvalue -c problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981, -c 141-152. -c -c *****parameter description: -c -c on input: -c -c ma,mb integer -c row dimensions of the arrays containing matrices -c a and b respectively, as declared in the main calling -c program dimension statement; -c -c n integer -c order of the matrices a and b; -c -c a real(ma,n) -c contains the a matrix of the generalized eigenproblem -c defined above; -c -c b real(mb,n) -c contains the b matrix of the generalized eigenproblem -c defined above; -c -c low integer -c specifies the beginning -1 for the rows and -c columns of a and b to be scaled; -c -c igh integer -c specifies the ending -1 for the rows and columns -c of a and b to be scaled; -c -c cperm real(n) -c work array. only locations low through igh are -c referenced and altered by this subroutine; -c -c wk real(n,6) -c work array that must contain at least 6*n locations. -c only locations low through igh, n+low through n+igh, -c ..., 5*n+low through 5*n+igh are referenced and -c altered by this subroutine. -c -c on output: -c -c a,b contain the scaled a and b matrices; -c -c cscale real(n) -c contains in its low through igh locations the integer -c exponents of 2 used for the column scaling factors. -c the other locations are not referenced; -c -c wk contains in its low through igh locations the integer -c exponents of 2 used for the row scaling factors. -c -c *****algorithm notes: -c none. -c -c *****history: -c written by r. c. ward....... -c modified 8/86 by bobby bodenheimer so that if -c sum = 0 (corresponding to the case where the matrix -c doesn't need to be scaled) the routine returns. -c -c --------------------------------------------------------------- -c - if (low .eq. igh) go to 410 - do 210 i = low,igh - wk(i,1) = 0.0d0 - wk(i,2) = 0.0d0 - wk(i,3) = 0.0d0 - wk(i,4) = 0.0d0 - wk(i,5) = 0.0d0 - wk(i,6) = 0.0d0 - cscale(i) = 0.0d0 - cperm(i) = 0.0d0 - 210 continue -c -c compute right side vector in resulting linear equations -c - basl = dlog10(2.0d0) - do 240 i = low,igh - do 240 j = low,igh - tb = b(i,j) - ta = a(i,j) - if (ta .eq. 0.0d0) go to 220 - ta = dlog10(dabs(ta)) / basl - 220 continue - if (tb .eq. 0.0d0) go to 230 - tb = dlog10(dabs(tb)) / basl - 230 continue - wk(i,5) = wk(i,5) - ta - tb - wk(j,6) = wk(j,6) - ta - tb - 240 continue - nr = igh-low+1 - coef = 1.0d0/float(2*nr) - coef2 = coef*coef - coef5 = 0.5d0*coef2 - nrp2 = nr+2 - beta = 0.0d0 - it = 1 -c -c start generalized conjugate gradient iteration -c - 250 continue - ew = 0.0d0 - ewc = 0.0d0 - gamma = 0.0d0 - do 260 i = low,igh - gamma = gamma + wk(i,5)*wk(i,5) + wk(i,6)*wk(i,6) - ew = ew + wk(i,5) - ewc = ewc + wk(i,6) - 260 continue - gamma = coef*gamma - coef2*(ew**2 + ewc**2) - + - coef5*(ew - ewc)**2 - if (it .ne. 1) beta = gamma / pgamma - t = coef5*(ewc - 3.0d0*ew) - tc = coef5*(ew - 3.0d0*ewc) - do 270 i = low,igh - wk(i,2) = beta*wk(i,2) + coef*wk(i,5) + t - cperm(i) = beta*cperm(i) + coef*wk(i,6) + tc - 270 continue -c -c apply matrix to vector -c - do 300 i = low,igh - kount = 0 - sum = 0.0d0 - do 290 j = low,igh - if (a(i,j) .eq. 0.0d0) go to 280 - kount = kount+1 - sum = sum + cperm(j) - 280 continue - if (b(i,j) .eq. 0.0d0) go to 290 - kount = kount+1 - sum = sum + cperm(j) - 290 continue - wk(i,3) = float(kount)*wk(i,2) + sum - 300 continue - do 330 j = low,igh - kount = 0 - sum = 0.0d0 - do 320 i = low,igh - if (a(i,j) .eq. 0.0d0) go to 310 - kount = kount+1 - sum = sum + wk(i,2) - 310 continue - if (b(i,j) .eq. 0.0d0) go to 320 - kount = kount+1 - sum = sum + wk(i,2) - 320 continue - wk(j,4) = float(kount)*cperm(j) + sum - 330 continue - sum = 0.0d0 - do 340 i = low,igh - sum = sum + wk(i,2)*wk(i,3) + cperm(i)*wk(i,4) - 340 continue - if(sum.eq.0.0d0) return - alpha = gamma / sum -c -c determine correction to current iterate -c - cmax = 0.0d0 - do 350 i = low,igh - cor = alpha * wk(i,2) - if (dabs(cor) .gt. cmax) cmax = dabs(cor) - wk(i,1) = wk(i,1) + cor - cor = alpha * cperm(i) - if (dabs(cor) .gt. cmax) cmax = dabs(cor) - cscale(i) = cscale(i) + cor - 350 continue - if (cmax .lt. 0.5d0) go to 370 - do 360 i = low,igh - wk(i,5) = wk(i,5) - alpha*wk(i,3) - wk(i,6) = wk(i,6) - alpha*wk(i,4) - 360 continue - pgamma = gamma - it = it+1 - if (it .le. nrp2) go to 250 -c -c end generalized conjugate gradient iteration -c - 370 continue - do 380 i = low,igh - ir = wk(i,1) + dsign(0.5d0,wk(i,1)) - wk(i,1) = ir - jc = cscale(i) + dsign(0.5d0,cscale(i)) - cscale(i) = jc - 380 continue -c -c scale a and b -c - do 400 i = 1,igh - ir = wk(i,1) - fi = 2.0d0**ir - if (i .lt. low) fi = 1.0d0 - do 400 j =low,n - jc = cscale(j) - fj = 2.0d0**jc - if (j .le. igh) go to 390 - if (i .lt. low) go to 400 - fj = 1.0d0 - 390 continue - a(i,j) = a(i,j)*fi*fj - b(i,j) = b(i,j)*fi*fj - 400 continue - 410 continue - return -c -c last line of scaleg -c - end
--- a/libcruft/eispack/Makefile.in Thu Sep 24 19:00:19 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,19 +0,0 @@ -# -# Makefile for octave's libcruft/eispack directory -# -# John W. Eaton -# jwe@bevo.che.wisc.edu -# University of Wisconsin-Madison -# Department of Chemical Engineering - -TOPDIR = ../.. - -srcdir = @srcdir@ -top_srcdir = @top_srcdir@ -VPATH = @srcdir@ - -EXTERNAL_DISTFILES = $(DISTFILES) - -include $(TOPDIR)/Makeconf - -include ../Makerules
--- a/libcruft/eispack/epslon.f Thu Sep 24 19:00:19 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,36 +0,0 @@ - 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
--- a/libcruft/eispack/qzhes.f Thu Sep 24 19:00:19 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,195 +0,0 @@ - subroutine qzhes(nm,n,a,b,matz,z) -c - integer i,j,k,l,n,lb,l1,nm,nk1,nm1,nm2 - double precision a(nm,n),b(nm,n),z(nm,n) - double precision r,s,t,u1,u2,v1,v2,rho - logical matz -c -c this subroutine is the first step of the qz algorithm -c for solving generalized matrix eigenvalue problems, -c siam j. numer. anal. 10, 241-256(1973) by moler and stewart. -c -c this subroutine accepts a pair of real general matrices and -c reduces one of them to upper hessenberg form and the other -c to upper triangular form using orthogonal transformations. -c it is usually followed by qzit, qzval and, possibly, qzvec. -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. -c -c a contains a real general matrix. -c -c b contains a real general matrix. -c -c matz should be set to .true. if the right hand transformations -c are to be accumulated for later use in computing -c eigenvectors, and to .false. otherwise. -c -c on output -c -c a has been reduced to upper hessenberg form. the elements -c below the first subdiagonal have been set to zero. -c -c b has been reduced to upper triangular form. the elements -c below the main diagonal have been set to zero. -c -c z contains the product of the right hand transformations if -c matz has been set to .true. otherwise, z is not referenced. -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 -c .......... initialize z .......... - if (.not. matz) go to 10 -c - do 3 j = 1, n -c - do 2 i = 1, n - z(i,j) = 0.0d0 - 2 continue -c - z(j,j) = 1.0d0 - 3 continue -c .......... reduce b to upper triangular form .......... - 10 if (n .le. 1) go to 170 - nm1 = n - 1 -c - do 100 l = 1, nm1 - l1 = l + 1 - s = 0.0d0 -c - do 20 i = l1, n - s = s + dabs(b(i,l)) - 20 continue -c - if (s .eq. 0.0d0) go to 100 - s = s + dabs(b(l,l)) - r = 0.0d0 -c - do 25 i = l, n - b(i,l) = b(i,l) / s - r = r + b(i,l)**2 - 25 continue -c - r = dsign(dsqrt(r),b(l,l)) - b(l,l) = b(l,l) + r - rho = r * b(l,l) -c - do 50 j = l1, n - t = 0.0d0 -c - do 30 i = l, n - t = t + b(i,l) * b(i,j) - 30 continue -c - t = -t / rho -c - do 40 i = l, n - b(i,j) = b(i,j) + t * b(i,l) - 40 continue -c - 50 continue -c - do 80 j = 1, n - t = 0.0d0 -c - do 60 i = l, n - t = t + b(i,l) * a(i,j) - 60 continue -c - t = -t / rho -c - do 70 i = l, n - a(i,j) = a(i,j) + t * b(i,l) - 70 continue -c - 80 continue -c - b(l,l) = -s * r -c - do 90 i = l1, n - b(i,l) = 0.0d0 - 90 continue -c - 100 continue -c .......... reduce a to upper hessenberg form, while -c keeping b triangular .......... - if (n .eq. 2) go to 170 - nm2 = n - 2 -c - do 160 k = 1, nm2 - nk1 = nm1 - k -c .......... for l=n-1 step -1 until k+1 do -- .......... - do 150 lb = 1, nk1 - l = n - lb - l1 = l + 1 -c .......... zero a(l+1,k) .......... - s = dabs(a(l,k)) + dabs(a(l1,k)) - if (s .eq. 0.0d0) go to 150 - u1 = a(l,k) / s - u2 = a(l1,k) / s - r = dsign(dsqrt(u1*u1+u2*u2),u1) - v1 = -(u1 + r) / r - v2 = -u2 / r - u2 = v2 / v1 -c - do 110 j = k, n - t = a(l,j) + u2 * a(l1,j) - a(l,j) = a(l,j) + t * v1 - a(l1,j) = a(l1,j) + t * v2 - 110 continue -c - a(l1,k) = 0.0d0 -c - do 120 j = l, n - t = b(l,j) + u2 * b(l1,j) - b(l,j) = b(l,j) + t * v1 - b(l1,j) = b(l1,j) + t * v2 - 120 continue -c .......... zero b(l+1,l) .......... - s = dabs(b(l1,l1)) + dabs(b(l1,l)) - if (s .eq. 0.0d0) go to 150 - u1 = b(l1,l1) / s - u2 = b(l1,l) / s - r = dsign(dsqrt(u1*u1+u2*u2),u1) - v1 = -(u1 + r) / r - v2 = -u2 / r - u2 = v2 / v1 -c - do 130 i = 1, l1 - t = b(i,l1) + u2 * b(i,l) - b(i,l1) = b(i,l1) + t * v1 - b(i,l) = b(i,l) + t * v2 - 130 continue -c - b(l1,l) = 0.0d0 -c - do 140 i = 1, n - t = a(i,l1) + u2 * a(i,l) - a(i,l1) = a(i,l1) + t * v1 - a(i,l) = a(i,l) + t * v2 - 140 continue -c - if (.not. matz) go to 150 -c - do 145 i = 1, n - t = z(i,l1) + u2 * z(i,l) - z(i,l1) = z(i,l1) + t * v1 - z(i,l) = z(i,l) + t * v2 - 145 continue -c - 150 continue -c - 160 continue -c - 170 return - end
--- a/libcruft/eispack/qzit.f Thu Sep 24 19:00:19 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,354 +0,0 @@ - subroutine qzit(nm,n,a,b,eps1,matz,z,ierr) -c - integer i,j,k,l,n,en,k1,k2,ld,ll,l1,na,nm,ish,itn,its,km1,lm1, - x enm2,ierr,lor1,enorn - double precision a(nm,n),b(nm,n),z(nm,n) - double precision r,s,t,a1,a2,a3,ep,sh,u1,u2,u3,v1,v2,v3,ani,a11, - x a12,a21,a22,a33,a34,a43,a44,bni,b11,b12,b22,b33,b34, - x b44,epsa,epsb,eps1,anorm,bnorm,epslon - logical matz,notlas -c -c this subroutine is the second step of the qz algorithm -c for solving generalized matrix eigenvalue problems, -c siam j. numer. anal. 10, 241-256(1973) by moler and stewart, -c as modified in technical note nasa tn d-7305(1973) by ward. -c -c this subroutine accepts a pair of real matrices, one of them -c in upper hessenberg form and the other in upper triangular form. -c it reduces the hessenberg matrix to quasi-triangular form using -c orthogonal transformations while maintaining the triangular form -c of the other matrix. it is usually preceded by qzhes and -c followed by qzval and, possibly, qzvec. -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. -c -c a contains a real upper hessenberg matrix. -c -c b contains a real upper triangular matrix. -c -c eps1 is a tolerance used to determine negligible elements. -c eps1 = 0.0 (or negative) may be input, in which case an -c element will be neglected only if it is less than roundoff -c error times the norm of its matrix. if the input eps1 is -c positive, then an element will be considered negligible -c if it is less than eps1 times the norm of its matrix. a -c positive value of eps1 may result in faster execution, -c but less accurate results. -c -c matz should be set to .true. if the right hand transformations -c are to be accumulated for later use in computing -c eigenvectors, and to .false. otherwise. -c -c z contains, if matz has been set to .true., the -c transformation matrix produced in the reduction -c by qzhes, if performed, or else the identity matrix. -c if matz has been set to .false., z is not referenced. -c -c on output -c -c a has been reduced to quasi-triangular form. the elements -c below the first subdiagonal are still zero and no two -c consecutive subdiagonal elements are nonzero. -c -c b is still in upper triangular form, although its elements -c have been altered. the location b(n,1) is used to store -c eps1 times the norm of b for later use by qzval and qzvec. -c -c z contains the product of the right hand transformations -c (for both steps) if matz has been set to .true.. -c -c ierr is set to -c zero for normal return, -c j if the limit of 30*n iterations is exhausted -c while the j-th eigenvalue is being sought. -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 -c .......... compute epsa,epsb .......... - anorm = 0.0d0 - bnorm = 0.0d0 -c - do 30 i = 1, n - ani = 0.0d0 - if (i .ne. 1) ani = dabs(a(i,i-1)) - bni = 0.0d0 -c - do 20 j = i, n - ani = ani + dabs(a(i,j)) - bni = bni + dabs(b(i,j)) - 20 continue -c - if (ani .gt. anorm) anorm = ani - if (bni .gt. bnorm) bnorm = bni - 30 continue -c - if (anorm .eq. 0.0d0) anorm = 1.0d0 - if (bnorm .eq. 0.0d0) bnorm = 1.0d0 - ep = eps1 - if (ep .gt. 0.0d0) go to 50 -c .......... use roundoff level if eps1 is zero .......... - ep = epslon(1.0d0) - 50 epsa = ep * anorm - epsb = ep * bnorm -c .......... reduce a to quasi-triangular form, while -c keeping b triangular .......... - lor1 = 1 - enorn = n - en = n - itn = 30*n -c .......... begin qz step .......... - 60 if (en .le. 2) go to 1001 - if (.not. matz) enorn = en - its = 0 - na = en - 1 - enm2 = na - 1 - 70 ish = 2 -c .......... check for convergence or reducibility. -c for l=en step -1 until 1 do -- .......... - do 80 ll = 1, en - lm1 = en - ll - l = lm1 + 1 - if (l .eq. 1) go to 95 - if (dabs(a(l,lm1)) .le. epsa) go to 90 - 80 continue -c - 90 a(l,lm1) = 0.0d0 - if (l .lt. na) go to 95 -c .......... 1-by-1 or 2-by-2 block isolated .......... - en = lm1 - go to 60 -c .......... check for small top of b .......... - 95 ld = l - 100 l1 = l + 1 - b11 = b(l,l) - if (dabs(b11) .gt. epsb) go to 120 - b(l,l) = 0.0d0 - s = dabs(a(l,l)) + dabs(a(l1,l)) - u1 = a(l,l) / s - u2 = a(l1,l) / s - r = dsign(dsqrt(u1*u1+u2*u2),u1) - v1 = -(u1 + r) / r - v2 = -u2 / r - u2 = v2 / v1 -c - do 110 j = l, enorn - t = a(l,j) + u2 * a(l1,j) - a(l,j) = a(l,j) + t * v1 - a(l1,j) = a(l1,j) + t * v2 - t = b(l,j) + u2 * b(l1,j) - b(l,j) = b(l,j) + t * v1 - b(l1,j) = b(l1,j) + t * v2 - 110 continue -c - if (l .ne. 1) a(l,lm1) = -a(l,lm1) - lm1 = l - l = l1 - go to 90 - 120 a11 = a(l,l) / b11 - a21 = a(l1,l) / b11 - if (ish .eq. 1) go to 140 -c .......... iteration strategy .......... - if (itn .eq. 0) go to 1000 - if (its .eq. 10) go to 155 -c .......... determine type of shift .......... - b22 = b(l1,l1) - if (dabs(b22) .lt. epsb) b22 = epsb - b33 = b(na,na) - if (dabs(b33) .lt. epsb) b33 = epsb - b44 = b(en,en) - if (dabs(b44) .lt. epsb) b44 = epsb - a33 = a(na,na) / b33 - a34 = a(na,en) / b44 - a43 = a(en,na) / b33 - a44 = a(en,en) / b44 - b34 = b(na,en) / b44 - t = 0.5d0 * (a43 * b34 - a33 - a44) - r = t * t + a34 * a43 - a33 * a44 - if (r .lt. 0.0d0) go to 150 -c .......... determine single shift zeroth column of a .......... - ish = 1 - r = dsqrt(r) - sh = -t + r - s = -t - r - if (dabs(s-a44) .lt. dabs(sh-a44)) sh = s -c .......... look for two consecutive small -c sub-diagonal elements of a. -c for l=en-2 step -1 until ld do -- .......... - do 130 ll = ld, enm2 - l = enm2 + ld - ll - if (l .eq. ld) go to 140 - lm1 = l - 1 - l1 = l + 1 - t = a(l,l) - if (dabs(b(l,l)) .gt. epsb) t = t - sh * b(l,l) - if (dabs(a(l,lm1)) .le. dabs(t/a(l1,l)) * epsa) go to 100 - 130 continue -c - 140 a1 = a11 - sh - a2 = a21 - if (l .ne. ld) a(l,lm1) = -a(l,lm1) - go to 160 -c .......... determine double shift zeroth column of a .......... - 150 a12 = a(l,l1) / b22 - a22 = a(l1,l1) / b22 - b12 = b(l,l1) / b22 - a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11) - x / a21 + a12 - a11 * b12 - a2 = (a22 - a11) - a21 * b12 - (a33 - a11) - (a44 - a11) - x + a43 * b34 - a3 = a(l1+1,l1) / b22 - go to 160 -c .......... ad hoc shift .......... - 155 a1 = 0.0d0 - a2 = 1.0d0 - a3 = 1.1605d0 - 160 its = its + 1 - itn = itn - 1 - if (.not. matz) lor1 = ld -c .......... main loop .......... - do 260 k = l, na - notlas = k .ne. na .and. ish .eq. 2 - k1 = k + 1 - k2 = k + 2 - km1 = max0(k-1,l) - ll = min0(en,k1+ish) - if (notlas) go to 190 -c .......... zero a(k+1,k-1) .......... - if (k .eq. l) go to 170 - a1 = a(k,km1) - a2 = a(k1,km1) - 170 s = dabs(a1) + dabs(a2) - if (s .eq. 0.0d0) go to 70 - u1 = a1 / s - u2 = a2 / s - r = dsign(dsqrt(u1*u1+u2*u2),u1) - v1 = -(u1 + r) / r - v2 = -u2 / r - u2 = v2 / v1 -c - do 180 j = km1, enorn - t = a(k,j) + u2 * a(k1,j) - a(k,j) = a(k,j) + t * v1 - a(k1,j) = a(k1,j) + t * v2 - t = b(k,j) + u2 * b(k1,j) - b(k,j) = b(k,j) + t * v1 - b(k1,j) = b(k1,j) + t * v2 - 180 continue -c - if (k .ne. l) a(k1,km1) = 0.0d0 - go to 240 -c .......... zero a(k+1,k-1) and a(k+2,k-1) .......... - 190 if (k .eq. l) go to 200 - a1 = a(k,km1) - a2 = a(k1,km1) - a3 = a(k2,km1) - 200 s = dabs(a1) + dabs(a2) + dabs(a3) - if (s .eq. 0.0d0) go to 260 - u1 = a1 / s - u2 = a2 / s - u3 = a3 / s - r = dsign(dsqrt(u1*u1+u2*u2+u3*u3),u1) - v1 = -(u1 + r) / r - v2 = -u2 / r - v3 = -u3 / r - u2 = v2 / v1 - u3 = v3 / v1 -c - do 210 j = km1, enorn - t = a(k,j) + u2 * a(k1,j) + u3 * a(k2,j) - a(k,j) = a(k,j) + t * v1 - a(k1,j) = a(k1,j) + t * v2 - a(k2,j) = a(k2,j) + t * v3 - t = b(k,j) + u2 * b(k1,j) + u3 * b(k2,j) - b(k,j) = b(k,j) + t * v1 - b(k1,j) = b(k1,j) + t * v2 - b(k2,j) = b(k2,j) + t * v3 - 210 continue -c - if (k .eq. l) go to 220 - a(k1,km1) = 0.0d0 - a(k2,km1) = 0.0d0 -c .......... zero b(k+2,k+1) and b(k+2,k) .......... - 220 s = dabs(b(k2,k2)) + dabs(b(k2,k1)) + dabs(b(k2,k)) - if (s .eq. 0.0d0) go to 240 - u1 = b(k2,k2) / s - u2 = b(k2,k1) / s - u3 = b(k2,k) / s - r = dsign(dsqrt(u1*u1+u2*u2+u3*u3),u1) - v1 = -(u1 + r) / r - v2 = -u2 / r - v3 = -u3 / r - u2 = v2 / v1 - u3 = v3 / v1 -c - do 230 i = lor1, ll - t = a(i,k2) + u2 * a(i,k1) + u3 * a(i,k) - a(i,k2) = a(i,k2) + t * v1 - a(i,k1) = a(i,k1) + t * v2 - a(i,k) = a(i,k) + t * v3 - t = b(i,k2) + u2 * b(i,k1) + u3 * b(i,k) - b(i,k2) = b(i,k2) + t * v1 - b(i,k1) = b(i,k1) + t * v2 - b(i,k) = b(i,k) + t * v3 - 230 continue -c - b(k2,k) = 0.0d0 - b(k2,k1) = 0.0d0 - if (.not. matz) go to 240 -c - do 235 i = 1, n - t = z(i,k2) + u2 * z(i,k1) + u3 * z(i,k) - z(i,k2) = z(i,k2) + t * v1 - z(i,k1) = z(i,k1) + t * v2 - z(i,k) = z(i,k) + t * v3 - 235 continue -c .......... zero b(k+1,k) .......... - 240 s = dabs(b(k1,k1)) + dabs(b(k1,k)) - if (s .eq. 0.0d0) go to 260 - u1 = b(k1,k1) / s - u2 = b(k1,k) / s - r = dsign(dsqrt(u1*u1+u2*u2),u1) - v1 = -(u1 + r) / r - v2 = -u2 / r - u2 = v2 / v1 -c - do 250 i = lor1, ll - t = a(i,k1) + u2 * a(i,k) - a(i,k1) = a(i,k1) + t * v1 - a(i,k) = a(i,k) + t * v2 - t = b(i,k1) + u2 * b(i,k) - b(i,k1) = b(i,k1) + t * v1 - b(i,k) = b(i,k) + t * v2 - 250 continue -c - b(k1,k) = 0.0d0 - if (.not. matz) go to 260 -c - do 255 i = 1, n - t = z(i,k1) + u2 * z(i,k) - z(i,k1) = z(i,k1) + t * v1 - z(i,k) = z(i,k) + t * v2 - 255 continue -c - 260 continue -c .......... end qz step .......... - go to 70 -c .......... set error -- all eigenvalues have not -c converged after 30*n iterations .......... - 1000 ierr = en -c .......... save epsb for use by qzval and qzvec .......... - 1001 if (n .gt. 1) b(n,1) = epsb - return - end
--- a/libcruft/eispack/qzval.f Thu Sep 24 19:00:19 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,277 +0,0 @@ - subroutine qzval(nm,n,a,b,alfr,alfi,beta,matz,z) -c - integer i,j,n,en,na,nm,nn,isw - double precision a(nm,n),b(nm,n),alfr(n),alfi(n),beta(n),z(nm,n) - double precision c,d,e,r,s,t,an,a1,a2,bn,cq,cz,di,dr,ei,ti,tr,u1, - x u2,v1,v2,a1i,a11,a12,a2i,a21,a22,b11,b12,b22,sqi,sqr, - x ssi,ssr,szi,szr,a11i,a11r,a12i,a12r,a22i,a22r,epsb - logical matz -c -c this subroutine is the third step of the qz algorithm -c for solving generalized matrix eigenvalue problems, -c siam j. numer. anal. 10, 241-256(1973) by moler and stewart. -c -c this subroutine accepts a pair of real matrices, one of them -c in quasi-triangular form and the other in upper triangular form. -c it reduces the quasi-triangular matrix further, so that any -c remaining 2-by-2 blocks correspond to pairs of complex -c eigenvalues, and returns quantities whose ratios give the -c generalized eigenvalues. it is usually preceded by qzhes -c and qzit and may be followed by qzvec. -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. -c -c a contains a real upper quasi-triangular matrix. -c -c b contains a real upper triangular matrix. in addition, -c location b(n,1) contains the tolerance quantity (epsb) -c computed and saved in qzit. -c -c matz should be set to .true. if the right hand transformations -c are to be accumulated for later use in computing -c eigenvectors, and to .false. otherwise. -c -c z contains, if matz has been set to .true., the -c transformation matrix produced in the reductions by qzhes -c and qzit, if performed, or else the identity matrix. -c if matz has been set to .false., z is not referenced. -c -c on output -c -c a has been reduced further to a quasi-triangular matrix -c in which all nonzero subdiagonal elements correspond to -c pairs of complex eigenvalues. -c -c b is still in upper triangular form, although its elements -c have been altered. b(n,1) is unaltered. -c -c alfr and alfi contain the real and imaginary parts of the -c diagonal elements of the triangular matrix that would be -c obtained if a were reduced completely to triangular form -c by unitary transformations. non-zero values of alfi occur -c in pairs, the first member positive and the second negative. -c -c beta contains the diagonal elements of the corresponding b, -c normalized to be real and non-negative. the generalized -c eigenvalues are then the ratios ((alfr+i*alfi)/beta). -c -c z contains the product of the right hand transformations -c (for all three steps) if matz has been set to .true. -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 - epsb = b(n,1) - isw = 1 -c .......... find eigenvalues of quasi-triangular matrices. -c for en=n step -1 until 1 do -- .......... - do 510 nn = 1, n - en = n + 1 - nn - na = en - 1 - if (isw .eq. 2) go to 505 - if (en .eq. 1) go to 410 - if (a(en,na) .ne. 0.0d0) go to 420 -c .......... 1-by-1 block, one real root .......... - 410 alfr(en) = a(en,en) - if (b(en,en) .lt. 0.0d0) alfr(en) = -alfr(en) - beta(en) = dabs(b(en,en)) - alfi(en) = 0.0d0 - go to 510 -c .......... 2-by-2 block .......... - 420 if (dabs(b(na,na)) .le. epsb) go to 455 - if (dabs(b(en,en)) .gt. epsb) go to 430 - a1 = a(en,en) - a2 = a(en,na) - bn = 0.0d0 - go to 435 - 430 an = dabs(a(na,na)) + dabs(a(na,en)) + dabs(a(en,na)) - x + dabs(a(en,en)) - bn = dabs(b(na,na)) + dabs(b(na,en)) + dabs(b(en,en)) - a11 = a(na,na) / an - a12 = a(na,en) / an - a21 = a(en,na) / an - a22 = a(en,en) / an - b11 = b(na,na) / bn - b12 = b(na,en) / bn - b22 = b(en,en) / bn - e = a11 / b11 - ei = a22 / b22 - s = a21 / (b11 * b22) - t = (a22 - e * b22) / b22 - if (dabs(e) .le. dabs(ei)) go to 431 - e = ei - t = (a11 - e * b11) / b11 - 431 c = 0.5d0 * (t - s * b12) - d = c * c + s * (a12 - e * b12) - if (d .lt. 0.0d0) go to 480 -c .......... two real roots. -c zero both a(en,na) and b(en,na) .......... - e = e + (c + dsign(dsqrt(d),c)) - a11 = a11 - e * b11 - a12 = a12 - e * b12 - a22 = a22 - e * b22 - if (dabs(a11) + dabs(a12) .lt. - x dabs(a21) + dabs(a22)) go to 432 - a1 = a12 - a2 = a11 - go to 435 - 432 a1 = a22 - a2 = a21 -c .......... choose and apply real z .......... - 435 s = dabs(a1) + dabs(a2) - u1 = a1 / s - u2 = a2 / s - r = dsign(dsqrt(u1*u1+u2*u2),u1) - v1 = -(u1 + r) / r - v2 = -u2 / r - u2 = v2 / v1 -c - do 440 i = 1, en - t = a(i,en) + u2 * a(i,na) - a(i,en) = a(i,en) + t * v1 - a(i,na) = a(i,na) + t * v2 - t = b(i,en) + u2 * b(i,na) - b(i,en) = b(i,en) + t * v1 - b(i,na) = b(i,na) + t * v2 - 440 continue -c - if (.not. matz) go to 450 -c - do 445 i = 1, n - t = z(i,en) + u2 * z(i,na) - z(i,en) = z(i,en) + t * v1 - z(i,na) = z(i,na) + t * v2 - 445 continue -c - 450 if (bn .eq. 0.0d0) go to 475 - if (an .lt. dabs(e) * bn) go to 455 - a1 = b(na,na) - a2 = b(en,na) - go to 460 - 455 a1 = a(na,na) - a2 = a(en,na) -c .......... choose and apply real q .......... - 460 s = dabs(a1) + dabs(a2) - if (s .eq. 0.0d0) go to 475 - u1 = a1 / s - u2 = a2 / s - r = dsign(dsqrt(u1*u1+u2*u2),u1) - v1 = -(u1 + r) / r - v2 = -u2 / r - u2 = v2 / v1 -c - do 470 j = na, n - t = a(na,j) + u2 * a(en,j) - a(na,j) = a(na,j) + t * v1 - a(en,j) = a(en,j) + t * v2 - t = b(na,j) + u2 * b(en,j) - b(na,j) = b(na,j) + t * v1 - b(en,j) = b(en,j) + t * v2 - 470 continue -c - 475 a(en,na) = 0.0d0 - b(en,na) = 0.0d0 - alfr(na) = a(na,na) - alfr(en) = a(en,en) - if (b(na,na) .lt. 0.0d0) alfr(na) = -alfr(na) - if (b(en,en) .lt. 0.0d0) alfr(en) = -alfr(en) - beta(na) = dabs(b(na,na)) - beta(en) = dabs(b(en,en)) - alfi(en) = 0.0d0 - alfi(na) = 0.0d0 - go to 505 -c .......... two complex roots .......... - 480 e = e + c - ei = dsqrt(-d) - a11r = a11 - e * b11 - a11i = ei * b11 - a12r = a12 - e * b12 - a12i = ei * b12 - a22r = a22 - e * b22 - a22i = ei * b22 - if (dabs(a11r) + dabs(a11i) + dabs(a12r) + dabs(a12i) .lt. - x dabs(a21) + dabs(a22r) + dabs(a22i)) go to 482 - a1 = a12r - a1i = a12i - a2 = -a11r - a2i = -a11i - go to 485 - 482 a1 = a22r - a1i = a22i - a2 = -a21 - a2i = 0.0d0 -c .......... choose complex z .......... - 485 cz = dsqrt(a1*a1+a1i*a1i) - if (cz .eq. 0.0d0) go to 487 - szr = (a1 * a2 + a1i * a2i) / cz - szi = (a1 * a2i - a1i * a2) / cz - r = dsqrt(cz*cz+szr*szr+szi*szi) - cz = cz / r - szr = szr / r - szi = szi / r - go to 490 - 487 szr = 1.0d0 - szi = 0.0d0 - 490 if (an .lt. (dabs(e) + ei) * bn) go to 492 - a1 = cz * b11 + szr * b12 - a1i = szi * b12 - a2 = szr * b22 - a2i = szi * b22 - go to 495 - 492 a1 = cz * a11 + szr * a12 - a1i = szi * a12 - a2 = cz * a21 + szr * a22 - a2i = szi * a22 -c .......... choose complex q .......... - 495 cq = dsqrt(a1*a1+a1i*a1i) - if (cq .eq. 0.0d0) go to 497 - sqr = (a1 * a2 + a1i * a2i) / cq - sqi = (a1 * a2i - a1i * a2) / cq - r = dsqrt(cq*cq+sqr*sqr+sqi*sqi) - cq = cq / r - sqr = sqr / r - sqi = sqi / r - go to 500 - 497 sqr = 1.0d0 - sqi = 0.0d0 -c .......... compute diagonal elements that would result -c if transformations were applied .......... - 500 ssr = sqr * szr + sqi * szi - ssi = sqr * szi - sqi * szr - i = 1 - tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21 - x + ssr * a22 - ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22 - dr = cq * cz * b11 + cq * szr * b12 + ssr * b22 - di = cq * szi * b12 + ssi * b22 - go to 503 - 502 i = 2 - tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21 - x + cq * cz * a22 - ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21 - dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22 - di = -ssi * b11 - sqi * cz * b12 - 503 t = ti * dr - tr * di - j = na - if (t .lt. 0.0d0) j = en - r = dsqrt(dr*dr+di*di) - beta(j) = bn * r - alfr(j) = an * (tr * dr + ti * di) / r - alfi(j) = an * t / r - if (i .eq. 1) go to 502 - 505 isw = 3 - isw - 510 continue - b(n,1) = epsb -c - return - end
--- a/src/DLD-FUNCTIONS/balance.cc Thu Sep 24 19:00:19 1998 +0000 +++ b/src/DLD-FUNCTIONS/balance.cc Fri Sep 25 02:53:39 1998 +0000 @@ -32,14 +32,32 @@ #include "CmplxAEPBAL.h" #include "dbleAEPBAL.h" #include "dbleAEPBAL.h" -#include "dbleGEPBAL.h" #include "defun-dld.h" #include "error.h" +#include "f77-fcn.h" #include "gripes.h" #include "oct-obj.h" #include "utils.h" +extern "C" +{ + int F77_FCN( dggbal, DGGBAL) (const char* JOB, const int& N, + double* A, const int& LDA, double* B, const int& LDB, + int& ILO, int & IHI, double* LSCALE, + double* RSCALE, double* WORK, int& INFO, long ); + + int F77_FCN( dggbak, DGGBAK) (const char* JOB, const char* SIDE, + const int& N, const int& ILO, const int& IHI, + double* LSCALE, double* RSCALE, int& M, + double* V, const int& LDV, int& INFO, long, long); + + int F77_FCN( zggbal, ZGGBAL) ( const char* JOB, const int& N, + Complex* A, const int& LDA, Complex* B, const int& LDB, + int& ILO, int & IHI, double* LSCALE, + double* RSCALE, double* WORK, int& INFO, long ); +} + DEFUN_DLD (balance, args, nargout, "AA = balance (A [, OPT]) or [[DD,] AA] = balance (A [, OPT])\n\ \n\ @@ -55,10 +73,11 @@ B: (default) permute and scale, in that order. Rows/columns\n\ of a (and b) that are isolated by permutation are not scaled\n\ \n\ -[DD, AA] = balance (A, OPT) returns aa = inv(dd)*a*dd,\n\ +[DD, AA] = balance (A, OPT) returns aa = dd*a*dd,\n\ \n\ [CC, DD, AA, BB] = balance (A, B, OPT) returns AA (BB) = CC*A*DD (CC*B*DD)") { + octave_value_list retval; int nargin = args.length (); @@ -69,212 +88,192 @@ return retval; } + // determine if it's AEP or GEP + int AEPcase = (nargin == 1 ? 1 : args(1).is_string() ); string bal_job; - int my_nargin; // # args w/o optional string arg - - // Determine if balancing option is listed. Set my_nargin to the - // number of matrix inputs. - if (args(nargin-1).is_string ()) - { - bal_job = args(nargin-1).string_value (); - my_nargin = nargin-1; - } - else - { - bal_job = "B"; - my_nargin = nargin; - } + // problem dimension + int nn = args(0).rows (); - octave_value arg_a = args(0); - - int a_nr = arg_a.rows (); - int a_nc = arg_a.columns (); - - // Check argument 1 dimensions. - - int arg_is_empty = empty_arg ("balance", a_nr, a_nc); + int arg_is_empty = empty_arg ("balance", nn, args(0).columns()); if (arg_is_empty < 0) return retval; if (arg_is_empty > 0) return octave_value_list (2, Matrix ()); - if (a_nr != a_nc) - { - gripe_square_matrix_required ("balance"); - return retval; - } + if (nn != args(0).columns()) + { + gripe_square_matrix_required ("balance"); + return retval; + } // Extract argument 1 parameter for both AEP and GEP. - Matrix aa; ComplexMatrix caa; - if (arg_a.is_complex_type ()) - caa = arg_a.complex_matrix_value (); - else - aa = arg_a.matrix_value (); - - if (error_state) - return retval; + if (args(0).is_complex_type ()) caa = args(0).complex_matrix_value (); + else aa = args(0).matrix_value (); + if (error_state) return retval; // Treat AEP/GEP cases. + if(AEPcase) + { + // Algebraic eigenvalue problem. + if(nargin == 1) + bal_job = "B"; + else if(args(1).is_string()) + bal_job = args(1).string_value(); + // the next line should never execute, but better safe than sorry. + else error("balance: AEP argument 2 must be a string"); - switch (my_nargin) + // balance the AEP + if (args(0).is_complex_type ()) + { + ComplexAEPBALANCE result (caa, bal_job); + + if (nargout == 0 || nargout == 1) + retval(0) = result.balanced_matrix (); + else + { + retval(1) = result.balanced_matrix (); + retval(0) = result.balancing_matrix (); + } + } + else { - case 1: - - // Algebraic eigenvalue problem. + AEPBALANCE result (aa, bal_job); - if (arg_a.is_complex_type ()) - { - ComplexAEPBALANCE result (caa, bal_job); + if (nargout == 0 || nargout == 1) + retval(0) = result.balanced_matrix (); + else + { + retval(1) = result.balanced_matrix (); + retval(0) = result.balancing_matrix (); + } + } + } + // + // end of AEP case, now do GEP case + else + { + // Generalized eigenvalue problem. + if(nargin == 2) + bal_job = "B"; + else if(args(2).is_string()) + bal_job = args(2).string_value(); + else error("balance: GEP argument 3 must be a string"); + + if( (nn != args(1).columns()) || (nn != args(1).rows() ) ) + { + gripe_nonconformant (); + return retval; + } + Matrix bb; + ComplexMatrix cbb; + if (args(1).is_complex_type ()) cbb = args(1).complex_matrix_value (); + else bb = args(1).matrix_value (); + if (error_state) return retval; - if (nargout == 0 || nargout == 1) - retval(0) = result.balanced_matrix (); - else - { - retval(1) = result.balanced_matrix (); - retval(0) = result.balancing_matrix (); - } - } + // + // Both matrices loaded, now let's check what kind of arithmetic: + // first, declare variables used in both the real and complex case + int ilo, ihi, info; + RowVector lscale(nn), rscale(nn), work(6*nn); + char job = bal_job[0]; + static int complex_case = (args(0).is_complex_type() + || args(1).is_complex_type()); + + // now balance + if (complex_case) + { + if (args(0).is_real_type ()) caa = aa; + if (args(1).is_real_type ()) cbb = bb; + + F77_XFCN( zggbal, ZGGBAL, ( &job, nn, caa.fortran_vec(), nn, + cbb.fortran_vec(), nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), work.fortran_vec(), info, 1L)); + } + else // real matrices case + { + F77_XFCN( dggbal, DGGBAL, (&job, nn, aa.fortran_vec(), + nn, bb.fortran_vec() , nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), work.fortran_vec(), info , 1L)); + + if(f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in balance GEP"); + } + + // Since we just want the balancing matrices, we can use dggbal + // for both the real and complex cases; + Matrix Pl(nn,nn), Pr(nn,nn); + for(int ii=0; ii < nn ; ii++) + for( int jj=0; jj < nn ; jj++) + Pl(ii,jj) = Pr(ii,jj) = (ii == jj ? 1.0 : 0.0); + + // left first + F77_XFCN( dggbak, DGGBAK, (&job, "L", + nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), nn, Pl.fortran_vec(), + nn, info, 1L, 1L)); + + if(f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in balance GEP(L)"); + + // then right + F77_XFCN(dggbak, DGGBAK, (&job, "R", + nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), nn, Pr.fortran_vec(), + nn, info, 1L, 1L)); + if(f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in balance GEP(R)"); + + switch (nargout) + { + case 0: + case 1: + warning ("balance: used GEP, should have two output arguments"); + if(complex_case) + retval(0) = caa; else - { - AEPBALANCE result (aa, bal_job); - - if (nargout == 0 || nargout == 1) - retval(0) = result.balanced_matrix (); - else - { - retval(1) = result.balanced_matrix (); - retval(0) = result.balancing_matrix (); - } - } + retval(0) = aa; break; case 2: + if(complex_case) { - // Generalized eigenvalue problem. - - // 1st we have to check argument 2 dimensions and type... - - octave_value arg_b = args(1); - - int b_nr = arg_b.rows (); - int b_nc = arg_b.columns (); - - // Check argument 2 dimensions -- must match arg 1. - - if (b_nr != b_nc || b_nr != a_nr) - { - gripe_nonconformant (); - return retval; - } - - // Now, extract the second matrix... - // Extract argument 1 parameter for both AEP and GEP. - - Matrix bb; - ComplexMatrix cbb; - if (arg_b.is_complex_type ()) - cbb = arg_b.complex_matrix_value (); - else - bb = arg_b.matrix_value (); - - if (error_state) - return retval; - - // Both matrices loaded, now let's check what kind of arithmetic: - - if (arg_a.is_complex_type () || arg_b.is_complex_type ()) - { - if (arg_a.is_real_type ()) - caa = aa; - - if (arg_b.is_real_type ()) - cbb = bb; - - // Compute magnitudes of elements for balancing purposes. - // Surely there's a function I can call someplace! - - for (int i = 0; i < a_nr; i++) - for (int j = 0; j < a_nc; j++) - { - aa (i, j) = abs (caa (i, j)); - bb (i, j) = abs (cbb (i, j)); - } - } - - GEPBALANCE result (aa, bb, bal_job); - - if (arg_a.is_complex_type () || arg_b.is_complex_type ()) - { - caa = result.left_balancing_matrix () * caa - * result.right_balancing_matrix (); - - cbb = result.left_balancing_matrix () * cbb - * result.right_balancing_matrix (); - - switch (nargout) - { - case 0: - case 1: - warning ("balance: should use two output arguments"); - retval(0) = caa; - break; - - case 2: - retval(1) = cbb; - retval(0) = caa; - break; - - case 4: - retval(3) = cbb; - retval(2) = caa; - retval(1) = result.right_balancing_matrix (); - retval(0) = result.left_balancing_matrix (); - break; - - default: - error ("balance: invalid number of output arguments"); - break; - } - } - else - { - switch (nargout) - { - case 0: - case 1: - warning ("balance: should use two output arguments"); - retval(0) = result.balanced_a_matrix (); - break; - - case 2: - retval(1) = result.balanced_b_matrix (); - retval(0) = result.balanced_a_matrix (); - break; - - case 4: - retval(3) = result.balanced_b_matrix (); - retval(2) = result.balanced_a_matrix (); - retval(1) = result.right_balancing_matrix (); - retval(0) = result.left_balancing_matrix (); - break; - - default: - error ("balance: invalid number of output arguments"); - break; - } - } + retval(1) = cbb; + retval(0) = caa; + } + else + { + retval(1) = bb; + retval(0) = aa; } break; + case 4: + if(complex_case) + { + retval(3) = cbb; + retval(2) = caa; + } + else + { + retval(3) = bb; + retval(2) = aa; + } + retval(1) = Pr; + retval(0) = Pl.transpose(); // so that aa_bal = cc*aa*dd, etc. + break; + default: - error ("balance requires one (AEP) or two (GEP) numeric arguments"); + error ("balance: invalid number of output arguments"); break; } - + } return retval; }
--- a/src/DLD-FUNCTIONS/qzval.cc Thu Sep 24 19:00:19 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,72 +0,0 @@ -/* - -Copyright (C) 1996, 1997 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 2, or (at your option) any -later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, write to the Free -Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - -*/ - -// Written by A. S. Hodel <scotte@eng.auburn.edu> - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <cfloat> - -#include "defun-dld.h" -#include "error.h" -#include "gripes.h" -#include "oct-obj.h" - -DEFUN_DLD (qzval, args, , - "X = qzval (A, B)\n\ -\n\ -compute generalized eigenvalues of the matrix pencil (A - lambda B).\n\ -A and B must be real matrices.") -{ - octave_value retval; - - int nargin = args.length (); - - if (nargin == 2) - { - octave_value arg_a = args(0); - octave_value arg_b = args(1); - - Matrix a = arg_a.matrix_value (); - Matrix b = arg_b.matrix_value (); - - if (! error_state) - { - ComplexColumnVector tmp = Qzval (a, b); - - if (! error_state) - retval = tmp; - } - } - else - print_usage ("qzval"); - - return retval; -} - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/