File:  [local] / rpl / lapack / blas / zhpr2.f
Revision 1.16: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:38:45 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    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, 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, 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, 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: *> \ingroup complex16_blas_level2
  128: *
  129: *> \par Further Details:
  130: *  =====================
  131: *>
  132: *> \verbatim
  133: *>
  134: *>  Level 2 Blas routine.
  135: *>
  136: *>  -- Written on 22-October-1986.
  137: *>     Jack Dongarra, Argonne National Lab.
  138: *>     Jeremy Du Croz, Nag Central Office.
  139: *>     Sven Hammarling, Nag Central Office.
  140: *>     Richard Hanson, Sandia National Labs.
  141: *> \endverbatim
  142: *>
  143: *  =====================================================================
  144:       SUBROUTINE ZHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
  145: *
  146: *  -- Reference BLAS level2 routine --
  147: *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
  148: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  149: *
  150: *     .. Scalar Arguments ..
  151:       COMPLEX*16 ALPHA
  152:       INTEGER INCX,INCY,N
  153:       CHARACTER UPLO
  154: *     ..
  155: *     .. Array Arguments ..
  156:       COMPLEX*16 AP(*),X(*),Y(*)
  157: *     ..
  158: *
  159: *  =====================================================================
  160: *
  161: *     .. Parameters ..
  162:       COMPLEX*16 ZERO
  163:       PARAMETER (ZERO= (0.0D+0,0.0D+0))
  164: *     ..
  165: *     .. Local Scalars ..
  166:       COMPLEX*16 TEMP1,TEMP2
  167:       INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
  168: *     ..
  169: *     .. External Functions ..
  170:       LOGICAL LSAME
  171:       EXTERNAL LSAME
  172: *     ..
  173: *     .. External Subroutines ..
  174:       EXTERNAL XERBLA
  175: *     ..
  176: *     .. Intrinsic Functions ..
  177:       INTRINSIC DBLE,DCONJG
  178: *     ..
  179: *
  180: *     Test the input parameters.
  181: *
  182:       INFO = 0
  183:       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  184:           INFO = 1
  185:       ELSE IF (N.LT.0) THEN
  186:           INFO = 2
  187:       ELSE IF (INCX.EQ.0) THEN
  188:           INFO = 5
  189:       ELSE IF (INCY.EQ.0) THEN
  190:           INFO = 7
  191:       END IF
  192:       IF (INFO.NE.0) THEN
  193:           CALL XERBLA('ZHPR2 ',INFO)
  194:           RETURN
  195:       END IF
  196: *
  197: *     Quick return if possible.
  198: *
  199:       IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
  200: *
  201: *     Set up the start points in X and Y if the increments are not both
  202: *     unity.
  203: *
  204:       IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
  205:           IF (INCX.GT.0) THEN
  206:               KX = 1
  207:           ELSE
  208:               KX = 1 - (N-1)*INCX
  209:           END IF
  210:           IF (INCY.GT.0) THEN
  211:               KY = 1
  212:           ELSE
  213:               KY = 1 - (N-1)*INCY
  214:           END IF
  215:           JX = KX
  216:           JY = KY
  217:       END IF
  218: *
  219: *     Start the operations. In this version the elements of the array AP
  220: *     are accessed sequentially with one pass through AP.
  221: *
  222:       KK = 1
  223:       IF (LSAME(UPLO,'U')) THEN
  224: *
  225: *        Form  A  when upper triangle is stored in AP.
  226: *
  227:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  228:               DO 20 J = 1,N
  229:                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
  230:                       TEMP1 = ALPHA*DCONJG(Y(J))
  231:                       TEMP2 = DCONJG(ALPHA*X(J))
  232:                       K = KK
  233:                       DO 10 I = 1,J - 1
  234:                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
  235:                           K = K + 1
  236:    10                 CONTINUE
  237:                       AP(KK+J-1) = DBLE(AP(KK+J-1)) +
  238:      +                             DBLE(X(J)*TEMP1+Y(J)*TEMP2)
  239:                   ELSE
  240:                       AP(KK+J-1) = DBLE(AP(KK+J-1))
  241:                   END IF
  242:                   KK = KK + J
  243:    20         CONTINUE
  244:           ELSE
  245:               DO 40 J = 1,N
  246:                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
  247:                       TEMP1 = ALPHA*DCONJG(Y(JY))
  248:                       TEMP2 = DCONJG(ALPHA*X(JX))
  249:                       IX = KX
  250:                       IY = KY
  251:                       DO 30 K = KK,KK + J - 2
  252:                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
  253:                           IX = IX + INCX
  254:                           IY = IY + INCY
  255:    30                 CONTINUE
  256:                       AP(KK+J-1) = DBLE(AP(KK+J-1)) +
  257:      +                             DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
  258:                   ELSE
  259:                       AP(KK+J-1) = DBLE(AP(KK+J-1))
  260:                   END IF
  261:                   JX = JX + INCX
  262:                   JY = JY + INCY
  263:                   KK = KK + J
  264:    40         CONTINUE
  265:           END IF
  266:       ELSE
  267: *
  268: *        Form  A  when lower triangle is stored in AP.
  269: *
  270:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  271:               DO 60 J = 1,N
  272:                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
  273:                       TEMP1 = ALPHA*DCONJG(Y(J))
  274:                       TEMP2 = DCONJG(ALPHA*X(J))
  275:                       AP(KK) = DBLE(AP(KK)) +
  276:      +                         DBLE(X(J)*TEMP1+Y(J)*TEMP2)
  277:                       K = KK + 1
  278:                       DO 50 I = J + 1,N
  279:                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
  280:                           K = K + 1
  281:    50                 CONTINUE
  282:                   ELSE
  283:                       AP(KK) = DBLE(AP(KK))
  284:                   END IF
  285:                   KK = KK + N - J + 1
  286:    60         CONTINUE
  287:           ELSE
  288:               DO 80 J = 1,N
  289:                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
  290:                       TEMP1 = ALPHA*DCONJG(Y(JY))
  291:                       TEMP2 = DCONJG(ALPHA*X(JX))
  292:                       AP(KK) = DBLE(AP(KK)) +
  293:      +                         DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
  294:                       IX = JX
  295:                       IY = JY
  296:                       DO 70 K = KK + 1,KK + N - J
  297:                           IX = IX + INCX
  298:                           IY = IY + INCY
  299:                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
  300:    70                 CONTINUE
  301:                   ELSE
  302:                       AP(KK) = DBLE(AP(KK))
  303:                   END IF
  304:                   JX = JX + INCX
  305:                   JY = JY + INCY
  306:                   KK = KK + N - J + 1
  307:    80         CONTINUE
  308:           END IF
  309:       END IF
  310: *
  311:       RETURN
  312: *
  313: *     End of ZHPR2
  314: *
  315:       END

CVSweb interface <joel.bertrand@systella.fr>