LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
chbtrd.f
Go to the documentation of this file.
1 *> \brief \b CHBTRD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CHBTRD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbtrd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbtrd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbtrd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
22 * WORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER UPLO, VECT
26 * INTEGER INFO, KD, LDAB, LDQ, N
27 * ..
28 * .. Array Arguments ..
29 * REAL D( * ), E( * )
30 * COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> CHBTRD reduces a complex Hermitian band matrix A to real symmetric
40 *> tridiagonal form T by a unitary similarity transformation:
41 *> Q**H * A * Q = T.
42 *> \endverbatim
43 *
44 * Arguments:
45 * ==========
46 *
47 *> \param[in] VECT
48 *> \verbatim
49 *> VECT is CHARACTER*1
50 *> = 'N': do not form Q;
51 *> = 'V': form Q;
52 *> = 'U': update a matrix X, by forming X*Q.
53 *> \endverbatim
54 *>
55 *> \param[in] UPLO
56 *> \verbatim
57 *> UPLO is CHARACTER*1
58 *> = 'U': Upper triangle of A is stored;
59 *> = 'L': Lower triangle of A is stored.
60 *> \endverbatim
61 *>
62 *> \param[in] N
63 *> \verbatim
64 *> N is INTEGER
65 *> The order of the matrix A. N >= 0.
66 *> \endverbatim
67 *>
68 *> \param[in] KD
69 *> \verbatim
70 *> KD is INTEGER
71 *> The number of superdiagonals of the matrix A if UPLO = 'U',
72 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
73 *> \endverbatim
74 *>
75 *> \param[in,out] AB
76 *> \verbatim
77 *> AB is COMPLEX array, dimension (LDAB,N)
78 *> On entry, the upper or lower triangle of the Hermitian band
79 *> matrix A, stored in the first KD+1 rows of the array. The
80 *> j-th column of A is stored in the j-th column of the array AB
81 *> as follows:
82 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
83 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
84 *> On exit, the diagonal elements of AB are overwritten by the
85 *> diagonal elements of the tridiagonal matrix T; if KD > 0, the
86 *> elements on the first superdiagonal (if UPLO = 'U') or the
87 *> first subdiagonal (if UPLO = 'L') are overwritten by the
88 *> off-diagonal elements of T; the rest of AB is overwritten by
89 *> values generated during the reduction.
90 *> \endverbatim
91 *>
92 *> \param[in] LDAB
93 *> \verbatim
94 *> LDAB is INTEGER
95 *> The leading dimension of the array AB. LDAB >= KD+1.
96 *> \endverbatim
97 *>
98 *> \param[out] D
99 *> \verbatim
100 *> D is REAL array, dimension (N)
101 *> The diagonal elements of the tridiagonal matrix T.
102 *> \endverbatim
103 *>
104 *> \param[out] E
105 *> \verbatim
106 *> E is REAL array, dimension (N-1)
107 *> The off-diagonal elements of the tridiagonal matrix T:
108 *> E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
109 *> \endverbatim
110 *>
111 *> \param[in,out] Q
112 *> \verbatim
113 *> Q is COMPLEX array, dimension (LDQ,N)
114 *> On entry, if VECT = 'U', then Q must contain an N-by-N
115 *> matrix X; if VECT = 'N' or 'V', then Q need not be set.
116 *>
117 *> On exit:
118 *> if VECT = 'V', Q contains the N-by-N unitary matrix Q;
119 *> if VECT = 'U', Q contains the product X*Q;
120 *> if VECT = 'N', the array Q is not referenced.
121 *> \endverbatim
122 *>
123 *> \param[in] LDQ
124 *> \verbatim
125 *> LDQ is INTEGER
126 *> The leading dimension of the array Q.
127 *> LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
128 *> \endverbatim
129 *>
130 *> \param[out] WORK
131 *> \verbatim
132 *> WORK is COMPLEX array, dimension (N)
133 *> \endverbatim
134 *>
135 *> \param[out] INFO
136 *> \verbatim
137 *> INFO is INTEGER
138 *> = 0: successful exit
139 *> < 0: if INFO = -i, the i-th argument had an illegal value
140 *> \endverbatim
141 *
142 * Authors:
143 * ========
144 *
145 *> \author Univ. of Tennessee
146 *> \author Univ. of California Berkeley
147 *> \author Univ. of Colorado Denver
148 *> \author NAG Ltd.
149 *
150 *> \date November 2011
151 *
152 *> \ingroup complexOTHERcomputational
153 *
154 *> \par Further Details:
155 * =====================
156 *>
157 *> \verbatim
158 *>
159 *> Modified by Linda Kaufman, Bell Labs.
160 *> \endverbatim
161 *>
162 * =====================================================================
163  SUBROUTINE chbtrd( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
164  $ work, info )
165 *
166 * -- LAPACK computational routine (version 3.4.0) --
167 * -- LAPACK is a software package provided by Univ. of Tennessee, --
168 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169 * November 2011
170 *
171 * .. Scalar Arguments ..
172  CHARACTER uplo, vect
173  INTEGER info, kd, ldab, ldq, n
174 * ..
175 * .. Array Arguments ..
176  REAL d( * ), e( * )
177  COMPLEX ab( ldab, * ), q( ldq, * ), work( * )
178 * ..
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183  REAL zero
184  parameter( zero = 0.0e+0 )
185  COMPLEX czero, cone
186  parameter( czero = ( 0.0e+0, 0.0e+0 ),
187  $ cone = ( 1.0e+0, 0.0e+0 ) )
188 * ..
189 * .. Local Scalars ..
190  LOGICAL initq, upper, wantq
191  INTEGER i, i2, ibl, inca, incx, iqaend, iqb, iqend, j,
192  $ j1, j1end, j1inc, j2, jend, jin, jinc, k, kd1,
193  $ kdm1, kdn, l, last, lend, nq, nr, nrt
194  REAL abst
195  COMPLEX t, temp
196 * ..
197 * .. External Subroutines ..
198  EXTERNAL clacgv, clar2v, clargv, clartg, clartv, claset,
199  $ crot, cscal, xerbla
200 * ..
201 * .. Intrinsic Functions ..
202  INTRINSIC abs, conjg, max, min, real
203 * ..
204 * .. External Functions ..
205  LOGICAL lsame
206  EXTERNAL lsame
207 * ..
208 * .. Executable Statements ..
209 *
210 * Test the input parameters
211 *
212  initq = lsame( vect, 'V' )
213  wantq = initq .OR. lsame( vect, 'U' )
214  upper = lsame( uplo, 'U' )
215  kd1 = kd + 1
216  kdm1 = kd - 1
217  incx = ldab - 1
218  iqend = 1
219 *
220  info = 0
221  IF( .NOT.wantq .AND. .NOT.lsame( vect, 'N' ) ) THEN
222  info = -1
223  ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
224  info = -2
225  ELSE IF( n.LT.0 ) THEN
226  info = -3
227  ELSE IF( kd.LT.0 ) THEN
228  info = -4
229  ELSE IF( ldab.LT.kd1 ) THEN
230  info = -6
231  ELSE IF( ldq.LT.max( 1, n ) .AND. wantq ) THEN
232  info = -10
233  END IF
234  IF( info.NE.0 ) THEN
235  CALL xerbla( 'CHBTRD', -info )
236  return
237  END IF
238 *
239 * Quick return if possible
240 *
241  IF( n.EQ.0 )
242  $ return
243 *
244 * Initialize Q to the unit matrix, if needed
245 *
246  IF( initq )
247  $ CALL claset( 'Full', n, n, czero, cone, q, ldq )
248 *
249 * Wherever possible, plane rotations are generated and applied in
250 * vector operations of length NR over the index set J1:J2:KD1.
251 *
252 * The real cosines and complex sines of the plane rotations are
253 * stored in the arrays D and WORK.
254 *
255  inca = kd1*ldab
256  kdn = min( n-1, kd )
257  IF( upper ) THEN
258 *
259  IF( kd.GT.1 ) THEN
260 *
261 * Reduce to complex Hermitian tridiagonal form, working with
262 * the upper triangle
263 *
264  nr = 0
265  j1 = kdn + 2
266  j2 = 1
267 *
268  ab( kd1, 1 ) = REAL( AB( KD1, 1 ) )
269  DO 90 i = 1, n - 2
270 *
271 * Reduce i-th row of matrix to tridiagonal form
272 *
273  DO 80 k = kdn + 1, 2, -1
274  j1 = j1 + kdn
275  j2 = j2 + kdn
276 *
277  IF( nr.GT.0 ) THEN
278 *
279 * generate plane rotations to annihilate nonzero
280 * elements which have been created outside the band
281 *
282  CALL clargv( nr, ab( 1, j1-1 ), inca, work( j1 ),
283  $ kd1, d( j1 ), kd1 )
284 *
285 * apply rotations from the right
286 *
287 *
288 * Dependent on the the number of diagonals either
289 * CLARTV or CROT is used
290 *
291  IF( nr.GE.2*kd-1 ) THEN
292  DO 10 l = 1, kd - 1
293  CALL clartv( nr, ab( l+1, j1-1 ), inca,
294  $ ab( l, j1 ), inca, d( j1 ),
295  $ work( j1 ), kd1 )
296  10 continue
297 *
298  ELSE
299  jend = j1 + ( nr-1 )*kd1
300  DO 20 jinc = j1, jend, kd1
301  CALL crot( kdm1, ab( 2, jinc-1 ), 1,
302  $ ab( 1, jinc ), 1, d( jinc ),
303  $ work( jinc ) )
304  20 continue
305  END IF
306  END IF
307 *
308 *
309  IF( k.GT.2 ) THEN
310  IF( k.LE.n-i+1 ) THEN
311 *
312 * generate plane rotation to annihilate a(i,i+k-1)
313 * within the band
314 *
315  CALL clartg( ab( kd-k+3, i+k-2 ),
316  $ ab( kd-k+2, i+k-1 ), d( i+k-1 ),
317  $ work( i+k-1 ), temp )
318  ab( kd-k+3, i+k-2 ) = temp
319 *
320 * apply rotation from the right
321 *
322  CALL crot( k-3, ab( kd-k+4, i+k-2 ), 1,
323  $ ab( kd-k+3, i+k-1 ), 1, d( i+k-1 ),
324  $ work( i+k-1 ) )
325  END IF
326  nr = nr + 1
327  j1 = j1 - kdn - 1
328  END IF
329 *
330 * apply plane rotations from both sides to diagonal
331 * blocks
332 *
333  IF( nr.GT.0 )
334  $ CALL clar2v( nr, ab( kd1, j1-1 ), ab( kd1, j1 ),
335  $ ab( kd, j1 ), inca, d( j1 ),
336  $ work( j1 ), kd1 )
337 *
338 * apply plane rotations from the left
339 *
340  IF( nr.GT.0 ) THEN
341  CALL clacgv( nr, work( j1 ), kd1 )
342  IF( 2*kd-1.LT.nr ) THEN
343 *
344 * Dependent on the the number of diagonals either
345 * CLARTV or CROT is used
346 *
347  DO 30 l = 1, kd - 1
348  IF( j2+l.GT.n ) THEN
349  nrt = nr - 1
350  ELSE
351  nrt = nr
352  END IF
353  IF( nrt.GT.0 )
354  $ CALL clartv( nrt, ab( kd-l, j1+l ), inca,
355  $ ab( kd-l+1, j1+l ), inca,
356  $ d( j1 ), work( j1 ), kd1 )
357  30 continue
358  ELSE
359  j1end = j1 + kd1*( nr-2 )
360  IF( j1end.GE.j1 ) THEN
361  DO 40 jin = j1, j1end, kd1
362  CALL crot( kd-1, ab( kd-1, jin+1 ), incx,
363  $ ab( kd, jin+1 ), incx,
364  $ d( jin ), work( jin ) )
365  40 continue
366  END IF
367  lend = min( kdm1, n-j2 )
368  last = j1end + kd1
369  IF( lend.GT.0 )
370  $ CALL crot( lend, ab( kd-1, last+1 ), incx,
371  $ ab( kd, last+1 ), incx, d( last ),
372  $ work( last ) )
373  END IF
374  END IF
375 *
376  IF( wantq ) THEN
377 *
378 * accumulate product of plane rotations in Q
379 *
380  IF( initq ) THEN
381 *
382 * take advantage of the fact that Q was
383 * initially the Identity matrix
384 *
385  iqend = max( iqend, j2 )
386  i2 = max( 0, k-3 )
387  iqaend = 1 + i*kd
388  IF( k.EQ.2 )
389  $ iqaend = iqaend + kd
390  iqaend = min( iqaend, iqend )
391  DO 50 j = j1, j2, kd1
392  ibl = i - i2 / kdm1
393  i2 = i2 + 1
394  iqb = max( 1, j-ibl )
395  nq = 1 + iqaend - iqb
396  iqaend = min( iqaend+kd, iqend )
397  CALL crot( nq, q( iqb, j-1 ), 1, q( iqb, 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 ), 1,
404  $ d( j ), conjg( work( j ) ) )
405  60 continue
406  END IF
407 *
408  END IF
409 *
410  IF( j2+kdn.GT.n ) THEN
411 *
412 * adjust J2 to keep within the bounds of the matrix
413 *
414  nr = nr - 1
415  j2 = j2 - kdn - 1
416  END IF
417 *
418  DO 70 j = j1, j2, kd1
419 *
420 * create nonzero element a(j-1,j+kd) outside the band
421 * and store it in WORK
422 *
423  work( j+kd ) = work( j )*ab( 1, j+kd )
424  ab( 1, j+kd ) = d( j )*ab( 1, j+kd )
425  70 continue
426  80 continue
427  90 continue
428  END IF
429 *
430  IF( kd.GT.0 ) THEN
431 *
432 * make off-diagonal elements real and copy them to E
433 *
434  DO 100 i = 1, n - 1
435  t = ab( kd, i+1 )
436  abst = abs( t )
437  ab( kd, i+1 ) = abst
438  e( i ) = abst
439  IF( abst.NE.zero ) THEN
440  t = t / abst
441  ELSE
442  t = cone
443  END IF
444  IF( i.LT.n-1 )
445  $ ab( kd, i+2 ) = ab( kd, i+2 )*t
446  IF( wantq ) THEN
447  CALL cscal( n, conjg( t ), q( 1, i+1 ), 1 )
448  END IF
449  100 continue
450  ELSE
451 *
452 * set E to zero if original matrix was diagonal
453 *
454  DO 110 i = 1, n - 1
455  e( i ) = zero
456  110 continue
457  END IF
458 *
459 * copy diagonal elements to D
460 *
461  DO 120 i = 1, n
462  d( i ) = ab( kd1, i )
463  120 continue
464 *
465  ELSE
466 *
467  IF( kd.GT.1 ) THEN
468 *
469 * Reduce to complex Hermitian tridiagonal form, working with
470 * the lower triangle
471 *
472  nr = 0
473  j1 = kdn + 2
474  j2 = 1
475 *
476  ab( 1, 1 ) = REAL( AB( 1, 1 ) )
477  DO 210 i = 1, n - 2
478 *
479 * Reduce i-th column of matrix to tridiagonal form
480 *
481  DO 200 k = kdn + 1, 2, -1
482  j1 = j1 + kdn
483  j2 = j2 + kdn
484 *
485  IF( nr.GT.0 ) THEN
486 *
487 * generate plane rotations to annihilate nonzero
488 * elements which have been created outside the band
489 *
490  CALL clargv( nr, ab( kd1, j1-kd1 ), inca,
491  $ work( j1 ), kd1, d( j1 ), kd1 )
492 *
493 * apply plane rotations from one side
494 *
495 *
496 * Dependent on the the number of diagonals either
497 * CLARTV or CROT is used
498 *
499  IF( nr.GT.2*kd-1 ) THEN
500  DO 130 l = 1, kd - 1
501  CALL clartv( nr, ab( kd1-l, j1-kd1+l ), inca,
502  $ ab( kd1-l+1, j1-kd1+l ), inca,
503  $ d( j1 ), work( j1 ), kd1 )
504  130 continue
505  ELSE
506  jend = j1 + kd1*( nr-1 )
507  DO 140 jinc = j1, jend, kd1
508  CALL crot( kdm1, ab( kd, jinc-kd ), incx,
509  $ ab( kd1, jinc-kd ), incx,
510  $ d( jinc ), work( jinc ) )
511  140 continue
512  END IF
513 *
514  END IF
515 *
516  IF( k.GT.2 ) THEN
517  IF( k.LE.n-i+1 ) THEN
518 *
519 * generate plane rotation to annihilate a(i+k-1,i)
520 * within the band
521 *
522  CALL clartg( ab( k-1, i ), ab( k, i ),
523  $ d( i+k-1 ), work( i+k-1 ), temp )
524  ab( k-1, i ) = temp
525 *
526 * apply rotation from the left
527 *
528  CALL crot( k-3, ab( k-2, i+1 ), ldab-1,
529  $ ab( k-1, i+1 ), ldab-1, d( i+k-1 ),
530  $ work( i+k-1 ) )
531  END IF
532  nr = nr + 1
533  j1 = j1 - kdn - 1
534  END IF
535 *
536 * apply plane rotations from both sides to diagonal
537 * blocks
538 *
539  IF( nr.GT.0 )
540  $ CALL clar2v( nr, ab( 1, j1-1 ), ab( 1, j1 ),
541  $ ab( 2, j1-1 ), inca, d( j1 ),
542  $ work( j1 ), kd1 )
543 *
544 * apply plane rotations from the right
545 *
546 *
547 * Dependent on the the number of diagonals either
548 * CLARTV or CROT is used
549 *
550  IF( nr.GT.0 ) THEN
551  CALL clacgv( nr, work( j1 ), kd1 )
552  IF( nr.GT.2*kd-1 ) THEN
553  DO 150 l = 1, kd - 1
554  IF( j2+l.GT.n ) THEN
555  nrt = nr - 1
556  ELSE
557  nrt = nr
558  END IF
559  IF( nrt.GT.0 )
560  $ CALL clartv( nrt, ab( l+2, j1-1 ), inca,
561  $ ab( l+1, j1 ), inca, d( j1 ),
562  $ work( j1 ), kd1 )
563  150 continue
564  ELSE
565  j1end = j1 + kd1*( nr-2 )
566  IF( j1end.GE.j1 ) THEN
567  DO 160 j1inc = j1, j1end, kd1
568  CALL crot( kdm1, ab( 3, j1inc-1 ), 1,
569  $ ab( 2, j1inc ), 1, d( j1inc ),
570  $ work( j1inc ) )
571  160 continue
572  END IF
573  lend = min( kdm1, n-j2 )
574  last = j1end + kd1
575  IF( lend.GT.0 )
576  $ CALL crot( lend, ab( 3, last-1 ), 1,
577  $ ab( 2, last ), 1, d( last ),
578  $ work( last ) )
579  END IF
580  END IF
581 *
582 *
583 *
584  IF( wantq ) THEN
585 *
586 * accumulate product of plane rotations in Q
587 *
588  IF( initq ) THEN
589 *
590 * take advantage of the fact that Q was
591 * initially the Identity matrix
592 *
593  iqend = max( iqend, j2 )
594  i2 = max( 0, k-3 )
595  iqaend = 1 + i*kd
596  IF( k.EQ.2 )
597  $ iqaend = iqaend + kd
598  iqaend = min( iqaend, iqend )
599  DO 170 j = j1, j2, kd1
600  ibl = i - i2 / kdm1
601  i2 = i2 + 1
602  iqb = max( 1, j-ibl )
603  nq = 1 + iqaend - iqb
604  iqaend = min( iqaend+kd, iqend )
605  CALL crot( nq, q( iqb, j-1 ), 1, q( iqb, j ),
606  $ 1, d( j ), work( j ) )
607  170 continue
608  ELSE
609 *
610  DO 180 j = j1, j2, kd1
611  CALL crot( n, q( 1, j-1 ), 1, q( 1, j ), 1,
612  $ d( j ), work( j ) )
613  180 continue
614  END IF
615  END IF
616 *
617  IF( j2+kdn.GT.n ) THEN
618 *
619 * adjust J2 to keep within the bounds of the matrix
620 *
621  nr = nr - 1
622  j2 = j2 - kdn - 1
623  END IF
624 *
625  DO 190 j = j1, j2, kd1
626 *
627 * create nonzero element a(j+kd,j-1) outside the
628 * band and store it in WORK
629 *
630  work( j+kd ) = work( j )*ab( kd1, j )
631  ab( kd1, j ) = d( j )*ab( kd1, j )
632  190 continue
633  200 continue
634  210 continue
635  END IF
636 *
637  IF( kd.GT.0 ) THEN
638 *
639 * make off-diagonal elements real and copy them to E
640 *
641  DO 220 i = 1, n - 1
642  t = ab( 2, i )
643  abst = abs( t )
644  ab( 2, i ) = abst
645  e( i ) = abst
646  IF( abst.NE.zero ) THEN
647  t = t / abst
648  ELSE
649  t = cone
650  END IF
651  IF( i.LT.n-1 )
652  $ ab( 2, i+1 ) = ab( 2, i+1 )*t
653  IF( wantq ) THEN
654  CALL cscal( n, t, q( 1, i+1 ), 1 )
655  END IF
656  220 continue
657  ELSE
658 *
659 * set E to zero if original matrix was diagonal
660 *
661  DO 230 i = 1, n - 1
662  e( i ) = zero
663  230 continue
664  END IF
665 *
666 * copy diagonal elements to D
667 *
668  DO 240 i = 1, n
669  d( i ) = ab( 1, i )
670  240 continue
671  END IF
672 *
673  return
674 *
675 * End of CHBTRD
676 *
677  END