LAPACK 3.12.0
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 161 of file chbtrd.f.

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