Return to dgelsx.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack.

    1: *> \brief <b> DGELSX solves overdetermined or underdetermined systems for GE matrices</b>
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DGELSX + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelsx.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelsx.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelsx.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
   22: *                          WORK, INFO )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       INTEGER            INFO, LDA, LDB, M, N, NRHS, RANK
   26: *       DOUBLE PRECISION   RCOND
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       INTEGER            JPVT( * )
   30: *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
   31: *       ..
   32: *  
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> This routine is deprecated and has been replaced by routine DGELSY.
   40: *>
   41: *> DGELSX computes the minimum-norm solution to a real linear least
   42: *> squares problem:
   43: *>     minimize || A * X - B ||
   44: *> using a complete orthogonal factorization of A.  A is an M-by-N
   45: *> matrix which may be rank-deficient.
   46: *>
   47: *> Several right hand side vectors b and solution vectors x can be
   48: *> handled in a single call; they are stored as the columns of the
   49: *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
   50: *> matrix X.
   51: *>
   52: *> The routine first computes a QR factorization with column pivoting:
   53: *>     A * P = Q * [ R11 R12 ]
   54: *>                 [  0  R22 ]
   55: *> with R11 defined as the largest leading submatrix whose estimated
   56: *> condition number is less than 1/RCOND.  The order of R11, RANK,
   57: *> is the effective rank of A.
   58: *>
   59: *> Then, R22 is considered to be negligible, and R12 is annihilated
   60: *> by orthogonal transformations from the right, arriving at the
   61: *> complete orthogonal factorization:
   62: *>    A * P = Q * [ T11 0 ] * Z
   63: *>                [  0  0 ]
   64: *> The minimum-norm solution is then
   65: *>    X = P * Z**T [ inv(T11)*Q1**T*B ]
   66: *>                 [        0         ]
   67: *> where Q1 consists of the first RANK columns of Q.
   68: *> \endverbatim
   69: *
   70: *  Arguments:
   71: *  ==========
   72: *
   73: *> \param[in] M
   74: *> \verbatim
   75: *>          M is INTEGER
   76: *>          The number of rows of the matrix A.  M >= 0.
   77: *> \endverbatim
   78: *>
   79: *> \param[in] N
   80: *> \verbatim
   81: *>          N is INTEGER
   82: *>          The number of columns of the matrix A.  N >= 0.
   83: *> \endverbatim
   84: *>
   85: *> \param[in] NRHS
   86: *> \verbatim
   87: *>          NRHS is INTEGER
   88: *>          The number of right hand sides, i.e., the number of
   89: *>          columns of matrices B and X. NRHS >= 0.
   90: *> \endverbatim
   91: *>
   92: *> \param[in,out] A
   93: *> \verbatim
   94: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   95: *>          On entry, the M-by-N matrix A.
   96: *>          On exit, A has been overwritten by details of its
   97: *>          complete orthogonal factorization.
   98: *> \endverbatim
   99: *>
  100: *> \param[in] LDA
  101: *> \verbatim
  102: *>          LDA is INTEGER
  103: *>          The leading dimension of the array A.  LDA >= max(1,M).
  104: *> \endverbatim
  105: *>
  106: *> \param[in,out] B
  107: *> \verbatim
  108: *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
  109: *>          On entry, the M-by-NRHS right hand side matrix B.
  110: *>          On exit, the N-by-NRHS solution matrix X.
  111: *>          If m >= n and RANK = n, the residual sum-of-squares for
  112: *>          the solution in the i-th column is given by the sum of
  113: *>          squares of elements N+1:M in that column.
  114: *> \endverbatim
  115: *>
  116: *> \param[in] LDB
  117: *> \verbatim
  118: *>          LDB is INTEGER
  119: *>          The leading dimension of the array B. LDB >= max(1,M,N).
  120: *> \endverbatim
  121: *>
  122: *> \param[in,out] JPVT
  123: *> \verbatim
  124: *>          JPVT is INTEGER array, dimension (N)
  125: *>          On entry, if JPVT(i) .ne. 0, the i-th column of A is an
  126: *>          initial column, otherwise it is a free column.  Before
  127: *>          the QR factorization of A, all initial columns are
  128: *>          permuted to the leading positions; only the remaining
  129: *>          free columns are moved as a result of column pivoting
  130: *>          during the factorization.
  131: *>          On exit, if JPVT(i) = k, then the i-th column of A*P
  132: *>          was the k-th column of A.
  133: *> \endverbatim
  134: *>
  135: *> \param[in] RCOND
  136: *> \verbatim
  137: *>          RCOND is DOUBLE PRECISION
  138: *>          RCOND is used to determine the effective rank of A, which
  139: *>          is defined as the order of the largest leading triangular
  140: *>          submatrix R11 in the QR factorization with pivoting of A,
  141: *>          whose estimated condition number < 1/RCOND.
  142: *> \endverbatim
  143: *>
  144: *> \param[out] RANK
  145: *> \verbatim
  146: *>          RANK is INTEGER
  147: *>          The effective rank of A, i.e., the order of the submatrix
  148: *>          R11.  This is the same as the order of the submatrix T11
  149: *>          in the complete orthogonal factorization of A.
  150: *> \endverbatim
  151: *>
  152: *> \param[out] WORK
  153: *> \verbatim
  154: *>          WORK is DOUBLE PRECISION array, dimension
  155: *>                      (max( min(M,N)+3*N, 2*min(M,N)+NRHS )),
  156: *> \endverbatim
  157: *>
  158: *> \param[out] INFO
  159: *> \verbatim
  160: *>          INFO is INTEGER
  161: *>          = 0:  successful exit
  162: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  163: *> \endverbatim
  164: *
  165: *  Authors:
  166: *  ========
  167: *
  168: *> \author Univ. of Tennessee 
  169: *> \author Univ. of California Berkeley 
  170: *> \author Univ. of Colorado Denver 
  171: *> \author NAG Ltd. 
  172: *
  173: *> \date November 2011
  174: *
  175: *> \ingroup doubleGEsolve
  176: *
  177: *  =====================================================================
  178:       SUBROUTINE DGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
  179:      $                   WORK, INFO )
  180: *
  181: *  -- LAPACK driver routine (version 3.4.0) --
  182: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  183: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  184: *     November 2011
  185: *
  186: *     .. Scalar Arguments ..
  187:       INTEGER            INFO, LDA, LDB, M, N, NRHS, RANK
  188:       DOUBLE PRECISION   RCOND
  189: *     ..
  190: *     .. Array Arguments ..
  191:       INTEGER            JPVT( * )
  192:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
  193: *     ..
  194: *
  195: *  =====================================================================
  196: *
  197: *     .. Parameters ..
  198:       INTEGER            IMAX, IMIN
  199:       PARAMETER          ( IMAX = 1, IMIN = 2 )
  200:       DOUBLE PRECISION   ZERO, ONE, DONE, NTDONE
  201:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, DONE = ZERO,
  202:      $                   NTDONE = ONE )
  203: *     ..
  204: *     .. Local Scalars ..
  205:       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
  206:       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
  207:      $                   SMAXPR, SMIN, SMINPR, SMLNUM, T1, T2
  208: *     ..
  209: *     .. External Functions ..
  210:       DOUBLE PRECISION   DLAMCH, DLANGE
  211:       EXTERNAL           DLAMCH, DLANGE
  212: *     ..
  213: *     .. External Subroutines ..
  214:       EXTERNAL           DGEQPF, DLAIC1, DLASCL, DLASET, DLATZM, DORM2R,
  215:      $                   DTRSM, DTZRQF, XERBLA
  216: *     ..
  217: *     .. Intrinsic Functions ..
  218:       INTRINSIC          ABS, MAX, MIN
  219: *     ..
  220: *     .. Executable Statements ..
  221: *
  222:       MN = MIN( M, N )
  223:       ISMIN = MN + 1
  224:       ISMAX = 2*MN + 1
  225: *
  226: *     Test the input arguments.
  227: *
  228:       INFO = 0
  229:       IF( M.LT.0 ) THEN
  230:          INFO = -1
  231:       ELSE IF( N.LT.0 ) THEN
  232:          INFO = -2
  233:       ELSE IF( NRHS.LT.0 ) THEN
  234:          INFO = -3
  235:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  236:          INFO = -5
  237:       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
  238:          INFO = -7
  239:       END IF
  240: *
  241:       IF( INFO.NE.0 ) THEN
  242:          CALL XERBLA( 'DGELSX', -INFO )
  243:          RETURN
  244:       END IF
  245: *
  246: *     Quick return if possible
  247: *
  248:       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
  249:          RANK = 0
  250:          RETURN
  251:       END IF
  252: *
  253: *     Get machine parameters
  254: *
  255:       SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
  256:       BIGNUM = ONE / SMLNUM
  257:       CALL DLABAD( SMLNUM, BIGNUM )
  258: *
  259: *     Scale A, B if max elements outside range [SMLNUM,BIGNUM]
  260: *
  261:       ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
  262:       IASCL = 0
  263:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  264: *
  265: *        Scale matrix norm up to SMLNUM
  266: *
  267:          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
  268:          IASCL = 1
  269:       ELSE IF( ANRM.GT.BIGNUM ) THEN
  270: *
  271: *        Scale matrix norm down to BIGNUM
  272: *
  273:          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
  274:          IASCL = 2
  275:       ELSE IF( ANRM.EQ.ZERO ) THEN
  276: *
  277: *        Matrix all zero. Return zero solution.
  278: *
  279:          CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
  280:          RANK = 0
  281:          GO TO 100
  282:       END IF
  283: *
  284:       BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
  285:       IBSCL = 0
  286:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
  287: *
  288: *        Scale matrix norm up to SMLNUM
  289: *
  290:          CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
  291:          IBSCL = 1
  292:       ELSE IF( BNRM.GT.BIGNUM ) THEN
  293: *
  294: *        Scale matrix norm down to BIGNUM
  295: *
  296:          CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
  297:          IBSCL = 2
  298:       END IF
  299: *
  300: *     Compute QR factorization with column pivoting of A:
  301: *        A * P = Q * R
  302: *
  303:       CALL DGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), INFO )
  304: *
  305: *     workspace 3*N. Details of Householder rotations stored
  306: *     in WORK(1:MN).
  307: *
  308: *     Determine RANK using incremental condition estimation
  309: *
  310:       WORK( ISMIN ) = ONE
  311:       WORK( ISMAX ) = ONE
  312:       SMAX = ABS( A( 1, 1 ) )
  313:       SMIN = SMAX
  314:       IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
  315:          RANK = 0
  316:          CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
  317:          GO TO 100
  318:       ELSE
  319:          RANK = 1
  320:       END IF
  321: *
  322:    10 CONTINUE
  323:       IF( RANK.LT.MN ) THEN
  324:          I = RANK + 1
  325:          CALL DLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
  326:      $                A( I, I ), SMINPR, S1, C1 )
  327:          CALL DLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
  328:      $                A( I, I ), SMAXPR, S2, C2 )
  329: *
  330:          IF( SMAXPR*RCOND.LE.SMINPR ) THEN
  331:             DO 20 I = 1, RANK
  332:                WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
  333:                WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
  334:    20       CONTINUE
  335:             WORK( ISMIN+RANK ) = C1
  336:             WORK( ISMAX+RANK ) = C2
  337:             SMIN = SMINPR
  338:             SMAX = SMAXPR
  339:             RANK = RANK + 1
  340:             GO TO 10
  341:          END IF
  342:       END IF
  343: *
  344: *     Logically partition R = [ R11 R12 ]
  345: *                             [  0  R22 ]
  346: *     where R11 = R(1:RANK,1:RANK)
  347: *
  348: *     [R11,R12] = [ T11, 0 ] * Y
  349: *
  350:       IF( RANK.LT.N )
  351:      $   CALL DTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO )
  352: *
  353: *     Details of Householder rotations stored in WORK(MN+1:2*MN)
  354: *
  355: *     B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
  356: *
  357:       CALL DORM2R( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ),
  358:      $             B, LDB, WORK( 2*MN+1 ), INFO )
  359: *
  360: *     workspace NRHS
  361: *
  362: *     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
  363: *
  364:       CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
  365:      $            NRHS, ONE, A, LDA, B, LDB )
  366: *
  367:       DO 40 I = RANK + 1, N
  368:          DO 30 J = 1, NRHS
  369:             B( I, J ) = ZERO
  370:    30    CONTINUE
  371:    40 CONTINUE
  372: *
  373: *     B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
  374: *
  375:       IF( RANK.LT.N ) THEN
  376:          DO 50 I = 1, RANK
  377:             CALL DLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA,
  378:      $                   WORK( MN+I ), B( I, 1 ), B( RANK+1, 1 ), LDB,
  379:      $                   WORK( 2*MN+1 ) )
  380:    50    CONTINUE
  381:       END IF
  382: *
  383: *     workspace NRHS
  384: *
  385: *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
  386: *
  387:       DO 90 J = 1, NRHS
  388:          DO 60 I = 1, N
  389:             WORK( 2*MN+I ) = NTDONE
  390:    60    CONTINUE
  391:          DO 80 I = 1, N
  392:             IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN
  393:                IF( JPVT( I ).NE.I ) THEN
  394:                   K = I
  395:                   T1 = B( K, J )
  396:                   T2 = B( JPVT( K ), J )
  397:    70             CONTINUE
  398:                   B( JPVT( K ), J ) = T1
  399:                   WORK( 2*MN+K ) = DONE
  400:                   T1 = T2
  401:                   K = JPVT( K )
  402:                   T2 = B( JPVT( K ), J )
  403:                   IF( JPVT( K ).NE.I )
  404:      $               GO TO 70
  405:                   B( I, J ) = T1
  406:                   WORK( 2*MN+K ) = DONE
  407:                END IF
  408:             END IF
  409:    80    CONTINUE
  410:    90 CONTINUE
  411: *
  412: *     Undo scaling
  413: *
  414:       IF( IASCL.EQ.1 ) THEN
  415:          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
  416:          CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
  417:      $                INFO )
  418:       ELSE IF( IASCL.EQ.2 ) THEN
  419:          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
  420:          CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
  421:      $                INFO )
  422:       END IF
  423:       IF( IBSCL.EQ.1 ) THEN
  424:          CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
  425:       ELSE IF( IBSCL.EQ.2 ) THEN
  426:          CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
  427:       END IF
  428: *
  429:   100 CONTINUE
  430: *
  431:       RETURN
  432: *
  433: *     End of DGELSX
  434: *
  435:       END