Return to dlauum.f CVS log

Up to [local] / rpl / lapack / lapack

Première mise à jour de lapack et blas.

    1: *> \brief \b DLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm).
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DLAUUM + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlauum.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlauum.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlauum.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLAUUM( UPLO, N, A, LDA, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       CHARACTER          UPLO
   25: *       INTEGER            INFO, LDA, N
   26: *       ..
   27: *       .. Array Arguments ..
   28: *       DOUBLE PRECISION   A( LDA, * )
   29: *       ..
   30: *
   31: *
   32: *> \par Purpose:
   33: *  =============
   34: *>
   35: *> \verbatim
   36: *>
   37: *> DLAUUM computes the product U * U**T or L**T * L, where the triangular
   38: *> factor U or L is stored in the upper or lower triangular part of
   39: *> the array A.
   40: *>
   41: *> If UPLO = 'U' or 'u' then the upper triangle of the result is stored,
   42: *> overwriting the factor U in A.
   43: *> If UPLO = 'L' or 'l' then the lower triangle of the result is stored,
   44: *> overwriting the factor L in A.
   45: *>
   46: *> This is the blocked form of the algorithm, calling Level 3 BLAS.
   47: *> \endverbatim
   48: *
   49: *  Arguments:
   50: *  ==========
   51: *
   52: *> \param[in] UPLO
   53: *> \verbatim
   54: *>          UPLO is CHARACTER*1
   55: *>          Specifies whether the triangular factor stored in the array A
   56: *>          is upper or lower triangular:
   57: *>          = 'U':  Upper triangular
   58: *>          = 'L':  Lower triangular
   59: *> \endverbatim
   60: *>
   61: *> \param[in] N
   62: *> \verbatim
   63: *>          N is INTEGER
   64: *>          The order of the triangular factor U or L.  N >= 0.
   65: *> \endverbatim
   66: *>
   67: *> \param[in,out] A
   68: *> \verbatim
   69: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   70: *>          On entry, the triangular factor U or L.
   71: *>          On exit, if UPLO = 'U', the upper triangle of A is
   72: *>          overwritten with the upper triangle of the product U * U**T;
   73: *>          if UPLO = 'L', the lower triangle of A is overwritten with
   74: *>          the lower triangle of the product L**T * L.
   75: *> \endverbatim
   76: *>
   77: *> \param[in] LDA
   78: *> \verbatim
   79: *>          LDA is INTEGER
   80: *>          The leading dimension of the array A.  LDA >= max(1,N).
   81: *> \endverbatim
   82: *>
   83: *> \param[out] INFO
   84: *> \verbatim
   85: *>          INFO is INTEGER
   86: *>          = 0: successful exit
   87: *>          < 0: if INFO = -k, the k-th argument had an illegal value
   88: *> \endverbatim
   89: *
   90: *  Authors:
   91: *  ========
   92: *
   93: *> \author Univ. of Tennessee
   94: *> \author Univ. of California Berkeley
   95: *> \author Univ. of Colorado Denver
   96: *> \author NAG Ltd.
   97: *
   98: *> \ingroup doubleOTHERauxiliary
   99: *
  100: *  =====================================================================
  101:       SUBROUTINE DLAUUM( UPLO, N, A, LDA, INFO )
  102: *
  103: *  -- LAPACK auxiliary routine --
  104: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  105: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  106: *
  107: *     .. Scalar Arguments ..
  108:       CHARACTER          UPLO
  109:       INTEGER            INFO, LDA, N
  110: *     ..
  111: *     .. Array Arguments ..
  112:       DOUBLE PRECISION   A( LDA, * )
  113: *     ..
  114: *
  115: *  =====================================================================
  116: *
  117: *     .. Parameters ..
  118:       DOUBLE PRECISION   ONE
  119:       PARAMETER          ( ONE = 1.0D+0 )
  120: *     ..
  121: *     .. Local Scalars ..
  122:       LOGICAL            UPPER
  123:       INTEGER            I, IB, NB
  124: *     ..
  125: *     .. External Functions ..
  126:       LOGICAL            LSAME
  127:       INTEGER            ILAENV
  128:       EXTERNAL           LSAME, ILAENV
  129: *     ..
  130: *     .. External Subroutines ..
  131:       EXTERNAL           DGEMM, DLAUU2, DSYRK, DTRMM, XERBLA
  132: *     ..
  133: *     .. Intrinsic Functions ..
  134:       INTRINSIC          MAX, MIN
  135: *     ..
  136: *     .. Executable Statements ..
  137: *
  138: *     Test the input parameters.
  139: *
  140:       INFO = 0
  141:       UPPER = LSAME( UPLO, 'U' )
  142:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  143:          INFO = -1
  144:       ELSE IF( N.LT.0 ) THEN
  145:          INFO = -2
  146:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  147:          INFO = -4
  148:       END IF
  149:       IF( INFO.NE.0 ) THEN
  150:          CALL XERBLA( 'DLAUUM', -INFO )
  151:          RETURN
  152:       END IF
  153: *
  154: *     Quick return if possible
  155: *
  156:       IF( N.EQ.0 )
  157:      $   RETURN
  158: *
  159: *     Determine the block size for this environment.
  160: *
  161:       NB = ILAENV( 1, 'DLAUUM', UPLO, N, -1, -1, -1 )
  162: *
  163:       IF( NB.LE.1 .OR. NB.GE.N ) THEN
  164: *
  165: *        Use unblocked code
  166: *
  167:          CALL DLAUU2( UPLO, N, A, LDA, INFO )
  168:       ELSE
  169: *
  170: *        Use blocked code
  171: *
  172:          IF( UPPER ) THEN
  173: *
  174: *           Compute the product U * U**T.
  175: *
  176:             DO 10 I = 1, N, NB
  177:                IB = MIN( NB, N-I+1 )
  178:                CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-unit',
  179:      $                     I-1, IB, ONE, A( I, I ), LDA, A( 1, I ),
  180:      $                     LDA )
  181:                CALL DLAUU2( 'Upper', IB, A( I, I ), LDA, INFO )
  182:                IF( I+IB.LE.N ) THEN
  183:                   CALL DGEMM( 'No transpose', 'Transpose', I-1, IB,
  184:      $                        N-I-IB+1, ONE, A( 1, I+IB ), LDA,
  185:      $                        A( I, I+IB ), LDA, ONE, A( 1, I ), LDA )
  186:                   CALL DSYRK( 'Upper', 'No transpose', IB, N-I-IB+1,
  187:      $                        ONE, A( I, I+IB ), LDA, ONE, A( I, I ),
  188:      $                        LDA )
  189:                END IF
  190:    10       CONTINUE
  191:          ELSE
  192: *
  193: *           Compute the product L**T * L.
  194: *
  195:             DO 20 I = 1, N, NB
  196:                IB = MIN( NB, N-I+1 )
  197:                CALL DTRMM( 'Left', 'Lower', 'Transpose', 'Non-unit', IB,
  198:      $                     I-1, ONE, A( I, I ), LDA, A( I, 1 ), LDA )
  199:                CALL DLAUU2( 'Lower', IB, A( I, I ), LDA, INFO )
  200:                IF( I+IB.LE.N ) THEN
  201:                   CALL DGEMM( 'Transpose', 'No transpose', IB, I-1,
  202:      $                        N-I-IB+1, ONE, A( I+IB, I ), LDA,
  203:      $                        A( I+IB, 1 ), LDA, ONE, A( I, 1 ), LDA )
  204:                   CALL DSYRK( 'Lower', 'Transpose', IB, N-I-IB+1, ONE,
  205:      $                        A( I+IB, I ), LDA, ONE, A( I, I ), LDA )
  206:                END IF
  207:    20       CONTINUE
  208:          END IF
  209:       END IF
  210: *
  211:       RETURN
  212: *
  213: *     End of DLAUUM
  214: *
  215:       END