Return to dlasd1.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
    2:      $                   IDXQ, IWORK, WORK, INFO )
    3: *
    4: *  -- LAPACK auxiliary routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       INTEGER            INFO, LDU, LDVT, NL, NR, SQRE
   11:       DOUBLE PRECISION   ALPHA, BETA
   12: *     ..
   13: *     .. Array Arguments ..
   14:       INTEGER            IDXQ( * ), IWORK( * )
   15:       DOUBLE PRECISION   D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
   16: *     ..
   17: *
   18: *  Purpose
   19: *  =======
   20: *
   21: *  DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
   22: *  where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
   23: *
   24: *  A related subroutine DLASD7 handles the case in which the singular
   25: *  values (and the singular vectors in factored form) are desired.
   26: *
   27: *  DLASD1 computes the SVD as follows:
   28: *
   29: *                ( D1(in)  0    0     0 )
   30: *    B = U(in) * (   Z1'   a   Z2'    b ) * VT(in)
   31: *                (   0     0   D2(in) 0 )
   32: *
   33: *      = U(out) * ( D(out) 0) * VT(out)
   34: *
   35: *  where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M
   36: *  with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
   37: *  elsewhere; and the entry b is empty if SQRE = 0.
   38: *
   39: *  The left singular vectors of the original matrix are stored in U, and
   40: *  the transpose of the right singular vectors are stored in VT, and the
   41: *  singular values are in D.  The algorithm consists of three stages:
   42: *
   43: *     The first stage consists of deflating the size of the problem
   44: *     when there are multiple singular values or when there are zeros in
   45: *     the Z vector.  For each such occurence the dimension of the
   46: *     secular equation problem is reduced by one.  This stage is
   47: *     performed by the routine DLASD2.
   48: *
   49: *     The second stage consists of calculating the updated
   50: *     singular values. This is done by finding the square roots of the
   51: *     roots of the secular equation via the routine DLASD4 (as called
   52: *     by DLASD3). This routine also calculates the singular vectors of
   53: *     the current problem.
   54: *
   55: *     The final stage consists of computing the updated singular vectors
   56: *     directly using the updated singular values.  The singular vectors
   57: *     for the current problem are multiplied with the singular vectors
   58: *     from the overall problem.
   59: *
   60: *  Arguments
   61: *  =========
   62: *
   63: *  NL     (input) INTEGER
   64: *         The row dimension of the upper block.  NL >= 1.
   65: *
   66: *  NR     (input) INTEGER
   67: *         The row dimension of the lower block.  NR >= 1.
   68: *
   69: *  SQRE   (input) INTEGER
   70: *         = 0: the lower block is an NR-by-NR square matrix.
   71: *         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
   72: *
   73: *         The bidiagonal matrix has row dimension N = NL + NR + 1,
   74: *         and column dimension M = N + SQRE.
   75: *
   76: *  D      (input/output) DOUBLE PRECISION array,
   77: *                        dimension (N = NL+NR+1).
   78: *         On entry D(1:NL,1:NL) contains the singular values of the
   79: *         upper block; and D(NL+2:N) contains the singular values of
   80: *         the lower block. On exit D(1:N) contains the singular values
   81: *         of the modified matrix.
   82: *
   83: *  ALPHA  (input/output) DOUBLE PRECISION
   84: *         Contains the diagonal element associated with the added row.
   85: *
   86: *  BETA   (input/output) DOUBLE PRECISION
   87: *         Contains the off-diagonal element associated with the added
   88: *         row.
   89: *
   90: *  U      (input/output) DOUBLE PRECISION array, dimension(LDU,N)
   91: *         On entry U(1:NL, 1:NL) contains the left singular vectors of
   92: *         the upper block; U(NL+2:N, NL+2:N) contains the left singular
   93: *         vectors of the lower block. On exit U contains the left
   94: *         singular vectors of the bidiagonal matrix.
   95: *
   96: *  LDU    (input) INTEGER
   97: *         The leading dimension of the array U.  LDU >= max( 1, N ).
   98: *
   99: *  VT     (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
  100: *         where M = N + SQRE.
  101: *         On entry VT(1:NL+1, 1:NL+1)' contains the right singular
  102: *         vectors of the upper block; VT(NL+2:M, NL+2:M)' contains
  103: *         the right singular vectors of the lower block. On exit
  104: *         VT' contains the right singular vectors of the
  105: *         bidiagonal matrix.
  106: *
  107: *  LDVT   (input) INTEGER
  108: *         The leading dimension of the array VT.  LDVT >= max( 1, M ).
  109: *
  110: *  IDXQ  (output) INTEGER array, dimension(N)
  111: *         This contains the permutation which will reintegrate the
  112: *         subproblem just solved back into sorted order, i.e.
  113: *         D( IDXQ( I = 1, N ) ) will be in ascending order.
  114: *
  115: *  IWORK  (workspace) INTEGER array, dimension( 4 * N )
  116: *
  117: *  WORK   (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
  118: *
  119: *  INFO   (output) INTEGER
  120: *          = 0:  successful exit.
  121: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  122: *          > 0:  if INFO = 1, an singular value did not converge
  123: *
  124: *  Further Details
  125: *  ===============
  126: *
  127: *  Based on contributions by
  128: *     Ming Gu and Huan Ren, Computer Science Division, University of
  129: *     California at Berkeley, USA
  130: *
  131: *  =====================================================================
  132: *
  133: *     .. Parameters ..
  134: *
  135:       DOUBLE PRECISION   ONE, ZERO
  136:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  137: *     ..
  138: *     .. Local Scalars ..
  139:       INTEGER            COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
  140:      $                   IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
  141:       DOUBLE PRECISION   ORGNRM
  142: *     ..
  143: *     .. External Subroutines ..
  144:       EXTERNAL           DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA
  145: *     ..
  146: *     .. Intrinsic Functions ..
  147:       INTRINSIC          ABS, MAX
  148: *     ..
  149: *     .. Executable Statements ..
  150: *
  151: *     Test the input parameters.
  152: *
  153:       INFO = 0
  154: *
  155:       IF( NL.LT.1 ) THEN
  156:          INFO = -1
  157:       ELSE IF( NR.LT.1 ) THEN
  158:          INFO = -2
  159:       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
  160:          INFO = -3
  161:       END IF
  162:       IF( INFO.NE.0 ) THEN
  163:          CALL XERBLA( 'DLASD1', -INFO )
  164:          RETURN
  165:       END IF
  166: *
  167:       N = NL + NR + 1
  168:       M = N + SQRE
  169: *
  170: *     The following values are for bookkeeping purposes only.  They are
  171: *     integer pointers which indicate the portion of the workspace
  172: *     used by a particular array in DLASD2 and DLASD3.
  173: *
  174:       LDU2 = N
  175:       LDVT2 = M
  176: *
  177:       IZ = 1
  178:       ISIGMA = IZ + M
  179:       IU2 = ISIGMA + N
  180:       IVT2 = IU2 + LDU2*N
  181:       IQ = IVT2 + LDVT2*M
  182: *
  183:       IDX = 1
  184:       IDXC = IDX + N
  185:       COLTYP = IDXC + N
  186:       IDXP = COLTYP + N
  187: *
  188: *     Scale.
  189: *
  190:       ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
  191:       D( NL+1 ) = ZERO
  192:       DO 10 I = 1, N
  193:          IF( ABS( D( I ) ).GT.ORGNRM ) THEN
  194:             ORGNRM = ABS( D( I ) )
  195:          END IF
  196:    10 CONTINUE
  197:       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
  198:       ALPHA = ALPHA / ORGNRM
  199:       BETA = BETA / ORGNRM
  200: *
  201: *     Deflate singular values.
  202: *
  203:       CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU,
  204:      $             VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2,
  205:      $             WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ),
  206:      $             IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO )
  207: *
  208: *     Solve Secular Equation and update singular vectors.
  209: *
  210:       LDQ = K
  211:       CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ),
  212:      $             U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ),
  213:      $             LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ),
  214:      $             INFO )
  215:       IF( INFO.NE.0 ) THEN
  216:          RETURN
  217:       END IF
  218: *
  219: *     Unscale.
  220: *
  221:       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
  222: *
  223: *     Prepare the IDXQ sorting permutation.
  224: *
  225:       N1 = K
  226:       N2 = N - K
  227:       CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
  228: *
  229:       RETURN
  230: *
  231: *     End of DLASD1
  232: *
  233:       END