File:  [local] / rpl / lapack / lapack / zlaqp2.f
Revision 1.20: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:30 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    1: *> \brief \b ZLAQP2 computes a QR factorization with column pivoting of the matrix block.
    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: *> \ingroup complex16OTHERauxiliary
  126: *
  127: *> \par Contributors:
  128: *  ==================
  129: *>
  130: *>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
  131: *>    X. Sun, Computer Science Dept., Duke University, USA
  132: *> \n
  133: *>  Partial column norm updating strategy modified on April 2011
  134: *>    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
  135: *>    University of Zagreb, Croatia.
  136: *
  137: *> \par References:
  138: *  ================
  139: *>
  140: *> LAPACK Working Note 176
  141: *
  142: *> \htmlonly
  143: *> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a>
  144: *> \endhtmlonly
  145: *
  146: *  =====================================================================
  147:       SUBROUTINE ZLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
  148:      $                   WORK )
  149: *
  150: *  -- LAPACK auxiliary routine --
  151: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  152: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  153: *
  154: *     .. Scalar Arguments ..
  155:       INTEGER            LDA, M, N, OFFSET
  156: *     ..
  157: *     .. Array Arguments ..
  158:       INTEGER            JPVT( * )
  159:       DOUBLE PRECISION   VN1( * ), VN2( * )
  160:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
  161: *     ..
  162: *
  163: *  =====================================================================
  164: *
  165: *     .. Parameters ..
  166:       DOUBLE PRECISION   ZERO, ONE
  167:       COMPLEX*16         CONE
  168:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
  169:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
  170: *     ..
  171: *     .. Local Scalars ..
  172:       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
  173:       DOUBLE PRECISION   TEMP, TEMP2, TOL3Z
  174:       COMPLEX*16         AII
  175: *     ..
  176: *     .. External Subroutines ..
  177:       EXTERNAL           ZLARF, ZLARFG, ZSWAP
  178: *     ..
  179: *     .. Intrinsic Functions ..
  180:       INTRINSIC          ABS, DCONJG, MAX, MIN, SQRT
  181: *     ..
  182: *     .. External Functions ..
  183:       INTEGER            IDAMAX
  184:       DOUBLE PRECISION   DLAMCH, DZNRM2
  185:       EXTERNAL           IDAMAX, DLAMCH, DZNRM2
  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 ZSWAP( 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 ZLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
  215:      $                   TAU( I ) )
  216:          ELSE
  217:             CALL ZLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
  218:          END IF
  219: *
  220:          IF( I.LT.N ) THEN
  221: *
  222: *           Apply H(i)**H to A(offset+i:m,i+1:n) from the left.
  223: *
  224:             AII = A( OFFPI, I )
  225:             A( OFFPI, I ) = CONE
  226:             CALL ZLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
  227:      $                  DCONJG( TAU( I ) ), A( OFFPI, I+1 ), LDA,
  228:      $                  WORK( 1 ) )
  229:             A( OFFPI, I ) = AII
  230:          END IF
  231: *
  232: *        Update partial column norms.
  233: *
  234:          DO 10 J = I + 1, N
  235:             IF( VN1( J ).NE.ZERO ) THEN
  236: *
  237: *              NOTE: The following 4 lines follow from the analysis in
  238: *              Lapack Working Note 176.
  239: *
  240:                TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
  241:                TEMP = MAX( TEMP, ZERO )
  242:                TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
  243:                IF( TEMP2 .LE. TOL3Z ) THEN
  244:                   IF( OFFPI.LT.M ) THEN
  245:                      VN1( J ) = DZNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
  246:                      VN2( J ) = VN1( J )
  247:                   ELSE
  248:                      VN1( J ) = ZERO
  249:                      VN2( J ) = ZERO
  250:                   END IF
  251:                ELSE
  252:                   VN1( J ) = VN1( J )*SQRT( TEMP )
  253:                END IF
  254:             END IF
  255:    10    CONTINUE
  256: *
  257:    20 CONTINUE
  258: *
  259:       RETURN
  260: *
  261: *     End of ZLAQP2
  262: *
  263:       END

CVSweb interface <joel.bertrand@systella.fr>