Return to dlaqp2.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1: *> \brief \b DLAQP2
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DLAQP2 + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqp2.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqp2.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqp2.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
   22: *                          WORK )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       INTEGER            LDA, M, N, OFFSET
   26: *       ..
   27: *       .. Array Arguments ..
   28: *       INTEGER            JPVT( * )
   29: *       DOUBLE PRECISION   A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
   30: *      $                   WORK( * )
   31: *       ..
   32: *  
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> DLAQP2 computes a QR factorization with column pivoting of
   40: *> the block A(OFFSET+1:M,1:N).
   41: *> The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
   42: *> \endverbatim
   43: *
   44: *  Arguments:
   45: *  ==========
   46: *
   47: *> \param[in] M
   48: *> \verbatim
   49: *>          M is INTEGER
   50: *>          The number of rows of the matrix A. M >= 0.
   51: *> \endverbatim
   52: *>
   53: *> \param[in] N
   54: *> \verbatim
   55: *>          N is INTEGER
   56: *>          The number of columns of the matrix A. N >= 0.
   57: *> \endverbatim
   58: *>
   59: *> \param[in] OFFSET
   60: *> \verbatim
   61: *>          OFFSET is INTEGER
   62: *>          The number of rows of the matrix A that must be pivoted
   63: *>          but no factorized. OFFSET >= 0.
   64: *> \endverbatim
   65: *>
   66: *> \param[in,out] A
   67: *> \verbatim
   68: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   69: *>          On entry, the M-by-N matrix A.
   70: *>          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is 
   71: *>          the triangular factor obtained; the elements in block
   72: *>          A(OFFSET+1:M,1:N) below the diagonal, together with the
   73: *>          array TAU, represent the orthogonal matrix Q as a product of
   74: *>          elementary reflectors. Block A(1:OFFSET,1:N) has been
   75: *>          accordingly pivoted, but no factorized.
   76: *> \endverbatim
   77: *>
   78: *> \param[in] LDA
   79: *> \verbatim
   80: *>          LDA is INTEGER
   81: *>          The leading dimension of the array A. LDA >= max(1,M).
   82: *> \endverbatim
   83: *>
   84: *> \param[in,out] JPVT
   85: *> \verbatim
   86: *>          JPVT is INTEGER array, dimension (N)
   87: *>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
   88: *>          to the front of A*P (a leading column); if JPVT(i) = 0,
   89: *>          the i-th column of A is a free column.
   90: *>          On exit, if JPVT(i) = k, then the i-th column of A*P
   91: *>          was the k-th column of A.
   92: *> \endverbatim
   93: *>
   94: *> \param[out] TAU
   95: *> \verbatim
   96: *>          TAU is DOUBLE PRECISION array, dimension (min(M,N))
   97: *>          The scalar factors of the elementary reflectors.
   98: *> \endverbatim
   99: *>
  100: *> \param[in,out] VN1
  101: *> \verbatim
  102: *>          VN1 is DOUBLE PRECISION array, dimension (N)
  103: *>          The vector with the partial column norms.
  104: *> \endverbatim
  105: *>
  106: *> \param[in,out] VN2
  107: *> \verbatim
  108: *>          VN2 is DOUBLE PRECISION array, dimension (N)
  109: *>          The vector with the exact column norms.
  110: *> \endverbatim
  111: *>
  112: *> \param[out] WORK
  113: *> \verbatim
  114: *>          WORK is DOUBLE PRECISION array, dimension (N)
  115: *> \endverbatim
  116: *
  117: *  Authors:
  118: *  ========
  119: *
  120: *> \author Univ. of Tennessee 
  121: *> \author Univ. of California Berkeley 
  122: *> \author Univ. of Colorado Denver 
  123: *> \author NAG Ltd. 
  124: *
  125: *> \date November 2011
  126: *
  127: *> \ingroup doubleOTHERauxiliary
  128: *
  129: *> \par Contributors:
  130: *  ==================
  131: *>
  132: *>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
  133: *>    X. Sun, Computer Science Dept., Duke University, USA
  134: *> \n
  135: *>  Partial column norm updating strategy modified on April 2011
  136: *>    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
  137: *>    University of Zagreb, Croatia.
  138: *
  139: *> \par References:
  140: *  ================
  141: *>
  142: *> LAPACK Working Note 176
  143: *
  144: *> \htmlonly
  145: *> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a> 
  146: *> \endhtmlonly 
  147: *
  148: *  =====================================================================
  149:       SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
  150:      $                   WORK )
  151: *
  152: *  -- LAPACK auxiliary routine (version 3.4.0) --
  153: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  154: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  155: *     November 2011
  156: *
  157: *     .. Scalar Arguments ..
  158:       INTEGER            LDA, M, N, OFFSET
  159: *     ..
  160: *     .. Array Arguments ..
  161:       INTEGER            JPVT( * )
  162:       DOUBLE PRECISION   A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
  163:      $                   WORK( * )
  164: *     ..
  165: *
  166: *  =====================================================================
  167: *
  168: *     .. Parameters ..
  169:       DOUBLE PRECISION   ZERO, ONE
  170:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  171: *     ..
  172: *     .. Local Scalars ..
  173:       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
  174:       DOUBLE PRECISION   AII, TEMP, TEMP2, TOL3Z
  175: *     ..
  176: *     .. External Subroutines ..
  177:       EXTERNAL           DLARF, DLARFG, DSWAP
  178: *     ..
  179: *     .. Intrinsic Functions ..
  180:       INTRINSIC          ABS, MAX, MIN, SQRT
  181: *     ..
  182: *     .. External Functions ..
  183:       INTEGER            IDAMAX
  184:       DOUBLE PRECISION   DLAMCH, DNRM2
  185:       EXTERNAL           IDAMAX, DLAMCH, DNRM2
  186: *     ..
  187: *     .. Executable Statements ..
  188: *
  189:       MN = MIN( M-OFFSET, N )
  190:       TOL3Z = SQRT(DLAMCH('Epsilon'))
  191: *
  192: *     Compute factorization.
  193: *
  194:       DO 20 I = 1, MN
  195: *
  196:          OFFPI = OFFSET + I
  197: *
  198: *        Determine ith pivot column and swap if necessary.
  199: *
  200:          PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 )
  201: *
  202:          IF( PVT.NE.I ) THEN
  203:             CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
  204:             ITEMP = JPVT( PVT )
  205:             JPVT( PVT ) = JPVT( I )
  206:             JPVT( I ) = ITEMP
  207:             VN1( PVT ) = VN1( I )
  208:             VN2( PVT ) = VN2( I )
  209:          END IF
  210: *
  211: *        Generate elementary reflector H(i).
  212: *
  213:          IF( OFFPI.LT.M ) THEN
  214:             CALL DLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
  215:      $                   TAU( I ) )
  216:          ELSE
  217:             CALL DLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
  218:          END IF
  219: *
  220:          IF( I.LE.N ) THEN
  221: *
  222: *           Apply H(i)**T to A(offset+i:m,i+1:n) from the left.
  223: *
  224:             AII = A( OFFPI, I )
  225:             A( OFFPI, I ) = ONE
  226:             CALL DLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
  227:      $                  TAU( I ), A( OFFPI, I+1 ), LDA, WORK( 1 ) )
  228:             A( OFFPI, I ) = AII
  229:          END IF
  230: *
  231: *        Update partial column norms.
  232: *
  233:          DO 10 J = I + 1, N
  234:             IF( VN1( J ).NE.ZERO ) THEN
  235: *
  236: *              NOTE: The following 4 lines follow from the analysis in
  237: *              Lapack Working Note 176.
  238: *
  239:                TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
  240:                TEMP = MAX( TEMP, ZERO )
  241:                TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
  242:                IF( TEMP2 .LE. TOL3Z ) THEN
  243:                   IF( OFFPI.LT.M ) THEN
  244:                      VN1( J ) = DNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
  245:                      VN2( J ) = VN1( J )
  246:                   ELSE
  247:                      VN1( J ) = ZERO
  248:                      VN2( J ) = ZERO
  249:                   END IF
  250:                ELSE
  251:                   VN1( J ) = VN1( J )*SQRT( TEMP )
  252:                END IF
  253:             END IF
  254:    10    CONTINUE
  255: *
  256:    20 CONTINUE
  257: *
  258:       RETURN
  259: *
  260: *     End of DLAQP2
  261: *
  262:       END