File:  [local] / rpl / lapack / lapack / zlaqp2.f
Revision 1.11: download - view: text, annotated - select for diffs - revision graph
Mon Nov 21 22:19:52 2011 UTC (12 years, 5 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_8, rpl-4_1_7, rpl-4_1_6, rpl-4_1_5, rpl-4_1_4, HEAD
Cohérence

    1: *> \brief \b ZLAQP2
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZLAQP2 + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqp2.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqp2.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqp2.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZLAQP2( 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   VN1( * ), VN2( * )
   30: *       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
   31: *       ..
   32: *  
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> ZLAQP2 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 complex16OTHERauxiliary
  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 ZLAQP2( 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   VN1( * ), VN2( * )
  163:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
  164: *     ..
  165: *
  166: *  =====================================================================
  167: *
  168: *     .. Parameters ..
  169:       DOUBLE PRECISION   ZERO, ONE
  170:       COMPLEX*16         CONE
  171:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
  172:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
  173: *     ..
  174: *     .. Local Scalars ..
  175:       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
  176:       DOUBLE PRECISION   TEMP, TEMP2, TOL3Z
  177:       COMPLEX*16         AII
  178: *     ..
  179: *     .. External Subroutines ..
  180:       EXTERNAL           ZLARF, ZLARFG, ZSWAP
  181: *     ..
  182: *     .. Intrinsic Functions ..
  183:       INTRINSIC          ABS, DCONJG, MAX, MIN, SQRT
  184: *     ..
  185: *     .. External Functions ..
  186:       INTEGER            IDAMAX
  187:       DOUBLE PRECISION   DLAMCH, DZNRM2
  188:       EXTERNAL           IDAMAX, DLAMCH, DZNRM2
  189: *     ..
  190: *     .. Executable Statements ..
  191: *
  192:       MN = MIN( M-OFFSET, N )
  193:       TOL3Z = SQRT(DLAMCH('Epsilon'))
  194: *
  195: *     Compute factorization.
  196: *
  197:       DO 20 I = 1, MN
  198: *
  199:          OFFPI = OFFSET + I
  200: *
  201: *        Determine ith pivot column and swap if necessary.
  202: *
  203:          PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 )
  204: *
  205:          IF( PVT.NE.I ) THEN
  206:             CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
  207:             ITEMP = JPVT( PVT )
  208:             JPVT( PVT ) = JPVT( I )
  209:             JPVT( I ) = ITEMP
  210:             VN1( PVT ) = VN1( I )
  211:             VN2( PVT ) = VN2( I )
  212:          END IF
  213: *
  214: *        Generate elementary reflector H(i).
  215: *
  216:          IF( OFFPI.LT.M ) THEN
  217:             CALL ZLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
  218:      $                   TAU( I ) )
  219:          ELSE
  220:             CALL ZLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
  221:          END IF
  222: *
  223:          IF( I.LT.N ) THEN
  224: *
  225: *           Apply H(i)**H to A(offset+i:m,i+1:n) from the left.
  226: *
  227:             AII = A( OFFPI, I )
  228:             A( OFFPI, I ) = CONE
  229:             CALL ZLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
  230:      $                  DCONJG( TAU( I ) ), A( OFFPI, I+1 ), LDA,
  231:      $                  WORK( 1 ) )
  232:             A( OFFPI, I ) = AII
  233:          END IF
  234: *
  235: *        Update partial column norms.
  236: *
  237:          DO 10 J = I + 1, N
  238:             IF( VN1( J ).NE.ZERO ) THEN
  239: *
  240: *              NOTE: The following 4 lines follow from the analysis in
  241: *              Lapack Working Note 176.
  242: *
  243:                TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
  244:                TEMP = MAX( TEMP, ZERO )
  245:                TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
  246:                IF( TEMP2 .LE. TOL3Z ) THEN
  247:                   IF( OFFPI.LT.M ) THEN
  248:                      VN1( J ) = DZNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
  249:                      VN2( J ) = VN1( J )
  250:                   ELSE
  251:                      VN1( J ) = ZERO
  252:                      VN2( J ) = ZERO
  253:                   END IF
  254:                ELSE
  255:                   VN1( J ) = VN1( J )*SQRT( TEMP )
  256:                END IF
  257:             END IF
  258:    10    CONTINUE
  259: *
  260:    20 CONTINUE
  261: *
  262:       RETURN
  263: *
  264: *     End of ZLAQP2
  265: *
  266:       END

CVSweb interface <joel.bertrand@systella.fr>