xref: /dports/math/xlife++/xlifepp-sources-v2.0.1-2018-05-09/src/term/TermMatrix.cpp

/*
XLiFE++ is an extended library of finite elements written in C++
    Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*!
  \file TermMatrix.cpp
  \author E. Lunéville
  \authors D. Martin, E. Lunéville
  \since 03 apr 2012
  \date 20 mar 2013

  \brief Implementation of xlifepp::TermMatrix class members
*/

#include "TermMatrix.hpp"

namespace xlifepp
{

//===============================================================================
//member functions of TermMatrix class
//===============================================================================
//constructors
TermMatrix::TermMatrix(const string_t& na)
{
  termType_ = _termMatrix;
  name_ = na;
  entries_p = 0;
  scalar_entries_p=0;
  rhs_matrix_p=0;
  constraints_u_p = 0;
  constraints_v_p = 0;
}

TermMatrix::TermMatrix(const BilinearForm& blf, const string_t& na, bool noass)
{
  initFromBlf(blf, na, noass);
  constraints_u_p = 0;
  constraints_v_p = 0;
}

//! constructor with no essential boundary conditions
TermMatrix::TermMatrix(const BilinearForm& blf, const string_t& na)
{
  std::vector<TermOption> opts;
  build(blf,0,0,ReductionMethod(_noReduction),opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, TermOption op1, const string_t& na)
{
  std::vector<TermOption> opts(1,op1);
  build(blf,0,0,ReductionMethod(_noReduction),opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, TermOption op1, TermOption op2, const string_t& na)
{
  std::vector<TermOption> opts(2,op1); opts[1]=op2;
  build(blf,0,0,ReductionMethod(_noReduction),opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, TermOption op1, TermOption op2, TermOption op3, const string_t& na)
{
  std::vector<TermOption> opts(3,op1); opts[1]=op2; opts[2]=op3;
  build(blf,0,0,ReductionMethod(_noReduction),opts,na);
}

//! constructor with one essential boundary conditions
TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ec, const string_t& na)
{
  std::vector<TermOption> opts;
  build(blf,&ec,0,ReductionMethod(_pseudoReduction),opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ec, TermOption op1, const string_t& na)
{
  std::vector<TermOption> opts(1,op1);
  build(blf,&ec,0,ReductionMethod(_pseudoReduction),opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ec, TermOption op1, TermOption op2, const string_t& na)
{
  std::vector<TermOption> opts(2,op1); opts[1]=op2;
  build(blf,&ec,0,ReductionMethod(_pseudoReduction),opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ec, TermOption op1, TermOption op2, TermOption op3, const string_t& na)
{
  std::vector<TermOption> opts(3,op1); opts[1]=op2; opts[2]=op3;
  build(blf,&ec,0,ReductionMethod(_pseudoReduction),opts,na);
}

//! constructor with essential boundary condition and explicit reduction method
TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ec, const ReductionMethod& rm, const string_t& na)
{
  std::vector<TermOption> opts;
  build(blf,&ec,0,rm,opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ec, const ReductionMethod& rm, TermOption op1, const string_t& na)
{
  std::vector<TermOption> opts(1,op1);
  build(blf,&ec,0,rm,opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ec, const ReductionMethod& rm, TermOption op1, TermOption op2, const string_t& na)
{
  std::vector<TermOption> opts(2,op1); opts[1]=op2;
  build(blf,&ec,0,rm,opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ec, const ReductionMethod& rm, TermOption op1, TermOption op2, TermOption op3, const string_t& na)
{
  std::vector<TermOption> opts(3,op1); opts[1]=op2; opts[2]=op3;
  build(blf,&ec,0,rm,opts,na);
}

//! constructor with two essential boundary conditions
TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ecu, const EssentialConditions& ecv, const string_t& na)
{
  std::vector<TermOption> opts;
  build(blf,&ecu,&ecv,ReductionMethod(_pseudoReduction),opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ecu, const EssentialConditions& ecv, TermOption op1, const string_t& na)
{
  std::vector<TermOption> opts(1,op1);
  build(blf,&ecu,&ecv,ReductionMethod(_pseudoReduction),opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ecu, const EssentialConditions& ecv, TermOption op1, TermOption op2, const string_t& na)
{
  std::vector<TermOption> opts(2,op1); opts[1]=op2;
  build(blf,&ecu,&ecv,ReductionMethod(_pseudoReduction),opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ecu, const EssentialConditions& ecv, TermOption op1, TermOption op2, TermOption op3, const string_t& na)
{
  std::vector<TermOption> opts(3,op1); opts[1]=op2; opts[2]=op3;
  build(blf,&ecu,&ecv,ReductionMethod(_pseudoReduction),opts,na);
}

//! constructor with two essential boundary conditions and explicit reduction method
TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ecu, const EssentialConditions& ecv, const ReductionMethod& rm, const string_t& na)
{
  std::vector<TermOption> opts;
  build(blf,&ecu,&ecv,rm,opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ecu, const EssentialConditions& ecv, const ReductionMethod& rm, TermOption op1, const string_t& na)
{
  std::vector<TermOption> opts(1,op1);
  build(blf,&ecu,&ecv,rm,opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ecu, const EssentialConditions& ecv, const ReductionMethod& rm,TermOption op1, TermOption op2, const string_t& na)
{
  std::vector<TermOption> opts(2,op1); opts[1]=op2;
  build(blf,&ecu,&ecv,rm,opts,na);
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ecu, const EssentialConditions& ecv, const ReductionMethod& rm, TermOption op1, TermOption op2, TermOption op3, const string_t& na)
{
  std::vector<TermOption> opts(3,op1); opts[1]=op2; opts[2]=op3;
  build(blf,&ecu,&ecv,rm,opts,na);
}


//! effective constructor
void TermMatrix::build(const BilinearForm& blf, const EssentialConditions* ecu, const EssentialConditions* ecv, const ReductionMethod& rm,
                       const std::vector<TermOption>& opts, const string_t& na)
{
  std::vector<TermOption>::const_iterator ito=opts.begin();
  computingInfo_.isComputed=false;
  computingInfo_.noAssembly=false;
  bool toCompute=true, toAssemble=true;
  StorageType st=_noStorage;  AccessType ac=_noAccess;
  SymType sym=_undefSymmetry;
  computingInfo_.reductionMethod = rm;
  for (; ito!=opts.end(); ++ito)
  {
    switch (*ito)
    {
      case _compute    : break;
      case _notCompute : toCompute=false; break;
      case _assembled   : break;
      case _unassembled : toAssemble=false; break;
      case _csRowStorage : st=_cs; ac=_row; break;
      case _csColStorage : st=_cs; ac=_col; break;
      case _csDualStorage: st=_cs; ac=_dual; break;
      case _csSymStorage : st=_cs; ac=_sym; break;
      case _denseRowStorage : st=_dense; ac=_row; break;
      case _denseColStorage : st=_dense; ac=_col; break;
      case _denseDualStorage: st=_dense; ac=_dual; break;
      case _skylineSymStorage : st=_skyline; ac=_sym; break;
      case _skylineDualStorage: st=_skyline; ac=_dual; break;
      case _nonSymmetricMatrix : sym=_noSymmetry; break;
      case _symmetricMatrix    : sym=_symmetric; break;
      case _selfAdjointMatrix  : sym=_selfAdjoint; break;
      case _skewSymmetricMatrix: sym=_skewSymmetric; break;
      case _skewAdjointMatrix  : sym=_skewAdjoint; break;
      case _pseudoReductionMethod : computingInfo_.reductionMethod.method=_pseudoReduction; break;
      case _realReductionMethod   : computingInfo_.reductionMethod.method=_realReduction; break;
      case _penalizationReductionMethod : computingInfo_.reductionMethod.method=_penalizationReduction; break;
      default:
        where("TermMatrix::build(...)");
        error("undef_option");
    }
  }

  initFromBlf(blf, na, !toAssemble);
  constraints_u_p = 0;
  if (ecu!=0) constraints_u_p = new SetOfConstraints(buildConstraints(*ecu));
  constraints_v_p = constraints_u_p;
  if (ecv!=0 && ecv!=ecu) constraints_v_p = new SetOfConstraints(buildConstraints(*ecv));
  if (sym!=_undefSymmetry)
  {
    where("TermMatrix::build(...)");
    error("symmetry_not_handled", words("symmetry", sym));
  }
  if (st!=_noStorage) setStorage(st,ac);
  if (toCompute) compute();
}

TermMatrix::TermMatrix(const BilinearForm& blf,  const EssentialConditions& ec, const string_t& na, bool noass)
{
  initFromBlf(blf, na, noass);
  constraints_u_p = new SetOfConstraints(buildConstraints(ec));
  constraints_v_p = constraints_u_p;
  computingInfo_.reductionMethod.method=_pseudoReduction;
}

TermMatrix::TermMatrix(const BilinearForm& blf, const EssentialConditions& ecu, const EssentialConditions& ecv, const string_t& na, bool noass)
{
  initFromBlf(blf, na, noass);
  constraints_u_p = new SetOfConstraints(buildConstraints(ecu));
  constraints_v_p = new SetOfConstraints(buildConstraints(ecv));
  computingInfo_.reductionMethod.method=_pseudoReduction;
}

TermMatrix::TermMatrix(const BilinearForm& blf, const SetOfConstraints& soc, const string_t& na, bool noass)
{
  initFromBlf(blf, na, noass);
  constraints_u_p = new SetOfConstraints(soc);   //full copy of soc !!!
  constraints_v_p = constraints_u_p;
  computingInfo_.reductionMethod.method=_pseudoReduction;
}

TermMatrix::TermMatrix(const BilinearForm& blf, const SetOfConstraints& socu, const SetOfConstraints& socv,
                       const string_t& na, bool noass)
{
  initFromBlf(blf, na, noass);
  constraints_u_p = new SetOfConstraints(socu);
  constraints_v_p = new SetOfConstraints(socv);
  computingInfo_.reductionMethod.method=_pseudoReduction;
}

//create suTermMatrix's
void TermMatrix::initFromBlf(const BilinearForm& blf, const string_t& na, bool noass)
{
  bilinForm_ = blf; //copy the linearform
  computingInfo_.noAssembly = noass;
  termType_ = _termMatrix;
  name_ = na;
  entries_p = 0;
  scalar_entries_p=0;
  rhs_matrix_p=0;
  for(it_mublc it = bilinForm_.begin(); it != bilinForm_.end(); it++)  //travel the map of single unknown bilinear forms
  {
    string_t nam = name_ + "_" + it->second.up()->name() + "_" + it->second.vp()->name();
    // string_t nam = name_ + "_sub";
    suTerms_[it->first] = new SuTermMatrix(&(it->second), nam, computingInfo_); //construct SuTermMatrix objects
  }
}

//constructor from a linear combination of TermVector's
TermMatrix::TermMatrix(const LcTerm<TermMatrix>& lc, const string_t& na)
{
  trace_p->push("TermMatrix::TermMatrix(LcTerm)");
  entries_p = 0;
  scalar_entries_p=0;
  rhs_matrix_p=0;
  constraints_u_p = 0;
  constraints_v_p = 0;
  termType_ = _termMatrix;
  compute(lc, na);
  trace_p->pop();
}

// destructor
TermMatrix::~TermMatrix()
{
  for(it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++) delete it->second;
  if(suTerms_.size()>1 && entries_p != 0) delete entries_p;
  if(scalar_entries_p != entries_p && scalar_entries_p != 0) delete scalar_entries_p;
  if(rhs_matrix_p != 0) delete rhs_matrix_p;
  if(constraints_u_p!=0) delete constraints_u_p;
  if(constraints_v_p!=0 && constraints_v_p!= constraints_u_p) delete constraints_v_p;
}

//copy constructor
TermMatrix::TermMatrix(const TermMatrix& mat, const string_t& na)
{
  computingInfo_.noAssembly = false;
  termType_ = _termMatrix;
  name_ = na;
  entries_p = 0;
  scalar_entries_p=0;
  rhs_matrix_p=0;
  constraints_u_p = 0;
  constraints_v_p = 0;
  copy(mat);
  if(na != "") name_ = na;
}

//! constructor from a SuTermMatrix, create a one block TermMatrix
TermMatrix::TermMatrix(const SuTermMatrix& sutm, const string_t& na)
{
  trace_p->push("TermMatrix::TermMatrix(SuTermMatrix)");
  computingInfo_.noAssembly = false;
  termType_ = _termMatrix;
  name_ = na;
  entries_p = 0;
  scalar_entries_p=0;
  rhs_matrix_p=0;
  constraints_u_p = 0;
  constraints_v_p = 0;
  suTerms_[uvPair(sutm.up(),sutm.vp())]=new SuTermMatrix(sutm);  //full copy
  computingInfo_.isComputed=sutm.computed();
  trace_p->pop();
}

TermMatrix::TermMatrix(SuTermMatrix* sutm, const string_t& na)
{
  trace_p->push("TermMatrix::TermMatrix(SuTermMatrix*)");
  computingInfo_.noAssembly = false;
  termType_ = _termMatrix;
  name_ = na;
  entries_p = 0;
  scalar_entries_p=0;
  rhs_matrix_p=0;
  constraints_u_p = 0;
  constraints_v_p = 0;
  suTerms_[uvPair(sutm->up(),sutm->vp())]=sutm;  //pointer copy
  computingInfo_.isComputed=sutm->computed();
  trace_p->pop();
}

// insert a SuTermMatrix  with full copy
void TermMatrix::insert(const SuTermMatrix& sutm)
{
  trace_p->push("TermMatrix::insert(SuTermMatrix)");
  suTerms_[uvPair(sutm.up(),sutm.vp())]=new SuTermMatrix(sutm);  //full copy
  if(!sutm.computed()) computingInfo_.isComputed= false;
  trace_p->pop();
}

// insert a SuTermMatrix with pointer copy
void TermMatrix::insert(SuTermMatrix* sutm)
{
  trace_p->push("TermMatrix::insert(SuTermMatrix*)");
  suTerms_[uvPair(sutm->up(),sutm->vp())]=sutm;       //pointer copy
  if(!sutm->computed()) computingInfo_.isComputed= false;
  trace_p->pop();
}


/*! constructor of a diagonal TermMatrix from a TermVector
    if TermVector is a multiple unknowns vector, the matrix will be a multiple unknowns TermMatrix, each block being diagonal !
*/
TermMatrix::TermMatrix(const TermVector& tv, const string_t& na)
{
  trace_p->push("TermMatrix::TermMatrix(TermVector)");
  computingInfo_.noAssembly = false;
  termType_ = _termMatrix;
  name_ = na;
  entries_p = 0;
  scalar_entries_p=0;
  rhs_matrix_p=0;
  constraints_u_p = 0;
  constraints_v_p = 0;
  for(cit_mustv it = tv.begin(); it != tv.end(); it++)
    {
      const Unknown* u=it->first;
      suTerms_[uvPair(u,u)]=new SuTermMatrix(*it->second,na+"_"+u->name()+"_"+u->name());
    }
  computingInfo_.isComputed=true;
  trace_p->pop();
}

/*! constructor of matrix opker(xi,yj) - Lagrange interpolation -
    opker is evaluated on dof coordinates linked to Unknown and GeomDomains
     ux, domx is associated to matrix rows and  uy, domy to matrix columns
     NOTE : do not use this constructor if kernel is singular on some couple of points
*/
TermMatrix::TermMatrix(const Unknown& ux, const GeomDomain& domx,
                       const Unknown& uy, const GeomDomain& domy, const OperatorOnKernel& opk, const string_t& na)
{
  trace_p->push("TermMatrix::TermMatrix(Unknown, GeomDomain, Unknown, GeomDomain, OperatorOnKernel)");
  computingInfo_.noAssembly = false;
  termType_ = _termMatrix;
  name_ = na;
  entries_p = 0;
  scalar_entries_p=0;
  rhs_matrix_p=0;
  constraints_u_p = 0;
  constraints_v_p = 0;
  suTerms_[uvPair(&uy,&ux)]=new SuTermMatrix(ux, domx, uy, domy, opk, na+"_"+uy.name()+"_"+ux.name());
  computingInfo_.isComputed=true;
  trace_p->pop();
}

/*! construct special matrix from a structure given by a TermMatrix
  up to now, only IdMatrix is available
*/
TermMatrix::TermMatrix(const TermMatrix& A, SpecialMatrix sm, const string_t& na)
{
  trace_p->push("TermMatrix(TermMatrix, SpecialMatrix)");
  if (sm!=_idMatrix) error("matrix_type_not_handled", words("matrix", sm));
  computingInfo_.noAssembly = false;
  termType_ = _termMatrix;
  if (na=="") name_="Id"; else name_ = na;
  entries_p = 0;
  scalar_entries_p=0;
  rhs_matrix_p=0;
  constraints_u_p = 0;
  constraints_v_p = 0;
  for (cit_mustm it = A.begin(); it != A.end(); it++)  //travel the map of single unknown linear forms
  {
    uvPair uv = it->first;
    if(uv.first==uv.second || uv.first->dual_p()==uv.second)  suTerms_[uv] = new SuTermMatrix(*it->second, sm);
  }
  computingInfo_ .isComputed=true;;
  if (A.scalar_entries_p != 0)
  {
    cdofs_c= A.cdofs_c;
    cdofs_r= A.cdofs_r;
  }
  trace_p->pop();
}

//assign operator
TermMatrix& TermMatrix::operator=(const TermMatrix& mat)
{
  if(this==&mat) return *this;
  clear();
  copy(mat);
  return *this;
}

// specify the storage of SuTermMatrix's of TermMatrix
// induce a storage conversion if terms are computed
void TermMatrix::setStorage(StorageType st, AccessType at)
{
  trace_p->push("TermMatrix::setStorage(...)");
  computingInfo_.storageType = st;
  computingInfo_.storageAccess = at;
  for(it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)
    it->second->setStorage(st, at);
  if(entries_p!=0)
    {
      if(entries_p->storageType()!=st || entries_p->accessType()!=at)
        {
          std::vector<std::vector<number_t> > colIndices = entries_p->storagep()->scalarColIndices();  //retry column indices
          MatrixStorage* mstorage = buildStorage(st, at, entries_p->nbOfRows(), entries_p->nbOfCols(), colIndices);
          entries_p->toStorage(mstorage);
        }
    }
  if(scalar_entries_p!=0)
    {
      if(scalar_entries_p->storageType()!=st || scalar_entries_p->accessType()!=at)
        {
          std::vector<std::vector<number_t> > colIndices = scalar_entries_p->storagep()->scalarColIndices();  //retry column indices
          MatrixStorage* mstorage = buildStorage(st, at, scalar_entries_p->nbOfRows(), scalar_entries_p->nbOfCols(), colIndices);
          scalar_entries_p->toStorage(mstorage);
        }
    }
  trace_p->pop();
}

// copy a TermMatrix into current TermMatrix
void TermMatrix::copy(const TermMatrix& mat)
{
  trace_p->push("TermMatrix::copy(TermMatrix)");
  bilinForm_ = mat.bilinForm_;
  //if(name_=="" || name_=="?") name_ = mat.name_ + "_copy";
  if(name_=="" || name_=="?") name_ = mat.name_ + "@";
  termType_ = _termMatrix;
  for(cit_mustm it = mat.begin(); it != mat.end(); it++)  //travel the map of single unknown linear forms
    {
      uvPair uv = it->first;
      SuTermMatrix* sutm = it->second;
      suTerms_[uv] = new SuTermMatrix(*sutm);
      suTerms_[uv]->sublfp() = bilinForm_.subLfp(uv);   //maintain links between bilinearform and sublinearforms
    }
  computingInfo_ = mat.computingInfo_;
  if(mat.entries_p != 0) entries_p = new MatrixEntry(*mat.entries_p);
  if(mat.scalar_entries_p != 0)
    {
      if(mat.scalar_entries_p!=mat.entries_p) scalar_entries_p = new MatrixEntry(*mat.scalar_entries_p);
      else scalar_entries_p = entries_p;
      cdofs_c= mat.cdofs_c;
      cdofs_r= mat.cdofs_r;
    }
  if(mat.rhs_matrix_p != 0) rhs_matrix_p = new MatrixEntry(*mat.rhs_matrix_p);
  if(mat.constraints_u_p !=0) constraints_u_p = new SetOfConstraints(*mat.constraints_u_p);
  if(mat.constraints_v_p !=0)
    {
      if(mat.constraints_v_p == mat.constraints_u_p) constraints_v_p = constraints_u_p;
      else constraints_v_p = new SetOfConstraints(*mat.constraints_v_p);
    }
  trace_p->pop();
}

//deallocate memory used by a TermMatrix
void TermMatrix::clear()
{
  trace_p->push("TermMatrix::clear");
  if(entries_p != 0) delete entries_p;
  if(scalar_entries_p!=entries_p && scalar_entries_p != 0) delete scalar_entries_p;
  if(rhs_matrix_p != 0) delete rhs_matrix_p;
  if(constraints_u_p!=0) delete constraints_u_p;
  if(constraints_v_p!=0 && constraints_v_p!= constraints_u_p) delete constraints_v_p;
  entries_p = 0;
  scalar_entries_p = 0;
  rhs_matrix_p =0;
  constraints_u_p = 0;
  constraints_v_p = 0;
  cdofs_r.clear();
  cdofs_c.clear();
  for(it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++) it->second->clear();
  //suTerms_.clear();
  computed() = false;
  trace_p->pop();
}

/*! compute a LcTerm (linear combination of TermMatrix's) to a TermMatrix
 all the TermMatrix's have to be computed before
 and the current TermMatrix is 'void'

 NOTE : this current version works only for "standard" representation (entries pointers) but NOT for scalar representation (scalar_entries pointers)
        it means that TermMatrix that have been converted using toScalar or toGlobal functions cannot be used in linear combinations.
        Such conversions occur when essential conditions have been applied to a TermMatrix
        To be improved in future !
*/
void TermMatrix::compute(const LcTerm<TermMatrix>& lc, const string_t& na)
{
  trace_p->push("TermMatrix::compute(LcTerm)");
  if(lc.size() == 0) error("is_void", "LcTerm");  //abnormal !!
  if(lc.size() == 1)
    {
      const TermMatrix* tM = dynamic_cast<const TermMatrix*>(lc.begin()->first);
      if(tM == 0) error("null_pointer", "TermMatrix");
      complex_t a = lc.begin()->second;
      for(cit_mustm it = tM->suTerms_.begin(); it != tM->suTerms_.end(); it++)
        {
          SuTermMatrix* sut = new SuTermMatrix(*it->second);   //copy SuTermMatrix
          sut->name() = tostring(a) + " * " + it->second->name();
          if(a != complex_t(1., 0)) *sut *= a;                    //product if a!=1
          suTerms_[it->first] = sut;
        }
      name_ = na;
      if(na == "")
        {
          if(a.imag() == 0) name_ = tostring(a.real());
          else name_ = tostring(a);
          name_ += " * ( " + tM->name() + " )" ;
        }
      computed() = true;
      trace_p->pop();
      return;
    }

  //case of a linear combination
  //create linear combinations of SuTermMatrix indexed by unknown pairs
  std::map<uvPair, LcTerm<SuTermMatrix> > lcs;
  std::map<uvPair, LcTerm<SuTermMatrix> >::iterator itm;
  name_ = na;
  string_t pl = "";
  for(LcTerm<TermMatrix>::const_iterator itt = lc.begin(); itt != lc.end(); itt++)
    {
      const TermMatrix* tM = dynamic_cast<const TermMatrix*>(itt->first);
      if(tM == 0) error("null_pointer", "TermMatrix");
      complex_t a = itt->second;
      for(cit_mustm it = tM->suTerms_.begin(); it != tM->suTerms_.end(); it++)
        {
          uvPair uv = it->first;
          SuTermMatrix* sut = it->second;
          itm = lcs.find(uv);
          if(itm == lcs.end()) lcs.insert(std::pair<uvPair, LcTerm<SuTermMatrix> >(uv, LcTerm<SuTermMatrix>(sut, a)));
          else
            {
              LcTerm<SuTermMatrix> lci = itm->second;
              lci.push_back(sut, a);
              itm->second = lci;
            }
        }
      if(na == "")
        {
          if(a.imag() == 0) name_ += pl + tostring(a.real()) ;
          else  name_ += pl + tostring(a) ;
          name_ += " * ( " + tM->name() + " )" ;
          pl = " + ";
        }
    }

  //compute linear combination of SuTermMatrix's
  for(itm = lcs.begin(); itm != lcs.end(); itm++)
    {
      SuTermMatrix* sut = new SuTermMatrix(itm->second);   //create SuTermMatrix
      suTerms_[itm->first] = sut;
    }

  computed() = true;
  trace_p->pop();
  return;
}

// assign operator from a linear combinations of TermVector's
TermMatrix& TermMatrix::operator=(const LcTerm<TermMatrix>& lc)
{
  trace_p->push("TermMatrix::operator =(LcTerm)");
  //prevent self adding case (for instance Tm=Tm+Tm2+...), Tm has to be copied before cleared !!!
  LcTerm<TermMatrix> lcc = lc;
  TermMatrix* clone = 0;
  for(LcTerm<TermMatrix>::iterator it = lcc.begin(); it != lcc.end(); it++)
    {
      if(it->first == const_cast<const TermMatrix*>(this))   //self adding
        {
          if(clone == 0) clone = new TermMatrix(*this);  //full copy !!!
          complex_t a = it->second;
          *it = std::pair<const TermMatrix*, complex_t>(clone, a); //modify linear combination by replacing current TermMatrix by its copy
        }
    }

  //clear current TermMatrix
  bilinForm_.clear();
  if(entries_p != 0)
    {
      delete entries_p;
      entries_p = 0;
    }
  for(it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++) delete it->second;  //delete SuTermMatrix's

  //compute linear combination
  compute(lcc);

  if(clone != 0) delete clone;
  trace_p->pop();
  return *this;
}

// value type (_real,_complex)
ValueType TermMatrix::valueType() const
{
  if(scalar_entries_p!=0) return scalar_entries_p->valueType_;
  if(entries_p!=0) return entries_p->valueType_;
  if(suTerms_.size() > 0)
    {
      for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)   //travel the map of single unknown matrices
        {
          if(it->second!=0 && it->second->valueType() == _complex) return _complex;
        }
    }
  return _real;
}

// is scalar if all SuTermMatrix's are scalar matrices
bool TermMatrix::isScalar() const
{
    if(suTerms_.size() > 0)
    {
      for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)   //travel the map of single unknown matrices
        {
          if(it->first.first->nbOfComponents()>1 || it->first.second->nbOfComponents()>1) return false;
        }
    }
  return true;
}

// TermMatrix and SutermMatrix's computed state = true or false
void TermMatrix::markAsComputed(bool c)
{
  computingInfo_.isComputed=c;
  for(it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)  //travel the map of single unknown matrices
    it->second->computed()=c;
}

// return factorization type if factorized else return _noFactorization
FactorizationType TermMatrix::factorization() const
{
  if(scalar_entries_p!=0) return scalar_entries_p->factorization();
  if(entries_p!=0) return entries_p->factorization();
  if(suTerms_.size()==1) return suTerms_.begin()->second->factorization();
  return _noFactorization;
}

// return set of row/col Unknowns (say v/u unknowns)
std::set<const Unknown*> TermMatrix::rowUnknowns() const
{
  std::set<const Unknown*> sun;
  for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)  //travel the map of single unknown matrices
    sun.insert(it->first.second);
  return sun;
}

std::set<const Unknown*> TermMatrix::colUnknowns() const
{
  std::set<const Unknown*> sun;
  for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)  //travel the map of single unknown matrices
    sun.insert(it->first.first);
  return sun;
}

//! list of involved unknown spaces
std::set<const Space*> TermMatrix::unknownSpaces() const
{
  std::set<const Space*> sps;
  for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)  //travel the map of single unknown matrices
    {
      sps.insert(it->first.first->space());
      sps.insert(it->first.second->space());
    }
  return sps;
}

/*!
  update the name of a TermMatrix and its SuTermMatrixs
  SuTermMatrix name is defined as follows: name(TermMatrix)_name(row_Unknown)_name(col_Unknown)
 */
void TermMatrix::name(const string_t& nam)
{
  name_=nam;
  for(it_mustm it = begin(); it != end(); it++)  //travel the map of single unknown linear forms
  {
    it->second->name() = nam + "_" + it->second->up()->name() + "_" + it->second->vp()->name();
  }
}

// return the actual pointer to entries (priority to scalar entry), return 0 is no allocated pointer
MatrixEntry* TermMatrix::actual_entries() const
{
  MatrixEntry* me = scalar_entries_p;
  if (me != 0) return me;
  if (suTerms_.size() == 1)
  {
    me = firstSut()->scalar_entries();
    if (me == 0) me = firstSut()->entries();
  }
  return me;
}

/*! return global symmetry property (_noSymmetry, _symmetric, _adjoint, _skewAdjoint, _skewSymmetric)
    that is same symmetry property of diagonal unknowns blocks it there are no non diagonal blocks
    if there are non diagonal blocks, their symmetry is not checked (to do later)
    NOTE : this function works only for allocated SCALAR entries
*/
SymType TermMatrix::symmetry() const
{
  SymType sy=_symmetric;
  if(suTerms_.size()==0) return sy;
  bool unset = true;  //symmetry property not set
  //identify symmetry property from diagonal block
  for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)  //travel the map of single unknown matrices
    {
      const Unknown* u=it->first.first, *v=it->first.second;
      SuTermMatrix* sut = it->second;
      if((u==v || u==v->dual_p()))   //block diagonal
        {
          if(sut->scalar_entries()!=0)  //priority to scalar entries
            {
              if(unset) {sy=sut->scalar_entries()->symmetry(); unset=false;}
              else if(sut->scalar_entries()->symmetry()!=sy) return _noSymmetry;
            }
          else if(sut->entries()!=0)  //matrix of matrices, use entries
            {
              if(unset) {sy=sut->entries()->symmetry(); unset=false;}
              else if(sut->entries()->symmetry()!=sy) return _noSymmetry;
            }
        }
      else return _noSymmetry;  //TermMatrix has non diagonal blocks (symmetry test not yet performed)
    }
  return sy;
}

// number of rows (counted in scalar dof)
number_t TermMatrix::numberOfRows() const
{
  if(scalar_entries_p!=0) return cdofs_r.size();
  if(suTerms_.size()==1) return begin()->second->numberOfRows();
  number_t n=0;
  std::map<const Unknown*,std::set<number_t> > rowdofs;
  std::map<const Unknown*,std::set<number_t> >::iterator itr;
  for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)
    {
      const Unknown* v=it->first.second;
      std::vector<number_t> dofs=it->second->space_vp()->dofIds();
      itr=rowdofs.find(v);
      if(itr==rowdofs.end()) rowdofs[v]=std::set<number_t>(dofs.begin(),dofs.end());
      else  rowdofs[v].insert(dofs.begin(),dofs.end());
    }
  for(itr=rowdofs.begin(); itr!=rowdofs.end(); itr++)
    n+=itr->second.size() * itr->first->nbOfComponents();
  return n;
}

// number of columns (counted in scalar dof)
number_t TermMatrix::numberOfCols() const
{
  if(scalar_entries_p!=0) return cdofs_c.size();
  if(suTerms_.size()==1) return begin()->second->numberOfCols();
  number_t n=0;
  std::map<const Unknown*,std::set<number_t> > coldofs;
  std::map<const Unknown*,std::set<number_t> >::iterator itc;
  for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)
    {
      const Unknown* u=it->first.first;
      std::vector<number_t> dofs=it->second->space_up()->dofIds();
      itc=coldofs.find(u);
      if(itc==coldofs.end()) coldofs[u]=std::set<number_t>(dofs.begin(),dofs.end());
      else  coldofs[u].insert(dofs.begin(),dofs.end());
    }
  for(itc=coldofs.begin(); itc!=coldofs.end(); itc++)
    n+=itc->second.size() * itc->first->nbOfComponents();
  return n;
}

number_t TermMatrix::nnz() const
{
  if (scalar_entries_p != 0) return scalar_entries_p->size();
  if (entries_p != 0) return entries_p->scalarSize();
  number_t n=0;
  for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)
    n+=it->second->nnz();
  return n;
}

// -----------------------------------------------------------------------------------------------------------
// access to matrix element
// -----------------------------------------------------------------------------------------------------------

//access to single unknown term matrix (internal use)
SuTermMatrix& TermMatrix::subMatrix(const Unknown* up, const Unknown* vp)
{
  if (up == 0)
  {
    where("TermMatrix::subMatrix(Unknown*, Unknown*)");
    error("termmatrix_submatrix_not_found", "up=0");
  }
  if (vp == 0)
  {
    where("TermMatrix::subMatrix(Unknown*, Unknown*)");
    error("termmatrix_submatrix_not_found", "vp=0");
  }
  uvPair uv(up, vp);
  it_mustm it = suTerms_.find(uv);
  if (it == suTerms_.end())
  {
    where("TermMatrix::subMatrix(Unknown*, Unknown*)");
    error("termmatrix_submatrix_not_found", "(" + up->name() + " " + vp->name() + ")");
  }
  return *(it->second);
}

const SuTermMatrix& TermMatrix::subMatrix(const Unknown* up, const Unknown* vp) const
{
  if (up == 0)
  {
    where("TermMatrix::subMatrix(Unknown*, Unknown*) const");
    error("termmatrix_submatrix_not_found", "up=0");
  }
  if (vp == 0)
  {

    where("TermMatrix::subMatrix(Unknown*, Unknown*) const");
    error("termmatrix_submatrix_not_found", "vp=0");
  }
  uvPair uv(up, vp);
  cit_mustm it = suTerms_.find(uv);
  if (it == suTerms_.end())
  {
    where("TermMatrix::subMatrix(Unknown*, Unknown*) const");
    error("termmatrix_submatrix_not_found", "(" + up->name() + " " + vp->name() + ")");
  }
  return *(it->second);
}

//access to single unknown term matrix (internal use)
SuTermMatrix* TermMatrix::subMatrix_p(const Unknown* up, const Unknown* vp)
{
  if (up == 0)
  {
    where("TermMatrix::subMatrix_p(Unknown*, Unknown*)");
    error("termmatrix_submatrix_not_found", "up=0");
  }
  if (vp == 0)
  {
    where("TermMatrix::subMatrix_p(Unknown*, Unknown*)");
    error("termmatrix_submatrix_not_found", "vp=0");
  }
  uvPair uv(up, vp);
  it_mustm it = suTerms_.find(uv);
  if(it != suTerms_.end()) return it->second;
  return 0;
}
const SuTermMatrix* TermMatrix::subMatrix_p(const Unknown* up, const Unknown* vp) const
{
  if (up == 0)
  {
    where("TermMatrix::subMatrix_p(Unknown*, Unknown*) const");
    error("termmatrix_submatrix_not_found", "up=0");
  }
  if (vp == 0)
  {
    where("TermMatrix::subMatrix_p(Unknown*, Unknown*) const");
    error("termmatrix_submatrix_not_found", "vp=0");
  }
  uvPair uv(up, vp);
  cit_mustm it = suTerms_.find(uv);
  if(it != suTerms_.end()) return it->second;
  return 0;
}

/*!return single unknown term matrix as TermMatrix (copy constructor like),
   useful to extract single unknown term from multiple unknowns term
   induce a full copy !!!
*/
TermMatrix TermMatrix::operator()(const Unknown& u, const Unknown& v) const
{
  const SuTermMatrix& sut = this->subMatrix(&u, &v);
  //create a new single unknown TermVector
  string_t na = name_ + "_" + u.name() + "_" + v.name();
  BilinearForm blf(*sut.sublf_p);                         //linear form from single unknown linear form
  TermMatrix tm(blf, na, sut.computingInfo_.noAssembly);  //create new TermMatrix
  SuTermMatrix& tmsu = *tm.suTerms_.begin()->second;
  tmsu.entries_p=0;
  if(sut.entries_p!=0)   //suTermMatrix is computed
    {
      tmsu.entries_p = new MatrixEntry(*sut.entries_p);       //copy MatrixEntries
      tmsu.computed()=true;
      tm.computed()=true;
    }
  tm.entries_p = tmsu.entries_p;                          //update MatrixEntry pointer of TermMatrix
  return tm;
}

// return the i-th row or col as a TermVector (r=1,...)
// restrict to single unknown TermMatrix or multiple unknown TermMatrix having global representation
TermVector& TermMatrix::getRowCol(number_t i, AccessType at, TermVector& R) const
{
  bool single=isSingleUnknown();
  if (!single && scalar_entries_p==0 && entries_p==0)
  {
    where("TermMatrix::getRowCol(Nuber, AccessType, TermVector");
    error("term_not_global_rep", name());
  }

  //initTermVector(R,valueType(),at==_col);
  initTermVector(R,valueType(),at==_row);  // init using row unknown if at=col and col unknown if at=row
  //retry matrix pointer
  MatrixEntry* mat=scalar_entries_p;
  bool scal = true;
  if (single && suTerms_.size()==1)
  {
    mat = suTerms_.begin()->second->scalar_entries();
    if (mat==0) {mat = suTerms_.begin()->second->entries(); scal=false;}

  }
  else if (mat==0) {mat=entries_p; scal=false;}
  if (mat==0)
  {
    where("TermMatrix::getRowCol(Nuber, AccessType, TermVector");
    error("term_no_entries");
  }

  number_t nbrc = mat->nbOfRows();
  if (at==_col) nbrc = mat->nbOfCols();
  if (i<1 || i>nbrc)
  {
    where("TermMatrix::getRowCol(Nuber, AccessType, TermVector");
    error("index_out_of_range",i, 1, nbrc);
  }

  VectorEntry vr = mat->getRowCol(i,at);
  if (single)
  {
    if (scal) R.firstSut()->scalar_entries() = new VectorEntry(vr);
    else R.firstSut()->entries() = new VectorEntry(vr);
  }
  else
  {
    if (scal) R.scalar_entries() = new VectorEntry(vr);
    else R.entries() = new VectorEntry(vr);
  }
  return R;
}

TermVector TermMatrix::row(number_t r) const
{
  TermVector R;
  return getRowCol(r,_row, R);
}

TermVector TermMatrix::column(number_t c) const
{
  TermVector C;
  return getRowCol(c,_col, C);
}

//initialize a TermVector consistent with matrix column unknown (col=true), with matrix row unknown (col=false)
void TermMatrix::initTermVector(TermVector& tv, ValueType vt, bool col) const
{
  string_t na = tv.name(); //preserving only name
  tv.clear();
  tv.name() = na;
  //make list of row unknowns
  std::map<const Unknown*, SuTermMatrix*> terms;
  for(cit_mustm itA = begin(); itA != end(); itA++)
    {
      const Unknown* u;
      if(col) u = itA->first.first;
      else    u = itA->first.second;
      SuTermMatrix* Auv = itA->second;
      if(terms.find(u) == terms.end()) terms[u] = Auv;
      else //update terms if u-space is larger than previous one
        {
          if(terms[u]->space_up()->dimSpace() < Auv->space_up()->dimSpace()) terms[u] = Auv;
          //to be improved : determine the largest space
        }
    }
  std::map<const Unknown*, SuTermMatrix*>::iterator it;
  for(it = terms.begin(); it != terms.end(); it++)
    {
      const Unknown* u = it->first;
      SuTermMatrix* sut = it->second;
      dimen_t nv = u->nbOfComponents();
      Space* sp = sut->space_up();
      if(!col) sp = sut->space_vp();
      number_t n = sp->dimSpace();
      SuTermVector* sutX = new SuTermVector("", u, sp, vt, n, nv);
      sutX->computed()=true;
      tv.insert(u, sutX);
    }
  tv.computed() = true;
}

// -----------------------------------------------------------------------------------------------------------
// compute the matrix term from bilinear form
// -----------------------------------------------------------------------------------------------------------
void TermMatrix::compute()
{
  if (computed()) return;  //already computed
  trace_p->push("TermMatrix::compute()");
  //compute each block
  for (it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)  //travel the map of single unknown linear forms
    { it->second->compute(); }
  /* merge blocks
    to avoid complexity and efficiency lost in elementary computation, some blocks have to be merged after
    this is the case when a block has a component of unknown (say u_1) and an other block refers to u
    we have to merge them using the mergeBlocks function
  */
  mergeBlocks();

  if (constraints_u_p==0 && constraints_v_p==0)
  {
    computed() = true;
    if(theVerboseLevel>0) printSummary(std::cout);
    trace_p->pop();
    return;
  }

  /* deal with essential conditions (constraints data)
     first, we have to move the matrix into scalar representation if there are some vector unknowns
     and may be to global scalar representation if there are some coupled unknowns condition
     then four methods are handled :
      - pseudo reduction : replace some matrix coefficients, see pseudoReduction function for details
      - real reduction : remove some matrix coefficients, see realReduction function for details (storage is modified)
      - penalisation : add large coefficients to some matrix coefficients, see penalizationReduction function for details
      - duality : introduce new dual unknown in matrix, see dualReduction function for details
  */
  if (computingInfo_.reductionMethod.method!=_pseudoReduction)
    error("reduction_method_not_handled", words("reduction method",computingInfo_.reductionMethod.method));
  bool global=false;

  //basic consistancy tests
  if (constraints_u_p!=0)
  {
    if (constraints_u_p->size()==0) error("is_void","SetOfConstraints");
    if (constraints_u_p->begin()->second==0) error("null_pointer","constraints_u_p");
    global=constraints_u_p->isGlobal();
  }
  if (constraints_v_p!=0)
  {
    if (constraints_v_p->size()==0) error("is_void","SetOfConstraints");
    if (constraints_v_p->begin()->second==0) error("null_pointer","constraints_v_p");
    if (!global) global=constraints_u_p->isGlobal();
  }

  //convert matrix to scalar representation and to global if constraints are global
  if (global) toGlobal(storageType(), storageAccess(), _noSymmetry, true);   //go to global scalar representation
  else toScalar();        //go to scalar representation

  switch (computingInfo_.reductionMethod.method)
  {
    case _noReduction              :  break;
    case _pseudoReduction          :  pseudoReduction(); break;
    // case _realReduction         : realReduction(); trace_p->pop(); break;
    // case _penalizationReduction : penalizationReduction(); trace_p->pop(); break;
    // case _dualReduction         : dualReduction(); trace_p->pop(); break;
    default :
      error("reduction_method_not_handled", words("reduction method",computingInfo_.reductionMethod.method));
  }

  computed() = true;
  if (theVerboseLevel>0)
   {
      printSummary(std::cout);
      if (theCout.traceOnFile)  printSummary(theCout.printStream->currentStream());
   }
  trace_p->pop();
}

// -----------------------------------------------------------------------------------------------------------
// merge blocks referring to same vector unknown
//   - when a block has a component of unknown as row or col unknown, it is convert to its vector representation
//   - when some blocks have different components of same unknown as row or col unknown, the blocks are really merged
//     into a block with vector unknowns
// see the function mergeSuTermMatrix
// -----------------------------------------------------------------------------------------------------------
void TermMatrix::mergeBlocks()
{
  trace_p->push("TermMatrix::mergeBlocks()");
  std::map<uvPair, std::list<SuTermMatrix*> > suts; //to collect block by vector unknowns
  std::map<uvPair, std::list<SuTermMatrix*> >::iterator itm;
  bool mergeFlag=false;
  for(it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)  //travel the map of single unknown SuTermMatrix
    {
      SuTermMatrix* sut=it->second;
      const Unknown* up=sut->up(), *vp=sut->vp();
      if(up->isComponent()) {up=up->parent(); mergeFlag=true;}
      if(vp->isComponent()) {vp=vp->parent(); mergeFlag=true;}
      itm=suts.find(uvPair(up,vp));
      if(itm!=suts.end()) itm->second.push_back(sut);
      else suts[uvPair(up,vp)]=std::list<SuTermMatrix*>(1,sut);
    }
  if(!mergeFlag) {trace_p->pop(); return;}   //nothing to merge

  //now merge each collected block
  for(itm=suts.begin(); itm!=suts.end(); itm++)
    {
      std::list<SuTermMatrix*>::iterator itl;
      std::vector<uvPair> uvel(itm->second.size());
      std::vector<uvPair>::iterator ituv=uvel.begin();
      for(itl=itm->second.begin(); itl!=itm->second.end(); itl++,ituv++)  //memorize uv pairs
        *ituv= uvPair((*itl)->up(),(*itl)->vp());
      SuTermMatrix* msut=mergeSuTermMatrix(itm->second);   //merge blocks, NO DELETION
      if(msut!=0)  //real merging
        {
          std::list<SuTermMatrix*>::iterator itl;
          ituv=uvel.begin();
          for(itl=itm->second.begin(); itl!=itm->second.end(); itl++, ituv++)  //delete old SuTermMatrix's
            {
              suTerms_.erase(*ituv);
              delete *itl;
            }
          suTerms_[itm->first]=msut;  //assign the merged SuTermMatrix
        }
    }

  trace_p->pop();
}

// -----------------------------------------------------------------------------------------------------------
// deal with essential condition (constraints)
// -----------------------------------------------------------------------------------------------------------

/*!
  perform pseudo reduction of reduced essential conditions in a matrix. Essential conditions have the following form  :
                          U_E + F*U_R = s   for column unknown U
                          V_E + G*V_R = 0   for row test function V (generally related to unknown U)
  where E are the eliminated unknown/test function indices and R are the reduced unknown/test function indices
  The pseudo reduction of matrix consists in
    - modifying column A.j for j in R by adding -Fjk*A.k for k in E and replacing column A.k for k in E by a 0 column
    - modifying row Ai. for i in R by adding -Gik*Ak. for k in E and replacing row Ak. for k in E by a 0 row
    - if eliminated v unknowns are dual of eliminated u unknowns, the VE rows are replaced by the equation (or a part of)
                                      U_E + F*U_R = s
  to delay the right hand side modification, the (A.k) columns (k in E) are stored in a new structure

  At the end of the process, the eliminated system looks like
                   U_E     U_R    U_S       RHS
                 ----------------------   ------
                 |      |       |     |   |    |
            V_E  |  Id  |   F   |  0  |   |  S |
                 |      |       |     |   |    |
                 ---------------------    -----
                 |      |       |     |   |    |
            V_R  |  0   |  ARR  | ARS |   | FR |
                 |      |       |     |   |    |
                 ----------------------   ------
                 |      |       |     |   |    |
            V_S  |  0   |  ASR  | ASS |   | FS |
                 |      |       |     |   |    |
                 ----------------------   ------

  In some cases (F=G=0 or local conditions u^n=g, ...) the storage is not modified but in other cases (transmission condition for instance)
  the storage is modified

*/
void TermMatrix::pseudoReduction()
{
  trace_p->push("TermMatrix::pseudoReduction()");
  if (constraints_u_p == 0 && constraints_v_p ==0)   //no essential condition
    {
      warning("free_warning", " in TermMatrix::pseudoReduction, no essential conditions to deal with");
      return;
    }

  bool global = false;
  complex_t alpha=computingInfo().reductionMethod.alpha;
  if(constraints_u_p != 0 && constraints_u_p->isGlobal()) global=true;
  else if(constraints_v_p != 0 && constraints_v_p->isGlobal()) global=true;
  Constraints* cons_u=0, *cons_v=0;

  if(global)  //global reduction
    {
      //modify matrix storage and symmetry
      if(constraints_u_p != 0) cons_u=(*constraints_u_p)(0);
      if(constraints_v_p != 0) cons_v=(*constraints_v_p)(0);
      if(scalar_entries_p->symmetry()!=_noSymmetry && constraints_u_p != constraints_v_p) scalar_entries_p->toUnsymmetric(_sym);
      extendStorage(scalar_entries_p, cdofs_r, cdofs_c, cons_u, cons_v);  //extend storage if required
      if(cons_u!=0)  cons_u->pseudoColReduction(scalar_entries_p, cdofs_r, cdofs_c, rhs_matrix_p);
      else error("global_sub_reduction_not_yet_handled", words("column"));
      //if (cons_v!=0 && scalar_entries_p->symmetry()==_noSymmetry) //row reduction
      if(cons_v!=0)  //row reduction
        cons_v->pseudoRowReduction(scalar_entries_p, cdofs_r, cdofs_c, alpha);
      else if(constraints_v_p!=0 && cons_v==0)
        error("global_sub_reduction_not_yet_handled", words("row"));
    }
  else //sub reduction (i.e on SuTermMatrix's)
    {
      for(it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)
        {
          SuTermMatrix* sutm=it->second;
          const Unknown* u=it->first.first, *v=it->first.second;
          Constraints* cons_u = (*constraints_u_p)(u), *cons_v = 0;
          if(constraints_v_p ==0) cons_v = (*constraints_v_p)(v->dual_p());
          if(cons_v == 0) cons_v = (*constraints_v_p)(v->dual_p());  //try with dual unknown
          MatrixEntry* submat=sutm->scalar_entries_p;
          if(submat->symmetry()!=_noSymmetry && cons_u != cons_v)  submat->toUnsymmetric(_sym);  //goto non symmetric matrix stored with a symmetric storage
          extendStorage(submat, sutm->cdofs_v, sutm->cdofs_u, cons_u, cons_v);                  //extend storage if required
          if(cons_u!=0)   //col reduction
            cons_u->pseudoColReduction(submat, sutm->cdofs_v, sutm->cdofs_u, sutm->rhs_matrix_p);
          //if (cons_v!=0 && submat->symmetry()==_noSymmetry)  //row reduction
          if(cons_v!=0)   //row reduction
            cons_v->pseudoRowReduction(submat, sutm->cdofs_v, sutm->cdofs_u, alpha);
        }
    }
  trace_p->pop();
}

// -----------------------------------------------------------------------------------------------------------
// change representation of matrix (entries and storage)
// -----------------------------------------------------------------------------------------------------------
/*! regarding a multiple unknowns TermMatrix, suggest a storage according to the following rules:
      if dense part > 10*(cs part+skyline part) -> dense storage  else -> cs storage
      if umfpack installed -> column access else dual row/col access
      if matrix has symmetry -> symmetric access */
std::pair<StorageType,AccessType> TermMatrix::findGlobalStorageType() const
{
  number_t densePart=0, csPart=0, skyPart=0;
  number_t denseFull = numberOfRows()*numberOfCols();
  StorageType st=_noStorage; AccessType at=_noAccess;
  for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)
    {
      StorageType sts=it->second->storageType();
      number_t s = it->second->nnz();
      if(sts==_dense) densePart+=s;
      if(sts==_cs) csPart+=s;
      if(sts==_skyline) skyPart+=s;
    }
  at=_dual;
  //if(densePart > 10*(csPart+skyPart)) {st=_dense; at=_row;}
  if(4*densePart > denseFull) {st=_dense; at=_row;}  //dense part represents 25% of the matrix, choose dense
  else
    {
      st=_cs;
#ifdef XLIFEPP_WITH_UMFPACK
      at=_col;     //well adapted to umfpack
#endif
    }

  //if(symmetry()!=_noSymmetry) at=_sym;     //raises a bug in toglobal, sym -> sym
  return std::make_pair(st,at);
}

void TermMatrix::toSkyline()
{
  trace_p->push("TermMatrix::toSkyline()");
  if(!computed()) error("not_computed_term", name());
  if(isSingleUnknown())
    {
      SuTermMatrix* sutm = begin()->second;
      if(sutm->entries()->storageType() == _skyline)
        {
          warning("free_warning", "TermMatrix" + name() + " is already a skyline matrix, nothing to do!");
          trace_p->pop();
          return;
        }
      if(sutm->entries()->storageType() == _dense)
        error("storage_not_implemented", "TermMatrix::toSkyline", words("storage type",_dense));
      if(sutm->entries()->storageType() == _cs) sutm->entries()->toSkyline();  //convert entries to skyline storage
    }
  else
    {
      error("block_conversion_not_yet_implemented", words("storage type",_skyline));
    }
  trace_p->pop();
}

/*! create scalar representation of TermMatrix
    for each SuTermMatrix (if computed, entries_!=0), if not exists create its scalar representation (set scalar_entries_p and cdofs vectors)
      if SuTermMatrix is already scalar, scalar_entries_p = entries_p
    NOTE : scalar_entries_p of TermMatrix is not update by this function (see toGlobal function)

*/
void TermMatrix::toScalar(bool keepEntries)
{
  trace_p->push("TermMatrix::toScalar()");
  for(it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++) it->second->toScalar(keepEntries);
  trace_p->pop();
}

/*! create scalar global representation of TermMatrix
    make scalar representations of SuTermMatrix's if not exist and merge it in a scalar global representation by allocating
    scalar_entries_p and setting the cdofs vectors

    st, at : storage and access type of global matrix, if st = _nosStorage try to identify "best" storage
    symt   : to force symmetry of matrix (default = _noSymmetry)
    keepSuterms : if false (default value), SutermMatrix entries are deleted
*/
void TermMatrix::toGlobal(StorageType st, AccessType at, SymType symt, bool keepSuTerms)
{
  if (scalar_entries_p!=0) return;  //already done, delete it outside if you want to rebuild
  trace_p->push("TermMatrix::toGlobal()");

  //create scalar representation
  toScalar();   //nothing is done if SuTermMatrix's already scalar !

  //cpuTime("toGlobal:toScalar");
  if (suTerms_.size()==1) {trace_p->pop(); return;}  //nothing to merge

  //cpuTime("toGlobal:enter to");
  //merge all SutermMatrix's
  // compare row and col unknowns
  bool sameRowColUnknowns = true;
  std::set<const Unknown*> vs=rowUnknowns(), us=colUnknowns();
  std::set<const Unknown*>::iterator itvs=vs.begin(), itus;
  if (us.size()!=vs.size()) sameRowColUnknowns = false;
  else
    for (; itvs!=vs.end() && sameRowColUnknowns; itvs++)
    {
      itus=us.find(*itvs);
      if (itus==us.end()) itus=us.find((*itvs)->dual_p());
      if (itus==us.end()) sameRowColUnknowns=false;
    }

  if (!sameRowColUnknowns && symt!=_noSymmetry) error("symmetry_not_handled", words("symmetry", symt));

  //try to find storage type if not set (has to be improved)
  if (st==_noStorage)
  {
    std::pair<StorageType,AccessType> stat=findGlobalStorageType();
    st=stat.first; at=stat.second;
  }

  //create numbering of blocks in global numbering (row and size)
  std::map<const Unknown*, std::list<SuTermMatrix*> > rowsut, colsut;  //list of SuTermMatrix by row/col unknown
  std::map<SuTermMatrix*, std::vector<number_t> > rownum, colnum;      //block row/col numbering map for each SuTermMatrix
  ValueType vt=_real;
  for (it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)
  {
    SuTermMatrix* sut= it->second;
    rowsut[it->first.second].push_back(sut);
    colsut[it->first.first].push_back(sut);
    if(sut->scalar_entries_p->valueType_ ==_complex) vt=_complex;
  }
  mergeNumberingFast(rowsut, rownum, cdofs_r,_row);
  mergeNumberingFast(colsut, colnum, cdofs_c,_col);
  //cpuTime("toGlobal:mergeNumbering");

  //buiding row indices of global matrix and calling buildstorage
  MatrixStorage* ms=0;
  if (st!=_dense)   //update storage pointer
  {
    std::vector<std::set<number_t> > indices(cdofs_r.size());  //col indices by row
    for (it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)
    {
      SuTermMatrix* sut= it->second;
      MatrixStorage* msu=sut->scalarStoragep();
      std::vector<number_t>& rowmap=rownum[sut];
      std::vector<number_t>& colmap=colnum[sut];
      addIndices(indices, msu, rowmap, colmap);
    }
    //cpuTime("toGlobal:prepare buildstorage");
    ms=buildStorage(st,at,cdofs_r.size(), cdofs_c.size(),indices);
    //cpuTime("toGlobal:buildstorage");
  }
  else ms=buildStorage(st,at,cdofs_r.size(), cdofs_c.size());

  //cpuTime("toGlobal:create_storage");
  scalar_entries_p=new MatrixEntry(vt, _scalar, ms, dimPair(1,1), symt);
  //merge SutermMatrix's
  for (it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)
  {
    SuTermMatrix* sut= it->second;
    scalar_entries_p->assign(*sut->scalar_entries_p, rownum[sut], colnum[sut]);
  }
  //cpuTime("toGlobal:merge_submatrices "+tostring(suTerms_.size()));

  if (!keepSuTerms)   //clear SuTermMatrix's
    for (it_mustm it = suTerms_.begin(); it != suTerms_.end(); it++) it->second->clear();

  if (theVerboseLevel>0)
     theCout<<"multiple unknown TermMatrix "<<name()<<" moves to global representation : "
            << words("value",scalar_entries_p->valueType())<<" scalar matrix "<<numberOfRows()<< " X "<<numberOfCols()
            <<" in "<<scalar_entries_p->storagep()->name()<<" ("<< scalar_entries_p->storagep()->size()<<" coefficients)"<< eol;
  trace_p->pop();
}

// Turns to scalar representation in case of multiple unknowns and returns pointer to data held in the matrix
MatrixEntry* TermMatrix::matrixData() {
  MatrixEntry* matA=0;

  if (isSingleUnknown()) {
    SuTermMatrix* sutmA = begin()->second;
    if (sutmA->space_up() != sutmA->space_vp())
      error("term_inconsistent_unknown_spaces", name());
    matA = sutmA->entries();
    if (sutmA->scalar_entries()!=0) matA= sutmA->scalar_entries();
  }
  else {  //multiple unknowns (goto scalar global representation)
    SymType syA=symmetry();  //detect global symmetry of TermMatrix
    if (entries()==0 && scalar_entries()==0) toGlobal(storageType(), storageAccess(), syA, true);
    else if (scalar_entries()==0) toScalar(false);
    matA = scalar_entries();
  }
  return matA;
}

//toGlobal tool : add col indices of a matrix storage to current set of col indices with index mapping
void TermMatrix::addIndices(std::vector<std::set<number_t> >& indices, MatrixStorage* msu,
                            const std::vector<number_t>& rowmap, const std::vector<number_t>& colmap)
{
  std::vector<number_t>::const_iterator itrc;
  std::set<number_t>::iterator its;
  number_t k=1, shr=0, shc=0, nbc=msu->nbOfColumns(), nbr=msu->nbOfRows(), m;
  AccessType msuat=msu->accessType();
  if(rowmap.size()==1) shr=rowmap[0];
  if(colmap.size()==1) shc=colmap[0];

  if(msuat==_row || msuat==_dual || msuat==_sym)   //row in row
    {
      m = nbc;
      for(k=1; k<=nbr; ++k)
        {
          if(msuat!=_row) m=k;    //restrict to lower part
          std::set<number_t> cols=msu->getCols(k,1,m);      // retry col indices of row k in ms storage
          std::set<number_t>* indr=&indices[k+shr-1];
          if(rowmap.size()>1) indr=&indices[rowmap[k-1]-1];
          if(colmap.size()==1)
            {
              if(shc==0) indr->insert(cols.begin(), cols.end());
              else
                for(its=cols.begin(); its!=cols.end(); ++its)  indr->insert((*its)+shc);     // insert col numbering shift
            }
          else
            for(its=cols.begin(); its!=cols.end(); ++its)  indr->insert(colmap[(*its)-1]);  // insert col numbering map
        }
    }
  if(msuat==_col || msuat==_dual || msuat==_sym)   //col in row
    {
      m = nbr;
      number_t d=2;
      if(msuat==_col) d=1;
      for(k=d; k<=nbc; ++k)
        {
          number_t c=k+shc;
          if(colmap.size()>1) c=colmap[k-1];
          if(msuat!=_col) m=k-1;   // restrict to upper part, diag term has been already taken into account
          std::set<number_t> rows=msu->getRows(k,1,m);    // retry row indices of col k in msu storage
          if(rowmap.size()==1)
            for(its=rows.begin(); its!=rows.end(); ++its)  indices[(*its)+shr-1].insert(c);   // insert col numbering map
          else
            for(its=rows.begin(); its!=rows.end(); ++its)  indices[rowmap[(*its)-1]-1].insert(c);  // insert col numbering map
        }
    }
}

//tool for TermMatrix::toGlobal function
//  create the list of merged row/col dofs
//  and for each SuTermMatrix the map between is row/col numbering and the global row/col numbering
void TermMatrix::mergeNumbering(std::map<const Unknown*, std::list<SuTermMatrix*> >& rcsut,
                                std::map<SuTermMatrix*, std::vector<number_t > >& rcnum,
                                std::vector<DofComponent>& rcdof, AccessType rc)
{
  //order unknowns according to their ranks
  std::map<const Unknown*, std::list<SuTermMatrix*> >::iterator itmus;
  std::map<number_t,const Unknown*> order;
  for(itmus=rcsut.begin(); itmus!=rcsut.end(); itmus++)
    {
      const Unknown* u=itmus->first;
      number_t r=u->rank();
      if(order.find(r)==order.end()) order.insert(std::make_pair(r,u));
    }
  //merge global numbering
  std::vector<DofComponent>::iterator itd;
  std::map<number_t,const Unknown*>::iterator itmnu;
  for(itmnu=order.begin(); itmnu!=order.end(); itmnu++)
    {
      const Unknown* u=itmnu->second;
      std::list<SuTermMatrix* >::iterator itl=rcsut[u].begin(), itle=rcsut[u].end();
      if(rcsut[u].size()==1)  //simple case only one SutermMatrix
        {
          SuTermMatrix* sut = *itl;
          std::vector<DofComponent>* cd=&sut->cdofsu();
          if(rc==_row) cd=&sut->cdofsv();
          number_t nbr=cd->size();
          rcnum[sut]=trivialNumbering<number_t>(rcdof.size()+1,rcdof.size()+nbr);
          rcdof.insert(rcdof.end(),cd->begin(),cd->end());
        }
      else //merge local numbering
        {
          std::map<DofComponent, number_t> dofmap;
          std::map<DofComponent, number_t>::iterator itmc;
          number_t ld=rcdof.size();
          for(; itl!=itle; itl++)  //travel SuTermMatrix list
            {
              SuTermMatrix* sut = *itl;
              std::vector<DofComponent>* cd=&sut->cdofsu();
              if(rc==_row) cd=&sut->cdofsv();
              number_t nbcd=cd->size();
              std::vector<DofComponent> locd(nbcd);
              std::vector<number_t >& num = rcnum[sut];
              num.resize(nbcd);
              std::vector<number_t >::iterator itc=num.begin();
              number_t k=0;
              for(itd=cd->begin(); itd!=cd->end(); itd++, itc++)
                {
                  itmc=dofmap.find(*itd);
                  if(itmc==dofmap.end())
                    {
                      ld++;            //increase global dof index
                      dofmap[*itd]=ld; //update dofmap
                      *itc=ld;         //update local->global numbering
                      locd[k]=*itd;    //update local cdof
                      k++;
                    }
                  else *itc=itmc->second;
                }
              locd.resize(k);
              rcdof.insert(rcdof.end(),locd.begin(),locd.end());
            }
        }
    }
}

//tool for TermMatrix::toGlobal function, fast version
//  create the list of merged row/col dofs
//  and for each SuTermMatrix the map between is row/col numbering and the global row/col numbering
//  trivial list n->n+p are encoded only using the shift n : rcnum[sut]=[n]
void TermMatrix::mergeNumberingFast(std::map<const Unknown*, std::list<SuTermMatrix*> >& rcsut,
                                    std::map<SuTermMatrix*, std::vector<number_t > >& rcnum,
                                    std::vector<DofComponent>& rcdof, AccessType rc)
{
  //order unknowns according to their ranks
  std::map<const Unknown*, std::list<SuTermMatrix*> >::iterator itmus;
  std::map<number_t,const Unknown*> order;
  for(itmus=rcsut.begin(); itmus!=rcsut.end(); itmus++)
    {
      const Unknown* u=itmus->first;
      number_t r=u->rank();
      if(order.find(r)==order.end()) order.insert(std::make_pair(r,u));
    }
  //merge global numbering
  std::vector<DofComponent>::iterator itd;
  std::map<number_t,const Unknown*>::iterator itmnu;
  for(itmnu=order.begin(); itmnu!=order.end(); itmnu++)
    {
      const Unknown* u=itmnu->second;
      //order u-SuTermMatrix's by their row/col sizes largest to smallest
      std::vector<SuTermMatrix*> ordsut(rcsut[u].begin(),rcsut[u].end());
      //std::sort(ordsut.begin(),ordsut.end(),compNnzSutermMatrix);     //rearrange from the largest to the smallest SuTermMatrix
      if(rc==_row) std::sort(ordsut.begin(),ordsut.end(),compRowSize);  //rearrange from the largest to the smallest row sizes
      else          std::sort(ordsut.begin(),ordsut.end(),compColSize); //rearrange from the largest to the smallest col sizes

      //deal with first SuTermMatrix
      std::vector<SuTermMatrix* >::iterator itl=ordsut.begin(), itle=ordsut.end();
      SuTermMatrix* sut = *itl;
      std::vector<DofComponent>* cd=&sut->cdofsu();
      if(rc==_row) cd=&sut->cdofsv();
      number_t ld=rcdof.size()+1;
      rcnum[sut]=std::vector<number_t >(1,rcdof.size());   //trivial numbering
      rcdof.insert(rcdof.end(),cd->begin(),cd->end());
      std::map<DofComponent, number_t> dofmap;
      std::map<DofComponent, number_t>::iterator itmc;
      if(ordsut.size()>1)
        for(itd=cd->begin(); itd!=cd->end(); ++itd, ++ld) dofmap[*itd]=ld;
      itl++;
      ld=rcdof.size();

      for(; itl!=itle; itl++)  //travel SuTermMatrix list
        {
          bool trivial=true;
          SuTermMatrix* sut = *itl;
          std::vector<DofComponent>* cd=&sut->cdofsu();
          if(rc==_row) cd=&sut->cdofsv();
          number_t nbcd=cd->size();
          std::vector<DofComponent> locd(nbcd);
          std::vector<number_t >& num = rcnum[sut];
          num.resize(nbcd);
          std::vector<number_t >::iterator itc=num.begin();
          number_t k=0;
          for(itd=cd->begin(); itd!=cd->end(); itd++, itc++)
            {
              itmc=dofmap.find(*itd);
              if(itmc==dofmap.end())
                {
                  ld++;            //increase global dof index
                  dofmap[*itd]=ld; //update dofmap
                  *itc=ld;         //update local->global numbering
                  locd[k]=*itd;    //update local cdof
                  k++;
                }
              else {*itc=itmc->second; trivial=false;}
            }
          locd.resize(k);
          if(trivial) rcnum[sut]= std::vector<number_t >(1,rcdof.size());    //trivial numbering
          rcdof.insert(rcdof.end(),locd.begin(),locd.end());
        }
    }
}

// -----------------------------------------------------------------------------------------------------------
// Algebraic operations on TermMatrix (+=, -=)
// -----------------------------------------------------------------------------------------------------------

/*! Algebraic operation += and -=
    the result of these operations is a TermMatrix with all unknown blocks of current TermMatrix and added Termatrix
        | Avu Avp    |    | Mvu        |         | Avu+Mvu Avp        |
        | Aqu        | += | Mqu Mqp    |   ==>   | Aqu+Mqu Mqp        |
        |     Atp Ats|    | Mtu     Mts|         | Mtu     Atp Ats+Mts|
    in other words, the matrices are merged and summation is a real one on common unknowns blocks
*/
TermMatrix& TermMatrix::operator+=(const TermMatrix& tM)
{
  bilinForm_ += tM.bilinForm_;
  if(! computed()) return *this;      //no hard copy

  for(it_mustm it = begin(); it != end(); it++)  //travel the blocks of current matrix
    {
      const Unknown* up = it->first.first, *vp = it->first.second; //unknowns of current block
      SuTermMatrix* sut = it->second;                              //current block
      const SuTermMatrix* stM = tM.subMatrix_p(up, vp);            //find (up,vp) block in tM
      if(stM != 0)(*sut) += *stM;                                   //add block to current block
    }
  for(cit_mustm itM = tM.begin(); itM != tM.end(); itM++)  //travel the blocks of added matrix
    {
      const Unknown* up = itM->first.first, *vp = itM->first.second; //unknowns of current block
      SuTermMatrix* sutM = itM->second;                              //current block
      const SuTermMatrix* sut = subMatrix_p(up, vp);                 //find (up,vp) block in current matrix
      if(sut == 0) suTerms_[uvPair(up, vp)] = new SuTermMatrix(*sutM);  //create new block in current matrix
    }
  return *this;
}

TermMatrix& TermMatrix::operator-=(const TermMatrix& tM)
{
  bilinForm_ -= tM.bilinForm_;
  if(! computed()) return *this;      //no hard copy

  for(it_mustm it = begin(); it != end(); it++)  //travel the blocks of current matrix
    {
      const Unknown* up = it->first.first, *vp = it->first.second; //unknowns of current block
      SuTermMatrix* sut = it->second;                          //current block
      const SuTermMatrix* stM = tM.subMatrix_p(up, vp);        //find (up,vp) block in tM
      if(stM != 0)(*sut) -= *stM;                               //add block to current block
    }
  for(cit_mustm itM = tM.begin(); itM != tM.end(); itM++)  //travel the blocks of added matrix
    {
      const Unknown* up = itM->first.first, *vp = itM->first.second; //unknowns of current block
      SuTermMatrix* sutM = itM->second;                           //current block
      const SuTermMatrix* sut = subMatrix_p(up, vp);              //find (up,vp) block in current matrix
      if(sut == 0)  //create new block in current matrix
        {
          SuTermMatrix* nsut = new SuTermMatrix(*sutM);
          *nsut *= -1.;
          suTerms_[uvPair(up, vp)] = nsut;
        }
    }
  return *this;
}

// -----------------------------------------------------------------------------------------------------------
// get real/imaginary part or move from real/complex to complex/real representation
// -----------------------------------------------------------------------------------------------------------
TermMatrix& TermMatrix::toReal()
{
  if(valueType()==_real) return *this;
  for(it_mustm it = begin(); it != end(); it++) it->second->toReal();
  return *this;
}

TermMatrix real(const TermMatrix& tm)
{
  TermMatrix atm=tm;
  atm.name()="real("+tm.name()+")";
  if(atm.valueType()==_complex) atm.toReal();
  return atm;
}

TermMatrix& TermMatrix::toComplex()
{
  for(it_mustm it = begin(); it != end(); it++) it->second->toComplex();
  return *this;
}

TermMatrix toComplex(const TermMatrix& tm)
{
  TermMatrix atm=tm;
  atm.name()="complex("+tm.name()+")";
  atm.toComplex();
  return atm;
}

TermMatrix& TermMatrix::toImag()
{
  for(it_mustm it = begin(); it != end(); it++) it->second->toImag();
  return *this;
}

TermMatrix imag(const TermMatrix& tm)
{
  TermMatrix atm=tm;
  atm.name()="imag("+tm.name()+")";
  atm.toImag();
  return atm;
}

TermMatrix& TermMatrix::toConj()
{
  for(it_mustm it = begin(); it != end(); it++) it->second->toConj();
  return *this;
}

TermMatrix conj(const TermMatrix& tm)
{
  TermMatrix atm=tm;
  if(tm.valueType()==_real) warning("free_warning"," taking the conjugate of the REAL TermMatrix "+tm.name()+" is useless");
  else
  {
      atm.name()="conj("+tm.name()+")";
      atm.toConj();
  }
  return atm;
}



TermMatrix& TermMatrix::roundToZero(real_t aszero)
{
  for(it_mustm it = begin(); it != end(); it++) it->second->roundToZero(aszero);
  return *this;
}

TermMatrix roundToZero(const TermMatrix& tm, real_t aszero = 10*theEpsilon )  //! return the 0 rounded TermMatrix)
{
   TermMatrix atm=tm;
   atm.name()="roundTo0("+tm.name()+")";
   atm.roundToZero(aszero);
   return atm;
}

// Frobenius norm : sqrt(sum|a_ij|^2)
real_t TermMatrix::norm2() const
{
    if(scalar_entries_p!=0) return scalar_entries_p->norm2();
    real_t n=0.;
    for(cit_mustm it = begin(); it != end(); it++)  //travel the blocks of current matrix
    {
        real_t nb=it->second->norm2();
        n+=nb*nb;
    }
    return std::sqrt(n);
}

// infinite norm : sup|a_ij|
real_t TermMatrix::norminfty() const
{
    if(scalar_entries_p!=0) return scalar_entries_p->norminfty();
    real_t n=0.;
    for(cit_mustm it = begin(); it != end(); it++)  //travel the blocks of current matrix
        n=std::max(n,it->second->norminfty());
    return n;
}

// -----------------------------------------------------------------------------------------------------------
// direct solver tools
// -----------------------------------------------------------------------------------------------------------

// member functions
// solve linear system using LU factorization (current matrix is preserved)
void TermMatrix::luSolve(const TermVector& B, TermVector& X)
{
  TermMatrix Af;
  factorize((*this), Af, _lu);
  X=factSolve(Af, B);
}

// solve linear system using LDLt factorization (current matrix is preserved)
void TermMatrix::ldltSolve(const TermVector& B, TermVector& X)
{
  TermMatrix Af;
  factorize((*this), Af, _ldlt);
  X=factSolve(Af, B);
}

// solve linear system using LDL* factorization (current matrix is preserved)
void TermMatrix::ldlstarSolve(const TermVector& B, TermVector& X)
{
  TermMatrix Af;
  factorize((*this), Af, _ldlstar);
  X=factSolve(Af, B);
}

// solve linear system using umfpack(current matrix is preserved)
void TermMatrix::umfpackSolve(const TermVector& B, TermVector& X)
{
  X = xlifepp::umfpackSolve(*this, B);
}

//------------------------------------------------------------------------------------------------------------
// particular solver used by SOR and SSOR iterative method
//------------------------------------------------------------------------------------------------------------
void TermMatrix::sorDiagonalMatrixVector(const TermVector& B, TermVector& X, const real_t w)
{
  sorSolve(B, X, w, _matrixVectorSorS);
}
void TermMatrix::sorUpperSolve(const TermVector& B, TermVector& X, const real_t w)
{
  sorSolve(B, X, w, _upperSorS);
}
void TermMatrix::sorLowerSolve(const TermVector& B, TermVector& X, const real_t w)
{
  sorSolve(B, X, w, _lowerSorS);
}
void TermMatrix::sorDiagonalSolve(const TermVector& B, TermVector& X, const real_t w)
{
  sorSolve(B, X, w, _diagSorS);
}

void TermMatrix::sorSolve(const TermVector& B, TermVector& X, const real_t w, SorSolverType sType)
{
  trace_p->push("sorSolver(TermMatrix&, const TermVector&, real_t , SorSolvertype)");
  //create TermVector X
  ValueType vtX = _real;
  if (this->valueType() == _complex || B.valueType() == _complex) vtX = _complex;
  this->initTermVector(X, vtX, true);

  //solve linear system
  if (this->isSingleUnknown())
  {
    const Unknown* v = this->begin()->first.second, *u=this->begin()->first.first;
    SuTermMatrix* sutm = this->begin()->second;
    //check (*this), B unknown consistancy
    const SuTermVector* vB = B.subVector_p(u);
    if (vB == 0)
    {
      vB=B.subVector_p(u->dual_p());
      if (vB==0) error("term_inconsistent_unknowns");
    }
    //prepare matrix, right hand side and solution
    if (sutm->strucType()==_matrix)  sutm->toScalar(true);
    SuTermVector cvB(*vB);    //copy B entries
    SuTermVector* vX = X.subVector_p(u);
    if (vX == 0) {vX = X.subVector_p(u->dual_p()); }
    MatrixEntry* mat=0;
    VectorEntry* b=0, *x=0;
    if (sutm->scalar_entries()==0)   //scalar system
    {
      mat=sutm->entries();
      b=cvB.entries();
      std::vector<number_t> renum = renumber(sutm->space_up(), vB->spacep());
      if (renum.size() != 0) b->extendEntries(renum,sutm->space_up()->dimSpace());    //extend vector to matrix column numbering
      x=vX->entries();
    }
    else //vector system in scalar form
    {
      mat=sutm->scalar_entries();
      b=cvB.scalar_entries();
      if (b==0) {cvB.toScalar(); b=cvB.scalar_entries();}
      std::vector<number_t> renum;
      if (cvB.up() == v) renum = renumber(sutm->cdofsv(), cvB.cdofs());
      else               renum = renumber(sutm->cdofsv(), dualDofComponents(cvB.cdofs()));
      if (renum.size() != 0) b->extendEntries(renum,sutm->cdofsv().size());    //extend vector to matrix column numbering
      vX->toScalar();
      x=vX->scalar_entries();
    }
    if ((*this).hasConstraints())  // case of constraints on (*this)
    {
      mat=sutm->scalar_entries();
      if (vX->scalar_entries()==0) x->toScalar();
      else x=vX->scalar_entries();
      cvB.toScalar();         //force scalar representation
      b=cvB.scalar_entries();
      Constraints* cons_u = 0, *cons_v = 0;
      if (this->constraints_u_p !=0) cons_u = (*(this->constraints_u_p))(u);
      if (this->constraints_v_p !=0)
      {
        cons_v = (*(*this).constraints_v_p)(v);
        if(cons_v == 0) cons_v = (*(*this).constraints_v_p)(v->dual_p());          //try with dual unknown
      }
      appliedRhsCorrectorTo(b, cvB.cdofs(), sutm->rhs_matrix(), cons_u, cons_v, this->computingInfo_.reductionMethod); //correct rhs to take into account constraints
    }
    //goto solvers
    switch(sType)
    {
      case _diagSorS    :
        mat->sorDiagonalSolve(*b, *x, w);
        break;
      case _upperSorS :
        mat->sorUpperSolve(*b, *x, w);
        break;
      case _lowerSorS     :
        mat->sorLowerSolve(*b, *x, w);
        break;
      case _matrixVectorSorS     :
        mat->sorDiagonalMatrixVector(*b, *x, w);
        break;
      default :
        error("undef_option", "SorSolverType");
    }
    if (vX->up() == sutm->up())  vX->cdofs()=sutm->cdofsu();         //update cdofs
    else vX->cdofs()=dualDofComponents(sutm->cdofsu());
    if (vX->up()->nbOfComponents() > 1)  vX->toVector(true);  //return to vector representation
    vX->computed() = true;
    X.computed() = true;
  }
  else
  {
    TermVector Bc(B);      //copy of B
    if (this->hasConstraints())  // case of constraints on (*this), applied rhs correction
    {
      Constraints* cons_u = 0, *cons_v = 0;
      if (this->constraints_u_p !=0) cons_u = (*(this->constraints_u_p))(0);
      if (this->constraints_v_p !=0) cons_v = (*(this->constraints_v_p))(0);
      if (cons_u!=0 || cons_v!=0)   //global constraint
      {
        this->toGlobal(_noStorage, _noAccess, _noSymmetry, false);
        Bc.toGlobal(false);                            //go to scalar global representation
        Bc.adjustScalarEntries(this->cdofsr());            //adjust vector to matrix column scalar numbering
        appliedRhsCorrectorTo(Bc.scalar_entries(), Bc.cdofs(), this->rhs_matrix(),
                              cons_u, cons_v, this->computingInfo_.reductionMethod);
      }
      else //local constraints, correction of subsystems
      {
        Bc.toScalar();
        for (it_mustm it = this->begin(); it != this->end(); it++)
        {
          const Unknown* u=it->first.first, *v=it->first.second;
          SuTermMatrix* sutm= it->second;
          SuTermVector* sutv= Bc.subVector_p(u);
          if (sutv==0) sutv= Bc.subVector_p(u->dual_p());   //try with dual
          if (this->constraints_u_p !=0) cons_u = (*(this->constraints_u_p))(u);
          if (this->constraints_v_p !=0)
          {
            cons_v = (*(this->constraints_v_p))(v);
            if (cons_v == 0) cons_v = (*(this->constraints_v_p))(v->dual_p());   //try with dual unknown
          }
          if (cons_u!=0 && cons_v!=0 && sutv!=0)
            appliedRhsCorrectorTo(sutv->scalar_entries(), sutv->cdofs(), sutm->rhs_matrix(),
                                  cons_u, cons_v, this->computingInfo_.reductionMethod);
        }
        Bc.toGlobal();
        this->toGlobal(_noStorage, _noAccess, _noSymmetry, false);
        Bc.adjustScalarEntries(this->cdofsr());             //adjust vector to matrix column scalar numbering
      }
    }
    else //no constraints to appply
    {
      this->toGlobal(_noStorage, _noAccess, _noSymmetry, false);
      Bc.toGlobal();
      Bc.adjustScalarEntries(this->cdofsr());             //adjust vector to matrix column scalar numbering
    }

    //goto solvers
    X.toGlobal(false);
    VectorEntry* x=X.scalar_entries();
    VectorEntry* b=Bc.scalar_entries();
    MatrixEntry* mat=this->scalar_entries();

    //goto solvers
    switch(sType)
    {
      case _diagSorS    :
        mat->sorDiagonalSolve(*b, *x, w);
        break;
      case _upperSorS :
        mat->sorUpperSolve(*b, *x, w);
        break;
      case _lowerSorS     :
        mat->sorLowerSolve(*b, *x, w);
        break;
      default :
        error("undef_option", "SorSolverType");
    }
    X.computed() = true;
    X.toLocal(false);
  }
  trace_p->pop();
}

// -----------------------------------------------------------------------------------------------------------
// in/out utilities
// -----------------------------------------------------------------------------------------------------------

//print storage on stream
void TermMatrix::viewStorage(std::ostream& out) const
{
  for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)
  {
    out << "SuTermMatrix " << it->second->name() << "\n";
    it->second->viewStorage(out);
  }
}

//print TermMatrix
void TermMatrix::print(std::ostream& out) const
{
  if(theVerboseLevel == 0) return;
  out << "TermMatrix " << name() << " : ";
  if(suTerms_.size() == 0)
    {
      out << " " << words("void") << " ";
      return;
    }
  cit_mustm it = suTerms_.begin();
  if(suTerms_.size() == 1)
    {
      out << words("unknown") << " \"" << it->first.first->name() << "\" ";
      out << words("test function") << " \"" << it->first.second->name() << "\"\n";
    }
  else
    {
      out << "("<<words("unknowns")<<", " <<  words("test functions")<<") : ";
      for(it = suTerms_.begin(); it != suTerms_.end(); it++)
        out << "(" << it->first.first->name() << ", "<< it->first.second->name() << ") ";
      out << eol;
    }
  if(theVerboseLevel < 2) return;

  //print SuTermMatrix's
  number_t vb = theVerboseLevel;
  if(theVerboseLevel < 5) verboseLevel(2);
  bool header = (suTerms_.size() > 1);
  for(it = suTerms_.begin(); it != suTerms_.end(); it++) it->second->print(out, header);
  verboseLevel(vb);

  //print global entries
  if(entries_p != 0)  out << " global entries : "<<*entries_p<<eol;
  if(scalar_entries_p != 0)
    {
      out << " global scalar entries : "<<eol;
      out << words("firste") << " " << words("row") << " cdofs : ";
      number_t nbd=cdofs_r.size();
      number_t m = std::min(number_t(theVerboseLevel),nbd);
      std::vector<DofComponent>::const_iterator itd=cdofs_r.begin();
      for(number_t i = 0; i < m; i++, itd++)  out  << *itd << "  ";
      if(m < nbd) out<<"...";
      out<<eol;
      out << words("firste") << " " << words("column") << " cdofs : ";
      nbd=cdofs_c.size();
      m = std::min(number_t(theVerboseLevel),nbd);
      itd=cdofs_c.begin();
      for(number_t i = 0; i < m; i++, itd++)  out  << *itd << "  ";
      if(m < nbd) out<<"...";
      out<<eol;
      if(scalar_entries_p==entries_p) out<<" same values as global entries"<<eol;
      else  out << (*scalar_entries_p)<<eol;
    }

  //print SetOfConstraints
  if(constraints_u_p!=0)
    out << " constraints on unknown (column) : "<<*constraints_u_p;
  if(constraints_v_p!=0)
    {
      out << "  constraints on test fuction (row)  : ";
      if(constraints_v_p == constraints_u_p) out<<" same as constraints on column"<<eol;
      else out <<*constraints_v_p<<eol;
    }
  if(constraints_u_p!=0 || constraints_v_p!=0) out<<"  "<<words("reduction method", computingInfo_.reductionMethod.method)<<eol;
}

//print brief computation report
void TermMatrix::printSummary(std::ostream& out) const
{
  out << "TermMatrix " << name() << " computed, size "<<numberOfRows()<< " X "<<numberOfCols()<<" : ";
  if(!isSingleUnknown()) out<<eol;
  for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++) it->second->printSummary(out);
  if(scalar_entries_p!=0) out<<eol<<"    has global representation in "<< scalar_entries_p->storagep()->name();
  if(rhs_matrix_p!=0) out<<eol<<"    has an essential condition correction matrix "<< rhs_matrix_p->nbOfRows()<<" X "<<rhs_matrix_p->nbOfCols();
  out << eol;
}

// save TermMatrix to file (dense or Matlab format)
void TermMatrix::saveToFile(const string_t& fn, StorageType st, bool encode) const
{
  if(!computed() || suTerms_.size() == 0)
    {
      warning("free_warning", "matrix is not computed, no saving");
      return;
    }

  //save global representation if exists, priority to scalar representation
  if(scalar_entries_p!=0)
    {
      scalar_entries_p->saveToFile(fn,st,encode);
      return;
    }
  if(entries_p!=0)
    {
      entries_p->saveToFile(fn,st,encode);
      return;
    }

  //save SuTermmMatrix's
  if(suTerms_.size() == 0)
    {
      warning("free_warning", "matrix seems to be empty, no saving");
      return;
    }

  cit_mustm it = suTerms_.begin();
  if(isSingleUnknown())
    {
      it->second->saveToFile(fn, st, encode);
      return;
    }

  std::pair<string_t, string_t> rootext = fileRootExtension(fn);
  for(cit_mustm it = suTerms_.begin(); it != suTerms_.end(); it++)
    {
      string_t fnb = rootext.first + "_" + it->first.first->name() + "_" + it->first.second->name() + "."+rootext.second;
      it->second->saveToFile(fnb, st, encode);
    }
}

} // end of namespace xlifepp