[Trilinos-Users] Epetra_CrsMatrix multiplication with a complex vector

[Nico Schlömer]

Hi David,

when dealing with complex-valued vectors and matrices, Tpetra is your
friend: As opposed to Epetra, you can stuff basically any scalar type into
it, e.g, complex<double>.

As of now I don't believe it's possible to have two Tpetra objects
interact which don't exactly share their scalar types, e.g., double vs.
complex<double>. -- Trilinos devs, am I up-to-date here?

As a workaround, you could just create your Hamiltonian as a
complex<double> matrix, which admittedly wastes the memory that all the
imaginary part zeros occupy. Anyway, I attached a minimal example for
complex<double> Tpetra matrix-vector multiplication.

Cheers,
Nico



On Sun, 24 Jan 2010 13:23:41 +0100, "David Hochstuhl"
<Davidhochstuhl at web.de> wrote:
> Hello,
> 
> at first I want to thank the trilinos team for all those impressing
> program pakets.
> 
> Up to now, I only use the anasazi package to do Configuration
Interaction,
> i.e. to find the eigenvalues/-vectors of a quantum chemistry Hamiltonian
> matrix. 
> The solutions shall serve as initial states for a subsequent time
> evolution, in which an electromagnetic field is applied to the system.
> 
> Now, to solve this task with Runge-Kutta, one needs to apply the real
> Hamiltonian matrix to a complex vector.
> 
> So, is there any straightforward way to do this efficiently within
> trilinos ?
> 
> Thanks in advance,
> David
> 
> ______________________________________________________________________
> Haiti-Nothilfe! Helfen Sie per SMS: Sende UIHAITI an die Nummer 81190.
> Von 5 Euro je SMS (zzgl. SMS-Gebühr) gehen 4,83 Euro an UNICEF.
> 
> 
> 
> _______________________________________________
> Trilinos-Users mailing list
> Trilinos-Users at software.sandia.gov
> http://software.sandia.gov/mailman/listinfo/trilinos-users
-------------- next part --------------
#include <Teuchos_RCP.hpp>
#include <Teuchos_DefaultComm.hpp>
#include <Tpetra_Map.hpp>
#include <Tpetra_CrsMatrix.hpp>

using namespace std;

// date type of the matrix
typedef complex<double> matrixType;

typedef Tpetra::Vector<complex<double> > ComplexVector;
typedef Tpetra::CrsMatrix<matrixType>    Matrix;

int main () {

    // Initialize MPI
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
#endif

  // Create a communicator for Tpetra objects
  const Teuchos::RCP<const Teuchos::Comm<int> > Comm = Teuchos::DefaultComm<int>::getComm();

  int NumComplexUnknowns = 10;
  
  // create map
  Teuchos::RCP<const Tpetra::Map<int> >  ComplexMap
              = Teuchos::rcp(new Tpetra::Map<int>(NumComplexUnknowns, 0, Comm));
  
  Tpetra::Vector<complex<double> > z = ComplexVector(ComplexMap);
  z.randomize();

  // take a look at the vector
  Teuchos::ArrayRCP<const complex<double> > zView = z.get1dView();
  for (unsigned int k=0;k<ComplexMap->getGlobalNumElements(); k++ )
     cout << "z[" << k << "] = " << zView[k] << endl;

  // create the matrix as diagonal matrixType
  // ( 1, 1/2, 1/3, ... )
  Matrix H(ComplexMap,ComplexMap,0);
  Teuchos::Array<int>        cols(1);
  Teuchos::Array<matrixType> values(1);
  for (unsigned int row=0; row<ComplexMap->getNodeNumElements(); row++) {
    cols[0] = ComplexMap->getGlobalElement(row);
    values[0] = 1.0 / (cols[0]+1);
    H.insertLocalValues( row, cols, values );
  }
  H.fillComplete();

  // compute H*z and store the result in res
  Tpetra::Vector<complex<double> > res = ComplexVector(ComplexMap);
  H.multiply( z, res, Teuchos::NO_TRANS );

    // take a look at the result
  Teuchos::ArrayRCP<const complex<double> > resView = res.get1dView();
  for (unsigned int k=0;k<ComplexMap->getGlobalNumElements(); k++ )
     cout << "res[" << k << "] = " << resView[k] << 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=main.exe

default: $(EXE)