/*@
VecDot - Computes the vector dot product.
Collective on Vec
Input Parameters:
. x, y - the vectors
Output Parameter:
. val - the dot product
Performance Issues:
$ per-processor memory bandwidth
$ interprocessor latency
$ work load inbalance that causes certain processes to arrive much earlier
than others
Notes for Users of Complex Numbers:
For complex vectors, VecDot() computes
$ val = (x,y) = y^H x,
where y^H denotes the conjugate transpose of y.
Use VecTDot() for the indefinite form
$ val = (x,y) = y^T x,
where y^T denotes the transpose of y.
Level: intermediate
Concepts: inner product
Concepts: vector^inner product
.seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd()
/*@
VecTDot - Computes an indefinite vector dot product. That is, this
routine does NOT use the complex conjugate.
Collective on Vec
Input Parameters:
. x, y - the vectors
Output Parameter:
. val - the dot product
Notes for Users of Complex Numbers:
For complex vectors, VecTDot() computes the indefinite form
$ val = (x,y) = y^T x,
where y^T denotes the transpose of y.
Use VecDot() for the inner product
$ val = (x,y) = y^H x,
where y^H denotes the conjugate transpose of y.
Level: intermediate
Concepts: inner product^non-Hermitian
Concepts: vector^inner product
Concepts: non-Hermitian inner product
.seealso: VecDot(), VecMTDot()
@*/
But source code has
PetscErrorCode VecDot_Seq(Vec xin,Vec yin,PetscScalar *z)
{
const PetscScalar *ya,*xa;
PetscBLASInt one = 1,bn = PetscBLASIntCast(xin->map->n);
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = VecGetArrayRead(xin,&xa);CHKERRQ(ierr);
ierr = VecGetArrayRead(yin,&ya);CHKERRQ(ierr);
*z = BLASdot_(&bn,xa,&one,ya,&one);
ierr = VecRestoreArrayRead(xin,&xa);CHKERRQ(ierr);
ierr = VecRestoreArrayRead(yin,&ya);CHKERRQ(ierr);
if (xin->map->n > 0) {
ierr = PetscLogFlops(2.0*xin->map->n-1);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "VecTDot_Seq"
PetscErrorCode VecTDot_Seq(Vec xin,Vec yin,PetscScalar *z)
{
const PetscScalar *ya,*xa;
PetscBLASInt one = 1, bn = PetscBLASIntCast(xin->map->n);
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = VecGetArrayRead(xin,&xa);CHKERRQ(ierr);
ierr = VecGetArrayRead(yin,&ya);CHKERRQ(ierr);
*z = BLASdot_(&bn,xa,&one,ya,&one);
ierr = VecRestoreArrayRead(xin,&xa);CHKERRQ(ierr);
ierr = VecRestoreArrayRead(yin,&ya);CHKERRQ(ierr);
if (xin->map->n > 0) {
ierr = PetscLogFlops(2.0*xin->map->n-1);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
How can they both use BLASdot() when one is supposed to use the complex
conjugate of y the second argument.
Note also
DOUBLE COMPLEX FUNCTION ZDOTC(N,ZX,INCX,ZY,INCY)
* .. Scalar Arguments ..
INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
DOUBLE COMPLEX ZX(*),ZY(*)
* ..
*
* Purpose
* =======
*
* ZDOTC forms the dot product of a vector.
*
* Further Details
* ===============
*
* jack dongarra, 3/11/78.
* modified 12/3/93, array(1) declarations changed to array(*)
*
IX = 1
IY = 1
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
DO I = 1,N
ZTEMP = ZTEMP + DCONJG(ZX(IX))*ZY(IY)
IX = IX + INCX
IY = IY + INCY
END DO
END IF
ZDOTC = ZTEMP
uses the complex conjugate of the FIRST argument.
So it looks like everything is major broken for the complex case.
Barry
