From bokg@cs.umu.se Mon Apr 2 11:53:55 1990 Received: from zeus.cs.umu.se by cs.utk.edu with SMTP (5.61++/2.3-UTK) id AA08978; Mon, 2 Apr 90 11:08:25 -0400 Received: from ikaros.cs.umu.se by zeus.cs.umu.se (5.61+IDA/KTH/LTH/89-09-19) id AAzeus05652; Sun, 1 Apr 90 13:00:09 +0200 Received: by ikaros.cs.umu.se (5.61+IDA/KTH/LTH/89-09-19) id AAikaros01486; Sun, 1 Apr 90 12:59:59 +0200 Return-Path: Date: Sun, 1 Apr 90 12:59:59 +0200 From: bokg@cs.umu.se Message-Id: <9004011059.AAikaros01486@ikaros.cs.umu.se> To: dongarra@cs.utk.edu Subject: guptri_for_netlib Status: R #! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh readme <<'END_OF_readme' XThis package of routines for computing the generalized Schur decomposition Xof an arbitrary (singular) pencil A-zB consists of the following files Xcontaining F77 subroutines and functions belonging to this package. XThey are: X zblas.f X zbnd.f subroutines bound and evalbd described in software paper X zcmatmlr.f X zguptri.f subroutine guptri described in software paper X zlinpack.f X zlistr.f X zmiscl.f X zqz.f X zrcsvdc.f X zreorder.f subroutine reordr described in software paper X zrzstr.f X XAll these files start with a statement describing the contents of Xthe actual file. X X XEnclosed with these files are also X zgschurm.f example program X kcfin.c1 input file for zgschurm.f for example C1 in software paper X zgschur.c1 output file for example C1 in paper X X XA standard usage of the package is as follows: X X call guptri (...) Compute generalized Schur decomposition of singular A-zB. X call reordr (...) Reorder the eigenvalues in specified order. X call bound (...) Compute error bounds for selected eigenvalues X call evalbd (...) and reducing subspaces. X XThe following papers describe software,algorithms and error bounds Xused in the package: X XJ. Demmel and B. Kagstrom, " The generalized Schur decomposition X of an arbitrary pencil A - zB: robust software with error bounds X and applications", Report UMINF-170.90,Institute of Information X Processing, Univ. of Umea, S-901 87 UMEA, SWEDEN,January 1990 X (submitted to ACM TOMS) X XJ. Demmel and B. Kagstrom, "Accurate Solutions of Ill-posed Problems X in Control Theory", SIAM J. Matrix Anal Appl, Vol 9, 1988, pp 126-145 X XJ. Demmel and B. Kagstrom, "Stable Eigendecompositions of Matrix Pencils", X Linear Algebra Applic., Vol 88/89, 1987, pp 137-186 X XJ. Demmel and B. Kagstrom, "Stably Computing the Kronecker Structure X and Reducing Subspaces of Singular pencils A-zB for Uncertain Data", X in J. Cullum and R. Willoughby (eds), Large Scale Eigenvalue Problems, X North Holland, 1986, pp 283-323 X XB. Kagstrom, "RGSVD - An Algorithm for Computing the Kronecker Structure X and Reducing Subspaces of Singular Matrix Pencils", SIAM J. Sci. Stat. X Comp., Vol 7, 1986, pp 185-211 X XAny comments or questions should be sent to: X X Bo Kagstrom X Institute of Information Processing X University of Umea X S-901 87 Umea, Sweden X email: na.kagstrom@na-net.stanford.edu X or bokg@cs.umu.se X X or X X James Demmel X Courant Institute X New York University X 215 Mercer Street X New York, NY 10012, USA X email: na.demmel@na-net.stanford.edu X XNotices: X1. The main program in file zgschurm.f with input from kcfin.c1 X produced the output on file zgschur.c1 when run on X Sun 3/80 workstation. The output in file zgschur.c1 also X includes some information that is not explained in the software paper. X We refer to the source for more information. X X2. The current version of the package has been developed during a period X of 4-5 years. The current version of the routines does not make use X of level 2 or 3 BLAS. X X3. Before a production code of this package is produced we would like X to obtain and collect as much information from users as possible. X THANK YOU IN ADVANCE! X X X X END_OF_readme if test 3177 -ne `wc -c kcfin.c1 <<'END_OF_kcfin.c1' X 4 5 X(1., 0.) (-2., 0.) (0.,0.) ( 0., 0.) ( 0., 0.) X(1., 0.) (0., 0.) (-1., 0.) (0., 0.) (0., 0.) X(0., 0.) (0., 0.) (0., 0.) (1., 0.) (0., 0.) X(0., 0.) (0., 0.) (0., 0.) (0., 0.) (2., 0.) X(0., 0.) (1., 0.) (0., 0.) (0., 0.) (0., 0.) X(0., 0.) (0., 0.) (1., 0.) (0., 0.) (0., 0.) X(0., 0.) (0., 0.) (0., 0.) (1., 0.) (0., 0.) X(0., 0.) (0., 0.) (0., 0.) (0., 0.) (1., 0.) X11000000000000000100 X10000110 X1.d-8 1000. X1.d-10 X 1 1 3 X1.d-10 1.d-9 1.d-8 1.d-7 1.d-6 X1.d-5 1.d-4 1.d-3 END_OF_kcfin.c1 if test 527 -ne `wc -c zblas.f <<'END_OF_zblas.f' Xc On this file - blas routines: double precision complex Xc zaxpy, zswap, dcabs1, dznrm2, zcopy, zdotc, zdotu, zscal, Xc zrotg, zdrot, drotg Xc X subroutine zaxpy(n,za,zx,incx,zy,incy) Xc Xc constant times a vector plus a vector. Xc jack dongarra, 3/11/78. Xc X double complex zx(1),zy(1),za X double precision dcabs1 X if(n.le.0)return X if (dcabs1(za) .eq. 0.0d0) return X if (incx.eq.1.and.incy.eq.1)go to 20 Xc Xc code for unequal increments or equal increments Xc not equal to 1 Xc X ix = 1 X iy = 1 X if(incx.lt.0)ix = (-n+1)*incx + 1 X if(incy.lt.0)iy = (-n+1)*incy + 1 X do 10 i = 1,n X zy(iy) = zy(iy) + za*zx(ix) X ix = ix + incx X iy = iy + incy X 10 continue X return Xc Xc code for both increments equal to 1 Xc X 20 do 30 i = 1,n X zy(i) = zy(i) + za*zx(i) X 30 continue X return X end X X subroutine zswap (n,zx,incx,zy,incy) Xc Xc interchanges two vectors. Xc jack dongarra, 3/11/78. Xc X double complex zx(1),zy(1),ztemp Xc X if(n.le.0)return X if(incx.eq.1.and.incy.eq.1)go to 20 Xc Xc code for unequal increments or equal increments not equal Xc to 1 Xc X ix = 1 X iy = 1 X if(incx.lt.0)ix = (-n+1)*incx + 1 X if(incy.lt.0)iy = (-n+1)*incy + 1 X do 10 i = 1,n X ztemp = zx(ix) X zx(ix) = zy(iy) X zy(iy) = ztemp X ix = ix + incx X iy = iy + incy X 10 continue X return Xc Xc code for both increments equal to 1 X 20 do 30 i = 1,n X ztemp = zx(i) X zx(i) = zy(i) X zy(i) = ztemp X 30 continue X return X end X X double precision function dcabs1(z) X double complex z,zz X double precision t(2) X equivalence (zz,t(1)) X zz = z X dcabs1 = dabs(t(1)) + dabs(t(2)) X return X end X X X double precision function dznrm2( n, zx, incx) X logical imag, scale X integer next X double precision cutlo, cuthi, hitest, sum, xmax, absx, zero, one X double complex zx(1) X double precision dreal,dimag X double complex zdumr,zdumi X dreal(zdumr) = zdumr X dimag(zdumi) = (0.0d0,-1.0d0)*zdumi X data zero, one /0.0d0, 1.0d0/ Xc Xc unitary norm of the complex n-vector stored in zx() with storage Xc increment incx . Xc if n .le. 0 return with result = 0. Xc if n .ge. 1 then incx must be .ge. 1 Xc Xc c.l.lawson , 1978 jan 08 Xc Xc four phase method using two built-in constants that are Xc hopefully applicable to all machines. Xc cutlo = maximum of sqrt(u/eps) over all known machines. Xc cuthi = minimum of sqrt(v) over all known machines. Xc where Xc eps = smallest no. such that eps + 1. .gt. 1. Xc u = smallest positive no. (underflow limit) Xc v = largest no. (overflow limit) Xc Xc brief outline of algorithm.. Xc Xc phase 1 scans zero components. Xc move to phase 2 when a component is nonzero and .le. cutlo Xc move to phase 3 when a component is .gt. cutlo Xc move to phase 4 when a component is .ge. cuthi/m Xc where m = n for x() real and m = 2*n for complex. Xc Xc values for cutlo and cuthi.. Xc from the environmental parameters listed in the imsl converter Xc document the limiting values are as follows.. Xc cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are Xc univac and dec at 2**(-103) Xc thus cutlo = 2**(-51) = 4.44089e-16 Xc cuthi, s.p. v = 2**127 for univac, honeywell, and dec. Xc thus cuthi = 2**(63.5) = 1.30438e19 Xc cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. Xc thus cutlo = 2**(-33.5) = 8.23181d-11 Xc cuthi, d.p. same as s.p. cuthi = 1.30438d19 Xc data cutlo, cuthi / 8.232d-11, 1.304d19 / Xc data cutlo, cuthi / 4.441e-16, 1.304e19 / X data cutlo, cuthi / 8.232d-11, 1.304d19 / Xc X if(n .gt. 0) go to 10 X dznrm2 = zero X go to 300 Xc X 10 assign 30 to next X sum = zero X nn = n * incx Xc begin main loop X do 210 i=1,nn,incx X absx = dabs(dreal(zx(i))) X imag = .false. X go to next,(30, 50, 70, 90, 110) X 30 if( absx .gt. cutlo) go to 85 X assign 50 to next X scale = .false. Xc Xc phase 1. sum is zero Xc X 50 if( absx .eq. zero) go to 200 X if( absx .gt. cutlo) go to 85 Xc Xc prepare for phase 2. X assign 70 to next X go to 105 Xc Xc prepare for phase 4. Xc X 100 assign 110 to next X sum = (sum / absx) / absx X 105 scale = .true. X xmax = absx X go to 115 Xc Xc phase 2. sum is small. Xc scale to avoid destructive underflow. Xc X 70 if( absx .gt. cutlo ) go to 75 Xc Xc common code for phases 2 and 4. Xc in phase 4 sum is large. scale to avoid overflow. Xc X 110 if( absx .le. xmax ) go to 115 X sum = one + sum * (xmax / absx)**2 X xmax = absx X go to 200 Xc X 115 sum = sum + (absx/xmax)**2 X go to 200 Xc Xc Xc prepare for phase 3. Xc X 75 sum = (sum * xmax) * xmax Xc X 85 assign 90 to next X scale = .false. Xc Xc for real or d.p. set hitest = cuthi/n Xc for complex set hitest = cuthi/(2*n) Xc X hitest = cuthi/float( n ) Xc Xc phase 3. sum is mid-range. no scaling. Xc X 90 if(absx .ge. hitest) go to 100 X sum = sum + absx**2 X 200 continue Xc control selection of real and imaginary parts. Xc X if(imag) go to 210 X absx = dabs(dimag(zx(i))) X imag = .true. X go to next,( 50, 70, 90, 110 ) Xc X 210 continue Xc Xc end of main loop. Xc compute square root and adjust for scaling. Xc X dznrm2 = dsqrt(sum) X if(scale) dznrm2 = dznrm2 * xmax X 300 continue X return X end X X subroutine zcopy(n,zx,incx,zy,incy) Xc Xc copies a vector, x, to a vector, y. Xc jack dongarra, linpack, 4/11/78. Xc X double complex zx(1),zy(1) X integer i,incx,incy,ix,iy,n Xc X if(n.le.0)return X if(incx.eq.1.and.incy.eq.1)go to 20 Xc Xc code for unequal increments or equal increments Xc not equal to 1 Xc X ix = 1 X iy = 1 X if(incx.lt.0)ix = (-n+1)*incx + 1 X if(incy.lt.0)iy = (-n+1)*incy + 1 X do 10 i = 1,n X zy(iy) = zx(ix) X ix = ix + incx X iy = iy + incy X 10 continue X return Xc Xc code for both increments equal to 1 Xc X 20 do 30 i = 1,n X zy(i) = zx(i) X 30 continue X return X end X X double complex function zdotc(n,zx,incx,zy,incy) Xc Xc forms the dot product of a vector. Xc jack dongarra, 3/11/78. Xc X double complex zx(1),zy(1),ztemp X ztemp = (0.0d0,0.0d0) X zdotc = (0.0d0,0.0d0) X if(n.le.0)return X if(incx.eq.1.and.incy.eq.1)go to 20 Xc Xc code for unequal increments or equal increments Xc not equal to 1 Xc X ix = 1 X iy = 1 X if(incx.lt.0)ix = (-n+1)*incx + 1 X if(incy.lt.0)iy = (-n+1)*incy + 1 X do 10 i = 1,n X ztemp = ztemp + dconjg(zx(ix))*zy(iy) X ix = ix + incx X iy = iy + incy X 10 continue X zdotc = ztemp X return Xc Xc code for both increments equal to 1 Xc X 20 do 30 i = 1,n X ztemp = ztemp + dconjg(zx(i))*zy(i) X 30 continue X zdotc = ztemp X return X end X X double complex function zdotu(n,zx,incx,zy,incy) Xc Xc forms the dot product of a vector. Xc jack dongarra, 3/11/78. Xc X double complex zx(1),zy(1),ztemp X ztemp = (0.0d0,0.0d0) X zdotu = (0.0d0,0.0d0) X if(n.le.0)return X if(incx.eq.1.and.incy.eq.1)go to 20 Xc Xc code for unequal increments or equal increments Xc not equal to 1 Xc X ix = 1 X iy = 1 X if(incx.lt.0)ix = (-n+1)*incx + 1 X if(incy.lt.0)iy = (-n+1)*incy + 1 X do 10 i = 1,n X ztemp = ztemp + zx(ix)*zy(iy) X ix = ix + incx X iy = iy + incy X 10 continue X zdotu = ztemp X return Xc Xc code for both increments equal to 1 Xc X 20 do 30 i = 1,n X ztemp = ztemp + zx(i)*zy(i) X 30 continue X zdotu = ztemp X return X end X X subroutine zscal(n,za,zx,incx) Xc Xc scales a vector by a constant. Xc jack dongarra, 3/11/78. Xc X double complex za,zx(1) X if(n.le.0)return X if(incx.eq.1)go to 20 Xc Xc code for increments not equal to 1 Xc X ix = 1 X if(incx.lt.0)ix = (-n+1)*incx + 1 X do 10 i = 1,n X zx(ix) = za*zx(ix) X ix = ix + incx X 10 continue X return Xc Xc code for increments equal to 1 Xc X 20 do 30 i = 1,n X zx(i) = za*zx(i) X 30 continue X return X end X X subroutine zrotg(ca,cb,c,s) X double complex ca,cb,s X double precision c,dcabs1 X double precision norm,scale X double complex alpha X if (dcabs1(ca) .ne. 0.0d0) go to 10 X c = 0.0d0 X s = (1.0d0,0.0d0) X ca = cb X go to 20 X 10 continue X scale = dcabs1(ca) + dcabs1(cb) X norm = scale*dsqrt((dcabs1(ca/dcmplx(scale,0.0d0)))**2 + X * (dcabs1(cb/dcmplx(scale,0.0d0)))**2) X alpha = ca /dcabs1(ca) X c = dcabs1(ca) / norm X s = alpha * dconjg(cb) / norm X ca = alpha * norm X 20 continue X return X end X X subroutine zdrot (n,zx,incx,zy,incy,c,s) Xc Xc applies a plane rotation, where the cos and sin (c and s) are Xc double precision and the vectors zx and zy are double complex. Xc jack dongarra, linpack, 3/11/78. Xc X double complex zx(1),zy(1),ztemp X double precision c,s X integer i,incx,incy,ix,iy,n Xc X if(n.le.0)return X if(incx.eq.1.and.incy.eq.1)go to 20 Xc Xc code for unequal increments or equal increments not equal Xc to 1 Xc X ix = 1 X iy = 1 X if(incx.lt.0)ix = (-n+1)*incx + 1 X if(incy.lt.0)iy = (-n+1)*incy + 1 X do 10 i = 1,n X ztemp = c*zx(ix) + s*zy(iy) X zy(iy) = c*zy(iy) - s*zx(ix) X zx(ix) = ztemp X ix = ix + incx X iy = iy + incy X 10 continue X return Xc Xc code for both increments equal to 1 Xc X 20 do 30 i = 1,n X ztemp = c*zx(i) + s*zy(i) X zy(i) = c*zy(i) - s*zx(i) X zx(i) = ztemp X 30 continue X return X end X X subroutine drotg(da,db,c,s) Xc Xc construct givens plane rotation. Xc jack dongarra, linpack, 3/11/78. Xc X double precision da,db,c,s,roe,scale,r,z Xc X roe = db X if( dabs(da) .gt. dabs(db) ) roe = da X scale = dabs(da) + dabs(db) X if( scale .ne. 0.0d0 ) go to 10 X c = 1.0d0 X s = 0.0d0 X r = 0.0d0 X go to 20 X 10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2) X r = dsign(1.0d0,roe)*r X c = da/r X s = db/r X 20 z = 1.0d0 X if( dabs(da) .gt. dabs(db) ) z = s X if( dabs(db) .ge. dabs(da) .and. c .ne. 0.0d0 ) z = 1.0d0/c X da = r X db = z X return X end Xc Xc*** no more on this file X END_OF_zblas.f if test 11182 -ne `wc -c zbnd.f <<'END_OF_zbnd.f' Xc as of june 22, 1987 this file contains Xc bound, ebdreg, gvec, pbound, blddfl, blddfu, bldrhs, prml, Xc prmlct, svdiv, evalbd, bndwsp Xc X subroutine bound(a,b,ldab,m,n,irstrt,icstrt,dimreg, X + evala,evalb,edlmax,gvcond,pqnorm,ecase, X + sdlmax, difl, difu, qnorm, pnorm, scase, X + work, info) Xc Xc implicit none Xc Xc**** debug space X common /debug2/ idbg(20), outunit X integer idbg, outunit Xc Xc**** formal parameter declarations X integer ldab,m,n,irstrt,icstrt,dimreg,info,ecase,scase X complex*16 a(ldab,*), b(ldab,*), evala(*),evalb(*) X complex*16 work(*) X real*8 gvcond(*), edlmax, sdlmax, qnorm, pnorm, pqnorm X real*8 difl, difu Xc Xc******************************************************************** Xc Xc compute error bounds for selected eigenvalues of general pencil Xc and error bounds for left and right reducing subspaces Xc Xc this version requires all selected eigenvalues be simple Xc input pencil a - lambda b must be in guptri form Xc Xc theorems and corollaries referred to below appear in: Xc 'accurate solutions of ill-posed problems in control theory' Xc proc. of the 25th ieee conference on decision and control, Xc athens, greece, december 10-12, 1986, pp 558-563 Xc by j. demmel and b. kagstrom Xc Xc see also: Xc j. demmel and b. kagstrom, 'computing stable eigendecompositions Xc of matrix pencils', linear algebra and its applications, Xc vol 88/89, 1987, pp 139-185 Xc Xc inputs Xc Xc a(ldab,n), b(ldab,n) - complex*16 - input pencil in Xc guptri form Xc Xc lda - integer - leading dimension of a and b Xc Xc m,n - integer - row, column dimensions of a and b Xc Xc irstrt, icstrt - integer - starting row and column of selected Xc part of pencil for which eigenvalue bounds Xc are desired. reducing subspace bounds will be Xc supplied for right reducing subspace spanned Xc by leading icstrt-1 components and for left Xc reducing subspace spanned by leading icstrt-1 Xc components. Xc note: set icstrt=n+1 to make right reducing Xc subspace whole space Xc set irstrt=m+1 to make left reducing Xc subspace whole space Xc Xc dimreg - integer - number of selected eigenvalues; Xc if dimreg.eq.0 only subspace perturbation bounds will be Xc computed Xc (note - one can select a subset of the regular part only; Xc this gives generally different bounds for common eigenvalues Xc from a different selected subset; see paper above for Xc discussion) Xc Xc outputs Xc Xc evala(dimreg), evalb(dimreg) - complex*16 - Xc normalized selected eigenvalues; Xc evala(i)/evalb(i) is i-th eigenvalue and Xc abs(evala(i))**2 + abs(evalb(i))**2 = 1 Xc Xc edlmax - real*8 - maximum frobenius norm of perturbation for Xc which eigenvalue perturbation bounds hold. Xc if no maximum norm then edlmax=-1. Xc Xc gvcond(dimreg) - real*8 - condition numbers; suppose the pencil Xc is perturbed by amount delta .le. edlmax (if edlmax=-1. then Xc delta arbitrary) such that the conditions of theorem 5 or Xc corollary 1 hold (edlmax=-1. implies these conditions always Xc hold). then if c/s is a perturbed eigenvalue such that Xc abs(c)**2 + abs(s)**2 = 1, then for some i Xc abs(c*evalb(i)-s*evala(i)) .le. delta * gvcond(i) Xc Xc pqnorm - real*8 - overall condition number; under same Xc conditions as for gvcond, if areg - lambda breg is regular Xc part of unperturbed pencil in guptri form, then Xc sigma-min(c*breg - s*areg) .le. delta * pqnorm Xc (sigma-min is the smallest singular value) Xc Xc ecase - integer - which of 5 cases for eigenvalue bounds Xc the pencil falls depending on input dimensions; Xc the first four cases are for dimreg.gt.0, in which Xc case the description gives: Xc (part of KCF to above, left of selected part) and Xc (part of KCF to below, right of selected part) Xc ecase=1 - (right singular and/or regular part) and Xc (left singular and/or regular part) Xc ecase=2 - (right singular and/or regular part) and (nothing) Xc ecase=3 - (nothing) and (left singular and/or regular part) Xc ecase=4 - (nothing) and (nothing) Xc ecase=5 - dimreg.eq.0 (no eigenvalue bounds) Xc Xc sdlmax - real*8 - maximum frobenius norm of perturbation for Xc which reducing subspace perturbation bounds hold Xc (if scase=4 (see below) sdlmax=-1. to indicate that Xc this bound does not apply) Xc Xc difl, difu - real*8 - difl and difu functions (used to Xc compute sdlmax, see paper for details) Xc (if scase=4 (see below), both set to 0) Xc Xc qnorm, pnorm, - real*8 - norms of left and right projectors Xc (used in reducing subspace bounds) Xc (if scase=4 (see below), both set to 1) Xc Xc scase - integer - which of 4 cases for reducing subspace Xc bounds the pencil falls depending on input dimensions: Xc scase=1 - both left and right subspaces nontrivial Xc scase=2 - left space trivial (0) and right space nontrivial Xc scase=3 - left space nontrivial and right space trivial Xc (whole space) Xc scase=4 - both spaces trivial (either 0 or whole space) Xc Xc the reducing subspace bounds may be calculated from Xc scase, sdlmax, pnorm and qnorm as follows: Xc let delta be the distance in the frobenius norm from a Xc perturbed pencil with the same structure as a - lambda b Xc to a - lambda b (see the above paper by demmel and Xc kagstrom for more details). if delta.lt.sdlmax then the Xc following bounds apply, where relerr=delta/sdlmax : Xc Xc upper bound on angular perturbation in left reducing subspace Xc if scase=1 (theorem 4, case 1 in paper) Xc atan(relerr/(pnorm-relerr*sqrt(pnorm**2-1))) Xc if scase=2 Xc 0 (since left subspace trivial) Xc if scase=3 Xc atan(relerr/(1-relerr)) Xc if scase=4 Xc 0 (since left subspace trivial) Xc Xc upper bound on angular perturbation in right reducing subspace Xc if scase=1 (theorem 4, case 1 in paper) Xc atan(relerr/(qnorm-relerr*sqrt(qnorm**2-1))) Xc if scase=2 Xc atan(relerr/(1-relerr)) Xc if scase=3 Xc 0 (since right subspace trivial) Xc if scase=4 Xc 0 (since right subspace trivial) Xc Xc lower bound on angular perturbation in left reducing subspace Xc if scase=1 (theorem 4, case 2 in paper) Xc atan(1/(sqrt(2*min(irstrt-1,m-irstrt+1))*pnorm + Xc sqrt(pnorm**2-1))) Xc if scase=2 this bound does not apply Xc if scase=3 this bound does not apply Xc if scase=4 this bound does not apply Xc Xc lower bound on angular perturbation in right reducing subspace Xc if scase=1 (theorem 4, case 2 in paper) Xc atan(1/(sqrt(2*min(icstrt-1,n-icstrt+1))*qnorm + Xc sqrt(qnorm**2-1))) Xc if scase=2 this bound does not apply Xc if scase=3 this bound does not apply Xc if scase=4 this bound does not apply Xc Xc (note: given scase, sdlmax, pnorm, qnorm, m, n, icstrt, irstrt Xc and delta (the frobenius norm of a perturbation), subroutine Xc evalbd will compute the above upper and lower subspace bounds) Xc Xc info - integer - 0 if normal return Xc 1 if svd error in difu calculation in pbound Xc 2 if difu=0 in pbound Xc 3 if svd error in difl calculation in pbound Xc 4 if difl=0 in pbound Xc 5 if multiple eigenvalues Xc 6 if inconsistent input dimensions Xc Xc workspace Xc work(*) - complex*16 - exact amount is complicated function of Xc input dimensions and depends on ecase, and computed Xc as follows: Xc Xc irend=irstrt+dimreg-1; icend=icstrt+dimreg-1; Xc if ecase=1 - m11=irstrt-1; m21=m-m11; n11=icstrt-1; n21=n-n11; Xc m12=irend-irstrt+1; m22=m-irend; Xc n12=icend-icstrt+1; n22=n-icend; Xc workspace = max( (2*n21*m11*(n11*n21+m11*m21+ Xc 2*n21*m11+2)+n11*n21+m11*m21) , Xc (2*((m21*n11+1)*(n11*n21+ Xc m11*m21+1)-1)) , Xc (2*n22*m12*(n12*n22+m12*m22+ Xc 2*n22*m12+2)+n12*n22+m12*m22) , Xc (2*((m22*n12+1)*(n12*n22+ Xc m12*m22+1)-1)) ) Xc if ecase=2 or ecase=5 - Xc m11=irstrt-1; m21=m-m11; n11=icstrt-1; n21=n-n11; Xc workspace = max( (2*n21*m11*(n11*n21+m11*m21+ Xc 2*n21*m11+2)+n11*n21+m11*m21) , Xc (2*((m21*n11+1)*(n11*n21+ Xc m11*m21+1)-1)) ) Xc if ecase=3 - m11=irend; m21=m-m11; n11=icend; n21=n-icend; Xc workspace = max( (2*n21*m11*(n11*n21+m11*m21+ Xc 2*n21*m11+2)+n11*n21+m11*m21) , Xc (2*((m21*n11+1)*(n11*n21+ Xc m11*m21+1)-1)) ) Xc if ecase=4 - workspace = n*n Xc Xc the following simple expression bounds the workspace also, but Xc may occasionally be much too large (especially if ecase=4): Xc workspace .le. 2*m*n* (n*n + m*m + 2*n + m + 2) + n*n + m*m Xc********************************************************************* Xc Xc**** this version dated 16 june 1987 Xc authors: jim demmel and bo kagstrom Xc Xc addresses: Xc jim demmel, courant institute, 251 mercer str, Xc new york, new york 10012, usa Xc electronic address: demmel at nyu.edu or Xc na.demmel at score.stanford.edu Xc bo kagstrom, institute of information processing, Xc university of umea, s-90187 umea, sweden Xc electronic address: bokg at seumdc51.bitnet or Xc na.kagstrom at score.stanford.edu Xc Xc**** bound uses the following functions and subroutines Xc pbound, ebdreg, cmatpr (debug only), gvec, dznrm2 (blas), Xc blddfu, blddfl, bldrhs, prml, prmlct, svdiv, zsvdc (linpack) Xc Xc**** internal variables X integer irend,icend,idummy,i X real*8 rdummy, difu1, difu2, difl1, difl2, pnorm1, pnorm2 X real*8 qnorm1, qnorm2, pdelta1, pdelta2, delta Xc Xc test input dimensions for consistency X info = 0 X if (irstrt.gt.icstrt .or. irstrt.le.0 .or. X + n-icstrt-dimreg.gt.m-irstrt-dimreg .or. X + n-icstrt-dimreg+1.lt.0 .or. dimreg.lt.0) then Xc inconsistent input dimensions X info = 6 X return X endif X icend = icstrt+dimreg-1 X irend = irstrt+dimreg-1 X delta = 0. Xc X if (dimreg.gt.0) then Xc there are eigenvalue bounds to compute Xc Xc ecase 1 - in addition to selected regular part KCF has Xc (right singular part and/or regular part) and Xc (left singular part and/or regular part) X if (icstrt.ne.1 .and. irend.ne.m) then X ecase = 1 X if (irstrt.eq.1) then X scase = 2 X else X scase = 1 X endif Xc see corollary 1 for explanation of bounds X call pbound(a,b,ldab,m,n,irstrt-1,icstrt-1, X + delta,difl1,difu1,qnorm1,pnorm1, pdelta1, X + rdummy,rdummy,rdummy,rdummy,idummy,work,info) X if (info.ne.0) return X call pbound(a(irstrt,icstrt),b(irstrt,icstrt),ldab, X + m-irstrt+1,n-icstrt+1,irend-irstrt+1, X + icend-icstrt+1, X + delta,difl2,difu2,qnorm2,pnorm2,pdelta2, X + rdummy,rdummy,rdummy,rdummy,idummy,work,info) X if (info.ne.0) return X edlmax = min (pdelta1, pdelta2/(sqrt(2.)*qnorm1)) X pqnorm = 2.*pnorm2*qnorm1 Xc X sdlmax = pdelta1 X pnorm = pnorm1 X qnorm = qnorm1 X difl = difl1 X difu = difu1 X endif Xc Xc ecase 2 - in addition to selected regular part KCF has Xc (right singular part and/or regular part) and Xc (nothing) X if (icstrt.ne.1 .and. irend.eq.m) then X ecase=2 X if (irstrt.eq.1) then X scase = 2 X else X scase = 1 X endif Xc see part 1 of theorem 5 for explanation of bounds X call pbound(a,b,ldab,m,n,irstrt-1,icstrt-1,delta,difl1, X + difu1,qnorm1,pnorm1,pdelta1,rdummy,rdummy, X + rdummy,rdummy,idummy,work,info) X if (info.ne.0) return X edlmax= pdelta1 X pqnorm=1. X if (idummy.eq.1) pqnorm=sqrt(2.)*qnorm1 Xc X sdlmax = pdelta1 X pnorm = pnorm1 X qnorm = qnorm1 X difl = difl1 X difu = difu1 X endif Xc Xc ecase 3 - in addition to selected regular part KCF has Xc (nothing) and Xc (left singular part and/or regular part) X if (icstrt.eq.1 .and. irend.ne.m) then X ecase = 3 X scase = 4 Xc see part 2 of theorem 5 for explanation of bounds X call pbound(a,b,ldab,m,n,irend,icend, X + delta,difl2,difu2,qnorm2,pnorm2,pdelta2, X + rdummy,rdummy,rdummy,rdummy,idummy,work,info) X if (info.ne.0) return X edlmax = pdelta2 X pqnorm = 1. X if (idummy.eq.1) pqnorm = sqrt(2.)*pnorm2 X difl = 0. X difu = 0. X pnorm = 1. X qnorm = 1. X sdlmax = -1. X endif Xc Xc ecase 4 - pencil regular and entire spectrum selected X if (icstrt.eq.1 .and. irend.eq.m) then X ecase=4 X edlmax=-1. X pqnorm=1. Xc X scase = 4 X difl = 0. X difu = 0. X pnorm = 1. X qnorm = 1. X sdlmax = -1. X endif Xc X call ebdreg(a,b,ldab,irstrt,icstrt,dimreg, X + gvcond,evala,evalb,work,info) X if (info.ne.0) then X info = 5 X return X endif X if (pqnorm.ne.1.) then X do 1 i=1,dimreg X gvcond(i)=gvcond(i)*pqnorm X1 continue X endif Xc X else Xc dimreg.eq.0, so only compute subspace bounds X ecase = 5 X call pbound(a,b,ldab,m,n,irstrt-1,icstrt-1, X + delta,difl,difu,qnorm,pnorm,sdlmax, X + rdummy,rdummy,rdummy,rdummy,scase,work,info) X endif Xc X if (idbg(20).ne.0) then X write(outunit,100) ldab,m,n,irstrt,icstrt,dimreg,ecase,scase X100 format(' bound - ldab,m,n,irstrt,icstrt,dimreg,ecase,scase=', X + /,8i5) X if (ecase.ne.5) then X write(outunit,101) edlmax,pqnorm X101 format(' edlmax,pqnorm=',2d15.6,/,' gvcond=') X write(outunit,102) (gvcond(i),i=1,dimreg) X102 format(5d15.6) X call cmatpr(work,dimreg,dimreg,dimreg,'gvec') X endif X if (scase.ne.4) write(outunit,103) sdlmax,pnorm,qnorm X103 format(' sdlmax,pnorm,qnorm=',3d15.6) X endif X return X end Xc Xc X subroutine ebdreg(a,b,ldab,irstrt,icstrt,dimreg, X + gvcond,evala,evalb,work,info) Xc implicit none Xc**** formal parameter declarations X integer ldab, dimreg, irstrt, icstrt, info X complex*16 a(ldab,*), b(ldab,*), work(*), evala(*), evalb(*) X real*8 gvcond(*) Xc Xc***************************************************************** Xc Xc compute error bounds for eigenvalues of a regular pencil Xc requires all simple eigenvalues Xc Xc inputs: Xc a(ldab,*), b(ldab,*) - complex*16 - contain pencil Xc irstrt, icstrt - integer - starting row and column locations Xc of pencil within a and b Xc dimreg - integer - dimension of regular pencil Xc Xc outputs: Xc evala(dimreg), evalb(dimreg) - complex*16 - normalized Xc eigenvalues: Xc evala(i)/evalb(i) is i-th eigenvalue and Xc abs(evala(i))**2 + abs(evalb(i))**2 =1 Xc gvcond(dimreg) - real*8 - gvcond(i) is condition number of Xc i-th eigenvalue where if the pencil is perturbed by Xc frobenius norm delta and the perturbed eigenvalue Xc is c/s where Xc abs(c)**2 + abs(s)**2 = 1 then for some i Xc abs(c*evalb(i) - s*evala(i)) .le. delta * gvcond(i) Xc info - integer - returns 0 (normal) if no multiple eigenvalues, Xc else nonzero Xc Xc workspace: Xc work(dimreg**2) - complex*16 - work space Xc Xc*********************************************************************** Xc Xc**** this version dated 16 june 1987 Xc authors: jim demmel and bo kagstrom Xc Xc**** ebdreg uses the following functions and subroutines: Xc gvec Xc Xc**** internal variables Xc X real*8 scl X integer i Xc Xc compute eigenvectors X call gvec(a( irstrt , icstrt ), X + b( irstrt , icstrt ), ldab, X + dimreg , work, dimreg, gvcond, info) Xc Xc compute normalized eigenpairs X do 555 i=1,dimreg X scl=sqrt(abs(a(irstrt-1+i,icstrt-1+i))**2+ X + abs(b(irstrt-1+i,icstrt-1+i))**2) X evala(i) = a(irstrt-1+i,icstrt-1+i)/scl X evalb(i) = b(irstrt-1+i,icstrt-1+i)/scl X if (info.eq.0) gvcond(i)= dimreg * gvcond(i) / scl X555 continue Xc X return X end Xc Xc X subroutine gvec(a,b,ldab,n,vec,ldvec,gvcond,info) Xc Xc implicit none Xc**** debug space X common /debug2/ idbg(20),outunit X integer idbg,outunit X logical ldebug Xc**** formal parameter declarations X integer ldab, n, ldvec, info X complex*16 a(ldab,*), b(ldab,n), vec(ldvec,*) X real*8 gvcond(*) Xc Xc******************************************************************** Xc Xc compute the left and right eigenvectors of the upper triangular Xc regular pencil a - lambda b Xc compute condition numbers of eigenvalues Xc Xc inputs Xc a(ldab,n),b(ldab,n) - complex*16 - n by n matrices Xc ldab - integer - leading dimension of a, b Xc n - integer - dimension of a, b Xc ldvec - integer - leading dimension of vec Xc Xc idbg(10) - if idbg(10) ne 0, print debug output Xc Xc outputs Xc vec(ldvec,n) - complex*16 - matrix containing eigenvectors Xc vec(1:i,i) contains the right eigenvector of the i-th Xc eigenvalue, normalized so vec(i,i)=1. the other Xc components of the eigenvector are zero Xc vec(i:n,i) contains the left eigenvector of the i-th Xc eigenvalue, normalized so vec(i,i)=1. the other Xc components of the eigenvector are zero Xc gvcond(n) - real*8 - array of condition numbers of eigenvalues. Xc if right eigenvectors scaled by diagonal matrix d Xc to have unit norm, scale left eigenvectors by d**-1. Xc then condition number is norm of left eigenvector. Xc info - integer - 0 if pencil regular without multiple eigenvalues Xc nonzero index of a multiple or 0/0 eigenvalue otherwise X Xc*********************************************************************** Xc Xc**** this version dated 16 june 1987 Xc authors: jim demmel and bo kagstrom Xc Xc**** gvec uses the following external function: Xc dznrm2 (blas) X real*8 dznrm2 Xc**** internal variables X integer nm1, i, im1, im2, j, jp1, k, ip1, ip2, jm1 X complex*16 alpha, beta, diag, cmul, csum X real*8 ca, cb, dmax, dmin, d Xc X ldebug=(idbg(10).ne.0) X info=0 X nm1=n-1 X if (ldebug) write(outunit,99) X99 format(' entering gvec') X do 1 i=1,n Xc X if (ldebug) write(outunit,100) i X100 format(' i=',i4) X vec(i,i)=1. Xc Xc compute alpha, beta so that zz = beta*a - alpha*b is a Xc singular matrix whose left and right null spaces are the Xc left and right eigenspaces we seek X ca=abs(a(i,i)) X cb=abs(b(i,i)) X dmax=max(ca,cb) X if (ldebug) write(outunit,101) a(i,i),b(i,i),ca,cb,dmax X101 format(' a(i,i)=',2d20.5,/,' b(i,i)=',2d20.5,/,' ca=',d20.5,/, X + ' cb=',d20.5,/,' dmax=',d20.5) X if (dmax.eq.0.0) then Xc singular pencil X info=i X return X endif X dmin=min(ca,cb) X d=dmax*sqrt(1+(dmin/dmax)**2) X alpha = a(i,i)/d X beta = b(i,i)/d X if (ldebug) write(outunit,102) dmin,d,alpha,beta X102 format(' dmin=',d20.5,/,' d=',d20.5,/,' alpha=',2d20.5,/, X + ' beta=',2d20.5) Xc Xc compute right eigenvector X if (i.ne.1) then Xc Xc solve zz(1:i-1,1:i-1) * x = -zz(1:i-1,i) for Xc x = vec(1:i-1,i) X diag=beta*a(i-1,i-1) - alpha*b(i-1,i-1) X im1=i-1 X if (ldebug) write(outunit,103) im1,i,diag X103 format(' i,j,diag=',2i4,2d20.5) X if (abs(diag).eq.0.0) then Xc multiple eigenvalue X info=i-1 X return X endif X vec(i-1,i)=-(beta*a(i-1,i)-alpha*b(i-1,i))/diag X if (i.ne.2) then X im1=i-1 X im2=i-2 X do 2 j=im2,1,-1 X diag=beta*a(j,j)-alpha*b(j,j) X if (ldebug) write(outunit,103) j,i,diag X if (abs(diag).eq.0.0) then Xc multiple eigenvalue X info=j X return X endif X csum=-(beta*a(j,i)-alpha*b(j,i)) X jp1=j+1 X do 3 k=jp1,im1 X cmul=beta*a(j,k)-alpha*b(j,k) X csum=csum-cmul*vec(k,i) X3 continue X vec(j,i)=csum/diag X2 continue X endif X endif Xc Xc compute left eigenvector X if (i.ne.n) then Xc solve xt * zz(i+1:n,i+1:n) = -zz(i,i+1:n) for Xc x = vec(i+1:n,i) X diag=beta*a(i+1,i+1)-alpha*b(i+1,i+1) X ip1=i+1 X if (ldebug) write(outunit,103) i,ip1,diag X if (abs(diag).eq.0.0) then Xc multiple eigenvalue X info=i X return X endif X vec(i+1,i)=-(beta*a(i,i+1)-alpha*b(i,i+1))/diag X if (i.ne.nm1) then X ip1=i+1 X ip2=i+2 X do 4 j=ip2,n X diag=beta*a(j,j)-alpha*b(j,j) X if (ldebug) write(outunit,103) i,j,diag X if (abs(diag).eq.0.0) then Xc multiple eigenvalue X info=i X return X endif X csum=-(beta*a(i,j)-alpha*b(i,j)) X jm1=j-1 X do 5 k=ip1,jm1 X cmul=beta*a(k,j)-alpha*b(k,j) X csum=csum-cmul*vec(k,i) X5 continue X vec(j,i)=csum/diag X4 continue X endif X endif X1 continue Xc Xc compute condition numbers X do 6 i=1,n X gvcond(i)=dznrm2(i,vec(1,i),1)*dznrm2(n-i+1,vec(i,i),1) X6 continue X return X end Xc X subroutine pbound(a,b,ldab,m,n,rowred,colred,delta,difl,difu, X + qnorm,pnorm,pdelta,lbndup,rbndup,lbndlw,rbndlw,scase,work, X + ierr) Xc Xc implicit none Xc Xc**** formal parameter declarations X integer ldab,m,n,rowred,colred,ierr,scase X complex*16 a(ldab,*),b(ldab,*),work(*) X real*8 delta,difl,difu,qnorm,pnorm,pdelta,lbndup,rbndup X real*8 lbndlw, rbndlw Xc Xc******************************************************************* Xc Xc compute perturbation bounds for reducing subspaces of Xc singular pencil a - lambda b Xc assume a - lambda b has been reduced to generalized upper Xc triangular form by guptri Xc need rowred .le. colred and n-colred .le. m-rowred Xc as implied by generalized upper triangular form Xc Xc there are 4 cases, depending on dimension: Xc Xc case 1: 0 .lt. rowred and 0 .lt. n-colred so that Xc both left and right reducing subspaces nontrivial Xc Xc case 2: if rowred=0 and 0 .lt. colred .lt. n then left reducing Xc subspace 0 but right one nontrivial and bounds exist for it Xc Xc case 3: if colred=n and 0 .lt. rowred .lt. m then right reducing Xc subspace is entire space but left one nontrivial with bounds Xc Xc case 4: if ( (rowred=0 and colred=0) or Xc (rowred=0 and colred=n) or Xc (rowred=m and colred=n) ) then Xc both left and right subspaces trivial Xc Xc inputs: Xc Xc a(ldab,n),b(ldab,n) - complex*16 - m by n matrices Xc Xc ldab - integer - leading dimension of a and b Xc Xc m,n - integer - dimensions of a and b Xc Xc rowred,colred - integer - number of rows and columns in Xc (1,1) position of a,b. dimensions of desired left Xc and right reducing subspaces Xc Xc delta - real*8 - distance of perturbed pencil from a - lambda b Xc Xc idbg(9) - integer - if idbg(9) ne 0, print debug output Xc Xc outputs: (described in more detail in Xc 'accurate solutions of ill-posed problems in control theory' Xc 25th conference on decision and control, Xc j. demmel and b. kagstrom Xc Xc difl - real*8 - difl function (in case 4, difl=0) Xc Xc difu - real*8 - difu function (in case 4, difu=0) Xc Xc qnorm - real*8 - right projector norm ( sqrt(r0**2+1) ) Xc (in case 4, qnorm=1.) Xc Xc pnorm - real*8 - left projector norm ( sqrt(l0**2+1) ) Xc (in case 4, prnorm=1.) Xc Xc pdelta - real*8 - radius of ball around a - lambda b within Xc which perturbation bounds hold (in case 4, pdelta=-1. Xc to show pdelta does not apply). if delta.ge.pdelta, Xc the following bounds are set to -1. the following Xc outputs are given in terms of relerr = delta/pdelta Xc Xc lbndup - real*8 - upper bound on angular perturbation in left Xc reducing subspace (case 1 of theorem 4 of above paper) Xc in case 1: Xc lbndup=atan(relerr/(pnorm-relerr*sqrt(pnorm**2-1))) Xc in case 2: Xc lbndup=0 Xc in case 3: Xc lbndup=atan(relerr/(1-relerr)) Xc in case 4: Xc lbndup=0 Xc Xc rbndup - real*8 - upper bound on angular perturbation in right Xc reducing subspace (case 1 of theorem 4) Xc in case 1: Xc rbndup=atan(relerr/(qnorm-relerr*sqrt(qnorm**2-1))) Xc in case 2: Xc rbndup=atan(relerr/(1-relerr)) Xc in case 3: Xc rbndup=0 Xc in case 4: Xc rbndup=0 Xc Xc lbndlw - real*8 - lower bound on angular perturbation in left Xc reducing subspace (case 2 of theorem 4) Xc in case 1: Xc lbndlw=atan(1/(sqrt(2*min(rowred,m-rowred))*pnorm + Xc sqrt(pnorm**2-1))) Xc in case 2: lbndlw=-1 since this bound does not apply Xc in case 3: lbndlw=-1 since this bound does not apply Xc in case 4: lbndlw=-1 since this bound does not apply Xc Xc rbndlw - real*8 - lower bound on angular perturbation in right Xc reducing subspace (case 2 of theorem 4) Xc in case 1: Xc rbndlw=atan(1/(sqrt(2*min(colred,n-colred))*qnorm + Xc sqrt(qnorm**2-1))) Xc in case 2: rbndlw=-1 since this bound does not apply Xc in case 3: rbndlw=-1 since this bound does not apply Xc in case 4: rbndlw=-1 since this bound does not apply Xc Xc scase - integer - 1, 2, 3 or 4 as described above Xc Xc ierr - integer - error flag Xc 0 means no error (normal return) Xc 1 means error in svd of difu Xc 2 means difu = 0 Xc 3 means error in svd of difl Xc 4 means difl = 0 Xc 5 means bad rowred or colred Xc Xc work space Xc work - complex*16 - array of length at least Xc max ( rowdfu*coldfu+coldfu**2+2*coldfu+rowdfu , Xc rowdfl*coldfl+2*coldfl+rowdfl ) Xc where Xc rowdfu=coldfl=colred*(n-colred)+rowred*(m-rowred) Xc coldfu=2*(n-colred)*rowred Xc rowdfl=2*(m-rowred)*colred Xc Xc********************************************************************* Xc Xc**** this version dated 16 june 1987 Xc authors: jim demmel, courant institute, 251 mercer str, new york, Xc new york, 10012 Xc electronic address: demmel at nyu.edu Xc bo kagstrom, institute of information processing, Xc university of umea, s-90187 umea, sweden Xc electronic address: bokg at seumdc51.bitnet Xc Xc**** pbound uses the following subroutines and functions Xc dznrm2, blddfu, blddfl, bldrhs, prml, prmlct, svdiv, zsvdc Xc Xc**** internal variables Xc X complex*16 dummy X integer rowdfu,coldfu,sstrt,wstrt,estrt,rowdfl,coldfl,vstrt X integer isub, i, j, info, len X real*8 r0, l0, relerr, dznrm2 Xc X ierr=0 X if ((rowred.gt.colred).or.((n-colred).gt.(m-rowred))) then Xc inconsistent dimensions X ierr = 5 X elseif ((0.lt.rowred) .and. (0.lt.n-colred)) then Xc case 1 X scase = 1 Xc compute difu Xc build transposed difu matrix starting at work(1) Xc rowdfu = number of rows in difut X rowdfu=colred*(n-colred)+rowred*(m-rowred) Xc coldfu = number of columns in difut X coldfu=2*(n-colred)*rowred Xc X call blddfu(work,rowdfu,a,b,ldab,m,n,rowred,colred) Xc Xc setup workspace for svd Xc store left singular vectors u over difu starting at work(1) X sstrt=1+rowdfu*coldfu Xc store singular values starting at work(sstrt) X wstrt=sstrt+coldfu Xc store work array needed for svd starting at work(wstrt) X estrt=wstrt+rowdfu Xc store e array needed for svd starting at work(estrt) X vstrt=estrt+coldfu Xc store right singular vectors v starting at work(vstrt) Xc Xc compute svd X call zsvdc(work(1),rowdfu,rowdfu,coldfu,work(sstrt), X + work(estrt),work(1),rowdfu,work(vstrt),coldfu,work(wstrt), X + 21,info) Xc X if (info.eq.0) goto 10 X ierr=1 X return X10 continue Xc Xc extract difu X difu=dreal(work(sstrt-1+coldfu)) Xc X if (difu.gt.0.) goto 20 X ierr=2 X return X20 continue Xc Xc compute pnorm, qnorm Xc build rhs = (-col a12, -col b12) starting at work(wstrt) X call bldrhs(work(wstrt),a,b,ldab,m,n,rowred,colred) Xc Xc solve underdetermined least squares problem Xc premultiply rhs by v* storing result at work(estrt) X call prmlct(work(vstrt),coldfu,coldfu,coldfu, X + work(wstrt),work(estrt)) Xc Xc premultiply by inverted singular values X call svdiv(work(estrt),coldfu,work(sstrt)) Xc Xc premultiply by u storing result at work(wstrt) X call prml(work,rowdfu,rowdfu,coldfu,work(estrt),work(wstrt)) Xc X len=colred*(n-colred) Xc compute r0 = norm of leading len components X r0=dznrm2(len,work(wstrt),1) Xc Xc compute l0 = norm of remaining components X len=rowred*(m-rowred) X l0=dznrm2(len,work(wstrt+len),1) Xc compute pnorm, qnorm from l0, r0 X pnorm=sqrt(1+l0**2) X qnorm=sqrt(1+r0**2) Xc Xc compute difl Xc build difl matrix starting at work(1) Xc rowdfl = number of rows in difl X rowdfl=2*colred*(m-rowred) Xc coldfl=number of columns in difl X coldfl=rowred*(m-rowred)+colred*(n-colred) X call blddfl(work,rowdfl,a,b,ldab,m,n,rowred,colred) Xc Xc setup workspace for svd Xc do not compute any singular vectors X sstrt=1+rowdfl*coldfl Xc store singular values starting at work(sstrt) X wstrt=sstrt+coldfl Xc store work array needed by svd starting at work(wstrt) X estrt=wstrt+rowdfl Xc store e array needed by svd starting at work(estrt) Xc X call zsvdc(work(1),rowdfl,rowdfl,coldfl,work(sstrt), X + work(estrt),dummy,1,dummy,1,work(wstrt),0,info) Xc X if (info.eq.0) goto 30 X ierr=3 X return X30 continue Xc Xc extract difl X difl=dreal(work(sstrt-1+coldfl)) X if (difl.gt.0.) goto 40 X ierr=4 X return X40 continue Xc compute perturbation bounds X pdelta=min(difl,difu)/(sqrt(pnorm**2+qnorm**2)+ X + 2.*max(pnorm,qnorm)) X relerr=delta/pdelta X lbndup=-1. X rbndup=-1. X lbndlw=-1. X rbndlw=-1. X if (relerr.ge.1.) goto 50 X lbndup=atan(relerr/(pnorm-relerr*sqrt(pnorm**2-1.))) X rbndup=atan(relerr/(qnorm-relerr*sqrt(qnorm**2-1.))) X lbndlw=atan(1./(sqrt(2.*min(rowred,m-rowred))*pnorm+ X + sqrt(pnorm**2-1.))) X rbndlw=atan(1./(sqrt(2.*min(colred,n-colred))*qnorm+ X + sqrt(qnorm**2-1.))) X50 continue X elseif (rowred.eq.0.and.colred.gt.0.and.colred.lt.n) then Xc case 2 X scase = 2 Xc compute difl Xc build difl matrix ( (a**t b**t)**t ) starting at work(1) X isub = 0 X do 100 j=colred+1, n X do 101 i=1, m X isub = isub +1 X work(isub) = a(i,j) X101 continue X do 102 i=1,m X isub = isub +1 X work(isub) = b(i,j) X102 continue X100 continue Xc compute singular values X sstrt=1+isub X estrt=sstrt + n-colred X wstrt=estrt + n-colred X call zsvdc(work,2*m,2*m,n-colred,work(sstrt),work(estrt), X + dummy,1,dummy,1,work(wstrt),0,info) X if (info.ne.0) then X ierr=3 X return X endif Xc extract difl X difl = abs(work(sstrt+n-colred-1)) X difu=difl X if (difl.eq.0.) then X ierr=4 X return X endif X pdelta=difl X relerr=delta/pdelta X pnorm = 1. X qnorm = 1. X lbndlw = -1. X rbndlw = -1. X lbndup = -1. X rbndup = -1. X if (relerr.lt.1.) then X lbndup = 0. X rbndup = atan(relerr/(1.-relerr)) X endif X elseif (colred.eq.n.and.rowred.gt.0.and.rowred.lt.m) then Xc case 3 X scase = 3 Xc compute difu Xc build difu matrix (a,b) starting at work(1) X isub = 0 X do 104 j=1,n X do 105 i=1,rowred X isub = isub +1 X work(isub) = a(i,j) X105 continue X104 continue X do 106 j=1,n X do 107 i=1,rowred X isub = isub +1 X work(isub) = b(i,j) X107 continue X106 continue Xc compute singular values X sstrt=isub+1 X estrt=sstrt+rowred+1 X wstrt=estrt+2*n X call zsvdc(work,rowred,rowred,2*n,work(sstrt),work(estrt), X + dummy,1,dummy,1,work(wstrt),0,info) X if (info.ne.0) then X ierr = 1 X return X endif Xc extract difu X difu=abs(work(sstrt+rowred-1)) X difl = difu X if (difu.eq.0.0) then X ierr = 2 X return X endif X pdelta = difu X relerr = delta/pdelta X pnorm = 1. X qnorm = 1. X lbndup = -1. X rbndup = -1. X lbndlw = -1. X rbndlw = -1. X if ( relerr.lt.1.0) then X rbndup = 0. X lbndup = atan(relerr/(1.-relerr)) X endif X else Xc both left and right subspace trivial X scase = 4 X lbndup = 0. X rbndup = 0. X lbndlw = -1. X rbndlw = -1. X difl = 0. X difu = 0. X pdelta = -1. X pnorm = 1. X qnorm = 1. X endif X return X end Xc Xc X subroutine blddfl(work,wrow,a,b,ldab,m,n,rowred,colred) Xc implicit none Xc**** formal parameter declarations X integer ldab, m, n, rowred, colred, wrow X complex*16 work(wrow,*),a(ldab,*),b(ldab,*) Xc Xc*************************************************************** Xc Xc build difl matrix in work Xc in matlab notation Xc Xc difl matrix = < ; Xc > Xc Xc where a11 = a(1:rowred , 1:colred) Xc a22 = a(rowred+1 : m , colred+1 : n) Xc b11 = b(1:rowred , 1:colred) Xc b22 = b(rowred+1 : m , colred+1 : n) Xc Xc*************************************************************** Xc Xc**** this version dated 16 june 1987 Xc authors: jim demmel and bo kagstrom Xc Xc**** internal variables X integer wcol,rstrta,rstrtb,cstrt,cnt,i,j X integer row12,col1,col2,mmrwrd,nmclrd Xc Xc nmclrd = number of columns in (1,2), (2,2) blocks of a, b X nmclrd = n-colred Xc mmrwrd = number of rows in (2,1), (2,2) blocks of a, b X mmrwrd = m-rowred Xc row12 = numbers of rows in each subblock of difl matrix X row12 = colred*mmrwrd Xc col1 = number of columns in (1,1), (2,1) blocks of difl X col1 = rowred*mmrwrd Xc col2 = number of columns in (1,2), (2,2) blocks of difl X col2 = colred*nmclrd Xc wcol = total number of columns in difl X wcol = col1+col2 Xc Xc zero out difl X do 10 j=1,wcol X do 11 i=1,wrow X work(i,j)=0. X11 continue X10 continue Xc Xc fill in (1,1), (2,1) blocks of difl X rstrta=0 X rstrtb=row12 X cstrt=0 X do 1 j=1,colred X do 2 i=1,rowred X do 3 cnt=1,mmrwrd X work(cnt+rstrta,cnt+cstrt)=a(i,j) X work(cnt+rstrtb,cnt+cstrt)=b(i,j) X3 continue X cstrt=cstrt+mmrwrd X2 continue X cstrt=0 X rstrta=rstrta+mmrwrd X rstrtb=rstrta+row12 X1 continue Xc Xc fill in (1,2), (2,2) blocks of difl X rstrta=0 X cstrt=col1 X do 4 cnt=1,colred X rstrtb=rstrta+row12 X do 5 j=1,nmclrd X do 6 i=1,mmrwrd X work(rstrta+i,cstrt+j)=-a(i+rowred,j+colred) X work(rstrtb+i,cstrt+j)=-b(i+rowred,j+colred) X6 continue X5 continue X rstrta=rstrta+mmrwrd X cstrt=cstrt+nmclrd X4 continue X return X end Xc Xc X subroutine blddfu(work,wrow,a,b,ldab,m,n,rowred,colred) Xc implicit none Xc**** formal parameter declarations X integer ldab, m, n, rowred, colred, wrow X complex*16 work(wrow,*),a(ldab,*),b(ldab,*) Xc********************************************************************* Xc Xc build conjugate transpose difu matrix in work Xc in matlab notation Xc Xc (difu matrix)' = Xc Xc < < eye(n-colred) .*. a11' , eye(n-colred) .*. b11' >; Xc < -conj(a22) .*. eye(rowred) , -conj(b22) .*. eye(rowred) >> Xc Xc where a11 = a(1:rowred , 1:colred) Xc a22 = a(rowred+1 : m , colred+1 : n) Xc b11 = b(1:rowred , 1:colred) Xc b22 = b(rowred+1 : m , colred+1 : n) Xc Xc********************************************************************* Xc Xc**** this version dated 16 june 1987 Xc authors: jim demmel and bo kagstrom Xc Xc**** internal variables Xc X integer wcol,cstrta,cstrtb,rstrt,cnt,i,j X integer mmrwrd,nmclrd,rwrdp1,clrdp1 X integer row1, row2, col12 Xc Xc nmclrd = number of columns in (1,2), (2,2) entries of a, b X nmclrd=n-colred Xc col12 = number of columns in each subblock of difuct matrix X col12=rowred*nmclrd Xc mmrwrd = number of rows in (2,1), (2,2) entries of a, b X mmrwrd = m-rowred Xc row1 = number of rows in (1,1), (2,1) sublocks of difu X row1 = colred*nmclrd Xc row2 = number of rows in (1,2), (2,2) subblocks of difu X row2 = rowred*mmrwrd Xc wcol = total number of columns in difu matrix X wcol = 2*col12 Xc initialize difu to zero X do 1 j=1,wcol X do 2 i=1,wrow X work(i,j)=0. X2 continue X1 continue Xc Xc fill in (1,1), (1,2) positions of difu X cstrta=0 X rstrt=0 X do 3 cnt=1,nmclrd X cstrtb=cstrta+col12 X do 4 j=1,colred X do 5 i=1,rowred X work(rstrt+j,cstrta+i)=conjg(a(i,j)) X work(rstrt+j,cstrtb+i)=conjg(b(i,j)) X5 continue X4 continue X cstrta=cstrta+rowred X rstrt=rstrt+colred X3 continue Xc Xc fill in (2,1), (2,2) positions of difuct X rwrdp1=rowred+1 X clrdp1=colred+1 X cstrta=0 X cstrtb=col12 X rstrt=row1 X do 6 j=clrdp1,n X do 7 i=rwrdp1,m X do 8 cnt=1,rowred X work(cnt+rstrt,cnt+cstrta)=-conjg(a(i,j)) X work(cnt+rstrt,cnt+cstrtb)=-conjg(b(i,j)) X8 continue X rstrt=rstrt+rowred X7 continue X rstrt=row1 X cstrta=cstrta+rowred X cstrtb=cstrta+col12 X6 continue X return X end Xc Xc X subroutine bldrhs(work,a,b,ldab,m,n,rowred,colred) Xc implicit none Xc**** formal parameter declarations X integer ldab, m, n, rowred, colred X complex*16 work(*), a(ldab,*), b(ldab,*) Xc Xc********************************************************************* Xc Xc extract a12 = (1,2) block of a and b12 = (1,2) block of b Xc and store columnwise in work=(-col a12, -col b12) Xc Xc********************************************************************* Xc Xc**** this version dated 16 june 1987 Xc authors: jim demmel and bo kagstrom Xc Xc**** internal variables X integer clrdp1, j, i, loc Xc X clrdp1=colred+1 X loc=0 X do 1 j=clrdp1,n X do 2 i=1,rowred X loc=loc+1 X work(loc)=-a(i,j) X2 continue X1 continue X do 3 j=clrdp1,n X do 4 i=1,rowred X loc=loc+1 X work(loc)=-b(i,j) X4 continue X3 continue X return X end Xc Xc X subroutine prml(u,ldu,m,n,rhs,prod) Xc implicit none X integer ldu, m, n X complex*16 u(ldu,n),rhs(n),prod(m) Xc Xc********************************************************************* Xc compute prod = u * rhs Xc Xc**** this version dated 16 june 1987 Xc authors: jim demmel and bo kagstrom Xc X integer i, j Xc X do 1 j=1,m X prod(j)=rhs(1)*u(j,1) X1 continue X if (n.eq.1) return X do 2 i=2,n X call zaxpy(m,rhs(i),u(1,i),1,prod,1) X2 continue X return X end Xc Xc X subroutine prmlct(u,ldu,m,n,rhs,prod) Xc implicit none X integer ldu, m, n X complex*16 u(ldu,n),rhs(m),prod(n),zdotc Xc Xc********************************************************************* Xc compute prod = (conjugate transpose u) * rhs Xc Xc**** this version dated 16 june 1987 Xc authors: jim demmel and bo kagstrom Xc X integer j Xc X do 1 j=1,n X prod(j)=zdotc(m,u(1,j),1,rhs,1) X1 continue X return X end Xc Xc X subroutine svdiv(z,n,s) Xc implicit none X integer n X complex*16 z(n),s(n) Xc Xc********************************************************************* Xc divide one array by another Xc Xc**** this version dated 16 june 1987 Xc authors: jim demmel and bo kagstrom Xc X integer j Xc X do 1 j=1,n X z(j)=z(j)/s(j) X1 continue X return X end Xc X subroutine evalbd(delta, sdlmax, qnorm, pnorm, scase, X + m, n, irstrt, icstrt, X + lbndup, rbndup, lbndlw, rbndlw) Xc Xc implicit none Xc**** formal parameter declarations Xc X real*8 delta, sdlmax, qnorm, pnorm X real*8 lbndup, rbndup, lbndlw, rbndlw X integer scase, m, n, irstrt, icstrt Xc Xc****************************************************************** Xc Xc evaluate reducing subspace angular perturbation bounds computed Xc by subroutine bound for a perturbation of frobenius Xc norm delta. see documentation to subroutine bound for more details. Xc Xc inputs: Xc Xc sdlmax, qnorm, pnorm and scase are computed by bound. Xc m, n, irstrt and icstrt are dimensions also input to bound Xc in order to compute sdlmax, qnorm, pnorm and scase. Xc Xc outputs: Xc Xc lbndup - real*8 - upper bound on angular perturbation in Xc left reducing subspace Xc (0 if space trivial and -1 if inapplicable) Xc Xc rbndup - real*8 - upper bound on angular perturbation in Xc right reducing subspace Xc (0 if space trivial and -1 if inapplicable) Xc Xc lbndlw - real*8 - lower bound on angular perturbation in Xc left reducing subspace (-1 if inapplicable) Xc Xc rbndlw - real*8 - lower bound on angular perturbation in Xc right reducing subspace (-1 if inapplicable) Xc Xc************************************************************************ Xc Xc**** this version dated 16 june 87 Xc authors: jim demmel and bo kagstrom Xc Xc**** internal variables X real*8 relerr Xc X if (scase .ne. 4) relerr = delta/sdlmax X if (scase.eq.1) then X lbndup = atan(relerr/(pnorm-relerr*sqrt(pnorm**2-1.))) X rbndup = atan(relerr/(qnorm-relerr*sqrt(qnorm**2-1.))) X lbndlw = atan(1./(sqrt(2.*min(irstrt-1,m-irstrt+1))*pnorm + X + sqrt(pnorm**2-1.))) X rbndlw = atan(1./(sqrt(2.*min(icstrt-1,n-icstrt+1))*qnorm + X + sqrt(qnorm**2-1.))) X elseif (scase.eq.2) then X lbndup = 0. X rbndup = atan(relerr/(1.-relerr)) X lbndlw = -1. X rbndlw = -1. X elseif (scase.eq.3) then X lbndup = atan(relerr/(1.-relerr)) X rbndup = 0. X lbndlw = -1. X rbndlw = -1. X elseif (scase.eq.4) then X lbndup = 0. X rbndup = 0. X lbndlw = -1. X rbndlw = -1. X endif X return X end Xc X subroutine bndwsp(m,n,irstrt,icstrt,dimreg,ecase,space,info) Xc Xc implicit none Xc Xc**** debug space X common /debug2/ idbg(20), outunit X integer idbg, outunit Xc Xc**** formal parameter declarations X integer m,n,irstrt,icstrt,dimreg,info,ecase,space Xc Xc******************************************************************** Xc Xc compute work space needed by subroutine bound Xc Xc inputs Xc Xc m,n - integer - row, column dimensions of a and b Xc Xc irstrt, icstrt - integer - starting row and column of selected Xc part of pencil for which eigenvalue bounds Xc are desired. reducing subspace bounds will be Xc supplied for right reducing subspace spanned Xc by leading icstrt-1 components and for left Xc reducing subspace spanned by leading icstrt-1 Xc components. Xc note: set icstrt=n+1 to make right reducing Xc subspace whole space Xc set irstrt=m+1 to make left reducing Xc subspace whole space Xc Xc dimreg - integer - number of selected eigenvalues; Xc if dimreg.eq.0 only subspace perturbation bounds will be Xc computed Xc (note - one can select a subset of the regular part only; Xc this gives generally different bounds for common eigenvalues Xc from a different selected subset; see paper above for Xc discussion) Xc Xc outputs Xc Xc ecase - integer - which of 5 cases for eigenvalue bounds Xc the pencil falls depending on input dimensions; Xc the first four cases are for dimreg.gt.0, in which Xc case the description gives: Xc (part of KCF to above, left of selected part) and Xc (part of KCF to below, right of selected part) Xc ecase=1 - (right singular and/or regular part) and Xc (left singular and/or regular part) Xc ecase=2 - (right singular and/or regular part) and (nothing) Xc ecase=3 - (nothing) and (left singular and/or regular part) Xc ecase=4 - (nothing) and (nothing) Xc ecase=5 - dimreg.eq.0 (no eigenvalue bounds) Xc Xc space - integer - amount of workspace (double precision complex Xc words) needed by subroutine bound Xc (the following simple expression bounds the workspace also, but Xc may occasionally be much too large (especially if ecase=4): Xc workspace .le. 2*m*n* (n*n + m*m + 2*n + m + 2) + n*n + m*m) Xc Xc info - integer - 0 if normal return Xc 1 if inconsistent input dimensions Xc Xc************************************************************************* Xc Xc**** this version dated 22 june 1987 Xc authors: jim demmel, courant institute, 251 mercer str, Xc new york, new york, 10012 Xc electronic address: demmel at nyu.edu Xc bo kagstrom, institute of information processing, Xc university of umea, s-90187 umea, sweden Xc electronic address: bokg at seumdc51.bitnet Xc Xc**** internal variables X integer irend,icend,m11,m21,m12,m22,n11,n12,n21,n22 Xc Xc test input dimensions for consistency X info = 0 X icend = icstrt+dimreg-1 X irend = irstrt+dimreg-1 X if (irstrt.gt.icstrt .or. irstrt.le.0 .or. X + n-icstrt-dimreg.gt.m-irstrt-dimreg .or. X + n-icstrt-dimreg+1.lt.0 .or. dimreg.lt.0) then Xc inconsistent input dimensions X info = 1 X else X if (dimreg.gt.0) then Xc there are eigenvalue bounds to compute Xc Xc ecase 1 - in addition to selected regular part KCF has Xc (right singular part and/or regular part) and Xc (left singular part and/or regular part) X if (icstrt.ne.1 .and. irend.ne.m) then X ecase = 1 X endif Xc Xc ecase 2 - in addition to selected regular part KCF has Xc (right singular part and/or regular part) and Xc (nothing) X if (icstrt.ne.1 .and. irend.eq.m) then X ecase=2 X endif Xc Xc ecase 3 - in addition to selected regular part KCF has Xc (nothing) and Xc (left singular part and/or regular part) X if (icstrt.eq.1 .and. irend.ne.m) then X ecase = 3 X endif Xc Xc ecase 4 - pencil regular and entire spectrum selected X if (icstrt.eq.1 .and. irend.eq.m) then X ecase=4 X endif Xc X else Xc dimreg.eq.0, so only compute subspace bounds X ecase = 5 X endif Xc X if (ecase .eq. 1) then X m11=irstrt-1 X m21=m-m11 X n11=icstrt-1 X n21=n-n11 X m12=irend-irstrt+1 X m22=m-irend X n12=icend-icstrt+1 X n22=n-icend X space = max( (2*n21*m11*(n11*n21+m11*m21+ X + 2*n21*m11+2)+n11*n21+m11*m21) , X + (2*((m21*n11+1)*(n11*n21+ X + m11*m21+1)-1)) , X + (2*n22*m12*(n12*n22+m12*m22+ X + 2*n22*m12+2)+n12*n22+m12*m22) , X + (2*((m22*n12+1)*(n12*n22+ X + m12*m22+1)-1)) ) X elseif (ecase .eq. 2 .or. ecase .eq. 5) then X m11=irstrt-1 X m21=m-m11 X n11=icstrt-1 X n21=n-n11 X space = max( (2*n21*m11*(n11*n21+m11*m21+ X + 2*n21*m11+2)+n11*n21+m11*m21) , X + (2*((m21*n11+1)*(n11*n21+ X + m11*m21+1)-1)) ) X elseif (ecase .eq. 3) then X m11=irend X m21=m-m11 X n11=icend X n21=n-icend X space = max( (2*n21*m11*(n11*n21+m11*m21+ X + 2*n21*m11+2)+n11*n21+m11*m21) , X + (2*((m21*n11+1)*(n11*n21+ X + m11*m21+1)-1)) ) X elseif (ecase .eq. 4) then X space = n*n X endif X endif Xc X if (idbg(19).ne.0) then X write(outunit,100) m,n,irstrt,icstrt,dimreg,ecase, X + space,info X100 format(' bndwsp - m,n,irstrt,icstrt,dimreg' X + ',ecase,space,info=',/,8i5) X endif X return X end END_OF_zbnd.f if test 52873 -ne `wc -c zcmatmlr.f <<'END_OF_zcmatmlr.f' X Xc on this file june 7, 1987: cmatml, cmatmr Xc X subroutine cmatml(a,lda,rowa,cola,b,ldb,rowb,c,ldc,work,job) Xc Xc implicit none X integer lda,rowa,cola,ldb,rowb,ldc,job X complex*16 a(lda,lda),b(ldb,ldb),c(ldc,ldc),work(*) X complex*16 zdotu,zdotc Xc Xc*********************************************************************** Xc Xc cmatml performs a complex (left) matrix multiplication b * a, Xc or b' * a (' = transpose ,conjugate) where a is rowa * cola, Xc b is rowb * rowa. the result is stored in c or overwritten in a. Xc note the extra restrictions on dimensions of b when job = 3 or 4. Xc Xc on entry Xc Xc a complex(lda,cola), where lda>=rowa. Xc Xc lda integer Xc lda is the leading dimension of the array a. Xc Xc rowa integer Xc rowa is the number of rows of a, which is also Xc the number of columns of b. Xc cola integer Xc cola is the number of columns of a, which is also Xc the number columns of the resulting matrix. Xc Xc b complex(ldb,rowa), ldb>=rowb. Xc Xc ldb integer Xc ldb is the leading dimension of the array b. Xc Xc rowb integer Xc rowb is the number of rows of the array, which Xc is also the number of rows of the resulting matrix. Xc Xc ldc integer Xc ldc is the leading dimension of the array c Xc Xc work complex(rowa) Xc work is a scratch array. Xc Xc job integer Xc job controls the matrix multiplication, and has Xc the following meaning Xc job=1 a = b * a Xc job=2 c = b * a Xc job=3 a = b' * a Xc job=4 c = b' * a Xc Xc on return Xc Xc c complex(ldc,cola), where ldc>=rowb. Xc c is the matrix product of a and b. if rowa (=colb) Xc = rowb then it is possible to call cmatml with c Xc equals to a, and the result is overwritten in a. Xc Xc********************************************************************* Xc Xc this version dated june 7, 1987 Xc authors: jim demmel and bo kagstrom Xc Xc***** internal variables Xc X integer i,j Xc Xc***** cmatml uses the following functions and subroutines Xc Xc blas zcopy, zdotc, zdotu Xc Xc***** determine what is to be computed via nested if-then -else's Xc X do 20 j = 1, cola X do 10 i = 1, rowb X if (job .eq. 1) then X work(i) = zdotu(rowa,b(i,1),ldb,a(1,j),1) X elseif (job .eq. 2) then X c(i,j) = zdotu(rowa,b(i,1),ldb,a(1,j),1) X elseif (job .eq. 3) then X work(i) = zdotc(rowa,b(1,i),1,a(1,j),1) X else Xc (job .eq. 4) X c(i,j) = zdotc(rowa,b(1,i),1,a(1,j),1) X endif X 10 continue X if (job .eq. 1 .or. job .eq. 3) then X call zcopy(rowa,work,1,a(1,j),1) X endif X 20 continue X return X end X X X subroutine cmatmr(a,lda,rowa,cola,b,ldb,colb,c,ldc,work,job) Xc Xc implicit none X integer lda,rowa,cola,ldb,colb,ldc,job X complex*16 a(lda,lda),b(ldb,ldb),c(ldc,ldc),work(*) X complex*16 zdotu,zdotc Xc Xc*********************************************************************** Xc Xc cmatmr performs a complex (right) matrix multiplication a * b, Xc or a * b' ,(' = transpose ,conjugate), where a is rowa * cola, Xc b is cola * colb. the result is stored in c or overwritten in a. Xc note the extra restrictions in dimension of b when job = 3 or 4. Xc Xc on entry Xc Xc a complex(lda,cola), where lda>=rowa. Xc Xc lda integer Xc lda is the leading dimension of the array a. Xc Xc rowa integer Xc rowa is the number of rows of a, which is also Xc the number of rows in the resulting matrix. Xc cola integer Xc cola is the number of columns of a, which is also Xc the number of rows of b. Xc Xc b complex(ldb,colb), ldb>=cola. Xc Xc ldb integer Xc ldb is the leading dimension of the array b. Xc Xc colb integer Xc colb is the number of columns of b, which is Xc also the number of columns of the resulting matrix Xc Xc ldc integer Xc ldc is the leading dimension of the array c Xc Xc work complex(cola) Xc work is a scratch array. Xc Xc job integer Xc job controls the matrix multiplication, and has Xc the following meaning Xc job=1 a = a * b Xc job=2 c = a * b Xc job=3 a = a * b' Xc job=4 c = a * b' Xc Xc on return Xc Xc c complex(ldc,colb), where ldc>=rowa. Xc c is the matrix product of a and b. if cola(=rowb) Xc = colb then it is possible to call cmatmr with c Xc equals to a, and the result is overwritten in a. Xc Xc********************************************************************* Xc Xc this version dated june 7, 1987 Xc authors: jim demmel and bo kagstrom Xc Xc***** internal variables Xc X integer i,j Xc Xc***** cmatmr uses the following functions and subroutines Xc Xc blas zcopy, zdotc, zdotu Xc Xc***** determine what is to be computed via nested if-then -else's Xc X do 20 i = 1, rowa X do 10 j = 1, colb X if (job .eq. 1) then X work(j) = zdotu(cola,a(i,1),lda,b(1,j),1) X else if (job .eq. 2) then X c(i,j) = zdotu(cola,a(i,1),lda,b(1,j),1) X else if (job .eq. 3) then X work(j) = zdotc(cola,b(j,1),ldb,a(i,1),lda) X else Xc (job .eq. 4) X c(i,j) = zdotc(cola,b(j,1),ldb,a(i,1),lda) X end if X 10 continue X if (job .eq. 1 .or. job .eq. 3) then X call zcopy(cola,work,1,a(i,1),lda) X end if X 20 continue X return X end X END_OF_zcmatmlr.f if test 6590 -ne `wc -c zftest1.f <<'END_OF_zftest1.f' X integer function ftest(alpha,beta) Xc Xc implicit none X complex*16 alpha, beta Xc Xc**** fout checks if the complex root alpha/beta lies outside Xc the unit disc Xcc if (abs(beta) .eq. 0.0 ) then Xcc ftest = 1 Xcc elseif (abs(alpha/beta) .lt. 1.0) then X ftest = -1 Xcc else Xcc ftest = 1 Xcc endif X return X end Xc X integer function ftestp(alpha,beta) Xc Xc implicit none X complex*16 alpha, beta Xc Xc**** fout checks if the complex root alpha/beta lies outside Xc the unit disc Xcc if (abs(beta) .eq. 0.0 ) then Xcc ftestp = 1 Xcc elseif (abs(alpha/beta) .lt. 1.0) then X ftestp = -1 Xcc else Xcc ftestp = 1 Xcc endif X return X end END_OF_zftest1.f if test 752 -ne `wc -c zgschur.c1 <<'END_OF_zgschur.c1' XTestrun identification: April 1, 1990 Example C1 - zgschur X X X lda= 30 m= 4 n= 5 X final version of a input X ---------------------------------------------------------------------- X 0.10000000000000000D+01-0.20000000000000000D+01 0.00000000000000000D+00 X 0.10000000000000000D+01 0.00000000000000000D+00-0.10000000000000000D+01 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X X X 0.00000000000000000D+00 0.00000000000000000D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 X 0.10000000000000000D+01 0.00000000000000000D+00 X 0.00000000000000000D+00 0.20000000000000000D+01 X X X lda= 30 m= 4 n= 5 X final version of b input X ---------------------------------------------------------------------- X 0.00000000000000000D+00 0.10000000000000000D+01 0.00000000000000000D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 0.10000000000000000D+01 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X X X 0.00000000000000000D+00 0.00000000000000000D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 X 0.10000000000000000D+01 0.00000000000000000D+00 X 0.00000000000000000D+00 0.10000000000000000D+01 X X X debug controls - X 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 X 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 X input: epsu= 0.10000D-07 X gap= 0.10000D+04 X epsper = 0.10000D-09 numex = 1 numtst= 1 jobper = 3 X epsbnd = X 0.100D-09 X zero=T X nostat=T X epsper= 0.10000D-09 X norm(a,e)= 0.34641D+01 norm(b,e)= 0.20000D+01 X start guptri X guptri - workspace for rzstr - 4 5 4 5 1 X 6 26 31 61 36 41 86 46 51 56 X X Xguptri - m,n,epsu= 4 5 0.100000D-07 X Xreduction 1 X Xkstr, last= 3 X 1 2 3 X 1 1 1 X 1 1 0 Xaccumulated perturbations in a,b = 0.000000D+00 0.000000D+00 X Xreduction 4, kfirst= 5 X 1 2 3 4 5 X 1 1 1 -1 0 X 1 1 0 -1 0 Xaccumulated perturbations in a,b = 0.000000D+00 0.000000D+00 X X Xfinal kstr= X 1 2 3 4 5 6 7 X 1 1 1 -1 -1 2 -1 X 1 1 0 -1 -1 2 -1 X nsumrz= 3 X rsumrz= 2 X djordz= 0 X nsumli= 0 X rsumli= 0 X djordi= 0 X dimreg= 2 X X Xfinal pstruc= 3 3 3 3 Xfinal struc = X 1 1 1 X rtce= 3 zrce= 3 fnce= 5 ince= 5 X rtre= 2 zrre= 2 fnre= 4 inre= 4 X computed eigenvalues X eigenvalue= 0.10000D+01 0.00000D+00 X eigenvalue= 0.20000D+01 0.00000D+00 X eigenvalues before reordering X eigenvalue= 0.10000D+01 0.00000D+00 X eigenvalue= 0.20000D+01 0.00000D+00 X eigenvalues after reorder and X computed eigenvalues X eigenvalue= 0.10000D+01 0.00000D+00 X eigenvalue= 0.20000D+01 0.00000D+00 X results from guptri and reorder X rtce= 3 zrce= 3 fnce= 5 ince= 5 X rtre= 2 zrre= 2 fnre= 4 inre= 4 Xpstruc = 3 3 3 3 Xstruc = X 1 1 1 Xnsumrz= 3 Xrsumrz= 2 Xdjordz= 0 Xnsumli= 0 Xrsumli= 0 Xdjordi= 0 Xdimreg= 2 Xndim= 0 X Relative perturbation in a= 0.000000D+00 X Relative perturbation in b= 0.000000D+00 X Frobeniusnorm of deleted singular vaules=0.000000D+00 X kstr, step= 7 X 1 1 1 -1 -1 2 -1 X 1 1 0 -1 -1 2 -1 X lda= 30 m= 4 n= 5 X Transformed matrix A X ---------------------------------------------------------------------- X 0.00000000000000000D+00-0.15811388300841895D+01 0.94868329805051310D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 0.18973665961010271D+01 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X X X 0.00000000000000000D+00 0.35138622112200610D-15 X 0.00000000000000000D+00 0.37144952276981167D-15 X 0.10000000000000000D+01 0.00000000000000000D+00 X 0.00000000000000000D+00 0.19999999999999998D+01 X X X lda= 30 m= 4 n= 5 X Transformed matrix B X ---------------------------------------------------------------------- X -0.74535599249992979D+00 0.63245553203367566D+00-0.21081851067789167D+00 X 0.00000000000000000D+00-0.31622776601683794D+00-0.94868329805051332D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X X X 0.00000000000000000D+00-0.10753675309148582D-15 X 0.00000000000000000D+00-0.15259246943748575D-15 X 0.10000000000000000D+01 0.00000000000000000D+00 X 0.00000000000000000D+00 0.99999999999999989D+00 X X X cond(PP)=0.100000D+01 X cond(QQ)=0.100000D+01 X abs(a-acopy)=0.777156D-15 X relative dif.=0.224346D-15 X abs(b-bcopy)=0.416334D-15 X relative dif.=0.208167D-15 X fro(a-pp" * acopy * qq)=0.368219D-15 X relative fro for a-part=0.106296D-15 X fro(b-pp" * bcopy * qq)=0.209550D-15 X relative fro for b-part=0.104775D-15 X colrs= 3 X rowrs= 2 X len= 2 X bndwsp X ecase= 2 X space= 170 X info= 0 X icase= 1 X ecase= 2 X ierr= 0 X delmax= 0.33861D+00 X pdelta= 0.33861D+00 X difl= 0.11561D+01 X difu= 0.13095D+01 X qnorm= 0.10000D+01 X pnorm= 0.10000D+01 X pqnorm= 0.14142D+01 X dsvd= 0.54641D-07 X X results from pbound X difl= 0.11561D+01 X difu= 0.13095D+01 X qnorm= 0.10000D+01 X pnorm= 0.10000D+01 X delta= 0.54641D-07 X pdelta= 0.33861D+00 X lbndup= 0.16137D-06 X rbndup= 0.16137D-06 X lbndlw= 0.46365D+00 X rbndlw= 0.46365D+00 X ierr= 0 X eigenvalue bounds X delmax(capital delta for eigenv)= 0.338615D+00 X eigenvalue= 0.100000000000000D+01 0.000000000000000D+00 X aii= 0.70711D+00 0.00000D+00 bii= 0.70711D+00 0.00000D+00 k= 0.20000D+01 X eigenvalue= 0.200000000000000D+01 0.000000000000000D+00 X aii= 0.89443D+00 0.00000D+00 bii= 0.44721D+00 0.00000D+00 k= 0.12649D+01 X epsbnd= 0.10000D-09 X norm(aper,e)= 0.34641D+01 norm(bper,e)= 0.20000D+01 X start guptri for perturbed pair no. X iper= 1 X itst= 1 X guptri - workspace for rzstr - 4 5 4 5 1 X 6 26 31 61 36 41 86 46 51 56 X X Xguptri - m,n,epsu= 4 5 0.100000D-07 X Xreduction 1 X Xkstr, last= 3 X 1 2 3 X 1 1 1 X 1 1 0 Xaccumulated perturbations in a,b = 0.000000D+00 0.323534D-08 X Xreduction 4, kfirst= 5 X 1 2 3 4 5 X 1 1 1 -1 0 X 1 1 0 -1 0 Xaccumulated perturbations in a,b = 0.000000D+00 0.323534D-08 X X Xfinal kstr= X 1 2 3 4 5 6 7 X 1 1 1 -1 -1 2 -1 X 1 1 0 -1 -1 2 -1 X nsumrz= 3 X rsumrz= 2 X djordz= 0 X nsumli= 0 X rsumli= 0 X djordi= 0 X dimreg= 2 X X Xfinal pstruc= 3 3 3 3 Xfinal struc = X 1 1 1 X rtce= 3 zrce= 3 fnce= 5 ince= 5 X rtre= 2 zrre= 2 fnre= 4 inre= 4 X computed eigenvalues X eigenvalue= 0.10000D+01 0.00000D+00 X eigenvalue= 0.20000D+01 0.00000D+00 X eigenvalues before reordering X eigenvalue= 0.10000D+01 0.00000D+00 X eigenvalue= 0.20000D+01 0.00000D+00 X eigenvalues after reorder and X computed eigenvalues X eigenvalue= 0.10000D+01 0.00000D+00 X eigenvalue= 0.20000D+01 0.00000D+00 X results from guptri and reorder, iper= 1 X rtce= 3 zrce= 3 fnce= 5 ince= 5 X rtre= 2 zrre= 2 fnre= 4 inre= 4 Xpstruc = 3 3 3 3 Xstruc = X 1 1 1 Xnsumrz= 3 Xrsumrz= 2 Xdjordz= 0 Xnsumli= 0 Xrsumli= 0 Xdjordi= 0 Xdimreg= 2 Xndim= 0 X Relative perturbation in a= 0.000000D+00 X Relative perturbation in b= 0.161767D-08 X Frobeniusnorm of deleted singular values=0.323534D-08 X kstr, step= 7 X 1 1 1 -1 -1 2 -1 X 1 1 0 -1 -1 2 -1 X lda= 30 m= 4 n= 5 X Transformed matrix A X ---------------------------------------------------------------------- X 0.00000000000000000D+00 0.15811388300841898D+01 0.94868329800765749D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 0.18973665960867427D+01 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X X X -0.23334701776492633D-08 0.35986913438820403D-09 X -0.30001759273521061D-08-0.30392879514122899D-15 X 0.10000000000677638D+01 0.11389717383782053D-15 X 0.00000000000000000D+00 0.20000000000677636D+01 X X X lda= 30 m= 4 n= 5 X Transformed matrix B X ---------------------------------------------------------------------- X -0.74535599255605001D+00-0.63245553197415183D+00-0.21081851065805099D+00 X 0.00000000000000000D+00 0.31622776601683827D+00-0.94868329805051377D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X 0.00000000000000000D+00 0.00000000000000000D+00 0.00000000000000000D+00 X X X 0.66670575164805837D-09-0.14394773176978720D-09 X 0.30001760174657973D-08-0.64776443177036579D-09 X 0.10000000000000002D+01 0.16868641798082849D-15 X 0.00000000000000000D+00 0.10000000000000000D+01 X X X cond(PPper)=0.100000D+01 X cond(QQper)=0.100000D+01 X abs(a-acopy)=0.217928D-14 X relative dif.=0.629104D-15 X abs(b-bcopy)=0.384527D-08 X relative dif.=0.192263D-08 X fro(a-pp" * acopy * qq)=0.850587D-15 X relative fro for a-part=0.245543D-15 X fro(b-pp" * bcopy * qq)=0.323534D-08 X relative fro for b-part=0.161767D-08 X X X perturbation results for iper= 1 itst= 1 epsbnd = 0.10000D-09 X dist =0.324950D-08 X distup =0.547367D-07 X pcolrs = 3 X prowrs = 2 X pdelta =0.338615D+00 X case 1 of theorem holds X rbndlw =0.463648D+00 X lbndlw =0.463648D+00 X rbdupp =0.959645D-08 X lbdupp =0.959645D-08 X thetar=0.216924D-08 X thetal =0.112663D-08 X X Xnew eigenbound test for iper= 1 X compare eigenvalues 1 X unperturbed eigenvalue = 0.100000000000000D+01 0.000000000000000D+00 X perturbed eigenvalue = 0.100000000006776D+01 0.000000000000000D+00 X eigenbound holds with ebnd= 0.45955D-08 edif= 0.47916D-10 X compare eigenvalues 2 X unperturbed eigenvalue = 0.200000000000000D+01 0.000000000000000D+00 X perturbed eigenvalue = 0.200000000006776D+01 0.000000000000000D+00 X eigenbound holds with ebnd= 0.45955D-08 edif= 0.30305D-10 X X Xtest eigenbounds for iper= 1 X compare eigenvalues 1 X unperturbed eigenvalue = 0.100000000000000D+01 0.000000000000000D+00 X perturbed eigenvalue = 0.100000000006776D+01 0.000000000000000D+00 X eigenbound holds with ebnd= 0.64990D-08 edif= 0.33882D-10 X compare eigenvalues 2 X unperturbed eigenvalue = 0.200000000000000D+01 0.000000000000000D+00 X perturbed eigenvalue = 0.200000000006776D+01 0.000000000000000D+00 X eigenbound holds with ebnd= 0.41103D-08 edif= 0.13553D-10 X Summary of statistics: X ===================== X X Number of bad svds and qzs = ninfo = 0 X Number of inapplicable eigenbounds = badeig = 0 X X Distance between pencils on the surface X divided by the true distance between perturbed X and unperturbed input pencils X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X min = 10.722715194391 X average = 10.722715194391 X max = 10.722715194391 X X Distance between pencils on the surface X divided by the size of the perturbation (epsbnd) X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X min = 32.494973305147 X average = 32.494973305147 X max = 32.494973305147 X Reducing subspaces: X Different cases: X 0 0 0 0 0 0 0 0 0 0 0 0 X X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X Case 1: right upper bounds X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X min = 4.4238853094134 X average = 4.4238853094134 X max = 4.4238853094134 X Case 1: left upper bounds X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X min = 8.5178173535157 X average = 8.5178173535157 X max = 8.5178173535157 X Case 2: right lower bounds X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X min = 0. X average = 0. X max = 0. X Case 2: left lower bounds X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X min = 0. X average = 0. X max = 0. X Eigenvalues: number of them= 2 X Different cases (Gerschgorin type bounds): X Eigv. no. 1 X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X Eigv. no. 2 X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X Eigenvalue bounds (upper) X Eigv. no. 1 X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X min = 191.81443960941 X average = 191.81443960941 X max = 191.81443960941 X Eigv. no. 2 X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X min = 303.28426484587 X average = 303.28426484587 X max = 303.28426484587 X Different cases( new bounds from LAA87): X Eigv. no. 1 X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X Eigv. no. 2 X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X Eigenvalue bounds (upper) X Eigv. no. 1 X 0 0 0 0 0 0 0 0 0 0 0 0 X X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X min = 95.907059974105 X average = 95.907059974105 X max = 95.907059974105 X Eigv. no. 2 X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 1 0 0 0 0 0 0 0 0 0 1 100 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 X X min = 151.64229361397 X average = 151.64229361397 X max = 151.64229361397 XMaximum values of radife, rbdife 2.4554337601543D-16 1.6176676252446D-09 END_OF_zgschur.c1 if test 17039 -ne `wc -c zgschurm.f <<'END_OF_zgschurm.f' Xc On this file March 1990: zgschurm, edist Xc X program zgschurm Xc implicit none Xc**** debug space Xc the common-block declarations assume that the dimension of the Xc input matrix pencil a - lambda b is not larger than 30. Xc the debug space is used for producing debug outputs (optional, Xc see below) Xc X integer abdim, wdim, abdim6 Xc abdim6 = abdim + 6 X parameter ( abdim = 30, wdim = 20000, abdim6 = 36) X common /debug1/ acopy(abdim, abdim),bcopy(abdim, abdim), X * atest(abdim, abdim), btest(abdim, abdim), swap X common /debug2/ idbg(20), outunit X complex*16 acopy, bcopy, atest, btest X integer idbg, outunit X logical swap Xc***** This version of zgschurm computes pairs of reducing Xc subspaces associated with different subspaces of a Xc (generalized) state space system. Further, it collects statistics Xc for random examples Xc Xc Revision: 900323 (this version goes with final versions of Xc guptri and bounds) Xc Xc*+*+*+ Xc The program starts by asking for input and output Xc files (infile and outfile) where Xc infile contains A and B of dimension M by N and Xc debug and control inputs (see below) Xc outfile contains output from the program Xc Then it asks for a textstring identifying the run (.le. 80 chars) Xc*+*+*+ Xc Xc***** debug flags (20i1) Xc idbg(1) ne 0 - turn on debug output for kcfmain Xc idbg(2) ne 0 - turn on debug output for guptri Xc idbg(3) ne 0 - turn on debug output for krnstr Xc idbg(4) ne 0 - turn on debug output for rzstr Xc idbg(5) ne 0 - turn on debug output for listr Xc idbg(6) ne 0 - turn on debug output for rcsvdc Xc idbg(7) ne 0 - turn on debug output for reordr Xc idbg(8) ne 0 - turn on debug output for exchng Xc idbg(9) ne 0 - turn on debug output for pbound (no debug) Xc idbg(10) ne 0 - turn on debug output for gvec Xc+*+ idbg(11) ne 0 - turn on debug output for pertb1 860729 Xc idbg(12) ne 0 - turn on debug output for qz Xc idbg(19) ne 0 - turn on debug output for bndwsp Xc idbg(20) ne 0 - turn on debug output for bound Xc Xc***** control inputs (2i1,i4,i1) Xc izero ne 0 - zero out nonzero singular values during reduction Xc itrpose ne 0 - transpose input matrices a and b Xc job (4th digit) ne 0 - pre, postmultiply a by random nonsingular Xc matrices p, q, called wanta in output Xc job (3rd digit) ne 0 - pre, postmultiple a by random nonsingular Xc matrices p, q, called wantb in output Xc job (2nd digit) ne 0 - add random noise of size machep to a, b, Xc called pertur in output Xc job (1st digit) ne 0 - print block structured input a, b, and Xc final input a,b if different, Xc called prints in output Xc exprin ne 0 - print outs for each example and statistics Xc eq 0 - only print outs of statistics Xc*+*+*+ 860731 Xc epsu (2d10.0) user specified uncertainty in the input Xc A and B (used for deleting small singular Xc values) Xc gap gap between small singular values Xc epsper (d10.0) size of perturbation to A and B on input Xc (only used if job (2nd digit ne 0) Xc numex (3i5) number of values of epsbnd's Xc numtst number of times we shall add noise of Xc size epsbnd(iper) to A and B Xc jobper structure of the perturbations added Xc to A and B Xc epsbnd(numex) (5d10.0) size of perturbation that we add to A and B Xc*+*+*+ Xc Statistics are collected from numex*numtst random examples. Xc Starting from a nongeneric pencil and a rule(epsu,gap) for Xc choosing a particular set of reducing subspaces we add random noise Xc to get perturbed pencils as input for GUPTRI Xc Xc X complex*16 a(abdim,abdim),b(abdim,abdim), X * work(wdim), pp(abdim,abdim), qq(abdim,abdim) X complex*16 zat, zbt, aorig(abdim,abdim), borig(abdim,abdim), X * aprim(abdim,abdim), bprim(abdim,abdim), X * ppper(abdim,abdim), X * qqper(abdim,abdim) X complex*16 aortr(abdim,abdim),bortr(abdim,abdim) X integer rtre, rtce, zrre, zrce, fnre, fnce, inre, ince X integer pstruc(4), struc(abdim), space X integer nsumrz,rsumrz,nsumli,rsumli X integer djordz,djordi,dimreg Xc 06/16/87 X integer rowb, colb, rowe, cole Xc X integer kstr(4,abdim6), step, allreg, krstrt, kcstrt, icase X integer three, ithree, ecase X integer ndim, rindx(abdim6), ftest, colrs, rowrs, pcolrs, prowrs X integer sstrt,estrt,wstrt,ninfo X integer fout, fin, folhp, fcrhp X external fout, fin, folhp, fcrhp, ftest, ftestp X logical zero, ldebug X logical trpose, pbndok, nostat Xc*+*+ demmel, 7/3/86 X logical ebndok Xc*+*+ X complex*16 evala(abdim),evalb(abdim) Xc*+*+ X complex*16 evalap(abdim),evalbp(abdim) Xc*+*+ X real*8 lbndup, lbndlw, gvcond(abdim), anormf, bnormf, X * lbdupp, rbdupp, relerr, rbndup, rbndlw X real*8 scl, epsu, gap X real*8 epsper, epsbnd(20), delmax Xc*+*+*+ 860731 X integer numex, jobper, numtst, iper, itst, exprin X character*80 infile, outfile, ident Xc*** data for statistics /rev 870526 and 870626 X integer statrs(6,12), stateg(3,12,10), stateg1(2,12,10) X integer sdstqt(7,12), sdstqe(7,12) X integer srqtup(6,12), slqtup(6,12) X integer srqtlw(6,12), slqtlw(6,12), segqt(6,12,10) X integer segqt1(6,12,10), badeig X real*8 rqtup, lqtup, rqtlw, lqtlw, egqt Xc Xc variables for min, average and max computations /870526 X real*8 maxrda, maxrdb X real*8 minqt, avrqt, maxqt X real*8 minqe, avrqe, maxqe X real*8 minrup, avrrup, maxrup X real*8 minlup, avrlup, maxlup X real*8 minrlw, avrrlw, maxrlw X real*8 minllw, avrllw, maxllw X real*8 minegq(10), avregq(10), maxegq(10) X real*8 minegq1(10), avregq1(10), maxegq1(10) Xc end of new variables for statistics/ 870526 and 870626 X logical infnt, infntp Xc*+*+*+ X complex*16 dum, dummy X integer i, izero, itrpos, job, lda, ldb, m, n, ldab X integer ldqq, ierr, info, len, k, jjj, ieig, j, ldpp X real*8 cpp, cqq, difa, difb, adife, bdife, anore, bnore X real*8 dsvd, adsvd, bdsvd X real*8 rdifa, rdifb, radife, rbdife, difu, difl X real*8 pqnorm, qnorm, pnorm X real*8 pdelta, dist, dstqt, dstqe, dstpu X real*8 dsvdp, adsvdp, bdsvdp X real*8 distup X real*8 thetal, thetar, ebnd, edif X real*8 cond, cnorm, cdife Xc**** generate a singular matrix pencil Xc Xc data lda/20/, ldb/20/, ldpp/20/, ldqq/20/, ldab/20/ X lda = abdim X ldb = abdim X ldpp = abdim X ldqq = abdim X ldab = abdim Xc Xc*+*+*+ 860731 X write(*,*) 'Give infile and outfile:' X read(*,7034) infile X read(*,7034) outfile X 7034 format(A) Xc X write(*,*) 'Identify this testrun:' X read (*,7034) ident X open(5, file = infile, status = 'old') X outunit = 6 X open(6, file = outfile, status = 'new') X write(6,7035) 'Testrun identification: ',ident X 7035 format(A,A//) Xc*+*+*+ Xc Xc read in matrix dimensions and matrices a and b X read(5,6543) m,n X do 7010 i = 1, m X read(5, *) (a(i,j), j = 1, n) X7010 continue X do 7015 i = 1, m X read(5,*) (b(i,j), j = 1,n) X7015 continue X6543 format(2i5) Xc Xc copy a and b to acopy and bcopy, respectively Xc X call cmcopy(b,ldb,m,n,bcopy) X call cmcopy(a,lda,m,n,acopy) Xc X call cmatpr(a,lda,m,n,'final version of a input') X call cmatpr(b,lda,m,n,'final version of b input') Xc read in debug controls X read(5,1235) (idbg(i),i=1,20) X1235 format(20i1) X write (6,1236) (j,j=1,20), (idbg(j),j=1,20) X1236 format(' debug controls -',/,1x,20i3,/,1x,20i3) Xc*+*+*+ Xc read in job controls X read(5,1234) izero,itrpos,job, exprin X1234 format(2i1,i4,i1) Xc Xc read epsu (relative error in input matrices) and gap Xc (for nullity testing) X read(5, 202) epsu, gap X 202 format(2d10.0) X write(6, 203) 'input: epsu=', epsu, 'gap=', gap X 203 format(t5,a,d15.5) X read (5,204) epsper X 204 format(d10.0) Xc*+*+*+ X read(5,205) numex, numtst, jobper X 205 format(3i5) Xc*+*+*+ X if (numex .gt. 0) read(5,206) ( epsbnd(i), i= 1, numex) X 206 format (5d10.0) Xc X write(6,207) ' epsper =', epsper, ' numex =', numex, X * ' numtst=', numtst, ' jobper =', jobper X 207 format(t5, a, d15.5, 3 (a, i5)) X if (numex .gt. 0) then X write(6,207) ' epsbnd =' X write(6,208) ( epsbnd(i), i= 1, numex) X 208 format(t5, 5d12.3) X endif Xc*+*+*+ start (trpose never used in this code!) X trpose=.false. X if (itrpos.ne.0) trpose=.true. X zero=.false. X if (izero.ne.0) zero=.true. X write(6,201) 'zero=', zero X201 format(t5,a,l1) Xc X nostat = .true. X if (exprin .eq. 0) nostat = .false. X write(6,201) 'nostat=', nostat Xc*+*+*+ end Xc copy a and b to aorig and borig for later perturbing Xc aorig and borig should never be changed!!!!! X call cmcopy(a, ldab, m, n, aorig) X call cmcopy(b, ldab, m, n, borig) X anormf = cnorm(a, ldab, m, n, 0, work) X bnormf = cnorm(b, ldab, m, n, 0, work) Xc X write(6, 350) 'epsper=', epsper X write(6, 350) 'norm(a,e)=', anormf, 'norm(b,e)=', bnormf X 350 format(t5,a,d12.5,tr5,a,d12.5,tr5,d12.5) Xc X 200 format(t5,a,d12.6) X X write(6, 100) 'start guptri' Xc Xc**** 6/16/87 Xc X call guptri(a ,b , ldab, m, n, epsu, gap, zero, X * pp, ldpp, qq, ldqq, adsvd, bdsvd, X * rtre, rtce, zrre, zrce, fnre, fnce, inre, ince, X * pstruc, struc, work, kstr, info) Xc Xc*** 6/18/87 X if (info .ne. 0) then X write (6,2000) 'after first guptri, info=', info X endif X dsvd = sqrt ( (anormf*adsvd)**2 + (bnormf*bdsvd)**2 ) Xc Xc**** 6/16/87 Xc compute step by searching through kstr X three = 0 X do 61687 ithree = 1, 20 X if ( three .eq. 3) go to 61688 X if( kstr(1, ithree) .eq. -1) then X three = three + 1 X endif X61687 continue Xc X if ( three .lt. 3) then X write(*,*) 'ERROR in kstr (computing step in driver)' X stop X endif Xc X61688 continue X step = ithree - 1 Xc*** end of computing step Xc Xc**** 6/15/87 Xc compute ome structure infortmation (not parameters to guptri any more) X nsumrz = zrce X rsumrz = zrre X nsumli = n - fnce X rsumli = m - fnre X djordz = zrre - rtre X djordi = inre - fnre X dimreg = fnre - zrre X ndim = 0 Xc Xc*+*+ added 06/16/87 Xc**** reorder the eigenvalues according to the user specified Xc integer function ftest Xc set debug flag for guptri so we can compare with old version of Xc driver X ldebug = idbg(2) X allreg = dimreg + djordz + djordi X rowb = rsumrz - djordz + 1 X colb = nsumrz - djordz + 1 X rowe = rowb + allreg - 1 X cole = colb + allreg - 1 X if (ldebug) then X write(outunit, 2005) 'eigenvalues before reordering' X do 70 i = rowb, rowe X j = colb + i - rowb X if (abs(b(i ,j)) .eq. 0. ) then X write(outunit, 2005) 'infinite eigenvalue',a(i,j), b(i,j) X 2005 format(t5,a,4d15.5) X else X write(outunit, 2005) 'eigenvalue=', a(i,j)/b(i,j) X endif X 70 continue X endif X if (allreg .ge. 1) then X call reordr(a, b, ldab, m, n, rowb, colb, rowe, cole, X * ftest, ndim, rindx, pp, ldpp, qq, ldqq) Xc X if (idbg(2) .gt. 1) then X call cmatpr(qq,ldqq,n,n,'qq after reordr') X call cmatpr(pp,ldpp,m,m,'pp after reordr') X endif X endif Xc X if (ldebug) then X write(outunit, 2005) 'eigenvalues after reorder and' X write(outunit, 2005) 'computed eigenvalues' X do 75 i = rowb, rowe X j = colb + i - rowb X if (abs(b(i ,j)) .eq. 0. ) then X write(outunit, 2005) 'infinite eigenvalue',a(i,j), b(i,j) X else X write(outunit, 2005) 'eigenvalue=', a(i,j)/b(i,j) X endif X 75 continue X endif Xc Xc+*+ end add reorder* Xc save transformed original a and b in aprim, bprim,aortr,bortr Xc X call cmcopy(a, ldab, m, n, aprim) X call cmcopy(b, ldab, m, n, bprim) Xc X call cmcopy(a, ldab, m, n, aortr) X call cmcopy(b, ldab, m, n, bortr) Xc Xc compute aprim = pp * aprim * qq**H and Xc bprim = pp * bprim * qq**H X call cmatml(aprim,ldab,m,n,pp,ldpp,m,aprim,ldab,work,1) X call cmatmr(aprim,ldab,m,n,qq,ldqq,n,aprim,ldab,work,3) X call cmatml(bprim,ldab,m,n,pp,ldpp,m,bprim,ldab,work,1) X call cmatmr(bprim,ldab,m,n,qq,ldqq,n,bprim,ldab,work,3) X if (idbg(1) .ge. 2) then X call cmatpr(aprim,ldab,m,n,'final aprim') X call cmatpr(bprim,ldab,m,n,'final bprim') X endif X write(6, 100) 'results from guptri and reorder' X 100 format (t5, a, i4) Xc Xc**** 6/15/87 X write(6,7357) 'rtce=',rtce,'zrce=',zrce,'fnce=',fnce, X + 'ince=',ince,'rtre=',rtre,'zrre=',zrre, X + 'fnre=',fnre,'inre=',inre X write (6,7355) (pstruc(j),j=1,4) X if (pstruc(4).gt.0) write (6,7356)(struc(j),j=1,pstruc(4)) X 7355 format('pstruc = ',4i4,/,'struc =') X 7356 format(15i4) X 7357 format(4(3x,a,i4),/,4(3x,a,i4)) X write (6,123) nsumrz,rsumrz,djordz,nsumli,rsumli,djordi, X * dimreg, ndim X 123 format('nsumrz=',i5,/,'rsumrz=',i5,/,'djordz=',i5,/, X * 'nsumli=',i5,/,'rsumli=',i5,/,'djordi=',i5,/, X * 'dimreg=',i5,/,'ndim= ',i5) X write(6,200) 'Relative perturbation in a= ', adsvd X write(6,200) 'Relative perturbation in b= ', bdsvd X write(6,200) 'Frobeniusnorm of deleted singular vaules=', X * dsvd X write(6, 100) 'kstr, step=',step X do 10 i = 1, 2 X write(6, 300) (kstr(i,j), j = 1, step) X 10 continue X 300 format(t5, 20i3) X if(idbg(1).ge.1)call cmatpr(a,lda,m,n,'Transformed matrix A') X if(idbg(1).ge.1)call cmatpr(b,ldb,m,n,'Transformed matrix B') X if (idbg(1).ge.2) call cmatpr(pp, ldpp, m, m, 'PP') X if (idbg(1).ge.2) call cmatpr(qq, ldqq, n, n, 'QQ') X cpp=cond(pp,ldpp,m,m,work) X write(6, 105) 'cond(PP)=', cpp X 105 format(t5, a, d12.6) X cqq=cond(qq,ldqq,n,n,work) X write(6, 105) 'cond(QQ)=', cqq Xc X call cmcopy(acopy, ldab, m, n, atest) X call cmcopy(bcopy, ldab, m, n, btest) X call cmatml(atest,lda,m,n,pp,ldpp,m,atest,lda,work,3) X call cmatmr(atest,lda,m,n,qq,ldqq,n,atest,lda,work,1) X if(idbg(1).ge.2) call cmatpr(atest,lda,m,n,'pp'' * a * qq') X call cmatml(btest,ldb,m,n,pp,ldpp,m,btest,ldb,work,3) X call cmatmr(btest,ldb,m,n,qq,ldqq,n,btest,ldb,work,1) X if(idbg(1).ge.2) call cmatpr(btest,ldb,m,n,'pp'' * b * qq') X difa=0 X difb=0 X do 20 i=1,m X do 30 j=1,n X difa = difa+abs(a(i,j)-atest(i,j)) X difb = difb+abs(b(i,j)-btest(i,j)) X 30 continue X 20 continue Xc X adife = cdife(a,atest,ldab,m,n) X bdife = cdife(b,btest,ldab,m,n) X anore = cnorm(a,ldab,m,n,0,work) X bnore = cnorm(b,ldab,m,n,0,work) Xc X rdifa = difa X if (anore .gt. 0.) rdifa = difa / anore X rdifb = difb X if (bnore .gt. 0.) rdifb = difb / bnore X radife = adife X if (anore .gt. 0.) radife = adife / anore X rbdife = bdife X if (bnore .gt. 0.) rbdife = bdife / bnore X maxrda = radife X maxrdb = rbdife Xc X write(6,105) 'abs(a-acopy)=',difa,'relative dif.=',rdifa X write(6,105) 'abs(b-bcopy)=',difb,'relative dif.=',rdifb X write(6,105) 'fro(a-pp" * acopy * qq)=', adife X write (6,105) 'relative fro for a-part=', radife X write (6,105) 'fro(b-pp" * bcopy * qq)=', bdife X write (6,105) 'relative fro for b-part=', rbdife Xc Xc**** compute error bounds for reducing subspaces Xc containing right singular part and eigenvalues Xc specified by ftest Xc Xc skip if right or left reducing subspace is zero or full dimensional Xc colrs = dimension of right reducing subspace Xc rowrs = dimension of left reducing subspace Xc allreg = dimension of the whole regular part 06/16/87 X colrs = nsumrz - djordz + ndim X rowrs = rsumrz - djordz + ndim X len = allreg - ndim Xc X write(6, 2000) 'colrs=', colrs, 'rowrs=', rowrs, 'len=', len X 2000 format(t5,a,i5) Xc Xc**** 6/22/87, compute workspace, stop if insufficient X call bndwsp(m,n,rowrs+1,colrs+1,len,ecase,space,info) X write(6,2000) ' bndwsp' X write(6,2000) 'ecase=',ecase,'space=',space,'info=',info X if (info.eq.1 .or. space.gt.wdim) stop Xc Xc**** 6/21/87 stop if no tests desired X if (numtst .eq. 0) stop Xc Xc deleted singular values cannot be less than epsu*(norma+normb) X dsvd = max(dsvd, epsu*( anormf + bnormf )) Xc Xc*** 06/16/87 new version of bounds Xc compute difl, difu, qnorm,pnorm, etc and eigenvalue bounds X call bound(a, b, ldab, m ,n, rowrs+1, colrs+1, len, X * evala, evalb, delmax, gvcond, pqnorm, ecase, X * pdelta, difl, difu, qnorm, pnorm, icase, X * work, ierr) Xc X write(6, 2000) ' icase= ', icase, ' ecase= ', ecase, X + ' ierr= ', ierr X write(6,203) 'delmax=',delmax,'pdelta=',pdelta,'difl=',difl, X + 'difu=',difu,'qnorm=',qnorm,'pnorm=',pnorm, X + 'pqnorm=',pqnorm,'dsvd=',dsvd Xc Xc*** 6/18/87 bounds for trivial spaces handled by bound, icase Xc pbndok = colrs .gt. 0 .and. rowrs .lt. m X pbndok = .true. Xc X if (pbndok) then Xc**** evaluate space - bounds X call evalbd( dsvd, pdelta, qnorm, pnorm, icase, m, n, X * rowrs+1, colrs+1, lbndup, rbndup, lbndlw, rbndlw) Xc X write(6,106) difl,difu,qnorm,pnorm,dsvd,pdelta,lbndup,rbndup, X + lbndlw,rbndlw,ierr X106 format(/,' results from pbound',/,' difl= ',d20.5, X + /,' difu= ',d20.5,/,' qnorm= ',d20.5,/,' pnorm= ',d20.5, X + /,' delta= ',d20.5, X + /,' pdelta=',d20.5,/,' lbndup=',d20.5,/,' rbndup=',d20.5, X + /,' lbndlw=',d20.5,/,' rbndlw=',d20.5,/,' ierr= ',i3) X endif Xc Xc**** compute error bounds for remaining eigenvalues Xc only if there are any (allreg.gt.ndim) and Xc no left (Kronecker) indices Xc ( rsumli .eq. nsumli ) Xc note: the case with no right (Kronecker) indices Xc and a regular part can be handled by transposing the Xc output from guptri (!!??) Xc note: includes perturbation theory for regular pencils Xc ( rsumli .eq. nsumli .and. rsumrz .eq. nsumrz) Xc Xc allreg = dimreg + djordz + djordi Xc len = allreg - ndim Xc*+*+ Xc ebndok = allreg .gt. ndim .and. rsumli .eq.nsumli X ebndok = len .gt. 0 Xc**** changed by demmel, 6/30/86 X if ( ebndok ) then Xc X krstrt = rsumrz - djordz + ndim + 1 X kcstrt = nsumrz - djordz + ndim + 1 X info = ierr X if (info .eq. 0) then Xc no multiple eigenvalues X write(6, 184) 'eigenvalue bounds' X 184 format(t5,a,2d23.15) X write(6, 105) 'delmax(capital delta for eigenv)= ', X * delmax X do 183 i = 1, len Xc X zat= evala(i) X zbt = evalb(i) X if (abs(zbt) .eq. 0.) then X write(6, 184) 'infinite eigenvalue' X else X write(6,184) 'eigenvalue= ', zat / zbt X endif X write(6,108) zat, zbt, gvcond(i) X 108 format(' aii=',2d13.5,' bii=',2d13.5,' k=',d13.5) X 183 continue X else Xc there are multiple eigenvalues X write(6,184) 'multiple eigenvalues' Xc 061387 changed X ebndok = .false. X do 185 i=1,len X zat=evala(i) X zbt=evalb(i) X if (abs(zbt).eq.0.) then X write(6,184) 'infinite eigenvalue' X else X write(6,184) 'eigenvalue=',zat/zbt X endif X185 continue X endif Xc X endif Xc for doing perturbation theory for eigenvalues Xc Xc***** compute GUPTRI forms for perturbed pencils Xc Xc prepare for statistics X do 8020 i = 1, 12 X do 8010 j = 1, 7 X sdstqt(j,i) = 0 X sdstqe(j,i) = 0 X 8010 continue X do 8012 j = 1, 3 X do 8011 k = 1, 10 X stateg(j,i,k) = 0 X if ( j .le. 2) stateg1(j,i,k) = 0 X 8011 continue X 8012 continue X do 8015 j = 1, 6 X statrs(j,i) = 0 X srqtup(j,i) = 0 X slqtup(j,i) = 0 X srqtlw(j,i) = 0 X slqtlw(j,i) = 0 X do 8013 k = 1, 10 X segqt(j,i,k) = 0 X segqt1(j,i,k) = 0 X 8013 continue X 8015 continue X 8020 continue Xc write(6,*) 'statrs before 7000' Xc write(6,9500) ((statrs(i,j), j=1,11), i=1,6) Xc X badeig = 0 X ninfo = 0 X if ( numex. gt. 0 .and. numtst .gt. 0) then X do 7000 iper = 1, numex Xc X do 6900 itst = 1, numtst Xc perturb a and b ( copies in acopy, and bcopy) Xc*+*+*+ start change 860729 X call pertb1( aorig, borig, a, b, ldab, m, n, epsbnd(iper), X * work,jobper,nostat) Xc*+*+*+ end X anormf = cnorm(a, ldab, m, n, 0, work) X bnormf = cnorm(b, ldab, m, n, 0, work) Xc**** compute the Kronecker structure Xc X X if (nostat) then X write(6, 100) 'start guptri for perturbed pair no.' X write(6, 100) 'iper= ', iper, 'itst= ', itst X endif Xc Xc**** 6/16/87 X call guptri(a ,b , ldab, m, n, epsu, gap, zero, X * ppper, ldpp, qqper, ldqq, X * adsvdp, bdsvdp, X * rtre, rtce, zrre, zrce, fnre, fnce, inre, ince, X * pstruc, struc, work, kstr, info) Xc Xc**** 6/18/87 Xc if (info.ne.0) write(6,2000) 'after guptri, info=',info X if (info.ne.0) then X ninfo = ninfo +1 X write(6,2000) 'after guptri, info=',info Xc goto next perturbed pair X goto 6900 X endif X dsvdp = sqrt ( (anormf*adsvdp)**2 + (bnormf*bdsvdp)**2 ) Xc Xc Xc**** 6/16/87 Xc compute step by searching through kstr X three = 0 X do 61689 ithree = 1, 20 X if ( three .eq. 3) go to 61690 X if( kstr(1, ithree) .eq. -1) then X three = three + 1 X endif X61689 continue Xc X if ( three .lt. 3) then X write(*,*) 'ERROR in kstr (computing step in driver)' X stop X endif Xc X61690 continue X step = ithree - 1 Xc*** end of computing step Xc Xc Xc**** 6/15/87 Xc compute these (not parameters to guptri any more) X nsumrz = zrce X rsumrz = zrre X nsumli = n - fnce X rsumli = m - fnre X djordz = zrre - rtre X djordi = inre - fnre X dimreg = fnre - zrre X ndim = 0 Xc Xc*+*+ added 06/16/87 Xc**** reorder the eigenvalues according to the user specified Xc integer function ftest Xc X allreg = dimreg + djordz + djordi X rowb = rsumrz - djordz + 1 X colb = nsumrz - djordz + 1 X rowe = rowb + allreg - 1 X cole = colb + allreg - 1 X if (ldebug) then X write(outunit, 2005) 'eigenvalues before reordering' X do 770 i = rowb, rowe X j = colb + i - rowb X if (abs(b(i ,j)) .eq. 0. ) then X write(outunit, 2005) 'infinite eigenvalue',a(i,j), b(i,j) X else X write(outunit, 2005) 'eigenvalue=', a(i,j)/b(i,j) X endif X 770 continue X endif X if (allreg .ge. 1) then X call reordr(a, b, ldab, m, n, rowb, colb, rowe, cole, X * ftest, ndim, rindx, ppper, ldpp, qqper, ldqq) Xc X if (idbg(2) .gt. 1) then X call cmatpr(qqper,ldqq,n,n,'qqper after reordr') X call cmatpr(ppper,ldpp,m,m,'ppper after reordr') X endif X endif Xc X if (ldebug) then X write(outunit, 2005) 'eigenvalues after reorder and' X write(outunit, 2005) 'computed eigenvalues' X do 775 i = rowb, rowe X j = colb + i - rowb X if (abs(b(i ,j)) .eq. 0. ) then X write(outunit, 2005) 'infinite eigenvalue',a(i,j), b(i,j) X else X write(outunit, 2005) 'eigenvalue=', a(i,j)/b(i,j) X endif X 775 continue X endif Xc Xc+*+ end add reorder* Xc X if (nostat) then X write(6, 100) 'results from guptri and reorder, iper= ', X * iper Xc Xc**** 6/15/87 X write(6,7357) 'rtce=',rtce,'zrce=',zrce,'fnce=',fnce, X + 'ince=',ince,'rtre=',rtre,'zrre=',zrre, X + 'fnre=',fnre,'inre=',inre X write (6,7355) (pstruc(j),j=1,4) X if (pstruc(4).gt.0) write (6,7356)(struc(j),j=1,pstruc(4)) Xc X write (6,123) nsumrz,rsumrz,djordz,nsumli,rsumli,djordi, X * dimreg, ndim X write(6,200) 'Relative perturbation in a= ', adsvdp X write(6,200) 'Relative perturbation in b= ', bdsvdp X write(6,200) 'Frobeniusnorm of deleted singular values=', X * dsvdp X write(6, 100) 'kstr, step=',step X do 710 i = 1, 2 X write(6, 300) (kstr(i,j), j = 1, step) X 710 continue X if (idbg(1).ge.1) X * call cmatpr(a,lda,m,n,'Transformed matrix A') X if (idbg(1).ge.1) X * call cmatpr(b,ldb,m,n,'Transformed matrix B') X if(idbg(1).ge.2) call cmatpr(ppper, ldpp, m, m, 'PPper') X if(idbg(1).ge.2) call cmatpr(qqper, ldqq, n, n, 'QQper') X cpp=cond(ppper,ldpp,m,m,work) X write(6, 105) 'cond(PPper)=', cpp X cqq=cond(qqper,ldqq,n,n,work) X write(6, 105) 'cond(QQper)=', cqq X endif Xc X call cmcopy(acopy, ldab, m, n, atest) X call cmcopy(bcopy, ldab, m, n, btest) X call cmatml(atest,lda,m,n,ppper,ldpp,m,atest,lda,work,3) X call cmatmr(atest,lda,m,n,qqper,ldqq,n,atest,lda,work,1) X if (idbg(1).ge.2) X * call cmatpr(atest,lda,m,n,'ppper'' * aper * qqper') X call cmatml(btest,ldb,m,n,ppper,ldpp,m,btest,ldb,work,3) X call cmatmr(btest,ldb,m,n,qqper,ldqq,n,btest,ldb,work,1) X if (idbg(1).ge.2) X * call cmatpr(btest,ldb,m,n,'pperp'' * bper * qqper') X difa=0 X difb=0 X do 720 i=1,m X do 730 j=1,n X difa=difa+abs(a(i,j)-atest(i,j)) X difb=difb+abs(b(i,j)-btest(i,j)) X 730 continue X 720 continue Xc X adife = cdife(a,atest,ldab,m,n) X bdife = cdife(b,btest,ldab,m,n) X anore = cnorm(a,ldab,m,n,0,work) X bnore = cnorm(b,ldab,m,n,0,work) Xc X rdifa = difa X radife = adife X if (anore .gt. 0.) then X rdifa = difa / anore X radife = adife / anore X endif X rbdife = bdife X rdifb = difb X if (bnore .gt. 0.) then X rbdife = bdife / bnore X rdifb = difb / bnore X endif Xc Xc collect maximum values of radife, rbdife X maxrda = max(maxrda, radife) X maxrdb = max(maxrdb, rbdife) Xc X if (nostat) then X write(6,105)'abs(a-acopy)=',difa,'relative dif.=',rdifa X write(6,105)'abs(b-bcopy)=',difb,'relative dif.=',rdifb X write(6,105)'fro(a-pp" * acopy * qq)=', adife X write (6,105)'relative fro for a-part=', radife X write (6,105) 'fro(b-pp" * bcopy * qq)=', bdife X write (6,105) 'relative fro for b-part=', rbdife X endif Xc Xc Xc compute the dimensions of perturbed reducing subspaces X pcolrs = nsumrz - djordz + ndim X prowrs = rsumrz - djordz + ndim Xc Xc*+*+ Xc save eigenvalues for later use X do 223 jjj=pcolrs+1,n X evalap(jjj-pcolrs)=a(jjj-pcolrs+prowrs,jjj) X evalbp(jjj-pcolrs)=b(jjj-pcolrs+prowrs,jjj) X223 continue Xc compute the distance between the matrix pairs on Xc the (nongeneric) surface Xc Xc compute a = ppper * a * qqper**H and Xc b = ppper * b * qqper**H X call cmatml(a,ldab,m,n,ppper,ldpp,m,a,ldab,work,1) X call cmatmr(a,ldab,m,n,qqper,ldqq,n,a,ldab,work,3) X call cmatml(b,ldab,m,n,ppper,ldpp,m,b,ldab,work,1) X call cmatmr(b,ldab,m,n,qqper,ldqq,n,b,ldab,work,3) X if (idbg(1) .ge. 2) then X call cmatpr(a,ldab,m,n,'final aprimprim') X call cmatpr(b,ldab,m,n,'final bprimprim') X endif Xc Xc compute dist = distance between pencils on manifold X dist = sqrt( cdife(aprim, a, ldab, m, n) ** 2 + X * cdife(bprim, b, ldab, m, n) ** 2 ) X dstqt = dist/epsbnd(iper) Xc870526 seps1 = 1.0 / sqrt(epsbnd(iper)) Xc seps2 = 1.0 / (epsbnd(iper) ** 0.75) X if (dstqt .le. 1.0) then X sdstqt(1,iper) = sdstqt(1,iper) + 1 X elseif (dstqt .le. 10.0) then X sdstqt(2,iper) = sdstqt(2,iper) + 1 X elseif (dstqt .le. 100.0) then X sdstqt(3,iper) = sdstqt(3,iper) + 1 X elseif (dstqt .le. 1000.0) then X sdstqt(4,iper) = sdstqt(4,iper) + 1 X elseif (dstqt .le. 10000.0) then X sdstqt(5,iper) = sdstqt(5,iper) + 1 X elseif (dstqt .le. 100000.0) then X sdstqt(6,iper) = sdstqt(6,iper) + 1 X else X sdstqt(7,iper) = sdstqt(7,iper) + 1 X endif Xc X if (iper .eq. 1 .and. itst .eq. 1 ) then X minqt = dstqt X avrqt = dstqt X maxqt = dstqt X else X minqt = min(minqt,dstqt) X avrqt = avrqt + dstqt X maxqt = max(maxqt,dstqt) X endif Xc Xc compute the true distance between perturbed and unperturbed Xc input pencils X dstpu = sqrt( cdife(acopy, aorig, ldab, m, n)**2 X * + cdife(bcopy, borig, ldab, m, n)**2 ) X if (dstpu.eq.0.) dstpu = 1. X dstqe = dist / dstpu X if (dstqe .le. 1.0) then X sdstqe(1,iper) = sdstqe(1,iper) + 1 X elseif (dstqe .le. 10.0) then X sdstqe(2,iper) = sdstqe(2,iper) + 1 X elseif (dstqe .le. 100.0) then X sdstqe(3,iper) = sdstqe(3,iper) + 1 X elseif (dstqe .le. 1000.0) then X sdstqe(4,iper) = sdstqe(4,iper) + 1 X elseif (dstqe .le. 10000.0) then X sdstqe(5,iper) = sdstqe(5,iper) + 1 X elseif (dstqe .le. 100000.0) then X sdstqe(6,iper) = sdstqe(6,iper) + 1 X else X sdstqe(7,iper) = sdstqe(7,iper) + 1 X endif Xc X if (iper .eq. 1 .and. itst .eq. 1 ) then X