Return to dorhr_col.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack.

    1: *> \brief \b DORHR_COL
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DORHR_COL + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorhr_col.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorhr_col.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorhr_col.f">
   15: *> [TXT]</a>
   16: *>
   17: *  Definition:
   18: *  ===========
   19: *
   20: *       SUBROUTINE DORHR_COL( M, N, NB, A, LDA, T, LDT, D, INFO )
   21: *
   22: *       .. Scalar Arguments ..
   23: *       INTEGER           INFO, LDA, LDT, M, N, NB
   24: *       ..
   25: *       .. Array Arguments ..
   26: *       DOUBLE PRECISION  A( LDA, * ), D( * ), T( LDT, * )
   27: *       ..
   28: *
   29: *> \par Purpose:
   30: *  =============
   31: *>
   32: *> \verbatim
   33: *>
   34: *>  DORHR_COL takes an M-by-N real matrix Q_in with orthonormal columns
   35: *>  as input, stored in A, and performs Householder Reconstruction (HR),
   36: *>  i.e. reconstructs Householder vectors V(i) implicitly representing
   37: *>  another M-by-N matrix Q_out, with the property that Q_in = Q_out*S,
   38: *>  where S is an N-by-N diagonal matrix with diagonal entries
   39: *>  equal to +1 or -1. The Householder vectors (columns V(i) of V) are
   40: *>  stored in A on output, and the diagonal entries of S are stored in D.
   41: *>  Block reflectors are also returned in T
   42: *>  (same output format as DGEQRT).
   43: *> \endverbatim
   44: *
   45: *  Arguments:
   46: *  ==========
   47: *
   48: *> \param[in] M
   49: *> \verbatim
   50: *>          M is INTEGER
   51: *>          The number of rows of the matrix A. M >= 0.
   52: *> \endverbatim
   53: *>
   54: *> \param[in] N
   55: *> \verbatim
   56: *>          N is INTEGER
   57: *>          The number of columns of the matrix A. M >= N >= 0.
   58: *> \endverbatim
   59: *>
   60: *> \param[in] NB
   61: *> \verbatim
   62: *>          NB is INTEGER
   63: *>          The column block size to be used in the reconstruction
   64: *>          of Householder column vector blocks in the array A and
   65: *>          corresponding block reflectors in the array T. NB >= 1.
   66: *>          (Note that if NB > N, then N is used instead of NB
   67: *>          as the column block size.)
   68: *> \endverbatim
   69: *>
   70: *> \param[in,out] A
   71: *> \verbatim
   72: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   73: *>
   74: *>          On entry:
   75: *>
   76: *>             The array A contains an M-by-N orthonormal matrix Q_in,
   77: *>             i.e the columns of A are orthogonal unit vectors.
   78: *>
   79: *>          On exit:
   80: *>
   81: *>             The elements below the diagonal of A represent the unit
   82: *>             lower-trapezoidal matrix V of Householder column vectors
   83: *>             V(i). The unit diagonal entries of V are not stored
   84: *>             (same format as the output below the diagonal in A from
   85: *>             DGEQRT). The matrix T and the matrix V stored on output
   86: *>             in A implicitly define Q_out.
   87: *>
   88: *>             The elements above the diagonal contain the factor U
   89: *>             of the "modified" LU-decomposition:
   90: *>                Q_in - ( S ) = V * U
   91: *>                       ( 0 )
   92: *>             where 0 is a (M-N)-by-(M-N) zero matrix.
   93: *> \endverbatim
   94: *>
   95: *> \param[in] LDA
   96: *> \verbatim
   97: *>          LDA is INTEGER
   98: *>          The leading dimension of the array A.  LDA >= max(1,M).
   99: *> \endverbatim
  100: *>
  101: *> \param[out] T
  102: *> \verbatim
  103: *>          T is DOUBLE PRECISION array,
  104: *>          dimension (LDT, N)
  105: *>
  106: *>          Let NOCB = Number_of_output_col_blocks
  107: *>                   = CEIL(N/NB)
  108: *>
  109: *>          On exit, T(1:NB, 1:N) contains NOCB upper-triangular
  110: *>          block reflectors used to define Q_out stored in compact
  111: *>          form as a sequence of upper-triangular NB-by-NB column
  112: *>          blocks (same format as the output T in DGEQRT).
  113: *>          The matrix T and the matrix V stored on output in A
  114: *>          implicitly define Q_out. NOTE: The lower triangles
  115: *>          below the upper-triangular blcoks will be filled with
  116: *>          zeros. See Further Details.
  117: *> \endverbatim
  118: *>
  119: *> \param[in] LDT
  120: *> \verbatim
  121: *>          LDT is INTEGER
  122: *>          The leading dimension of the array T.
  123: *>          LDT >= max(1,min(NB,N)).
  124: *> \endverbatim
  125: *>
  126: *> \param[out] D
  127: *> \verbatim
  128: *>          D is DOUBLE PRECISION array, dimension min(M,N).
  129: *>          The elements can be only plus or minus one.
  130: *>
  131: *>          D(i) is constructed as D(i) = -SIGN(Q_in_i(i,i)), where
  132: *>          1 <= i <= min(M,N), and Q_in_i is Q_in after performing
  133: *>          i-1 steps of “modified” Gaussian elimination.
  134: *>          See Further Details.
  135: *> \endverbatim
  136: *>
  137: *> \param[out] INFO
  138: *> \verbatim
  139: *>          INFO is INTEGER
  140: *>          = 0:  successful exit
  141: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  142: *> \endverbatim
  143: *>
  144: *> \par Further Details:
  145: *  =====================
  146: *>
  147: *> \verbatim
  148: *>
  149: *> The computed M-by-M orthogonal factor Q_out is defined implicitly as
  150: *> a product of orthogonal matrices Q_out(i). Each Q_out(i) is stored in
  151: *> the compact WY-representation format in the corresponding blocks of
  152: *> matrices V (stored in A) and T.
  153: *>
  154: *> The M-by-N unit lower-trapezoidal matrix V stored in the M-by-N
  155: *> matrix A contains the column vectors V(i) in NB-size column
  156: *> blocks VB(j). For example, VB(1) contains the columns
  157: *> V(1), V(2), ... V(NB). NOTE: The unit entries on
  158: *> the diagonal of Y are not stored in A.
  159: *>
  160: *> The number of column blocks is
  161: *>
  162: *>     NOCB = Number_of_output_col_blocks = CEIL(N/NB)
  163: *>
  164: *> where each block is of order NB except for the last block, which
  165: *> is of order LAST_NB = N - (NOCB-1)*NB.
  166: *>
  167: *> For example, if M=6,  N=5 and NB=2, the matrix V is
  168: *>
  169: *>
  170: *>     V = (    VB(1),   VB(2), VB(3) ) =
  171: *>
  172: *>       = (   1                      )
  173: *>         ( v21    1                 )
  174: *>         ( v31  v32    1            )
  175: *>         ( v41  v42  v43   1        )
  176: *>         ( v51  v52  v53  v54    1  )
  177: *>         ( v61  v62  v63  v54   v65 )
  178: *>
  179: *>
  180: *> For each of the column blocks VB(i), an upper-triangular block
  181: *> reflector TB(i) is computed. These blocks are stored as
  182: *> a sequence of upper-triangular column blocks in the NB-by-N
  183: *> matrix T. The size of each TB(i) block is NB-by-NB, except
  184: *> for the last block, whose size is LAST_NB-by-LAST_NB.
  185: *>
  186: *> For example, if M=6,  N=5 and NB=2, the matrix T is
  187: *>
  188: *>     T  = (    TB(1),    TB(2), TB(3) ) =
  189: *>
  190: *>        = ( t11  t12  t13  t14   t15  )
  191: *>          (      t22       t24        )
  192: *>
  193: *>
  194: *> The M-by-M factor Q_out is given as a product of NOCB
  195: *> orthogonal M-by-M matrices Q_out(i).
  196: *>
  197: *>     Q_out = Q_out(1) * Q_out(2) * ... * Q_out(NOCB),
  198: *>
  199: *> where each matrix Q_out(i) is given by the WY-representation
  200: *> using corresponding blocks from the matrices V and T:
  201: *>
  202: *>     Q_out(i) = I - VB(i) * TB(i) * (VB(i))**T,
  203: *>
  204: *> where I is the identity matrix. Here is the formula with matrix
  205: *> dimensions:
  206: *>
  207: *>  Q(i){M-by-M} = I{M-by-M} -
  208: *>    VB(i){M-by-INB} * TB(i){INB-by-INB} * (VB(i))**T {INB-by-M},
  209: *>
  210: *> where INB = NB, except for the last block NOCB
  211: *> for which INB=LAST_NB.
  212: *>
  213: *> =====
  214: *> NOTE:
  215: *> =====
  216: *>
  217: *> If Q_in is the result of doing a QR factorization
  218: *> B = Q_in * R_in, then:
  219: *>
  220: *> B = (Q_out*S) * R_in = Q_out * (S * R_in) = O_out * R_out.
  221: *>
  222: *> So if one wants to interpret Q_out as the result
  223: *> of the QR factorization of B, then corresponding R_out
  224: *> should be obtained by R_out = S * R_in, i.e. some rows of R_in
  225: *> should be multiplied by -1.
  226: *>
  227: *> For the details of the algorithm, see [1].
  228: *>
  229: *> [1] "Reconstructing Householder vectors from tall-skinny QR",
  230: *>     G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
  231: *>     E. Solomonik, J. Parallel Distrib. Comput.,
  232: *>     vol. 85, pp. 3-31, 2015.
  233: *> \endverbatim
  234: *>
  235: *  Authors:
  236: *  ========
  237: *
  238: *> \author Univ. of Tennessee
  239: *> \author Univ. of California Berkeley
  240: *> \author Univ. of Colorado Denver
  241: *> \author NAG Ltd.
  242: *
  243: *> \date November 2019
  244: *
  245: *> \ingroup doubleOTHERcomputational
  246: *
  247: *> \par Contributors:
  248: *  ==================
  249: *>
  250: *> \verbatim
  251: *>
  252: *> November   2019, Igor Kozachenko,
  253: *>            Computer Science Division,
  254: *>            University of California, Berkeley
  255: *>
  256: *> \endverbatim
  257: *
  258: *  =====================================================================
  259:       SUBROUTINE DORHR_COL( M, N, NB, A, LDA, T, LDT, D, INFO )
  260:       IMPLICIT NONE
  261: *
  262: *  -- LAPACK computational routine (version 3.9.0) --
  263: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  264: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  265: *     November 2019
  266: *
  267: *     .. Scalar Arguments ..
  268:       INTEGER           INFO, LDA, LDT, M, N, NB
  269: *     ..
  270: *     .. Array Arguments ..
  271:       DOUBLE PRECISION  A( LDA, * ), D( * ), T( LDT, * )
  272: *     ..
  273: *
  274: *  =====================================================================
  275: *
  276: *     .. Parameters ..
  277:       DOUBLE PRECISION   ONE, ZERO
  278:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  279: *     ..
  280: *     .. Local Scalars ..
  281:       INTEGER            I, IINFO, J, JB, JBTEMP1, JBTEMP2, JNB,
  282:      $                   NPLUSONE
  283: *     ..
  284: *     .. External Subroutines ..
  285:       EXTERNAL           DCOPY, DLAORHR_COL_GETRFNP, DSCAL, DTRSM,
  286:      $                   XERBLA
  287: *     ..
  288: *     .. Intrinsic Functions ..
  289:       INTRINSIC          MAX, MIN
  290: *     ..
  291: *     .. Executable Statements ..
  292: *
  293: *     Test the input parameters
  294: *
  295:       INFO = 0
  296:       IF( M.LT.0 ) THEN
  297:          INFO = -1
  298:       ELSE IF( N.LT.0 .OR. N.GT.M ) THEN
  299:          INFO = -2
  300:       ELSE IF( NB.LT.1 ) THEN
  301:          INFO = -3
  302:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  303:          INFO = -5
  304:       ELSE IF( LDT.LT.MAX( 1, MIN( NB, N ) ) ) THEN
  305:          INFO = -7
  306:       END IF
  307: *
  308: *     Handle error in the input parameters.
  309: *
  310:       IF( INFO.NE.0 ) THEN
  311:          CALL XERBLA( 'DORHR_COL', -INFO )
  312:          RETURN
  313:       END IF
  314: *
  315: *     Quick return if possible
  316: *
  317:       IF( MIN( M, N ).EQ.0 ) THEN
  318:          RETURN
  319:       END IF
  320: *
  321: *     On input, the M-by-N matrix A contains the orthogonal
  322: *     M-by-N matrix Q_in.
  323: *
  324: *     (1) Compute the unit lower-trapezoidal V (ones on the diagonal
  325: *     are not stored) by performing the "modified" LU-decomposition.
  326: *
  327: *     Q_in - ( S ) = V * U = ( V1 ) * U,
  328: *            ( 0 )           ( V2 )
  329: *
  330: *     where 0 is an (M-N)-by-N zero matrix.
  331: *
  332: *     (1-1) Factor V1 and U.
  333: 
  334:       CALL DLAORHR_COL_GETRFNP( N, N, A, LDA, D, IINFO )
  335: *
  336: *     (1-2) Solve for V2.
  337: *
  338:       IF( M.GT.N ) THEN
  339:          CALL DTRSM( 'R', 'U', 'N', 'N', M-N, N, ONE, A, LDA,
  340:      $               A( N+1, 1 ), LDA )
  341:       END IF
  342: *
  343: *     (2) Reconstruct the block reflector T stored in T(1:NB, 1:N)
  344: *     as a sequence of upper-triangular blocks with NB-size column
  345: *     blocking.
  346: *
  347: *     Loop over the column blocks of size NB of the array A(1:M,1:N)
  348: *     and the array T(1:NB,1:N), JB is the column index of a column
  349: *     block, JNB is the column block size at each step JB.
  350: *
  351:       NPLUSONE = N + 1
  352:       DO JB = 1, N, NB
  353: *
  354: *        (2-0) Determine the column block size JNB.
  355: *
  356:          JNB = MIN( NPLUSONE-JB, NB )
  357: *
  358: *        (2-1) Copy the upper-triangular part of the current JNB-by-JNB
  359: *        diagonal block U(JB) (of the N-by-N matrix U) stored
  360: *        in A(JB:JB+JNB-1,JB:JB+JNB-1) into the upper-triangular part
  361: *        of the current JNB-by-JNB block T(1:JNB,JB:JB+JNB-1)
  362: *        column-by-column, total JNB*(JNB+1)/2 elements.
  363: *
  364:          JBTEMP1 = JB - 1
  365:          DO J = JB, JB+JNB-1
  366:             CALL DCOPY( J-JBTEMP1, A( JB, J ), 1, T( 1, J ), 1 )
  367:          END DO
  368: *
  369: *        (2-2) Perform on the upper-triangular part of the current
  370: *        JNB-by-JNB diagonal block U(JB) (of the N-by-N matrix U) stored
  371: *        in T(1:JNB,JB:JB+JNB-1) the following operation in place:
  372: *        (-1)*U(JB)*S(JB), i.e the result will be stored in the upper-
  373: *        triangular part of T(1:JNB,JB:JB+JNB-1). This multiplication
  374: *        of the JNB-by-JNB diagonal block U(JB) by the JNB-by-JNB
  375: *        diagonal block S(JB) of the N-by-N sign matrix S from the
  376: *        right means changing the sign of each J-th column of the block
  377: *        U(JB) according to the sign of the diagonal element of the block
  378: *        S(JB), i.e. S(J,J) that is stored in the array element D(J).
  379: *
  380:          DO J = JB, JB+JNB-1
  381:             IF( D( J ).EQ.ONE ) THEN
  382:                CALL DSCAL( J-JBTEMP1, -ONE, T( 1, J ), 1 )
  383:             END IF
  384:          END DO
  385: *
  386: *        (2-3) Perform the triangular solve for the current block
  387: *        matrix X(JB):
  388: *
  389: *               X(JB) * (A(JB)**T) = B(JB), where:
  390: *
  391: *               A(JB)**T  is a JNB-by-JNB unit upper-triangular
  392: *                         coefficient block, and A(JB)=V1(JB), which
  393: *                         is a JNB-by-JNB unit lower-triangular block
  394: *                         stored in A(JB:JB+JNB-1,JB:JB+JNB-1).
  395: *                         The N-by-N matrix V1 is the upper part
  396: *                         of the M-by-N lower-trapezoidal matrix V
  397: *                         stored in A(1:M,1:N);
  398: *
  399: *               B(JB)     is a JNB-by-JNB  upper-triangular right-hand
  400: *                         side block, B(JB) = (-1)*U(JB)*S(JB), and
  401: *                         B(JB) is stored in T(1:JNB,JB:JB+JNB-1);
  402: *
  403: *               X(JB)     is a JNB-by-JNB upper-triangular solution
  404: *                         block, X(JB) is the upper-triangular block
  405: *                         reflector T(JB), and X(JB) is stored
  406: *                         in T(1:JNB,JB:JB+JNB-1).
  407: *
  408: *             In other words, we perform the triangular solve for the
  409: *             upper-triangular block T(JB):
  410: *
  411: *               T(JB) * (V1(JB)**T) = (-1)*U(JB)*S(JB).
  412: *
  413: *             Even though the blocks X(JB) and B(JB) are upper-
  414: *             triangular, the routine DTRSM will access all JNB**2
  415: *             elements of the square T(1:JNB,JB:JB+JNB-1). Therefore,
  416: *             we need to set to zero the elements of the block
  417: *             T(1:JNB,JB:JB+JNB-1) below the diagonal before the call
  418: *             to DTRSM.
  419: *
  420: *        (2-3a) Set the elements to zero.
  421: *
  422:          JBTEMP2 = JB - 2
  423:          DO J = JB, JB+JNB-2
  424:             DO I = J-JBTEMP2, NB
  425:                T( I, J ) = ZERO
  426:             END DO
  427:          END DO
  428: *
  429: *        (2-3b) Perform the triangular solve.
  430: *
  431:          CALL DTRSM( 'R', 'L', 'T', 'U', JNB, JNB, ONE,
  432:      $               A( JB, JB ), LDA, T( 1, JB ), LDT )
  433: *
  434:       END DO
  435: *
  436:       RETURN
  437: *
  438: *     End of DORHR_COL
  439: *
  440:       END