[Trilinos-Users] Tpetra::dot for SCALAR=std::complex

[Nico Schlömer]

Hi,

I'm calculating dot-products of Tpetra::Vectors with complex scalar types,
and I was just trying to see whether Tpetra is smart enough to translate

    alpha = z0.dot( z1 ) // source code

into

    \alpha = z_0^H z_1 // maths

or whether it does

    \alpha = z_0^T z_1.

Well -- turns out it's neither of the two. Tpetra::dot apparently treats
std::real(z0) and std::imag(z0) as separate vector components, calculates
the |R-dot product of the 2*N vector, and effectively returns

    \alpha = \Re( z_0^T z_1 )

This is with the latest Git. A minimal example is attached to this mail.

Anyway, this is certainly not what I would expect a dot-product to do. I
suggest that either the documentation is updated, or rather the code fixed.
How about a test for this, too?

Cheers,
Nico
-------------- next part --------------
#include <Tpetra_Vector.hpp>
#include <Teuchos_DefaultComm.hpp>

typedef Tpetra::Vector<std::complex<double>,int> ComplexVector;

using namespace std;
using namespace Teuchos;

// ============================================================================
int
main( int argc, char* argv[] )
{
    // Initialize MPI
#ifdef HAVE_MPI
    MPI_Init ( &argc,&argv );
#endif

  const RCP<const Comm<int> > Comm = DefaultComm<int>::getComm();

  int N = 1;
  const RCP<Tpetra::Map<int> > complexMap = rcp ( new Tpetra::Map<int> ( N, 0, Comm ) );

  ComplexVector v0( complexMap );
  ComplexVector v1( complexMap );

  v0.putScalar( complex<double>(1.0,0.0) );
  v1.putScalar( complex<double>(0.0,1.0) );

  // expected:: (0,1), (0,-1)
  cout << v0.dot( v1 ) << endl;
  cout << v1.dot( v0 ) << endl;

#ifdef HAVE_MPI
    MPI_Finalize();
#endif

  return 0;
}
// ============================================================================
-------------- next part --------------
TRILINOS_ROOT=/opt/trilinos/dev/gcc/4.3.4/
include ${TRILINOS_ROOT}/include/Makefile.client.Tpetra

EXE=tpetraTest.exe

GFORTRAN_LIB=/usr/lib/gcc/x86_64-pc-linux-gnu/4.3.4/libgfortran.a

default: ${EXE}

%.exe : %.o $(EXTRA_OBJS)
	$(CXX) -o $@ $(LDFLAGS) $(CXXFLAGS) $< $(EXTRA_OBJS) -L$(LIB_PATH) $(LIBS) $(GFORTRAN_LIB)

clean:
	@rm -f *~ ${EXE}