Return to dgbbrd.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    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