computation/interfaces/direct/Mumps.cxx

// Copyright (C) 2003-2009 Marc Duruflé
//
// This file is part of the linear-algebra library Seldon,
// http://seldon.sourceforge.net/.
//
// Seldon is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 2.1 of the License, or (at your option)
// any later version.
//
// Seldon 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 Lesser General Public License for
// more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with Seldon. If not, see http://www.gnu.org/licenses/.


#ifndef SELDON_FILE_MUMPS_CXX

#include "Mumps.hxx"

namespace Seldon
{

  template<>
  inline void MatrixMumps<double>::CallMumps()
  {
    dmumps_c(&struct_mumps);
  }


  template<>
  inline void MatrixMumps<complex<double> >::CallMumps()
  {
    zmumps_c(&struct_mumps);
  }


  template<class T>
  inline MatrixMumps<T>::MatrixMumps()
  {
    // struct_mumps.comm_fortran = MPI_Comm_c2f(MPI_COMM_WORLD);
    struct_mumps.comm_fortran = -987654;

    // parameters for mumps
    struct_mumps.job = -1;
    struct_mumps.par = 1;
    struct_mumps.sym = 0; // 0 -> unsymmetric matrix

    // mumps is called
    CallMumps();

    // other parameters
    struct_mumps.n = 0;
    type_ordering = 7; // default : we let Mumps choose the ordering
    print_level = -1;
    out_of_core = false;
    new_communicator = false;
  }


  template<class T> template<class MatrixSparse>
  inline void MatrixMumps<T>
  ::InitMatrix(const MatrixSparse& A, bool distributed)
  {
    // we clear previous factorization
    Clear();

#ifdef SELDON_WITH_MPI
    // MPI initialization for parallel version
    if (distributed)
      {
        // for distributed matrix, every processor is assumed to be involved
        struct_mumps.comm_fortran = -987654;
      }
    else
      {
        // centralized matrix => a linear system per processor
        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
        MPI_Comm_group(MPI_COMM_WORLD, &single_group);
        MPI_Group_incl(single_group, 1, &rank, &single_group);
        MPI_Comm_create(MPI_COMM_WORLD, single_group, &single_comm);
        struct_mumps.comm_fortran = MPI_Comm_c2f(single_comm);
        new_communicator = true;
      }
#endif

    // symmetry is specified during the initialization stage
    struct_mumps.job = -1;
    if (IsSymmetricMatrix(A))
      struct_mumps.sym = 2; // general symmetric matrix
    else
      struct_mumps.sym = 0; // unsymmetric matrix

    // mumps is called
    CallMumps();

    struct_mumps.icntl[13] = 20;
    struct_mumps.icntl[6] = type_ordering;
    // setting out of core parameters
    if (out_of_core)
      struct_mumps.icntl[21] = 1;
    else
      struct_mumps.icntl[21] = 0;

    struct_mumps.icntl[17] = 0;

    // the print level is set in mumps
    if (print_level >= 0)
      {
        struct_mumps.icntl[0] = 6;
        struct_mumps.icntl[1] = 0;
        struct_mumps.icntl[2] = 6;
        struct_mumps.icntl[3] = 2;
      }
    else
      {
        struct_mumps.icntl[0] = -1;
        struct_mumps.icntl[1] = -1;
        struct_mumps.icntl[2] = -1;
        struct_mumps.icntl[3] = 0;
      }
  }


  template<class T>
  inline void MatrixMumps<T>::SelectOrdering(int num_ordering)
  {
    type_ordering = num_ordering;
  }


  template<class T>
  MatrixMumps<T>::~MatrixMumps()
  {
    Clear();
  }


  template<class T>
  inline void MatrixMumps<T>::Clear()
  {
    if (struct_mumps.n > 0)
      {
        struct_mumps.job = -2;
        CallMumps(); /* Terminate instance */

        num_row_glob.Clear(); num_col_glob.Clear();
        struct_mumps.n = 0;

      }

#ifdef SELDON_WITH_MPI
    if (new_communicator)
      {
        MPI_Comm_free(&single_comm);
        new_communicator = false;
      }
#endif

  }


  template<class T>
  inline void MatrixMumps<T>::HideMessages()
  {
    print_level = -1;

    struct_mumps.icntl[0] = -1;
    struct_mumps.icntl[1] = -1;
    struct_mumps.icntl[2] = -1;
    struct_mumps.icntl[3] = 0;

  }


  template<class T>
  inline void MatrixMumps<T>::ShowMessages()
  {
    print_level = 0;

    struct_mumps.icntl[0] = 6;
    struct_mumps.icntl[1] = 0;
    struct_mumps.icntl[2] = 6;
    struct_mumps.icntl[3] = 2;

  }


  template<class T>
  inline void MatrixMumps<T>::EnableOutOfCore()
  {
    out_of_core = true;
  }


  template<class T>
  inline void MatrixMumps<T>::DisableOutOfCore()
  {
    out_of_core = false;
  }



  template<class T> template<class Prop,class Storage,class Allocator>
  void MatrixMumps<T>::FindOrdering(Matrix<T, Prop, Storage, Allocator> & mat,
                                    IVect& numbers, bool keep_matrix)
  {
    InitMatrix(mat);

    int n = mat.GetM(), nnz = mat.GetNonZeros();
    // conversion in coordinate format
    IVect num_row, num_col; Vector<T, VectFull, Allocator> values;
    ConvertMatrix_to_Coordinates(mat, num_row, num_col, values, 1);
    // no values needed to renumber
    values.Clear();
    if (!keep_matrix)
      mat.Clear();

    struct_mumps.n = n; struct_mumps.nz = nnz;
    struct_mumps.irn = num_row.GetData();
    struct_mumps.jcn = num_col.GetData();

    /* Call the MUMPS package. */
    struct_mumps.job = 1; // we analyse the system
    CallMumps();

    numbers.Reallocate(n);
    for (int i = 0; i < n; i++)
      numbers(i) = struct_mumps.sym_perm[i]-1;
  }



  template<class T> template<class Prop, class Storage, class Allocator>
  void MatrixMumps<T>::FactorizeMatrix(Matrix<T,Prop,Storage,Allocator> & mat,
                                       bool keep_matrix)
  {
    InitMatrix(mat);

    int n = mat.GetM(), nnz = mat.GetNonZeros();
    // conversion in coordinate format with fortran convention (1-index)
    IVect num_row, num_col; Vector<T, VectFull, Allocator> values;
    ConvertMatrix_to_Coordinates(mat, num_row, num_col, values, 1);
    if (!keep_matrix)
      mat.Clear();

    struct_mumps.n = n; struct_mumps.nz = nnz;
    struct_mumps.irn = num_row.GetData();
    struct_mumps.jcn = num_col.GetData();
    struct_mumps.a = reinterpret_cast<pointer>(values.GetData());

    /* Call the MUMPS package. */
    struct_mumps.job = 4; // we analyse and factorize the system
    CallMumps();
  }


  template<class T> template<class Prop, class Storage, class Allocator>
  void MatrixMumps<T>
  ::PerformAnalysis(Matrix<T, Prop, Storage, Allocator> & mat)
  {
    InitMatrix(mat);

    int n = mat.GetM(), nnz = mat.GetNonZeros();
    // conversion in coordinate format with fortran convention (1-index)
    Vector<T, VectFull, Allocator> values;
    ConvertMatrix_to_Coordinates(mat, num_row_glob, num_col_glob, values, 1);

    struct_mumps.n = n; struct_mumps.nz = nnz;
    struct_mumps.irn = num_row_glob.GetData();
    struct_mumps.jcn = num_col_glob.GetData();
    struct_mumps.a = reinterpret_cast<pointer>(values.GetData());

    /* Call the MUMPS package. */
    struct_mumps.job = 1; // we analyse the system
    CallMumps();
  }



  template<class T> template<class Prop, class Storage, class Allocator>
  void MatrixMumps<T>
  ::PerformFactorization(Matrix<T, Prop, Storage, Allocator> & mat)
  {
    // we consider that the values are corresponding
    // to the row/column numbers given for the analysis
    struct_mumps.a = reinterpret_cast<pointer>(mat.GetData());

    /* Call the MUMPS package. */
    struct_mumps.job = 2; // we factorize the system
    CallMumps();
  }


  template<class T>
  int MatrixMumps<T>::GetInfoFactorization() const
  {
    return struct_mumps.info[0];
  }



  template<class T> template<class Prop1, class Storage1, class Allocator,
                             class Prop2, class Storage2, class Allocator2>
  void MatrixMumps<T>::
  GetSchurMatrix(Matrix<T, Prop1, Storage1, Allocator>& mat, const IVect& num,
                 Matrix<T, Prop2, Storage2, Allocator2> & mat_schur,
                 bool keep_matrix)
  {
    InitMatrix(mat);

    int n_schur = num.GetM(), n = mat.GetM();
    // Subscripts are changed to respect fortran convention
    IVect index_schur(n_schur);
    for (int i = 0; i < n_schur; i++)
      index_schur(i) = num(i)+1;

    // array that will contain values of Schur matrix
    Vector<T, VectFull, Allocator2> vec_schur(n_schur*n_schur);

    struct_mumps.icntl[18] = n_schur;
    struct_mumps.size_schur = n_schur;
    struct_mumps.listvar_schur = index_schur.GetData();
    struct_mumps.schur = reinterpret_cast<pointer>(vec_schur.GetData());

    // factorization of the matrix
    FactorizeMatrix(mat, keep_matrix);

    // resetting parameters related to Schur complement
    struct_mumps.icntl[18] = 0;
    struct_mumps.size_schur = 0;
    struct_mumps.listvar_schur = NULL;
    struct_mumps.schur = NULL;

    // schur complement stored by rows
    int nb = 0;
    mat_schur.Reallocate(n,n);
    for (int i = 0; i < n; i++)
      for (int j = 0; j < n ;j++)
        mat_schur(i,j) = vec_schur(nb++);

    vec_schur.Clear(); index_schur.Clear();
  }



  template<class T> template<class Allocator2, class Transpose_status>
  void MatrixMumps<T>::Solve(const Transpose_status& TransA,
                             Vector<T, VectFull, Allocator2>& x)
  {
#ifdef SELDON_CHECK_DIMENSIONS
    if (x.GetM() != struct_mumps.n)
      throw WrongDim("Mumps::Solve(TransA, c)",
                     string("The length of x is equal to ")
                     + to_str(x.GetM())
                     + " while the size of the matrix is equal to "
                     + to_str(struct_mumps.n) + ".");
#endif

    if (TransA.Trans())
      struct_mumps.icntl[8] = 0;
    else
      struct_mumps.icntl[8] = 1;

    struct_mumps.rhs = reinterpret_cast<pointer>(x.GetData());
    struct_mumps.job = 3; // we solve system
    CallMumps();
  }



  template<class T> template<class Allocator2>
  void MatrixMumps<T>::Solve(Vector<T, VectFull, Allocator2>& x)
  {
    Solve(SeldonNoTrans, x);
  }



  template<class T>
  template<class Allocator2, class Transpose_status, class Prop>
  void MatrixMumps<T>::Solve(const Transpose_status& TransA,
                             Matrix<T, Prop, ColMajor, Allocator2>& x)
  {

#ifdef SELDON_CHECK_BOUNDS
    if (x.GetM() != struct_mumps.n)
      throw WrongDim("Mumps::Solve", string("Row size of x is equal to ")
                     + to_str(x.GetM()) + " while size of matrix is equal to "
                     + to_str(struct_mumps.n));
#endif

    if (TransA.Trans())
      struct_mumps.icntl[8] = 0;
    else
      struct_mumps.icntl[8] = 1;

    struct_mumps.nrhs = x.GetN();
    struct_mumps.lrhs = x.GetM();
    struct_mumps.rhs = reinterpret_cast<pointer>(x.GetData());
    struct_mumps.job = 3; // we solve system
    CallMumps();
  }


#ifdef SELDON_WITH_MPI


  template<class T> template<class Prop, class Allocator>
  void MatrixMumps<T>::
  FactorizeDistributedMatrix(Matrix<T, General, ColSparse, Allocator> & mat,
                             const Prop& sym, const IVect& glob_number,
                             bool keep_matrix)
  {
    // initialization depending on symmetric of the matrix
    Matrix<T, Prop, RowSparse, Allocator> Atest;
    InitMatrix(Atest, true);

    // distributed matrix
    struct_mumps.icntl[17] = 3;

    // global number of rows : mat.GetM()
    int N = mat.GetM();
    int nnz = mat.GetNonZeros();
    // conversion in coordinate format with C-convention (0-index)
    IVect num_row, num_col; Vector<T, VectFull, Allocator> values;
    ConvertMatrix_to_Coordinates(mat, num_row,
                                 num_col, values, 0);

    // we replace num_col with global numbers
    for (int i = 0; i < num_row.GetM(); i++)
      {
        num_row(i)++;
        num_col(i) = glob_number(num_col(i)) + 1;
      }

    if (!keep_matrix)
      mat.Clear();

    /* Define the problem on the host */
    struct_mumps.n = N; struct_mumps.nz_loc = nnz;
    struct_mumps.irn_loc = num_row.GetData();
    struct_mumps.jcn_loc = num_col.GetData();
    struct_mumps.a_loc = reinterpret_cast<pointer>(values.GetData());

    /* Call the MUMPS package. */
    struct_mumps.job = 4; // we analyse and factorize the system
    CallMumps();
    cout<<"Factorization completed"<<endl;
  }


  template<class Alloc1, class Alloc2, class Alloc3, class Tint>
  void FactorizeDistributedMatrix(Vector<int, VectFull, Alloc1>&,
                                  Vector<int, VectFull, Alloc2>&,
                                  Vector<T, VectFull, Alloc3>&,
                                  const Vector<Tint>& glob_number,
                                  bool sym, bool keep_matrix)
  {
    throw Undefined("FactorizeDistributedMatrix(Vector<int>&, Vector<int>&, "
                    "Vector<T>&, Vector<Tint>, bool, bool)");
  }



  template<class T> template<class Allocator2, class Transpose_status>
  void MatrixMumps<T>::SolveDistributed(const Transpose_status& TransA,
                                        Vector<T, VectFull, Allocator2>& x,
                                        const IVect& glob_num)
  {
    Vector<T, VectFull, Allocator2> rhs;
    int cplx = sizeof(T)/8;
    // allocating the global right hand side
    rhs.Reallocate(struct_mumps.n); rhs.Zero();

    if (rank == 0)
      {
        // on the host, we retrieve datas of all the other processors
        int nb_procs; MPI_Status status;
        MPI_Comm_size(MPI_COMM_WORLD, &nb_procs);
        if (nb_procs > 1)
          {
            // assembling the right hand side
            Vector<T, VectFull, Allocator2> xp;
            IVect nump;
            for (int i = 0; i < nb_procs; i++)
              {

                if (i != 0)
                  {
                    int nb_dof;
                    MPI_Recv(&nb_dof, 1, MPI_INT, i, 34,
                             MPI_COMM_WORLD, &status);
                    xp.Reallocate(nb_dof);
                    nump.Reallocate(nb_dof);
                    MPI_Recv(xp.GetDataVoid(), cplx*nb_dof,
                             MPI_DOUBLE, i, 35, MPI_COMM_WORLD, &status);
                    MPI_Recv(nump.GetData(), nb_dof, MPI_INT, i, 36,
                             MPI_COMM_WORLD, &status);
                  }
                else
                  {
                    xp = x; nump = glob_num;
                  }

                for (int j = 0; j < nump.GetM(); j++)
                  rhs(nump(j)) = xp(j);
              }
          }
        else
          Copy(x, rhs);

        struct_mumps.rhs = reinterpret_cast<pointer>(rhs.GetData());
      }
    else
      {
        // on other processors, we send solution
        int nb = x.GetM();
        MPI_Send(&nb, 1, MPI_INT, 0, 34, MPI_COMM_WORLD);
        MPI_Send(x.GetDataVoid(), cplx*nb, MPI_DOUBLE, 0, 35, MPI_COMM_WORLD);
        MPI_Send(glob_num.GetData(), nb, MPI_INT, 0, 36, MPI_COMM_WORLD);
      }

    // we solve system
    if (TransA.Trans())
      struct_mumps.icntl[8] = 0;
    else
      struct_mumps.icntl[8] = 1;

    struct_mumps.job = 3;
    CallMumps();

    // we distribute solution on all the processors
    MPI_Bcast(rhs.GetDataVoid(), cplx*rhs.GetM(),
              MPI_DOUBLE, 0, MPI_COMM_WORLD);

    // and we extract the solution on provided numbers
    for (int i = 0; i < x.GetM(); i++)
      x(i) = rhs(glob_num(i));
  }


  template<class T> template<class Allocator2>
  void MatrixMumps<T>::SolveDistributed(Vector<T, VectFull, Allocator2>& x,
                                        const IVect& glob_num)
  {
    SolveDistributed(SeldonNoTrans, x, glob_num);
  }
#endif


  template<class T, class Storage, class Allocator>
  void GetLU(Matrix<T,Symmetric,Storage,Allocator>& A, MatrixMumps<T>& mat_lu,
             bool keep_matrix = false)
  {
    mat_lu.FactorizeMatrix(A, keep_matrix);
  }


  template<class T, class Storage, class Allocator>
  void GetLU(Matrix<T,General,Storage,Allocator>& A, MatrixMumps<T>& mat_lu,
             bool keep_matrix = false)
  {
    mat_lu.FactorizeMatrix(A, keep_matrix);
  }


  template<class T, class Storage, class Allocator, class MatrixFull>
  void GetSchurMatrix(Matrix<T, Symmetric, Storage, Allocator>& A,
                      MatrixMumps<T>& mat_lu, const IVect& num,
                      MatrixFull& schur_matrix, bool keep_matrix = false)
  {
    mat_lu.GetSchurMatrix(A, num, schur_matrix, keep_matrix);
  }


  template<class T, class Storage, class Allocator, class MatrixFull>
  void GetSchurMatrix(Matrix<T, General, Storage, Allocator>& A,
                      MatrixMumps<T>& mat_lu, const IVect& num,
                      MatrixFull& schur_matrix, bool keep_matrix = false)
  {
    mat_lu.GetSchurMatrix(A, num, schur_matrix, keep_matrix);
  }


  template<class T, class Allocator>
  void SolveLU(MatrixMumps<T>& mat_lu, Vector<T, VectFull, Allocator>& x)
  {
    mat_lu.Solve(x);
  }


  template<class T, class Allocator, class Transpose_status>
  void SolveLU(const Transpose_status& TransA,
               MatrixMumps<T>& mat_lu, Vector<T, VectFull, Allocator>& x)
  {
    mat_lu.Solve(TransA, x);
  }


  template<class T, class Prop, class Allocator>
  void SolveLU(MatrixMumps<T>& mat_lu,
               Matrix<T, Prop, ColMajor, Allocator>& x)
  {
    mat_lu.Solve(SeldonNoTrans, x);
  }


  template<class T, class Allocator, class Transpose_status>
  void SolveLU(const Transpose_status& TransA,
               MatrixMumps<T>& mat_lu, Matrix<T, ColMajor, Allocator>& x)
  {
    mat_lu.Solve(TransA, x);
  }

}

#define SELDON_FILE_MUMPS_CXX
#endif