SRC/cspr.f

   1 *> \brief \b CSPR performs the symmetrical rank-1 update of a complex symmetric packed matrix.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CSPR + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cspr.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cspr.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cspr.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CSPR( UPLO, N, ALPHA, X, INCX, AP )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INCX, N
  26 *       COMPLEX            ALPHA
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       COMPLEX            AP( * ), X( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> CSPR    performs the symmetric rank 1 operation
  39 *>
  40 *>    A := alpha*x*x**H + A,
  41 *>
  42 *> where alpha is a complex scalar, x is an n element vector and A is an
  43 *> n by n symmetric matrix, supplied in packed form.
  44 *> \endverbatim
  45 *
  46 *  Arguments:
  47 *  ==========
  48 *
  49 *> \param[in] UPLO
  50 *> \verbatim
  51 *>          UPLO is CHARACTER*1
  52 *>           On entry, UPLO specifies whether the upper or lower
  53 *>           triangular part of the matrix A is supplied in the packed
  54 *>           array AP as follows:
  55 *>
  56 *>              UPLO = 'U' or 'u'   The upper triangular part of A is
  57 *>                                  supplied in AP.
  58 *>
  59 *>              UPLO = 'L' or 'l'   The lower triangular part of A is
  60 *>                                  supplied in AP.
  61 *>
  62 *>           Unchanged on exit.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] N
  66 *> \verbatim
  67 *>          N is INTEGER
  68 *>           On entry, N specifies the order of the matrix A.
  69 *>           N must be at least zero.
  70 *>           Unchanged on exit.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] ALPHA
  74 *> \verbatim
  75 *>          ALPHA is COMPLEX
  76 *>           On entry, ALPHA specifies the scalar alpha.
  77 *>           Unchanged on exit.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] X
  81 *> \verbatim
  82 *>          X is COMPLEX array, dimension at least
  83 *>           ( 1 + ( N - 1 )*abs( INCX ) ).
  84 *>           Before entry, the incremented array X must contain the N-
  85 *>           element vector x.
  86 *>           Unchanged on exit.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] INCX
  90 *> \verbatim
  91 *>          INCX is INTEGER
  92 *>           On entry, INCX specifies the increment for the elements of
  93 *>           X. INCX must not be zero.
  94 *>           Unchanged on exit.
  95 *> \endverbatim
  96 *>
  97 *> \param[in,out] AP
  98 *> \verbatim
  99 *>          AP is COMPLEX array, dimension at least
 100 *>           ( ( N*( N + 1 ) )/2 ).
 101 *>           Before entry, with  UPLO = 'U' or 'u', the array AP must
 102 *>           contain the upper triangular part of the symmetric matrix
 103 *>           packed sequentially, column by column, so that AP( 1 )
 104 *>           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
 105 *>           and a( 2, 2 ) respectively, and so on. On exit, the array
 106 *>           AP is overwritten by the upper triangular part of the
 107 *>           updated matrix.
 108 *>           Before entry, with UPLO = 'L' or 'l', the array AP must
 109 *>           contain the lower triangular part of the symmetric matrix
 110 *>           packed sequentially, column by column, so that AP( 1 )
 111 *>           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
 112 *>           and a( 3, 1 ) respectively, and so on. On exit, the array
 113 *>           AP is overwritten by the lower triangular part of the
 114 *>           updated matrix.
 115 *>           Note that the imaginary parts of the diagonal elements need
 116 *>           not be set, they are assumed to be zero, and on exit they
 117 *>           are set to zero.
 118 *> \endverbatim
 119 *
 120 *  Authors:
 121 *  ========
 122 *
 123 *> \author Univ. of Tennessee
 124 *> \author Univ. of California Berkeley
 125 *> \author Univ. of Colorado Denver
 126 *> \author NAG Ltd.
 127 *
 128 *> \date September 2012
 129 *
 130 *> \ingroup complexOTHERauxiliary
 131 *
 132 *  =====================================================================
 133       SUBROUTINE CSPR( UPLO, N, ALPHA, X, INCX, AP )
 134 *
 135 *  -- LAPACK auxiliary routine (version 3.4.2) --
 136 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 137 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 138 *     September 2012
 139 *
 140 *     .. Scalar Arguments ..
 141       CHARACTER          UPLO
 142       INTEGER            INCX, N
 143       COMPLEX            ALPHA
 144 *     ..
 145 *     .. Array Arguments ..
 146       COMPLEX            AP( * ), X( * )
 147 *     ..
 148 *
 149 * =====================================================================
 150 *
 151 *     .. Parameters ..
 152       COMPLEX            ZERO
 153       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ) )
 154 *     ..
 155 *     .. Local Scalars ..
 156       INTEGER            I, INFO, IX, J, JX, K, KK, KX
 157       COMPLEX            TEMP
 158 *     ..
 159 *     .. External Functions ..
 160       LOGICAL            LSAME
 161       EXTERNAL           LSAME
 162 *     ..
 163 *     .. External Subroutines ..
 164       EXTERNAL           XERBLA
 165 *     ..
 166 *     .. Executable Statements ..
 167 *
 168 *     Test the input parameters.
 169 *
 170       INFO = 0
 171       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 172          INFO = 1
 173       ELSE IF( N.LT.0 ) THEN
 174          INFO = 2
 175       ELSE IF( INCX.EQ.0 ) THEN
 176          INFO = 5
 177       END IF
 178       IF( INFO.NE.0 ) THEN
 179          CALL XERBLA( 'CSPR  ', INFO )
 180          RETURN
 181       END IF
 182 *
 183 *     Quick return if possible.
 184 *
 185       IF( ( N.EQ.0 ) .OR. ( ALPHA.EQ.ZERO ) )
 186      $   RETURN
 187 *
 188 *     Set the start point in X if the increment is not unity.
 189 *
 190       IF( INCX.LE.0 ) THEN
 191          KX = 1 - ( N-1 )*INCX
 192       ELSE IF( INCX.NE.1 ) THEN
 193          KX = 1
 194       END IF
 195 *
 196 *     Start the operations. In this version the elements of the array AP
 197 *     are accessed sequentially with one pass through AP.
 198 *
 199       KK = 1
 200       IF( LSAME( UPLO, 'U' ) ) THEN
 201 *
 202 *        Form  A  when upper triangle is stored in AP.
 203 *
 204          IF( INCX.EQ.1 ) THEN
 205             DO 20 J = 1, N
 206                IF( X( J ).NE.ZERO ) THEN
 207                   TEMP = ALPHA*X( J )
 208                   K = KK
 209                   DO 10 I = 1, J - 1
 210                      AP( K ) = AP( K ) + X( I )*TEMP
 211                      K = K + 1
 212    10             CONTINUE
 213                   AP( KK+J-1 ) = AP( KK+J-1 ) + X( J )*TEMP
 214                ELSE
 215                   AP( KK+J-1 ) = AP( KK+J-1 )
 216                END IF
 217                KK = KK + J
 218    20       CONTINUE
 219          ELSE
 220             JX = KX
 221             DO 40 J = 1, N
 222                IF( X( JX ).NE.ZERO ) THEN
 223                   TEMP = ALPHA*X( JX )
 224                   IX = KX
 225                   DO 30 K = KK, KK + J - 2
 226                      AP( K ) = AP( K ) + X( IX )*TEMP
 227                      IX = IX + INCX
 228    30             CONTINUE
 229                   AP( KK+J-1 ) = AP( KK+J-1 ) + X( JX )*TEMP
 230                ELSE
 231                   AP( KK+J-1 ) = AP( KK+J-1 )
 232                END IF
 233                JX = JX + INCX
 234                KK = KK + J
 235    40       CONTINUE
 236          END IF
 237       ELSE
 238 *
 239 *        Form  A  when lower triangle is stored in AP.
 240 *
 241          IF( INCX.EQ.1 ) THEN
 242             DO 60 J = 1, N
 243                IF( X( J ).NE.ZERO ) THEN
 244                   TEMP = ALPHA*X( J )
 245                   AP( KK ) = AP( KK ) + TEMP*X( J )
 246                   K = KK + 1
 247                   DO 50 I = J + 1, N
 248                      AP( K ) = AP( K ) + X( I )*TEMP
 249                      K = K + 1
 250    50             CONTINUE
 251                ELSE
 252                   AP( KK ) = AP( KK )
 253                END IF
 254                KK = KK + N - J + 1
 255    60       CONTINUE
 256          ELSE
 257             JX = KX
 258             DO 80 J = 1, N
 259                IF( X( JX ).NE.ZERO ) THEN
 260                   TEMP = ALPHA*X( JX )
 261                   AP( KK ) = AP( KK ) + TEMP*X( JX )
 262                   IX = JX
 263                   DO 70 K = KK + 1, KK + N - J
 264                      IX = IX + INCX
 265                      AP( K ) = AP( K ) + X( IX )*TEMP
 266    70             CONTINUE
 267                ELSE
 268                   AP( KK ) = AP( KK )
 269                END IF
 270                JX = JX + INCX
 271                KK = KK + N - J + 1
 272    80       CONTINUE
 273          END IF
 274       END IF
 275 *
 276       RETURN
 277 *
 278 *     End of CSPR
 279 *
 280       END