[Trilinos-Users] SerialDenseSolver for std::complex<double>

Hoemmen, Mark mhoemme at sandia.gov
Sun Apr 8 18:01:57 MDT 2012


On Sun 08 Apr 2012 at 19:14:29 +0200, Holger Brandsmeier wrote:
> does SerialDenseSolver work for std::complex<double>? If I try to
> compile it I get the following error:
>
> In instantiation of member function 'Teuchos::SerialDenseSolver<int,
> ?std::complex<double> >::applyRefinement'
> 
> Teuchos_SerialDenseSolver.hpp:653:48: error: cannot initialize a
> parameter of type 'double *' with an rvalue of type 'int *'
> ? ? ? ? ? ? ?&FERR_[0], &BERR_[0], &WORK_[0], &IWORK_[0], &INFO_);
>
> (where he complains about &IWORK_[0])

I think I see what the problem is.  SerialDenseSolver inherits from
Teuchos::LAPACK.  The generic version of LAPACK (for any ScalarType
template parameter) has a GERFS() routine with an "int* IWORK"
argument.  However, the specialization of Teuchos::LAPACK for
ScalarType = std::complex<double> has a GERFS routine whose
second-to-last argument is "double* RWORK".  In fact, this correctly
imitates the arguments of the actual LAPACK routines DGERFS
resp. ZGERFS:

http://www.netlib.org/lapack/double/dgerfs.f
http://www.netlib.org/lapack/complex16/zgerfs.f

However, this means that SerialDenseSolver::applyRefinement() needs a
different specialization for real vs. complex ScalarType types.

> Similarly later I get the error:
> 
> Teuchos_SerialDenseSolver.hpp:835:65: error: cannot initialize a
> parameter of type 'double *' with an rvalue of type 'int *'
> ?this->GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0],
> &IWORK_[0], &INFO_);

This looks like the same problem as with GERFS above.

> Is this a bug, is this an error somewhere on my installation or is
> this not yet implemented?

This is definitely a bug in SerialDenseSolver.

> If this is not yet implemented, can you give me some hints how to
> implement it (I will send you the implementation if I succeed). I
> assume that it should be possible because lapack supports complex
> numbers. I am not very familiar to the specifics of lapack though,
> I always think lapack is too cryptic.

I think some LAPACK developers find LAPACK cryptic on occasion ;-)

We very much appreciate your offer to implement this, and we would
welcome your contribution.  Probably the easiest thing for now would
be to specialize the applyRefinement() and
reciprocalConditionEstimate() methods for ScalarType=std::complex<T>.
It's a little bit tricky because IWORK_ is a member datum of
SerialDenseSolver, but I think these two methods are the only methods
which use IWORK_ at all, and IWORK_ is just a workspace array and not
actually used for output.  Thus, I think the easiest thing to do for
now would be to declare and allocate an array just for these methods.
After the generic declaration of the routine, you could add the code
shown below (I haven't tested it but I plan on adding it once I get
the chance). 

You're welcome to submit a patch once you've tested it out, in case
there are any bugs in my code.  Thanks for your offer to help out!

mfh

// Make sure not to use std::complex unless building with complex
// arithmetic support.
#ifdef HAVE_TEUCHOS_COMPLEX

// Specialization for complex ScalarType types.
template<typename OrdinalType, typename T>
int 
SerialDenseSolver<OrdinalType, std::complex<T> >::applyRefinement()
{
  TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
    "SerialDenseSolver<T>::applyRefinement: Must have an existing solution!");
  TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
    "SerialDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");

  OrdinalType NRHS = RHS_->numCols();
  FERR_.resize( NRHS );
  BERR_.resize( NRHS );
  allocateWORK();
  allocateIWORK();
  
  INFO_ = 0;
  std::vector<T> RWORK (N_);
  this->GERFS(ETranspChar[TRANS_], N_, NRHS, A_, LDA_, AF_, LDAF_, &IPIV_[0],
	      RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(), 
              &FERR_[0], &BERR_[0], &WORK_[0], (N_ == 0 ? NULL : &RWORK[0]),
              &INFO_);
  solutionErrorsEstimated_ = true;
  reciprocalConditionEstimated_ = true;
  solutionRefined_ = true;
  
  return(INFO_);
}

// Specialization for complex ScalarType types.
template<typename OrdinalType, typename T>
int 
SerialDenseSolver<OrdinalType, std::complex<T> >::
reciprocalConditionEstimate (MagnitudeType& Value)
{
  if (reciprocalConditionEstimated()) {
    Value = RCOND_;
    return(0); // Already computed, just return it.
  }

  if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();

  int ierr = 0;
  if (!factored()) ierr = factor(); // Need matrix factored.
  if (ierr!=0) return(ierr);

  allocateWORK();
  allocateIWORK();

  std::vector<T> RWORK (N_);
  INFO_ = 0;
  // We will assume a one-norm condition number
  this->GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], 
               (N_ == 0 ? NULL : &RWORK[0]), &INFO_);
  reciprocalConditionEstimated_ = true;
  Value = RCOND_;

  return(INFO_);
}
#endif // HAVE_TEUCHOS_COMPLEX






More information about the Trilinos-Users mailing list