#include <iostream>
#include <cassert>
-
-std::ofstream logfile ("00-manipulate_petsc_scalar/output");
-
// These tests were used while building extensions for
// PetscSclar=complex. They check that std::complex<double> and
// PetscScalar complex are compatible without explicitly casting.
//
-// These tests are archaic and produce unusual verbose output for
-// historical reasons...
+// These tests look archaic - they are the original tests used by the
+// authors to figure out how PetscScalar complex works...
+
+// Multiply a PETSc complex number by a complex number.
+void multiply_petsc_complex_by_a_number ()
+{
+ deallog << "Check PetscScalar multiply" << std::endl;
+
+ // make alpha
+ const PetscScalar alpha = 1.0 + 2.0*PETSC_i;
+
+ // make beta
+ const PetscScalar beta = 2.0 + 1.0*PETSC_i;
+
+ // multiply them together
+ const PetscScalar gamma = alpha * beta;
+
+ // Check real access
+ Assert (PetscRealPart (gamma)==0., ExcInternalError());
+ Assert (PetscImaginaryPart (gamma) ==5., ExcInternalError());
+
+ // This is legal too!
+ Assert (gamma==std::complex<double> (0.,5), ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
// Divide a PETSc complex number by a real number.
void divide_petsc_complex_by_a_number ()
{
+ deallog << "Check PetscScalar division" << std::endl;
+
+ // make alpha
PetscScalar alpha = 1.0 + 2.0*PETSC_i;
+
+ // make beta <- alpha/2
const PetscScalar beta = alpha/2.;
- logfile << " divide a petsc complex by 2.: "
- << PetscRealPart (beta)
- << "+" << PetscImaginaryPart (beta)
- << "i"
- << " should be 0.5+1i"
- << std::endl;
+ // Check real access
+ Assert (PetscRealPart (alpha)==1., ExcInternalError());
+ Assert (PetscRealPart (beta) ==0.5, ExcInternalError());
+ // operate on alpha (now alpha==beta)
alpha /= 2.;
- logfile << " divide-equals a petsc complex by 2.: "
- << PetscRealPart (alpha)
- << "+" << PetscImaginaryPart (alpha)
- << "i"
- << " should be 0.5+1i"
- << std::endl;
-
+ // Check complex access
+ Assert (PetscImaginaryPart (alpha)==1., ExcInternalError());
+ Assert (PetscImaginaryPart (beta) ==1., ExcInternalError());
+
+ // This is legal too!
Assert (alpha==beta, ExcInternalError());
+
+ deallog << "OK" << std::endl;
}
// Initialize a std::complex number from an PETSc complex number
void make_std_complex_from_petsc_complex ()
{
+ deallog << "Check std initialised from PetscScalar" << std::endl;
const PetscScalar alpha = 1.0 + 2.0*PETSC_i;
const std::complex<double> beta = alpha;
- logfile << " make std::complex from petsc complex: "
- << std::real (beta)
- << "+" << std::imag (beta)
- << "i"
- << std::endl;
-
+ // These should be the same of course
Assert (alpha==beta, ExcInternalError());
+
+ deallog << "OK" << std::endl;
}
// Initialize a PETSc complex number from an std::complex number
void make_petsc_complex_from_std_complex ()
{
+ deallog << "Check PetscScalar initialised from std" << std::endl;
const std::complex<double> beta (1,2);
const PetscScalar alpha = beta;
-
- logfile << " make petsc complex from std::complex: "
- << PetscRealPart (alpha)
- << "+" << PetscImaginaryPart (alpha)
- << "i"
- << std::endl;
-
+
+ // These should be the same of course
Assert (alpha==beta, ExcInternalError());
+
+ deallog << "OK" << std::endl;
}
// Initialize a PETSc complex number directly.
void make_petsc_complex ()
{
+ deallog << "Check a nonzero PetscScalar" << std::endl;
+
+ // init
const PetscScalar alpha = 1.0 + 2.0*PETSC_i;
-
- logfile << " make petsc complex: "
- << PetscRealPart (alpha)
- << "+" << PetscImaginaryPart (alpha)
- << "i"
- << " should be 1+2i"
- << std::endl;
+ // Test if PetscScalar is an std::complex.
Assert (alpha==std::complex<double> (1.,2.), ExcInternalError());
+
+ deallog << "OK" << std::endl;
}
// initialised to 0+0i.
void init_petsc_complex ()
{
+ deallog << "Check PetscScalar initialisation" << std::endl;
+
+ // Initialise (no argument) to zero.
const PetscScalar alpha;
+ Assert (alpha==std::complex<double> (0.,0.), ExcInternalError());
- logfile << " make petsc complex: "
- << PetscRealPart (alpha)
- << "+" << PetscImaginaryPart (alpha)
- << "i"
- << " should be 0+0i"
- << std::endl;
-
+ // Initialise (real argument) to zero.
const PetscScalar beta = 0.;
-
- logfile << " make petsc complex: "
- << PetscRealPart (beta)
- << "+" << PetscImaginaryPart (beta)
- << "i"
- << " should be 0+0i"
- << std::endl;
+ Assert (beta==alpha, ExcInternalError());
- Assert (alpha==beta, ExcInternalError());
+ // Initialise (real+complex argument) to zero.
+ const PetscScalar gamma = 0. + 0.*PETSC_i;
+ Assert (gamma==beta, ExcInternalError());
+
+ // If alpha==beta==gamma, then:
+ Assert (alpha==gamma, ExcInternalError());
+
+ deallog << "OK" << std::endl;
}
int main (int argc, char **argv)
{
+ std::ofstream logfile ("00/output");
dealii::deallog.attach (logfile);
- dealii::deallog.depth_console (1);
+ dealii::deallog.depth_console (0);
+ deallog.threshold_double(1.e-10);
try
{
PetscInitialize (&argc, &argv, (char*) 0, (char*) 0);
{
- // some initialisations makes and so on...
+ // check initialisation (zero and nonzero)
init_petsc_complex ();
make_petsc_complex ();
- // first do conversions explicitly
+ // check initialisation from std::complex (and vice versa)
make_petsc_complex_from_std_complex ();
make_std_complex_from_petsc_complex ();
-
- // then try to use operators to do the same thing
- logfile << " make petsc complex from std::complex: ";
- const std::complex<double> number1 (1,2);
- PetscScalar alpha1;
- alpha1 += number1;
-
- logfile << PetscRealPart (alpha1)
- << "+" << PetscImaginaryPart (alpha1)
- << "i"
- << " should be 1+2i"
- << std::endl;
-
- // then try to use operators to do the same thing - except it
- // seems to work without(!)
- logfile << " make std::complex from petsc complex: ";
- const PetscScalar alpha2 = 1.0 + 2.0*PETSC_i;
- std::complex<double> number2;
- number2 += alpha2;
-
- logfile << std::real (number2)
- << "+" << std::imag (number2)
- << "i"
- << " should be 1+2i"
- << std::endl;
-
- logfile << " make std::complex from petsc complex: ";
- const PetscScalar alpha3 = 1.0 - 2.0*PETSC_i;
- number2 += alpha3;
-
- logfile << std::real (number2)
- << "+" << std::imag (number2)
- << "i"
- << " should be 2+0i"
- << std::endl;
-
- // Let's try some other things...
+
+ // Check some operators
divide_petsc_complex_by_a_number ();
-
+ multiply_petsc_complex_by_a_number ();
}
PetscFinalize ();
}