Return to zhpr2.f CVS log

Up to [local] / rpl / lapack / blas

Cohérence.

    1: *> \brief \b ZHPR2
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *  Definition:
    9: *  ===========
   10: *
   11: *       SUBROUTINE ZHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
   12: * 
   13: *       .. Scalar Arguments ..
   14: *       COMPLEX*16 ALPHA
   15: *       INTEGER INCX,INCY,N
   16: *       CHARACTER UPLO
   17: *       ..
   18: *       .. Array Arguments ..
   19: *       COMPLEX*16 AP(*),X(*),Y(*)
   20: *       ..
   21: *  
   22: *
   23: *> \par Purpose:
   24: *  =============
   25: *>
   26: *> \verbatim
   27: *>
   28: *> ZHPR2  performs the hermitian rank 2 operation
   29: *>
   30: *>    A := alpha*x*y**H + conjg( alpha )*y*x**H + A,
   31: *>
   32: *> where alpha is a scalar, x and y are n element vectors and A is an
   33: *> n by n hermitian matrix, supplied in packed form.
   34: *> \endverbatim
   35: *
   36: *  Arguments:
   37: *  ==========
   38: *
   39: *> \param[in] UPLO
   40: *> \verbatim
   41: *>          UPLO is CHARACTER*1
   42: *>           On entry, UPLO specifies whether the upper or lower
   43: *>           triangular part of the matrix A is supplied in the packed
   44: *>           array AP as follows:
   45: *>
   46: *>              UPLO = 'U' or 'u'   The upper triangular part of A is
   47: *>                                  supplied in AP.
   48: *>
   49: *>              UPLO = 'L' or 'l'   The lower triangular part of A is
   50: *>                                  supplied in AP.
   51: *> \endverbatim
   52: *>
   53: *> \param[in] N
   54: *> \verbatim
   55: *>          N is INTEGER
   56: *>           On entry, N specifies the order of the matrix A.
   57: *>           N must be at least zero.
   58: *> \endverbatim
   59: *>
   60: *> \param[in] ALPHA
   61: *> \verbatim
   62: *>          ALPHA is COMPLEX*16
   63: *>           On entry, ALPHA specifies the scalar alpha.
   64: *> \endverbatim
   65: *>
   66: *> \param[in] X
   67: *> \verbatim
   68: *>          X is COMPLEX*16 array of dimension at least
   69: *>           ( 1 + ( n - 1 )*abs( INCX ) ).
   70: *>           Before entry, the incremented array X must contain the n
   71: *>           element vector x.
   72: *> \endverbatim
   73: *>
   74: *> \param[in] INCX
   75: *> \verbatim
   76: *>          INCX is INTEGER
   77: *>           On entry, INCX specifies the increment for the elements of
   78: *>           X. INCX must not be zero.
   79: *> \endverbatim
   80: *>
   81: *> \param[in] Y
   82: *> \verbatim
   83: *>          Y is COMPLEX*16 array of dimension at least
   84: *>           ( 1 + ( n - 1 )*abs( INCY ) ).
   85: *>           Before entry, the incremented array Y must contain the n
   86: *>           element vector y.
   87: *> \endverbatim
   88: *>
   89: *> \param[in] INCY
   90: *> \verbatim
   91: *>          INCY is INTEGER
   92: *>           On entry, INCY specifies the increment for the elements of
   93: *>           Y. INCY must not be zero.
   94: *> \endverbatim
   95: *>
   96: *> \param[in,out] AP
   97: *> \verbatim
   98: *>          AP is COMPLEX*16 array of DIMENSION at least
   99: *>           ( ( n*( n + 1 ) )/2 ).
  100: *>           Before entry with  UPLO = 'U' or 'u', the array AP must
  101: *>           contain the upper triangular part of the hermitian matrix
  102: *>           packed sequentially, column by column, so that AP( 1 )
  103: *>           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
  104: *>           and a( 2, 2 ) respectively, and so on. On exit, the array
  105: *>           AP is overwritten by the upper triangular part of the
  106: *>           updated matrix.
  107: *>           Before entry with UPLO = 'L' or 'l', the array AP must
  108: *>           contain the lower triangular part of the hermitian matrix
  109: *>           packed sequentially, column by column, so that AP( 1 )
  110: *>           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  111: *>           and a( 3, 1 ) respectively, and so on. On exit, the array
  112: *>           AP is overwritten by the lower triangular part of the
  113: *>           updated matrix.
  114: *>           Note that the imaginary parts of the diagonal elements need
  115: *>           not be set, they are assumed to be zero, and on exit they
  116: *>           are set to zero.
  117: *> \endverbatim
  118: *
  119: *  Authors:
  120: *  ========
  121: *
  122: *> \author Univ. of Tennessee 
  123: *> \author Univ. of California Berkeley 
  124: *> \author Univ. of Colorado Denver 
  125: *> \author NAG Ltd. 
  126: *
  127: *> \date November 2011
  128: *
  129: *> \ingroup complex16_blas_level2
  130: *
  131: *> \par Further Details:
  132: *  =====================
  133: *>
  134: *> \verbatim
  135: *>
  136: *>  Level 2 Blas routine.
  137: *>
  138: *>  -- Written on 22-October-1986.
  139: *>     Jack Dongarra, Argonne National Lab.
  140: *>     Jeremy Du Croz, Nag Central Office.
  141: *>     Sven Hammarling, Nag Central Office.
  142: *>     Richard Hanson, Sandia National Labs.
  143: *> \endverbatim
  144: *>
  145: *  =====================================================================
  146:       SUBROUTINE ZHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
  147: *
  148: *  -- Reference BLAS level2 routine (version 3.4.0) --
  149: *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
  150: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  151: *     November 2011
  152: *
  153: *     .. Scalar Arguments ..
  154:       COMPLEX*16 ALPHA
  155:       INTEGER INCX,INCY,N
  156:       CHARACTER UPLO
  157: *     ..
  158: *     .. Array Arguments ..
  159:       COMPLEX*16 AP(*),X(*),Y(*)
  160: *     ..
  161: *
  162: *  =====================================================================
  163: *
  164: *     .. Parameters ..
  165:       COMPLEX*16 ZERO
  166:       PARAMETER (ZERO= (0.0D+0,0.0D+0))
  167: *     ..
  168: *     .. Local Scalars ..
  169:       COMPLEX*16 TEMP1,TEMP2
  170:       INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
  171: *     ..
  172: *     .. External Functions ..
  173:       LOGICAL LSAME
  174:       EXTERNAL LSAME
  175: *     ..
  176: *     .. External Subroutines ..
  177:       EXTERNAL XERBLA
  178: *     ..
  179: *     .. Intrinsic Functions ..
  180:       INTRINSIC DBLE,DCONJG
  181: *     ..
  182: *
  183: *     Test the input parameters.
  184: *
  185:       INFO = 0
  186:       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  187:           INFO = 1
  188:       ELSE IF (N.LT.0) THEN
  189:           INFO = 2
  190:       ELSE IF (INCX.EQ.0) THEN
  191:           INFO = 5
  192:       ELSE IF (INCY.EQ.0) THEN
  193:           INFO = 7
  194:       END IF
  195:       IF (INFO.NE.0) THEN
  196:           CALL XERBLA('ZHPR2 ',INFO)
  197:           RETURN
  198:       END IF
  199: *
  200: *     Quick return if possible.
  201: *
  202:       IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
  203: *
  204: *     Set up the start points in X and Y if the increments are not both
  205: *     unity.
  206: *
  207:       IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
  208:           IF (INCX.GT.0) THEN
  209:               KX = 1
  210:           ELSE
  211:               KX = 1 - (N-1)*INCX
  212:           END IF
  213:           IF (INCY.GT.0) THEN
  214:               KY = 1
  215:           ELSE
  216:               KY = 1 - (N-1)*INCY
  217:           END IF
  218:           JX = KX
  219:           JY = KY
  220:       END IF
  221: *
  222: *     Start the operations. In this version the elements of the array AP
  223: *     are accessed sequentially with one pass through AP.
  224: *
  225:       KK = 1
  226:       IF (LSAME(UPLO,'U')) THEN
  227: *
  228: *        Form  A  when upper triangle is stored in AP.
  229: *
  230:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  231:               DO 20 J = 1,N
  232:                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
  233:                       TEMP1 = ALPHA*DCONJG(Y(J))
  234:                       TEMP2 = DCONJG(ALPHA*X(J))
  235:                       K = KK
  236:                       DO 10 I = 1,J - 1
  237:                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
  238:                           K = K + 1
  239:    10                 CONTINUE
  240:                       AP(KK+J-1) = DBLE(AP(KK+J-1)) +
  241:      +                             DBLE(X(J)*TEMP1+Y(J)*TEMP2)
  242:                   ELSE
  243:                       AP(KK+J-1) = DBLE(AP(KK+J-1))
  244:                   END IF
  245:                   KK = KK + J
  246:    20         CONTINUE
  247:           ELSE
  248:               DO 40 J = 1,N
  249:                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
  250:                       TEMP1 = ALPHA*DCONJG(Y(JY))
  251:                       TEMP2 = DCONJG(ALPHA*X(JX))
  252:                       IX = KX
  253:                       IY = KY
  254:                       DO 30 K = KK,KK + J - 2
  255:                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
  256:                           IX = IX + INCX
  257:                           IY = IY + INCY
  258:    30                 CONTINUE
  259:                       AP(KK+J-1) = DBLE(AP(KK+J-1)) +
  260:      +                             DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
  261:                   ELSE
  262:                       AP(KK+J-1) = DBLE(AP(KK+J-1))
  263:                   END IF
  264:                   JX = JX + INCX
  265:                   JY = JY + INCY
  266:                   KK = KK + J
  267:    40         CONTINUE
  268:           END IF
  269:       ELSE
  270: *
  271: *        Form  A  when lower triangle is stored in AP.
  272: *
  273:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  274:               DO 60 J = 1,N
  275:                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
  276:                       TEMP1 = ALPHA*DCONJG(Y(J))
  277:                       TEMP2 = DCONJG(ALPHA*X(J))
  278:                       AP(KK) = DBLE(AP(KK)) +
  279:      +                         DBLE(X(J)*TEMP1+Y(J)*TEMP2)
  280:                       K = KK + 1
  281:                       DO 50 I = J + 1,N
  282:                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
  283:                           K = K + 1
  284:    50                 CONTINUE
  285:                   ELSE
  286:                       AP(KK) = DBLE(AP(KK))
  287:                   END IF
  288:                   KK = KK + N - J + 1
  289:    60         CONTINUE
  290:           ELSE
  291:               DO 80 J = 1,N
  292:                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
  293:                       TEMP1 = ALPHA*DCONJG(Y(JY))
  294:                       TEMP2 = DCONJG(ALPHA*X(JX))
  295:                       AP(KK) = DBLE(AP(KK)) +
  296:      +                         DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
  297:                       IX = JX
  298:                       IY = JY
  299:                       DO 70 K = KK + 1,KK + N - J
  300:                           IX = IX + INCX
  301:                           IY = IY + INCY
  302:                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
  303:    70                 CONTINUE
  304:                   ELSE
  305:                       AP(KK) = DBLE(AP(KK))
  306:                   END IF
  307:                   JX = JX + INCX
  308:                   JY = JY + INCY
  309:                   KK = KK + N - J + 1
  310:    80         CONTINUE
  311:           END IF
  312:       END IF
  313: *
  314:       RETURN
  315: *
  316: *     End of ZHPR2 .
  317: *
  318:       END