SRC/dlauum.f

   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 *> \date September 2012
  99 *
 100 *> \ingroup doubleOTHERauxiliary
 101 *
 102 *  =====================================================================
 103       SUBROUTINE DLAUUM( UPLO, N, A, LDA, INFO )
 104 *
 105 *  -- LAPACK auxiliary routine (version 3.4.2) --
 106 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 107 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 108 *     September 2012
 109 *
 110 *     .. Scalar Arguments ..
 111       CHARACTER          UPLO
 112       INTEGER            INFO, LDA, N
 113 *     ..
 114 *     .. Array Arguments ..
 115       DOUBLE PRECISION   A( LDA, * )
 116 *     ..
 117 *
 118 *  =====================================================================
 119 *
 120 *     .. Parameters ..
 121       DOUBLE PRECISION   ONE
 122       PARAMETER          ( ONE = 1.0D+0 )
 123 *     ..
 124 *     .. Local Scalars ..
 125       LOGICAL            UPPER
 126       INTEGER            I, IB, NB
 127 *     ..
 128 *     .. External Functions ..
 129       LOGICAL            LSAME
 130       INTEGER            ILAENV
 131       EXTERNAL           LSAME, ILAENV
 132 *     ..
 133 *     .. External Subroutines ..
 134       EXTERNAL           DGEMM, DLAUU2, DSYRK, DTRMM, XERBLA
 135 *     ..
 136 *     .. Intrinsic Functions ..
 137       INTRINSIC          MAX, MIN
 138 *     ..
 139 *     .. Executable Statements ..
 140 *
 141 *     Test the input parameters.
 142 *
 143       INFO = 0
 144       UPPER = LSAME( UPLO, 'U' )
 145       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 146          INFO = -1
 147       ELSE IF( N.LT.0 ) THEN
 148          INFO = -2
 149       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 150          INFO = -4
 151       END IF
 152       IF( INFO.NE.0 ) THEN
 153          CALL XERBLA( 'DLAUUM', -INFO )
 154          RETURN
 155       END IF
 156 *
 157 *     Quick return if possible
 158 *
 159       IF( N.EQ.0 )
 160      $   RETURN
 161 *
 162 *     Determine the block size for this environment.
 163 *
 164       NB = ILAENV( 1, 'DLAUUM', UPLO, N, -1, -1, -1 )
 165 *
 166       IF( NB.LE.1 .OR. NB.GE.N ) THEN
 167 *
 168 *        Use unblocked code
 169 *
 170          CALL DLAUU2( UPLO, N, A, LDA, INFO )
 171       ELSE
 172 *
 173 *        Use blocked code
 174 *
 175          IF( UPPER ) THEN
 176 *
 177 *           Compute the product U * U**T.
 178 *
 179             DO 10 I = 1, N, NB
 180                IB = MIN( NB, N-I+1 )
 181                CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-unit',
 182      $                     I-1, IB, ONE, A( I, I ), LDA, A( 1, I ),
 183      $                     LDA )
 184                CALL DLAUU2( 'Upper', IB, A( I, I ), LDA, INFO )
 185                IF( I+IB.LE.N ) THEN
 186                   CALL DGEMM( 'No transpose', 'Transpose', I-1, IB,
 187      $                        N-I-IB+1, ONE, A( 1, I+IB ), LDA,
 188      $                        A( I, I+IB ), LDA, ONE, A( 1, I ), LDA )
 189                   CALL DSYRK( 'Upper', 'No transpose', IB, N-I-IB+1,
 190      $                        ONE, A( I, I+IB ), LDA, ONE, A( I, I ),
 191      $                        LDA )
 192                END IF
 193    10       CONTINUE
 194          ELSE
 195 *
 196 *           Compute the product L**T * L.
 197 *
 198             DO 20 I = 1, N, NB
 199                IB = MIN( NB, N-I+1 )
 200                CALL DTRMM( 'Left', 'Lower', 'Transpose', 'Non-unit', IB,
 201      $                     I-1, ONE, A( I, I ), LDA, A( I, 1 ), LDA )
 202                CALL DLAUU2( 'Lower', IB, A( I, I ), LDA, INFO )
 203                IF( I+IB.LE.N ) THEN
 204                   CALL DGEMM( 'Transpose', 'No transpose', IB, I-1,
 205      $                        N-I-IB+1, ONE, A( I+IB, I ), LDA,
 206      $                        A( I+IB, 1 ), LDA, ONE, A( I, 1 ), LDA )
 207                   CALL DSYRK( 'Lower', 'Transpose', IB, N-I-IB+1, ONE,
 208      $                        A( I+IB, I ), LDA, ONE, A( I, I ), LDA )
 209                END IF
 210    20       CONTINUE
 211          END IF
 212       END IF
 213 *
 214       RETURN
 215 *
 216 *     End of DLAUUM
 217 *
 218       END