1
subroutine bdiag(lda,n,a,epsshr,rmax,er,ei,bs,x,xi,scale,
5
c dbdiag reduces a matrix a to block diagonal form by first
6
c reducing it to quasi-triangular form by hqror2 and then by
7
c solving the matrix equation -a11*p+p*a22=a12 to introduce zeros
9
c right transformation is factored : p*d*u*y ;where:
10
c p is a permutation matrix and d positive diagonal matrix
11
c u is orthogonal and y block upper triangular with identity
12
c blocks on the diagonal
16
c subroutine bdiag(lda,n,a,epsshr,rmax,er,ei,bs,x,xi,scale,
19
c integer lda, n, bs, job
20
c double precision a,er,ei,x,xi,rmax,epsshr,scale
21
c dimension a(lda,n),x(lda,n),xi(lda,n),er(n),ei(n),bs(n)
25
c starred parameters are altered by the subroutine
28
c *a an array that initially contains the m x n matrix
29
c to be reduced. on return, see job
31
c lda the leading dimension of array a. and array x,xi.
33
c n the order of the matrices a,x,xi
35
c epsshr the minimal conditionnement allowed for linear sytems
37
c rmax the maximum size allowed for any element of the
40
c *er a singly subscripted real array containing the real
41
c parts of the eigenvalues.
43
c *ei a singly subscripted real array containg the imaginary
44
c parts of the eigenvalues.
46
c *bs a singly subscripted integer array that contains block
47
c structure information. if there is a block of order
48
c k starting at a(l,l) in the output matrix a, then
49
c bs(l) contains the positive integer k, bs(l+1) contains
50
c -(k-1), bs(la+2) = -(k-2), ..., bs(l+k-1) = -1.
51
c thus a positive integer in the l-th entry of bs
52
c indicates a new block of order bs(l) starting at a(l,l).
54
c *x contains, either right reducing transformation u*y,
55
c either orthogonal tranformation u (see job)
57
c *xi xi contains the inverse reducing matrix transformation
58
c or y matrix (see job)
60
c *scale contains the scale factor and definitions of p size(n)
62
c job integer parametre specifying outputed transformations
63
c job=0 : a contains block diagonal form
64
c x right transformation
66
c job=1:like job=0 and xi contain x**-1
67
c job=2 a contains block diagonal form
68
c x contains u and xi contain y
70
c -block diagonal form in the diagonal blocks
71
c -a factorisation of y in the upper triangular
74
c *fail a logical variable which is false on normal return and
75
c true if there is any error in bdiag.
79
c orthes ortran (eispack)
80
c hqror2 exch split (eispack.extensions)
82
c real dble abs (fortran)
87
integer lda, n, bs, job
88
double precision a,er,ei,x,xi,rmax,epsshr,scale(n)
89
dimension a(lda,n),x(lda,n),xi(lda,n),er(n),ei(n),bs(n)
92
double precision c,cav,d,e1,e2,rav,temp,zero,one,mone,ddot,eps
93
double precision dlamch
94
integer da11,da22,i,j,k,km1,km2,l11,l22,l22m1,nk,ino
96
data zero, one, mone /0.0d+0,1.0d+0,-1.0d+0/
103
c compute eps the l1 norm of the a matrix
109
temp=temp+abs(a(i,j))
113
if (eps.eq.0.0d0) eps=1.0d0
118
c convert a to upper hessenberg form.
120
call balanc(lda, n, a, low, igh, scale)
121
call orthes(lda, n, low, igh, a, er)
122
call ortran(lda, n, low, igh, a, er, x)
124
c convert a to quasi-upper triangular form by qr method.
126
call hqror2(lda,n,1,n,a,er,ei,x,ierr,21)
128
c check to see if hqror2 failed in computing any eigenvalue
131
if(ierr.gt.1) goto 600
133
c reduce a to block diagonal form
136
c segment a into 4 matrices: a11, a da11 x da11 block
137
c whose (1,1)-element is at a(l11,l11)) a22, a da22 x da22
138
c block whose (1,1)-element is at a(l22,l22)) a12,
139
c a da11 x da22 block whose (1,1)-element is at a(l11,l22))
140
c and a21, a da22 x da11 block = 0 whose (1,1)-
141
c element is at a(l22,l11).
145
c this loop uses l11 as loop index and splits off a block
146
c starting at a(l11,l11).
151
if (l11.gt.n) go to 350
154
c this loop uses da11 as loop variable and attempts to split
155
c off a block of size da11 starting at a(l11,l11)
158
if (l22.ne.l11) go to 60
160
if(l11 .eq. n) go to 51
161
if(abs(a(l11+1,l11)) .gt.eps ) then
171
c compute the average of the eigenvalues in a11
177
cav = cav + abs(ei(i))
179
rav = rav/dble(real(da11) )
180
cav = cav/dble(real(da11) )
182
c loop on eigenvalues of a22 to find the one closest to the av
184
d = (rav-er(l22))**2 + (cav-ei(l22))**2
187
if(l22 .eq. n) go to 71
188
if(abs(a(l22+1,l22)) .gt. eps) l = l22 + 2
191
if (l.gt.n) go to 100
192
c = (rav-er(l))**2 + (cav-ei(l))**2
199
if (abs(a(l,l-1)).gt.eps) l=l+1
204
c loop to move the eigenvalue just located
205
c into first position of block a22.
208
if (abs(a(k+1,k)).gt.eps) go to 150
210
c the block we're moving to add to a11 is a 1 x 1
214
if (k.eq.l22) go to 230
216
if (abs(a(km1,k-2)).lt.eps) go to 140
218
c we're swapping the closest block with a 2 x 2
221
call exch(lda,n,a, x, km2, 2, 1)
223
c try to split this block into 2 real eigenvalues
225
call split(a, x, n, km1, e1, e2, lda, lda)
226
if (a(k,km1).eq.zero) go to 120
228
c block is still complex.
238
c complex block split into two real eigenvalues.
247
if (k.le.l22) go to 230
251
c we're swapping the closest block with a 1 x 1.
254
call exch(lda,n,a, x, km1, 1, 1)
259
if (k.le.l22) go to 230
262
c the block we're moving to add to a11 is a 2 x 2.
267
if (k.eq.l22) go to 230
269
if (abs(a(km1,k-2)).lt.eps) goto 190
271
c we're swapping the closest block with a 2 x 2 block.
274
call exch(lda,n,a, x, km2, 2, 2)
276
c try to split swapped block into two reals.
278
call split(a, x, n, k, e1, e2, lda, lda)
283
if (a(k+1,k).eq.zero) go to 170
285
c still complex block.
302
if (k.eq.l22) go to 210
305
c we're swapping the closest block with a 1 x 1.
308
call exch(lda,n,a, x, km1, 1, 2)
318
if (k.eq.l22) go to 210
321
c try to split relocated complex block.
324
call split(a, x, n, k, e1, e2, lda, lda)
325
if (a(k+1,k).eq.zero) go to 220
335
c split into two real eigenvalues.
348
if (l22.gt.n) go to 290
350
c attempt to split off a block of size da11.
354
c save a12 in its transpose form in block a21.
363
c convert a11 to lower quasi-triangular and multiply it by -1 and
364
c a12 appropriately (for solving -a11*p+p*a22=a12).
366
call dad(a, lda, l11, l22m1, l11, n, one, 0)
367
call dad(a, lda, l11, l22m1, l11, l22m1, mone, 1)
369
c solve -a11*p + p*a22 = a12.
371
call shrslv(a(l11,l11), a(l22,l22), a(l11,l22), da11,
372
* da22, lda, lda, lda, eps,epsshr,rmax, fails)
373
if (.not.fails) go to 290
375
c change a11 back to upper quasi-triangular.
377
call dad(a, lda, l11, l22m1, l11, l22m1, one, 1)
378
call dad(a, lda, l11, l22m1, l11, l22m1, mone, 0)
380
c was unable to solve for p - try again
383
c move saved a12 back into its correct position.
396
c change solution to p to proper form.
398
if (l22.gt.n) go to 300
399
call dad(a, lda, l11, l22m1, l11, n, one, 0)
400
call dad(a, lda, l11, l22m1, l11, l22m1, mone, 1)
403
c store block size in array bs.
407
if (j.eq.0) go to 320
410
bs(l11pi) = -(da11-i)
418
c set transformations matrices as required
422
c compute inverse transformation
423
if(job.ne.1) goto 450
431
if(l22.gt.n) goto 431
435
xi(i,j)=xi(i,j)-ddot(n-l22m1,a(i,l22),lda,xi(l22,j),1)
438
c in-lines back-tranfc in-lines right transformations of xi
440
if (igh .ne. low) then
442
temp=1.0d+00/scale(j)
450
if (i.ge.low .and. i.le.igh) goto 445
451
if (i.lt.low) i=low-ii
462
if(job.eq.2) goto 500
463
c compute right transformation
467
if(l22.gt.n) goto 480
470
x(i,j)=x(i,j)+ddot(l22-l11,x(i,l11),lda,a(l11,j),1)
475
call balbak( lda, n, low, igh, scale, n, x)
478
c extract non orthogonal tranformation from a
481
call dset(n,zero,xi(1,j),1)
483
call dset(n,one,xi(1,1),lda+1)
486
if(l11.gt.n) goto 550
490
xi(i,j)=xi(i,j)+ddot(l22-l11,xi(i,l11),lda,a(l11,j),1)
494
c set zeros in the matrix a
500
call dset(n-l22m1,zero,a(j,l22),lda)
501
call dset(n-l22m1,zero,a(l22,j),1)
added
removed
routines/control/bdiag.f
