SRC/dgbbrd.f

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