File:  [local] / rpl / lapack / lapack / dlaorhr_col_getrfnp.f
Revision 1.2: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:38:55 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 DLAORHR_COL_GETRFNP
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DLAORHR_COL_GETRFNP + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaorhr_col_getrfnp.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaorhr_col_getrfnp.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaorhr_col_getrfnp.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLAORHR_COL_GETRFNP( M, N, A, LDA, D, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            INFO, LDA, M, N
   25: *       ..
   26: *       .. Array Arguments ..
   27: *       DOUBLE PRECISION   A( LDA, * ), D( * )
   28: *       ..
   29: *
   30: *
   31: *> \par Purpose:
   32: *  =============
   33: *>
   34: *> \verbatim
   35: *>
   36: *> DLAORHR_COL_GETRFNP computes the modified LU factorization without
   37: *> pivoting of a real general M-by-N matrix A. The factorization has
   38: *> the form:
   39: *>
   40: *>     A - S = L * U,
   41: *>
   42: *> where:
   43: *>    S is a m-by-n diagonal sign matrix with the diagonal D, so that
   44: *>    D(i) = S(i,i), 1 <= i <= min(M,N). The diagonal D is constructed
   45: *>    as D(i)=-SIGN(A(i,i)), where A(i,i) is the value after performing
   46: *>    i-1 steps of Gaussian elimination. This means that the diagonal
   47: *>    element at each step of "modified" Gaussian elimination is
   48: *>    at least one in absolute value (so that division-by-zero not
   49: *>    not possible during the division by the diagonal element);
   50: *>
   51: *>    L is a M-by-N lower triangular matrix with unit diagonal elements
   52: *>    (lower trapezoidal if M > N);
   53: *>
   54: *>    and U is a M-by-N upper triangular matrix
   55: *>    (upper trapezoidal if M < N).
   56: *>
   57: *> This routine is an auxiliary routine used in the Householder
   58: *> reconstruction routine DORHR_COL. In DORHR_COL, this routine is
   59: *> applied to an M-by-N matrix A with orthonormal columns, where each
   60: *> element is bounded by one in absolute value. With the choice of
   61: *> the matrix S above, one can show that the diagonal element at each
   62: *> step of Gaussian elimination is the largest (in absolute value) in
   63: *> the column on or below the diagonal, so that no pivoting is required
   64: *> for numerical stability [1].
   65: *>
   66: *> For more details on the Householder reconstruction algorithm,
   67: *> including the modified LU factorization, see [1].
   68: *>
   69: *> This is the blocked right-looking version of the algorithm,
   70: *> calling Level 3 BLAS to update the submatrix. To factorize a block,
   71: *> this routine calls the recursive routine DLAORHR_COL_GETRFNP2.
   72: *>
   73: *> [1] "Reconstructing Householder vectors from tall-skinny QR",
   74: *>     G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
   75: *>     E. Solomonik, J. Parallel Distrib. Comput.,
   76: *>     vol. 85, pp. 3-31, 2015.
   77: *> \endverbatim
   78: *
   79: *  Arguments:
   80: *  ==========
   81: *
   82: *> \param[in] M
   83: *> \verbatim
   84: *>          M is INTEGER
   85: *>          The number of rows of the matrix A.  M >= 0.
   86: *> \endverbatim
   87: *>
   88: *> \param[in] N
   89: *> \verbatim
   90: *>          N is INTEGER
   91: *>          The number of columns of the matrix A.  N >= 0.
   92: *> \endverbatim
   93: *>
   94: *> \param[in,out] A
   95: *> \verbatim
   96: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   97: *>          On entry, the M-by-N matrix to be factored.
   98: *>          On exit, the factors L and U from the factorization
   99: *>          A-S=L*U; the unit diagonal elements of L are not stored.
  100: *> \endverbatim
  101: *>
  102: *> \param[in] LDA
  103: *> \verbatim
  104: *>          LDA is INTEGER
  105: *>          The leading dimension of the array A.  LDA >= max(1,M).
  106: *> \endverbatim
  107: *>
  108: *> \param[out] D
  109: *> \verbatim
  110: *>          D is DOUBLE PRECISION array, dimension min(M,N)
  111: *>          The diagonal elements of the diagonal M-by-N sign matrix S,
  112: *>          D(i) = S(i,i), where 1 <= i <= min(M,N). The elements can
  113: *>          be only plus or minus one.
  114: *> \endverbatim
  115: *>
  116: *> \param[out] INFO
  117: *> \verbatim
  118: *>          INFO is INTEGER
  119: *>          = 0:  successful exit
  120: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  121: *> \endverbatim
  122: *>
  123: *  Authors:
  124: *  ========
  125: *
  126: *> \author Univ. of Tennessee
  127: *> \author Univ. of California Berkeley
  128: *> \author Univ. of Colorado Denver
  129: *> \author NAG Ltd.
  130: *
  131: *> \ingroup doubleGEcomputational
  132: *
  133: *> \par Contributors:
  134: *  ==================
  135: *>
  136: *> \verbatim
  137: *>
  138: *> November 2019, Igor Kozachenko,
  139: *>                Computer Science Division,
  140: *>                University of California, Berkeley
  141: *>
  142: *> \endverbatim
  143: *
  144: *  =====================================================================
  145:       SUBROUTINE DLAORHR_COL_GETRFNP( M, N, A, LDA, D, INFO )
  146:       IMPLICIT NONE
  147: *
  148: *  -- LAPACK computational routine --
  149: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  150: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  151: *
  152: *     .. Scalar Arguments ..
  153:       INTEGER            INFO, LDA, M, N
  154: *     ..
  155: *     .. Array Arguments ..
  156:       DOUBLE PRECISION   A( LDA, * ), D( * )
  157: *     ..
  158: *
  159: *  =====================================================================
  160: *
  161: *     .. Parameters ..
  162:       DOUBLE PRECISION   ONE
  163:       PARAMETER          ( ONE = 1.0D+0 )
  164: *     ..
  165: *     .. Local Scalars ..
  166:       INTEGER            IINFO, J, JB, NB
  167: *     ..
  168: *     .. External Subroutines ..
  169:       EXTERNAL           DGEMM, DLAORHR_COL_GETRFNP2, DTRSM, XERBLA
  170: *     ..
  171: *     .. External Functions ..
  172:       INTEGER            ILAENV
  173:       EXTERNAL           ILAENV
  174: *     ..
  175: *     .. Intrinsic Functions ..
  176:       INTRINSIC          MAX, MIN
  177: *     ..
  178: *     .. Executable Statements ..
  179: *
  180: *     Test the input parameters.
  181: *
  182:       INFO = 0
  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:       END IF
  190:       IF( INFO.NE.0 ) THEN
  191:          CALL XERBLA( 'DLAORHR_COL_GETRFNP', -INFO )
  192:          RETURN
  193:       END IF
  194: *
  195: *     Quick return if possible
  196: *
  197:       IF( MIN( M, N ).EQ.0 )
  198:      $   RETURN
  199: *
  200: *     Determine the block size for this environment.
  201: *
  202: 
  203:       NB = ILAENV( 1, 'DLAORHR_COL_GETRFNP', ' ', M, N, -1, -1 )
  204: 
  205:       IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
  206: *
  207: *        Use unblocked code.
  208: *
  209:          CALL DLAORHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
  210:       ELSE
  211: *
  212: *        Use blocked code.
  213: *
  214:          DO J = 1, MIN( M, N ), NB
  215:             JB = MIN( MIN( M, N )-J+1, NB )
  216: *
  217: *           Factor diagonal and subdiagonal blocks.
  218: *
  219:             CALL DLAORHR_COL_GETRFNP2( M-J+1, JB, A( J, J ), LDA,
  220:      $                                 D( J ), IINFO )
  221: *
  222:             IF( J+JB.LE.N ) THEN
  223: *
  224: *              Compute block row of U.
  225: *
  226:                CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
  227:      $                     N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
  228:      $                     LDA )
  229:                IF( J+JB.LE.M ) THEN
  230: *
  231: *                 Update trailing submatrix.
  232: *
  233:                   CALL DGEMM( 'No transpose', 'No transpose', M-J-JB+1,
  234:      $                        N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA,
  235:      $                        A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ),
  236:      $                        LDA )
  237:                END IF
  238:             END IF
  239:          END DO
  240:       END IF
  241:       RETURN
  242: *
  243: *     End of DLAORHR_COL_GETRFNP
  244: *
  245:       END

CVSweb interface <joel.bertrand@systella.fr>