File:  [local] / rpl / lapack / lapack / zgeqrf.f
Revision 1.19: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:19 2023 UTC (9 months 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 ZGEQRF
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download ZGEQRF + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeqrf.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeqrf.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeqrf.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZGEQRF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            INFO, LDA, LWORK, M, N
   25: *       ..
   26: *       .. Array Arguments ..
   27: *       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
   28: *       ..
   29: *
   30: *
   31: *> \par Purpose:
   32: *  =============
   33: *>
   34: *> \verbatim
   35: *>
   36: *> ZGEQRF computes a QR factorization of a complex M-by-N matrix A:
   37: *>
   38: *>    A = Q * ( R ),
   39: *>            ( 0 )
   40: *>
   41: *> where:
   42: *>
   43: *>    Q is a M-by-M orthogonal matrix;
   44: *>    R is an upper-triangular N-by-N matrix;
   45: *>    0 is a (M-N)-by-N zero matrix, if M > N.
   46: *>
   47: *> \endverbatim
   48: *
   49: *  Arguments:
   50: *  ==========
   51: *
   52: *> \param[in] M
   53: *> \verbatim
   54: *>          M is INTEGER
   55: *>          The number of rows of the matrix A.  M >= 0.
   56: *> \endverbatim
   57: *>
   58: *> \param[in] N
   59: *> \verbatim
   60: *>          N is INTEGER
   61: *>          The number of columns of the matrix A.  N >= 0.
   62: *> \endverbatim
   63: *>
   64: *> \param[in,out] A
   65: *> \verbatim
   66: *>          A is COMPLEX*16 array, dimension (LDA,N)
   67: *>          On entry, the M-by-N matrix A.
   68: *>          On exit, the elements on and above the diagonal of the array
   69: *>          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
   70: *>          upper triangular if m >= n); the elements below the diagonal,
   71: *>          with the array TAU, represent the unitary matrix Q as a
   72: *>          product of min(m,n) elementary reflectors (see Further
   73: *>          Details).
   74: *> \endverbatim
   75: *>
   76: *> \param[in] LDA
   77: *> \verbatim
   78: *>          LDA is INTEGER
   79: *>          The leading dimension of the array A.  LDA >= max(1,M).
   80: *> \endverbatim
   81: *>
   82: *> \param[out] TAU
   83: *> \verbatim
   84: *>          TAU is COMPLEX*16 array, dimension (min(M,N))
   85: *>          The scalar factors of the elementary reflectors (see Further
   86: *>          Details).
   87: *> \endverbatim
   88: *>
   89: *> \param[out] WORK
   90: *> \verbatim
   91: *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
   92: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
   93: *> \endverbatim
   94: *>
   95: *> \param[in] LWORK
   96: *> \verbatim
   97: *>          LWORK is INTEGER
   98: *>          The dimension of the array WORK.
   99: *>          LWORK >= 1, if MIN(M,N) = 0, and LWORK >= N, otherwise.
  100: *>          For optimum performance LWORK >= N*NB, where NB is
  101: *>          the optimal blocksize.
  102: *>
  103: *>          If LWORK = -1, then a workspace query is assumed; the routine
  104: *>          only calculates the optimal size of the WORK array, returns
  105: *>          this value as the first entry of the WORK array, and no error
  106: *>          message related to LWORK is issued by XERBLA.
  107: *> \endverbatim
  108: *>
  109: *> \param[out] INFO
  110: *> \verbatim
  111: *>          INFO is INTEGER
  112: *>          = 0:  successful exit
  113: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  114: *> \endverbatim
  115: *
  116: *  Authors:
  117: *  ========
  118: *
  119: *> \author Univ. of Tennessee
  120: *> \author Univ. of California Berkeley
  121: *> \author Univ. of Colorado Denver
  122: *> \author NAG Ltd.
  123: *
  124: *> \ingroup complex16GEcomputational
  125: *
  126: *> \par Further Details:
  127: *  =====================
  128: *>
  129: *> \verbatim
  130: *>
  131: *>  The matrix Q is represented as a product of elementary reflectors
  132: *>
  133: *>     Q = H(1) H(2) . . . H(k), where k = min(m,n).
  134: *>
  135: *>  Each H(i) has the form
  136: *>
  137: *>     H(i) = I - tau * v * v**H
  138: *>
  139: *>  where tau is a complex scalar, and v is a complex vector with
  140: *>  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
  141: *>  and tau in TAU(i).
  142: *> \endverbatim
  143: *>
  144: *  =====================================================================
  145:       SUBROUTINE ZGEQRF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
  146: *
  147: *  -- LAPACK computational routine --
  148: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  149: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  150: *
  151: *     .. Scalar Arguments ..
  152:       INTEGER            INFO, LDA, LWORK, M, N
  153: *     ..
  154: *     .. Array Arguments ..
  155:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
  156: *     ..
  157: *
  158: *  =====================================================================
  159: *
  160: *     .. Local Scalars ..
  161:       LOGICAL            LQUERY
  162:       INTEGER            I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB,
  163:      $                   NBMIN, NX
  164: *     ..
  165: *     .. External Subroutines ..
  166:       EXTERNAL           XERBLA, ZGEQR2, ZLARFB, ZLARFT
  167: *     ..
  168: *     .. Intrinsic Functions ..
  169:       INTRINSIC          MAX, MIN
  170: *     ..
  171: *     .. External Functions ..
  172:       INTEGER            ILAENV
  173:       EXTERNAL           ILAENV
  174: *     ..
  175: *     .. Executable Statements ..
  176: *
  177: *     Test the input arguments
  178: *
  179:       K = MIN( M, N )
  180:       INFO = 0
  181:       NB = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
  182:       LQUERY = ( LWORK.EQ.-1 )
  183:       IF( M.LT.0 ) THEN
  184:          INFO = -1
  185:       ELSE IF( N.LT.0 ) THEN
  186:          INFO = -2
  187:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  188:          INFO = -4
  189:       ELSE IF( .NOT.LQUERY ) THEN
  190:          IF( LWORK.LE.0 .OR. ( M.GT.0 .AND. LWORK.LT.MAX( 1, N ) ) )
  191:      $      INFO = -7
  192:       END IF
  193:       IF( INFO.NE.0 ) THEN
  194:          CALL XERBLA( 'ZGEQRF', -INFO )
  195:          RETURN
  196:       ELSE IF( LQUERY ) THEN
  197:          IF( K.EQ.0 ) THEN
  198:             LWKOPT = 1
  199:          ELSE
  200:             LWKOPT = N*NB
  201:          END IF
  202:          WORK( 1 ) = LWKOPT
  203:          RETURN
  204:       END IF
  205: *
  206: *     Quick return if possible
  207: *
  208:       IF( K.EQ.0 ) THEN
  209:          WORK( 1 ) = 1
  210:          RETURN
  211:       END IF
  212: *
  213:       NBMIN = 2
  214:       NX = 0
  215:       IWS = N
  216:       IF( NB.GT.1 .AND. NB.LT.K ) THEN
  217: *
  218: *        Determine when to cross over from blocked to unblocked code.
  219: *
  220:          NX = MAX( 0, ILAENV( 3, 'ZGEQRF', ' ', M, N, -1, -1 ) )
  221:          IF( NX.LT.K ) THEN
  222: *
  223: *           Determine if workspace is large enough for blocked code.
  224: *
  225:             LDWORK = N
  226:             IWS = LDWORK*NB
  227:             IF( LWORK.LT.IWS ) THEN
  228: *
  229: *              Not enough workspace to use optimal NB:  reduce NB and
  230: *              determine the minimum value of NB.
  231: *
  232:                NB = LWORK / LDWORK
  233:                NBMIN = MAX( 2, ILAENV( 2, 'ZGEQRF', ' ', M, N, -1,
  234:      $                 -1 ) )
  235:             END IF
  236:          END IF
  237:       END IF
  238: *
  239:       IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
  240: *
  241: *        Use blocked code initially
  242: *
  243:          DO 10 I = 1, K - NX, NB
  244:             IB = MIN( K-I+1, NB )
  245: *
  246: *           Compute the QR factorization of the current block
  247: *           A(i:m,i:i+ib-1)
  248: *
  249:             CALL ZGEQR2( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK,
  250:      $                   IINFO )
  251:             IF( I+IB.LE.N ) THEN
  252: *
  253: *              Form the triangular factor of the block reflector
  254: *              H = H(i) H(i+1) . . . H(i+ib-1)
  255: *
  256:                CALL ZLARFT( 'Forward', 'Columnwise', M-I+1, IB,
  257:      $                      A( I, I ), LDA, TAU( I ), WORK, LDWORK )
  258: *
  259: *              Apply H**H to A(i:m,i+ib:n) from the left
  260: *
  261:                CALL ZLARFB( 'Left', 'Conjugate transpose', 'Forward',
  262:      $                      'Columnwise', M-I+1, N-I-IB+1, IB,
  263:      $                      A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ),
  264:      $                      LDA, WORK( IB+1 ), LDWORK )
  265:             END IF
  266:    10    CONTINUE
  267:       ELSE
  268:          I = 1
  269:       END IF
  270: *
  271: *     Use unblocked code to factor the last or only block.
  272: *
  273:       IF( I.LE.K )
  274:      $   CALL ZGEQR2( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK,
  275:      $                IINFO )
  276: *
  277:       WORK( 1 ) = IWS
  278:       RETURN
  279: *
  280: *     End of ZGEQRF
  281: *
  282:       END

CVSweb interface <joel.bertrand@systella.fr>