LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ chbtrd()

subroutine chbtrd ( character vect,
character uplo,
integer n,
integer kd,
complex, dimension( ldab, * ) ab,
integer ldab,
real, dimension( * ) d,
real, dimension( * ) e,
complex, dimension( ldq, * ) q,
integer ldq,
complex, dimension( * ) work,
integer info )

CHBTRD

Download CHBTRD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> CHBTRD reduces a complex Hermitian band matrix A to real symmetric
!> tridiagonal form T by a unitary similarity transformation:
!> Q**H * A * Q = T.
!>
Parameters
[in]VECT
!>          VECT is CHARACTER*1
!>          = 'N':  do not form Q;
!>          = 'V':  form Q;
!>          = 'U':  update a matrix X, by forming X*Q.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in,out]AB
!>          AB is COMPLEX array, dimension (LDAB,N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
!>          On exit, the diagonal elements of AB are overwritten by the
!>          diagonal elements of the tridiagonal matrix T; if KD > 0, the
!>          elements on the first superdiagonal (if UPLO = 'U') or the
!>          first subdiagonal (if UPLO = 'L') are overwritten by the
!>          off-diagonal elements of T; the rest of AB is overwritten by
!>          values generated during the reduction.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD+1.
!>
[out]D
!>          D is REAL array, dimension (N)
!>          The diagonal elements of the tridiagonal matrix T.
!>
[out]E
!>          E is REAL array, dimension (N-1)
!>          The off-diagonal elements of the tridiagonal matrix T:
!>          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
!>
[in,out]Q
!>          Q is COMPLEX array, dimension (LDQ,N)
!>          On entry, if VECT = 'U', then Q must contain an N-by-N
!>          matrix X; if VECT = 'N' or 'V', then Q need not be set.
!>
!>          On exit:
!>          if VECT = 'V', Q contains the N-by-N unitary matrix Q;
!>          if VECT = 'U', Q contains the product X*Q;
!>          if VECT = 'N', the array Q is not referenced.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q.
!>          LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  Modified by Linda Kaufman, Bell Labs.
!>

Definition at line 159 of file chbtrd.f.

161*
162* -- LAPACK computational routine --
163* -- LAPACK is a software package provided by Univ. of Tennessee, --
164* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165*
166* .. Scalar Arguments ..
167 CHARACTER UPLO, VECT
168 INTEGER INFO, KD, LDAB, LDQ, N
169* ..
170* .. Array Arguments ..
171 REAL D( * ), E( * )
172 COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * )
173* ..
174*
175* =====================================================================
176*
177* .. Parameters ..
178 REAL ZERO
179 parameter( zero = 0.0e+0 )
180 COMPLEX CZERO, CONE
181 parameter( czero = ( 0.0e+0, 0.0e+0 ),
182 $ cone = ( 1.0e+0, 0.0e+0 ) )
183* ..
184* .. Local Scalars ..
185 LOGICAL INITQ, UPPER, WANTQ
186 INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
187 $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
188 $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
189 REAL ABST
190 COMPLEX T, TEMP
191* ..
192* .. External Subroutines ..
193 EXTERNAL clacgv, clar2v, clargv, clartg, clartv,
194 $ claset,
195 $ crot, cscal, xerbla
196* ..
197* .. Intrinsic Functions ..
198 INTRINSIC abs, conjg, max, min, real
199* ..
200* .. External Functions ..
201 LOGICAL LSAME
202 EXTERNAL lsame
203* ..
204* .. Executable Statements ..
205*
206* Test the input parameters
207*
208 initq = lsame( vect, 'V' )
209 wantq = initq .OR. lsame( vect, 'U' )
210 upper = lsame( uplo, 'U' )
211 kd1 = kd + 1
212 kdm1 = kd - 1
213 incx = ldab - 1
214 iqend = 1
215*
216 info = 0
217 IF( .NOT.wantq .AND. .NOT.lsame( vect, 'N' ) ) THEN
218 info = -1
219 ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
220 info = -2
221 ELSE IF( n.LT.0 ) THEN
222 info = -3
223 ELSE IF( kd.LT.0 ) THEN
224 info = -4
225 ELSE IF( ldab.LT.kd1 ) THEN
226 info = -6
227 ELSE IF( ldq.LT.max( 1, n ) .AND. wantq ) THEN
228 info = -10
229 END IF
230 IF( info.NE.0 ) THEN
231 CALL xerbla( 'CHBTRD', -info )
232 RETURN
233 END IF
234*
235* Quick return if possible
236*
237 IF( n.EQ.0 )
238 $ RETURN
239*
240* Initialize Q to the unit matrix, if needed
241*
242 IF( initq )
243 $ CALL claset( 'Full', n, n, czero, cone, q, ldq )
244*
245* Wherever possible, plane rotations are generated and applied in
246* vector operations of length NR over the index set J1:J2:KD1.
247*
248* The real cosines and complex sines of the plane rotations are
249* stored in the arrays D and WORK.
250*
251 inca = kd1*ldab
252 kdn = min( n-1, kd )
253 IF( upper ) THEN
254*
255 IF( kd.GT.1 ) THEN
256*
257* Reduce to complex Hermitian tridiagonal form, working with
258* the upper triangle
259*
260 nr = 0
261 j1 = kdn + 2
262 j2 = 1
263*
264 ab( kd1, 1 ) = real( ab( kd1, 1 ) )
265 DO 90 i = 1, n - 2
266*
267* Reduce i-th row of matrix to tridiagonal form
268*
269 DO 80 k = kdn + 1, 2, -1
270 j1 = j1 + kdn
271 j2 = j2 + kdn
272*
273 IF( nr.GT.0 ) THEN
274*
275* generate plane rotations to annihilate nonzero
276* elements which have been created outside the band
277*
278 CALL clargv( nr, ab( 1, j1-1 ), inca,
279 $ work( j1 ),
280 $ kd1, d( j1 ), kd1 )
281*
282* apply rotations from the right
283*
284*
285* Dependent on the the number of diagonals either
286* CLARTV or CROT is used
287*
288 IF( nr.GE.2*kd-1 ) THEN
289 DO 10 l = 1, kd - 1
290 CALL clartv( nr, ab( l+1, j1-1 ), inca,
291 $ ab( l, j1 ), inca, d( j1 ),
292 $ work( j1 ), kd1 )
293 10 CONTINUE
294*
295 ELSE
296 jend = j1 + ( nr-1 )*kd1
297 DO 20 jinc = j1, jend, kd1
298 CALL crot( kdm1, ab( 2, jinc-1 ), 1,
299 $ ab( 1, jinc ), 1, d( jinc ),
300 $ work( jinc ) )
301 20 CONTINUE
302 END IF
303 END IF
304*
305*
306 IF( k.GT.2 ) THEN
307 IF( k.LE.n-i+1 ) THEN
308*
309* generate plane rotation to annihilate a(i,i+k-1)
310* within the band
311*
312 CALL clartg( ab( kd-k+3, i+k-2 ),
313 $ ab( kd-k+2, i+k-1 ), d( i+k-1 ),
314 $ work( i+k-1 ), temp )
315 ab( kd-k+3, i+k-2 ) = temp
316*
317* apply rotation from the right
318*
319 CALL crot( k-3, ab( kd-k+4, i+k-2 ), 1,
320 $ ab( kd-k+3, i+k-1 ), 1, d( i+k-1 ),
321 $ work( i+k-1 ) )
322 END IF
323 nr = nr + 1
324 j1 = j1 - kdn - 1
325 END IF
326*
327* apply plane rotations from both sides to diagonal
328* blocks
329*
330 IF( nr.GT.0 )
331 $ CALL clar2v( nr, ab( kd1, j1-1 ), ab( kd1, j1 ),
332 $ ab( kd, j1 ), inca, d( j1 ),
333 $ work( j1 ), kd1 )
334*
335* apply plane rotations from the left
336*
337 IF( nr.GT.0 ) THEN
338 CALL clacgv( nr, work( j1 ), kd1 )
339 IF( 2*kd-1.LT.nr ) THEN
340*
341* Dependent on the the number of diagonals either
342* CLARTV or CROT is used
343*
344 DO 30 l = 1, kd - 1
345 IF( j2+l.GT.n ) THEN
346 nrt = nr - 1
347 ELSE
348 nrt = nr
349 END IF
350 IF( nrt.GT.0 )
351 $ CALL clartv( nrt, ab( kd-l, j1+l ),
352 $ inca,
353 $ ab( kd-l+1, j1+l ), inca,
354 $ d( j1 ), work( j1 ), kd1 )
355 30 CONTINUE
356 ELSE
357 j1end = j1 + kd1*( nr-2 )
358 IF( j1end.GE.j1 ) THEN
359 DO 40 jin = j1, j1end, kd1
360 CALL crot( kd-1, ab( kd-1, jin+1 ),
361 $ incx,
362 $ ab( kd, jin+1 ), incx,
363 $ d( jin ), work( jin ) )
364 40 CONTINUE
365 END IF
366 lend = min( kdm1, n-j2 )
367 last = j1end + kd1
368 IF( lend.GT.0 )
369 $ CALL crot( lend, ab( kd-1, last+1 ), incx,
370 $ ab( kd, last+1 ), incx, d( last ),
371 $ work( last ) )
372 END IF
373 END IF
374*
375 IF( wantq ) THEN
376*
377* accumulate product of plane rotations in Q
378*
379 IF( initq ) THEN
380*
381* take advantage of the fact that Q was
382* initially the Identity matrix
383*
384 iqend = max( iqend, j2 )
385 i2 = max( 0, k-3 )
386 iqaend = 1 + i*kd
387 IF( k.EQ.2 )
388 $ iqaend = iqaend + kd
389 iqaend = min( iqaend, iqend )
390 DO 50 j = j1, j2, kd1
391 ibl = i - i2 / kdm1
392 i2 = i2 + 1
393 iqb = max( 1, j-ibl )
394 nq = 1 + iqaend - iqb
395 iqaend = min( iqaend+kd, iqend )
396 CALL crot( nq, q( iqb, j-1 ), 1, q( iqb,
397 $ j ),
398 $ 1, d( j ), conjg( work( j ) ) )
399 50 CONTINUE
400 ELSE
401*
402 DO 60 j = j1, j2, kd1
403 CALL crot( n, q( 1, j-1 ), 1, q( 1, j ),
404 $ 1,
405 $ d( j ), conjg( work( j ) ) )
406 60 CONTINUE
407 END IF
408*
409 END IF
410*
411 IF( j2+kdn.GT.n ) THEN
412*
413* adjust J2 to keep within the bounds of the matrix
414*
415 nr = nr - 1
416 j2 = j2 - kdn - 1
417 END IF
418*
419 DO 70 j = j1, j2, kd1
420*
421* create nonzero element a(j-1,j+kd) outside the band
422* and store it in WORK
423*
424 work( j+kd ) = work( j )*ab( 1, j+kd )
425 ab( 1, j+kd ) = d( j )*ab( 1, j+kd )
426 70 CONTINUE
427 80 CONTINUE
428 90 CONTINUE
429 END IF
430*
431 IF( kd.GT.0 ) THEN
432*
433* make off-diagonal elements real and copy them to E
434*
435 DO 100 i = 1, n - 1
436 t = ab( kd, i+1 )
437 abst = abs( t )
438 ab( kd, i+1 ) = abst
439 e( i ) = abst
440 IF( abst.NE.zero ) THEN
441 t = t / abst
442 ELSE
443 t = cone
444 END IF
445 IF( i.LT.n-1 )
446 $ ab( kd, i+2 ) = ab( kd, i+2 )*t
447 IF( wantq ) THEN
448 CALL cscal( n, conjg( t ), q( 1, i+1 ), 1 )
449 END IF
450 100 CONTINUE
451 ELSE
452*
453* set E to zero if original matrix was diagonal
454*
455 DO 110 i = 1, n - 1
456 e( i ) = zero
457 110 CONTINUE
458 END IF
459*
460* copy diagonal elements to D
461*
462 DO 120 i = 1, n
463 d( i ) = real( ab( kd1, i ) )
464 120 CONTINUE
465*
466 ELSE
467*
468 IF( kd.GT.1 ) THEN
469*
470* Reduce to complex Hermitian tridiagonal form, working with
471* the lower triangle
472*
473 nr = 0
474 j1 = kdn + 2
475 j2 = 1
476*
477 ab( 1, 1 ) = real( ab( 1, 1 ) )
478 DO 210 i = 1, n - 2
479*
480* Reduce i-th column of matrix to tridiagonal form
481*
482 DO 200 k = kdn + 1, 2, -1
483 j1 = j1 + kdn
484 j2 = j2 + kdn
485*
486 IF( nr.GT.0 ) THEN
487*
488* generate plane rotations to annihilate nonzero
489* elements which have been created outside the band
490*
491 CALL clargv( nr, ab( kd1, j1-kd1 ), inca,
492 $ work( j1 ), kd1, d( j1 ), kd1 )
493*
494* apply plane rotations from one side
495*
496*
497* Dependent on the the number of diagonals either
498* CLARTV or CROT is used
499*
500 IF( nr.GT.2*kd-1 ) THEN
501 DO 130 l = 1, kd - 1
502 CALL clartv( nr, ab( kd1-l, j1-kd1+l ),
503 $ inca,
504 $ ab( kd1-l+1, j1-kd1+l ), inca,
505 $ d( j1 ), work( j1 ), kd1 )
506 130 CONTINUE
507 ELSE
508 jend = j1 + kd1*( nr-1 )
509 DO 140 jinc = j1, jend, kd1
510 CALL crot( kdm1, ab( kd, jinc-kd ), incx,
511 $ ab( kd1, jinc-kd ), incx,
512 $ d( jinc ), work( jinc ) )
513 140 CONTINUE
514 END IF
515*
516 END IF
517*
518 IF( k.GT.2 ) THEN
519 IF( k.LE.n-i+1 ) THEN
520*
521* generate plane rotation to annihilate a(i+k-1,i)
522* within the band
523*
524 CALL clartg( ab( k-1, i ), ab( k, i ),
525 $ d( i+k-1 ), work( i+k-1 ), temp )
526 ab( k-1, i ) = temp
527*
528* apply rotation from the left
529*
530 CALL crot( k-3, ab( k-2, i+1 ), ldab-1,
531 $ ab( k-1, i+1 ), ldab-1, d( i+k-1 ),
532 $ work( i+k-1 ) )
533 END IF
534 nr = nr + 1
535 j1 = j1 - kdn - 1
536 END IF
537*
538* apply plane rotations from both sides to diagonal
539* blocks
540*
541 IF( nr.GT.0 )
542 $ CALL clar2v( nr, ab( 1, j1-1 ), ab( 1, j1 ),
543 $ ab( 2, j1-1 ), inca, d( j1 ),
544 $ work( j1 ), kd1 )
545*
546* apply plane rotations from the right
547*
548*
549* Dependent on the the number of diagonals either
550* CLARTV or CROT is used
551*
552 IF( nr.GT.0 ) THEN
553 CALL clacgv( nr, work( j1 ), kd1 )
554 IF( nr.GT.2*kd-1 ) THEN
555 DO 150 l = 1, kd - 1
556 IF( j2+l.GT.n ) THEN
557 nrt = nr - 1
558 ELSE
559 nrt = nr
560 END IF
561 IF( nrt.GT.0 )
562 $ CALL clartv( nrt, ab( l+2, j1-1 ),
563 $ inca,
564 $ ab( l+1, j1 ), inca, d( j1 ),
565 $ work( j1 ), kd1 )
566 150 CONTINUE
567 ELSE
568 j1end = j1 + kd1*( nr-2 )
569 IF( j1end.GE.j1 ) THEN
570 DO 160 j1inc = j1, j1end, kd1
571 CALL crot( kdm1, ab( 3, j1inc-1 ), 1,
572 $ ab( 2, j1inc ), 1, d( j1inc ),
573 $ work( j1inc ) )
574 160 CONTINUE
575 END IF
576 lend = min( kdm1, n-j2 )
577 last = j1end + kd1
578 IF( lend.GT.0 )
579 $ CALL crot( lend, ab( 3, last-1 ), 1,
580 $ ab( 2, last ), 1, d( last ),
581 $ work( last ) )
582 END IF
583 END IF
584*
585*
586*
587 IF( wantq ) THEN
588*
589* accumulate product of plane rotations in Q
590*
591 IF( initq ) THEN
592*
593* take advantage of the fact that Q was
594* initially the Identity matrix
595*
596 iqend = max( iqend, j2 )
597 i2 = max( 0, k-3 )
598 iqaend = 1 + i*kd
599 IF( k.EQ.2 )
600 $ iqaend = iqaend + kd
601 iqaend = min( iqaend, iqend )
602 DO 170 j = j1, j2, kd1
603 ibl = i - i2 / kdm1
604 i2 = i2 + 1
605 iqb = max( 1, j-ibl )
606 nq = 1 + iqaend - iqb
607 iqaend = min( iqaend+kd, iqend )
608 CALL crot( nq, q( iqb, j-1 ), 1, q( iqb,
609 $ j ),
610 $ 1, d( j ), work( j ) )
611 170 CONTINUE
612 ELSE
613*
614 DO 180 j = j1, j2, kd1
615 CALL crot( n, q( 1, j-1 ), 1, q( 1, j ),
616 $ 1,
617 $ d( j ), work( j ) )
618 180 CONTINUE
619 END IF
620 END IF
621*
622 IF( j2+kdn.GT.n ) THEN
623*
624* adjust J2 to keep within the bounds of the matrix
625*
626 nr = nr - 1
627 j2 = j2 - kdn - 1
628 END IF
629*
630 DO 190 j = j1, j2, kd1
631*
632* create nonzero element a(j+kd,j-1) outside the
633* band and store it in WORK
634*
635 work( j+kd ) = work( j )*ab( kd1, j )
636 ab( kd1, j ) = d( j )*ab( kd1, j )
637 190 CONTINUE
638 200 CONTINUE
639 210 CONTINUE
640 END IF
641*
642 IF( kd.GT.0 ) THEN
643*
644* make off-diagonal elements real and copy them to E
645*
646 DO 220 i = 1, n - 1
647 t = ab( 2, i )
648 abst = abs( t )
649 ab( 2, i ) = abst
650 e( i ) = abst
651 IF( abst.NE.zero ) THEN
652 t = t / abst
653 ELSE
654 t = cone
655 END IF
656 IF( i.LT.n-1 )
657 $ ab( 2, i+1 ) = ab( 2, i+1 )*t
658 IF( wantq ) THEN
659 CALL cscal( n, t, q( 1, i+1 ), 1 )
660 END IF
661 220 CONTINUE
662 ELSE
663*
664* set E to zero if original matrix was diagonal
665*
666 DO 230 i = 1, n - 1
667 e( i ) = zero
668 230 CONTINUE
669 END IF
670*
671* copy diagonal elements to D
672*
673 DO 240 i = 1, n
674 d( i ) = real( ab( 1, i ) )
675 240 CONTINUE
676 END IF
677*
678 RETURN
679*
680* End of CHBTRD
681*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
clacgv
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:72
clar2v
subroutine clar2v(n, x, y, z, incx, c, s, incc)
CLAR2V applies a vector of plane rotations with real cosines and complex sines from both sides to a s...
Definition clar2v.f:109
clargv
subroutine clargv(n, x, incx, y, incy, c, incc)
CLARGV generates a vector of plane rotations with real cosines and complex sines.
Definition clargv.f:120
clartg
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition clartg.f90:116
clartv
subroutine clartv(n, x, incx, y, incy, c, s, incc)
CLARTV applies a vector of plane rotations with real cosines and complex sines to the elements of a p...
Definition clartv.f:105
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:104
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
crot
subroutine crot(n, cx, incx, cy, incy, c, s)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition crot.f:101
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
Here is the call graph for this function:
Here is the caller graph for this function: