SRC/zgbbrd.f

   1 *> \brief \b ZGBBRD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZGBBRD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbbrd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbbrd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbbrd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
  22 *                          LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          VECT
  26 *       INTEGER            INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
  30 *       COMPLEX*16         AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
  31 *      $                   Q( LDQ, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> ZGBBRD reduces a complex general m-by-n band matrix A to real upper
  41 *> bidiagonal form B by a unitary transformation: Q**H * A * P = B.
  42 *>
  43 *> The routine computes B, and optionally forms Q or P**H, or computes
  44 *> Q**H*C for a given matrix C.
  45 *> \endverbatim
  46 *
  47 *  Arguments:
  48 *  ==========
  49 *
  50 *> \param[in] VECT
  51 *> \verbatim
  52 *>          VECT is CHARACTER*1
  53 *>          Specifies whether or not the matrices Q and P**H are to be
  54 *>          formed.
  55 *>          = 'N': do not form Q or P**H;
  56 *>          = 'Q': form Q only;
  57 *>          = 'P': form P**H only;
  58 *>          = 'B': form both.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] M
  62 *> \verbatim
  63 *>          M is INTEGER
  64 *>          The number of rows of the matrix A.  M >= 0.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] N
  68 *> \verbatim
  69 *>          N is INTEGER
  70 *>          The number of columns of the matrix A.  N >= 0.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] NCC
  74 *> \verbatim
  75 *>          NCC is INTEGER
  76 *>          The number of columns of the matrix C.  NCC >= 0.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] KL
  80 *> \verbatim
  81 *>          KL is INTEGER
  82 *>          The number of subdiagonals of the matrix A. KL >= 0.
  83 *> \endverbatim
  84 *>
  85 *> \param[in] KU
  86 *> \verbatim
  87 *>          KU is INTEGER
  88 *>          The number of superdiagonals of the matrix A. KU >= 0.
  89 *> \endverbatim
  90 *>
  91 *> \param[in,out] AB
  92 *> \verbatim
  93 *>          AB is COMPLEX*16 array, dimension (LDAB,N)
  94 *>          On entry, the m-by-n band matrix A, stored in rows 1 to
  95 *>          KL+KU+1. The j-th column of A is stored in the j-th column of
  96 *>          the array AB as follows:
  97 *>          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
  98 *>          On exit, A is overwritten by values generated during the
  99 *>          reduction.
 100 *> \endverbatim
 101 *>
 102 *> \param[in] LDAB
 103 *> \verbatim
 104 *>          LDAB is INTEGER
 105 *>          The leading dimension of the array A. LDAB >= KL+KU+1.
 106 *> \endverbatim
 107 *>
 108 *> \param[out] D
 109 *> \verbatim
 110 *>          D is DOUBLE PRECISION array, dimension (min(M,N))
 111 *>          The diagonal elements of the bidiagonal matrix B.
 112 *> \endverbatim
 113 *>
 114 *> \param[out] E
 115 *> \verbatim
 116 *>          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
 117 *>          The superdiagonal elements of the bidiagonal matrix B.
 118 *> \endverbatim
 119 *>
 120 *> \param[out] Q
 121 *> \verbatim
 122 *>          Q is COMPLEX*16 array, dimension (LDQ,M)
 123 *>          If VECT = 'Q' or 'B', the m-by-m unitary matrix Q.
 124 *>          If VECT = 'N' or 'P', the array Q is not referenced.
 125 *> \endverbatim
 126 *>
 127 *> \param[in] LDQ
 128 *> \verbatim
 129 *>          LDQ is INTEGER
 130 *>          The leading dimension of the array Q.
 131 *>          LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise.
 132 *> \endverbatim
 133 *>
 134 *> \param[out] PT
 135 *> \verbatim
 136 *>          PT is COMPLEX*16 array, dimension (LDPT,N)
 137 *>          If VECT = 'P' or 'B', the n-by-n unitary matrix P'.
 138 *>          If VECT = 'N' or 'Q', the array PT is not referenced.
 139 *> \endverbatim
 140 *>
 141 *> \param[in] LDPT
 142 *> \verbatim
 143 *>          LDPT is INTEGER
 144 *>          The leading dimension of the array PT.
 145 *>          LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise.
 146 *> \endverbatim
 147 *>
 148 *> \param[in,out] C
 149 *> \verbatim
 150 *>          C is COMPLEX*16 array, dimension (LDC,NCC)
 151 *>          On entry, an m-by-ncc matrix C.
 152 *>          On exit, C is overwritten by Q**H*C.
 153 *>          C is not referenced if NCC = 0.
 154 *> \endverbatim
 155 *>
 156 *> \param[in] LDC
 157 *> \verbatim
 158 *>          LDC is INTEGER
 159 *>          The leading dimension of the array C.
 160 *>          LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0.
 161 *> \endverbatim
 162 *>
 163 *> \param[out] WORK
 164 *> \verbatim
 165 *>          WORK is COMPLEX*16 array, dimension (max(M,N))
 166 *> \endverbatim
 167 *>
 168 *> \param[out] RWORK
 169 *> \verbatim
 170 *>          RWORK is DOUBLE PRECISION array, dimension (max(M,N))
 171 *> \endverbatim
 172 *>
 173 *> \param[out] INFO
 174 *> \verbatim
 175 *>          INFO is INTEGER
 176 *>          = 0:  successful exit.
 177 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 178 *> \endverbatim
 179 *
 180 *  Authors:
 181 *  ========
 182 *
 183 *> \author Univ. of Tennessee
 184 *> \author Univ. of California Berkeley
 185 *> \author Univ. of Colorado Denver
 186 *> \author NAG Ltd.
 187 *
 188 *> \date November 2011
 189 *
 190 *> \ingroup complex16GBcomputational
 191 *
 192 *  =====================================================================
 193       SUBROUTINE ZGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
 194      $                   LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO )
 195 *
 196 *  -- LAPACK computational routine (version 3.4.0) --
 197 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 198 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 199 *     November 2011
 200 *
 201 *     .. Scalar Arguments ..
 202       CHARACTER          VECT
 203       INTEGER            INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
 204 *     ..
 205 *     .. Array Arguments ..
 206       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
 207       COMPLEX*16         AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
 208      $                   Q( LDQ, * ), WORK( * )
 209 *     ..
 210 *
 211 *  =====================================================================
 212 *
 213 *     .. Parameters ..
 214       DOUBLE PRECISION   ZERO
 215       PARAMETER          ( ZERO = 0.0D+0 )
 216       COMPLEX*16         CZERO, CONE
 217       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
 218      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
 219 *     ..
 220 *     .. Local Scalars ..
 221       LOGICAL            WANTB, WANTC, WANTPT, WANTQ
 222       INTEGER            I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1,
 223      $                   KUN, L, MINMN, ML, ML0, MU, MU0, NR, NRT
 224       DOUBLE PRECISION   ABST, RC
 225       COMPLEX*16         RA, RB, RS, T
 226 *     ..
 227 *     .. External Subroutines ..
 228       EXTERNAL           XERBLA, ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT,
 229      $                   ZSCAL
 230 *     ..
 231 *     .. Intrinsic Functions ..
 232       INTRINSIC          ABS, DCONJG, MAX, MIN
 233 *     ..
 234 *     .. External Functions ..
 235       LOGICAL            LSAME
 236       EXTERNAL           LSAME
 237 *     ..
 238 *     .. Executable Statements ..
 239 *
 240 *     Test the input parameters
 241 *
 242       WANTB = LSAME( VECT, 'B' )
 243       WANTQ = LSAME( VECT, 'Q' ) .OR. WANTB
 244       WANTPT = LSAME( VECT, 'P' ) .OR. WANTB
 245       WANTC = NCC.GT.0
 246       KLU1 = KL + KU + 1
 247       INFO = 0
 248       IF( .NOT.WANTQ .AND. .NOT.WANTPT .AND. .NOT.LSAME( VECT, 'N' ) )
 249      $     THEN
 250          INFO = -1
 251       ELSE IF( M.LT.0 ) THEN
 252          INFO = -2
 253       ELSE IF( N.LT.0 ) THEN
 254          INFO = -3
 255       ELSE IF( NCC.LT.0 ) THEN
 256          INFO = -4
 257       ELSE IF( KL.LT.0 ) THEN
 258          INFO = -5
 259       ELSE IF( KU.LT.0 ) THEN
 260          INFO = -6
 261       ELSE IF( LDAB.LT.KLU1 ) THEN
 262          INFO = -8
 263       ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. LDQ.LT.MAX( 1, M ) ) THEN
 264          INFO = -12
 265       ELSE IF( LDPT.LT.1 .OR. WANTPT .AND. LDPT.LT.MAX( 1, N ) ) THEN
 266          INFO = -14
 267       ELSE IF( LDC.LT.1 .OR. WANTC .AND. LDC.LT.MAX( 1, M ) ) THEN
 268          INFO = -16
 269       END IF
 270       IF( INFO.NE.0 ) THEN
 271          CALL XERBLA( 'ZGBBRD', -INFO )
 272          RETURN
 273       END IF
 274 *
 275 *     Initialize Q and P**H to the unit matrix, if needed
 276 *
 277       IF( WANTQ )
 278      $   CALL ZLASET( 'Full', M, M, CZERO, CONE, Q, LDQ )
 279       IF( WANTPT )
 280      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, PT, LDPT )
 281 *
 282 *     Quick return if possible.
 283 *
 284       IF( M.EQ.0 .OR. N.EQ.0 )
 285      $   RETURN
 286 *
 287       MINMN = MIN( M, N )
 288 *
 289       IF( KL+KU.GT.1 ) THEN
 290 *
 291 *        Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
 292 *        first to lower bidiagonal form and then transform to upper
 293 *        bidiagonal
 294 *
 295          IF( KU.GT.0 ) THEN
 296             ML0 = 1
 297             MU0 = 2
 298          ELSE
 299             ML0 = 2
 300             MU0 = 1
 301          END IF
 302 *
 303 *        Wherever possible, plane rotations are generated and applied in
 304 *        vector operations of length NR over the index set J1:J2:KLU1.
 305 *
 306 *        The complex sines of the plane rotations are stored in WORK,
 307 *        and the real cosines in RWORK.
 308 *
 309          KLM = MIN( M-1, KL )
 310          KUN = MIN( N-1, KU )
 311          KB = KLM + KUN
 312          KB1 = KB + 1
 313          INCA = KB1*LDAB
 314          NR = 0
 315          J1 = KLM + 2
 316          J2 = 1 - KUN
 317 *
 318          DO 90 I = 1, MINMN
 319 *
 320 *           Reduce i-th column and i-th row of matrix to bidiagonal form
 321 *
 322             ML = KLM + 1
 323             MU = KUN + 1
 324             DO 80 KK = 1, KB
 325                J1 = J1 + KB
 326                J2 = J2 + KB
 327 *
 328 *              generate plane rotations to annihilate nonzero elements
 329 *              which have been created below the band
 330 *
 331                IF( NR.GT.0 )
 332      $            CALL ZLARGV( NR, AB( KLU1, J1-KLM-1 ), INCA,
 333      $                         WORK( J1 ), KB1, RWORK( J1 ), KB1 )
 334 *
 335 *              apply plane rotations from the left
 336 *
 337                DO 10 L = 1, KB
 338                   IF( J2-KLM+L-1.GT.N ) THEN
 339                      NRT = NR - 1
 340                   ELSE
 341                      NRT = NR
 342                   END IF
 343                   IF( NRT.GT.0 )
 344      $               CALL ZLARTV( NRT, AB( KLU1-L, J1-KLM+L-1 ), INCA,
 345      $                            AB( KLU1-L+1, J1-KLM+L-1 ), INCA,
 346      $                            RWORK( J1 ), WORK( J1 ), KB1 )
 347    10          CONTINUE
 348 *
 349                IF( ML.GT.ML0 ) THEN
 350                   IF( ML.LE.M-I+1 ) THEN
 351 *
 352 *                    generate plane rotation to annihilate a(i+ml-1,i)
 353 *                    within the band, and apply rotation from the left
 354 *
 355                      CALL ZLARTG( AB( KU+ML-1, I ), AB( KU+ML, I ),
 356      $                            RWORK( I+ML-1 ), WORK( I+ML-1 ), RA )
 357                      AB( KU+ML-1, I ) = RA
 358                      IF( I.LT.N )
 359      $                  CALL ZROT( MIN( KU+ML-2, N-I ),
 360      $                             AB( KU+ML-2, I+1 ), LDAB-1,
 361      $                             AB( KU+ML-1, I+1 ), LDAB-1,
 362      $                             RWORK( I+ML-1 ), WORK( I+ML-1 ) )
 363                   END IF
 364                   NR = NR + 1
 365                   J1 = J1 - KB1
 366                END IF
 367 *
 368                IF( WANTQ ) THEN
 369 *
 370 *                 accumulate product of plane rotations in Q
 371 *
 372                   DO 20 J = J1, J2, KB1
 373                      CALL ZROT( M, Q( 1, J-1 ), 1, Q( 1, J ), 1,
 374      $                          RWORK( J ), DCONJG( WORK( J ) ) )
 375    20             CONTINUE
 376                END IF
 377 *
 378                IF( WANTC ) THEN
 379 *
 380 *                 apply plane rotations to C
 381 *
 382                   DO 30 J = J1, J2, KB1
 383                      CALL ZROT( NCC, C( J-1, 1 ), LDC, C( J, 1 ), LDC,
 384      $                          RWORK( J ), WORK( J ) )
 385    30             CONTINUE
 386                END IF
 387 *
 388                IF( J2+KUN.GT.N ) THEN
 389 *
 390 *                 adjust J2 to keep within the bounds of the matrix
 391 *
 392                   NR = NR - 1
 393                   J2 = J2 - KB1
 394                END IF
 395 *
 396                DO 40 J = J1, J2, KB1
 397 *
 398 *                 create nonzero element a(j-1,j+ku) above the band
 399 *                 and store it in WORK(n+1:2*n)
 400 *
 401                   WORK( J+KUN ) = WORK( J )*AB( 1, J+KUN )
 402                   AB( 1, J+KUN ) = RWORK( J )*AB( 1, J+KUN )
 403    40          CONTINUE
 404 *
 405 *              generate plane rotations to annihilate nonzero elements
 406 *              which have been generated above the band
 407 *
 408                IF( NR.GT.0 )
 409      $            CALL ZLARGV( NR, AB( 1, J1+KUN-1 ), INCA,
 410      $                         WORK( J1+KUN ), KB1, RWORK( J1+KUN ),
 411      $                         KB1 )
 412 *
 413 *              apply plane rotations from the right
 414 *
 415                DO 50 L = 1, KB
 416                   IF( J2+L-1.GT.M ) THEN
 417                      NRT = NR - 1
 418                   ELSE
 419                      NRT = NR
 420                   END IF
 421                   IF( NRT.GT.0 )
 422      $               CALL ZLARTV( NRT, AB( L+1, J1+KUN-1 ), INCA,
 423      $                            AB( L, J1+KUN ), INCA,
 424      $                            RWORK( J1+KUN ), WORK( J1+KUN ), KB1 )
 425    50          CONTINUE
 426 *
 427                IF( ML.EQ.ML0 .AND. MU.GT.MU0 ) THEN
 428                   IF( MU.LE.N-I+1 ) THEN
 429 *
 430 *                    generate plane rotation to annihilate a(i,i+mu-1)
 431 *                    within the band, and apply rotation from the right
 432 *
 433                      CALL ZLARTG( AB( KU-MU+3, I+MU-2 ),
 434      $                            AB( KU-MU+2, I+MU-1 ),
 435      $                            RWORK( I+MU-1 ), WORK( I+MU-1 ), RA )
 436                      AB( KU-MU+3, I+MU-2 ) = RA
 437                      CALL ZROT( MIN( KL+MU-2, M-I ),
 438      $                          AB( KU-MU+4, I+MU-2 ), 1,
 439      $                          AB( KU-MU+3, I+MU-1 ), 1,
 440      $                          RWORK( I+MU-1 ), WORK( I+MU-1 ) )
 441                   END IF
 442                   NR = NR + 1
 443                   J1 = J1 - KB1
 444                END IF
 445 *
 446                IF( WANTPT ) THEN
 447 *
 448 *                 accumulate product of plane rotations in P**H
 449 *
 450                   DO 60 J = J1, J2, KB1
 451                      CALL ZROT( N, PT( J+KUN-1, 1 ), LDPT,
 452      $                          PT( J+KUN, 1 ), LDPT, RWORK( J+KUN ),
 453      $                          DCONJG( WORK( J+KUN ) ) )
 454    60             CONTINUE
 455                END IF
 456 *
 457                IF( J2+KB.GT.M ) THEN
 458 *
 459 *                 adjust J2 to keep within the bounds of the matrix
 460 *
 461                   NR = NR - 1
 462                   J2 = J2 - KB1
 463                END IF
 464 *
 465                DO 70 J = J1, J2, KB1
 466 *
 467 *                 create nonzero element a(j+kl+ku,j+ku-1) below the
 468 *                 band and store it in WORK(1:n)
 469 *
 470                   WORK( J+KB ) = WORK( J+KUN )*AB( KLU1, J+KUN )
 471                   AB( KLU1, J+KUN ) = RWORK( J+KUN )*AB( KLU1, J+KUN )
 472    70          CONTINUE
 473 *
 474                IF( ML.GT.ML0 ) THEN
 475                   ML = ML - 1
 476                ELSE
 477                   MU = MU - 1
 478                END IF
 479    80       CONTINUE
 480    90    CONTINUE
 481       END IF
 482 *
 483       IF( KU.EQ.0 .AND. KL.GT.0 ) THEN
 484 *
 485 *        A has been reduced to complex lower bidiagonal form
 486 *
 487 *        Transform lower bidiagonal form to upper bidiagonal by applying
 488 *        plane rotations from the left, overwriting superdiagonal
 489 *        elements on subdiagonal elements
 490 *
 491          DO 100 I = 1, MIN( M-1, N )
 492             CALL ZLARTG( AB( 1, I ), AB( 2, I ), RC, RS, RA )
 493             AB( 1, I ) = RA
 494             IF( I.LT.N ) THEN
 495                AB( 2, I ) = RS*AB( 1, I+1 )
 496                AB( 1, I+1 ) = RC*AB( 1, I+1 )
 497             END IF
 498             IF( WANTQ )
 499      $         CALL ZROT( M, Q( 1, I ), 1, Q( 1, I+1 ), 1, RC,
 500      $                    DCONJG( RS ) )
 501             IF( WANTC )
 502      $         CALL ZROT( NCC, C( I, 1 ), LDC, C( I+1, 1 ), LDC, RC,
 503      $                    RS )
 504   100    CONTINUE
 505       ELSE
 506 *
 507 *        A has been reduced to complex upper bidiagonal form or is
 508 *        diagonal
 509 *
 510          IF( KU.GT.0 .AND. M.LT.N ) THEN
 511 *
 512 *           Annihilate a(m,m+1) by applying plane rotations from the
 513 *           right
 514 *
 515             RB = AB( KU, M+1 )
 516             DO 110 I = M, 1, -1
 517                CALL ZLARTG( AB( KU+1, I ), RB, RC, RS, RA )
 518                AB( KU+1, I ) = RA
 519                IF( I.GT.1 ) THEN
 520                   RB = -DCONJG( RS )*AB( KU, I )
 521                   AB( KU, I ) = RC*AB( KU, I )
 522                END IF
 523                IF( WANTPT )
 524      $            CALL ZROT( N, PT( I, 1 ), LDPT, PT( M+1, 1 ), LDPT,
 525      $                       RC, DCONJG( RS ) )
 526   110       CONTINUE
 527          END IF
 528       END IF
 529 *
 530 *     Make diagonal and superdiagonal elements real, storing them in D
 531 *     and E
 532 *
 533       T = AB( KU+1, 1 )
 534       DO 120 I = 1, MINMN
 535          ABST = ABS( T )
 536          D( I ) = ABST
 537          IF( ABST.NE.ZERO ) THEN
 538             T = T / ABST
 539          ELSE
 540             T = CONE
 541          END IF
 542          IF( WANTQ )
 543      $      CALL ZSCAL( M, T, Q( 1, I ), 1 )
 544          IF( WANTC )
 545      $      CALL ZSCAL( NCC, DCONJG( T ), C( I, 1 ), LDC )
 546          IF( I.LT.MINMN ) THEN
 547             IF( KU.EQ.0 .AND. KL.EQ.0 ) THEN
 548                E( I ) = ZERO
 549                T = AB( 1, I+1 )
 550             ELSE
 551                IF( KU.EQ.0 ) THEN
 552                   T = AB( 2, I )*DCONJG( T )
 553                ELSE
 554                   T = AB( KU, I+1 )*DCONJG( T )
 555                END IF
 556                ABST = ABS( T )
 557                E( I ) = ABST
 558                IF( ABST.NE.ZERO ) THEN
 559                   T = T / ABST
 560                ELSE
 561                   T = CONE
 562                END IF
 563                IF( WANTPT )
 564      $            CALL ZSCAL( N, T, PT( I+1, 1 ), LDPT )
 565                T = AB( KU+1, I+1 )*DCONJG( T )
 566             END IF
 567          END IF
 568   120 CONTINUE
 569       RETURN
 570 *
 571 *     End of ZGBBRD
 572 *
 573       END