LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f90:118
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:74
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
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
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
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
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
Here is the call graph for this function:
Here is the caller graph for this function: