Reference documentation for deal.II version 9.0.0
scalapack.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_scalapack_h
17 #define dealii_scalapack_h
18 
19 #include <deal.II/base/config.h>
20 
21 #ifdef DEAL_II_WITH_SCALAPACK
22 
23 #include <deal.II/base/exceptions.h>
24 #include <deal.II/base/process_grid.h>
25 #include <deal.II/lac/full_matrix.h>
26 #include <deal.II/lac/lapack_support.h>
27 #include <deal.II/base/mpi.h>
28 #include <deal.II/base/thread_management.h>
29 #include <mpi.h>
30 
31 #include <memory>
32 
33 DEAL_II_NAMESPACE_OPEN
34 
104 template <typename NumberType>
105 class ScaLAPACKMatrix : protected TransposeTable<NumberType>
106 {
107 public:
108 
112  typedef unsigned int size_type;
113 
119  const size_type n_columns,
120  const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
121  const size_type row_block_size = 32,
122  const size_type column_block_size = 32,
124 
130  const std::shared_ptr<const Utilities::MPI::ProcessGrid> process_grid,
131  const size_type block_size = 32,
132  const LAPACKSupport::Property property = LAPACKSupport::Property::symmetric);
133 
137  ~ScaLAPACKMatrix() = default;
138 
143  void reinit(const size_type n_rows,
144  const size_type n_columns,
145  const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
146  const size_type row_block_size = 32,
147  const size_type column_block_size = 32,
149 
153  void reinit(const size_type size,
154  const std::shared_ptr<const Utilities::MPI::ProcessGrid> process_grid,
155  const size_type block_size = 32,
156  const LAPACKSupport::Property property = LAPACKSupport::Property::symmetric);
157 
162 
167 
172 
181 
188  void copy_to (FullMatrix<NumberType> &matrix) const;
189 
195  void copy_to (ScaLAPACKMatrix<NumberType> &dest) const;
196 
218  const std::pair<unsigned int,unsigned int> &offset_A,
219  const std::pair<unsigned int,unsigned int> &offset_B,
220  const std::pair<unsigned int,unsigned int> &submatrix_size) const;
221 
230 
241  void add(const ScaLAPACKMatrix<NumberType> &B,
242  const NumberType a=0.,
243  const NumberType b=1.,
244  const bool transpose_B=false);
245 
254  void add(const NumberType b,
255  const ScaLAPACKMatrix<NumberType> &B);
256 
265  void Tadd(const NumberType b,
266  const ScaLAPACKMatrix<NumberType> &B);
267 
285  void mult(const NumberType b,
287  const NumberType c,
289  const bool transpose_A=false,
290  const bool transpose_B=false) const;
291 
309  const bool adding=false) const;
310 
328  const bool adding=false) const;
329 
347  const bool adding=false) const;
348 
367  const bool adding=false) const;
368 
387  void save(const char *filename,
388  const std::pair<unsigned int,unsigned int> &chunk_size=std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int)) const;
389 
401  void load(const char *filename);
402 
408 
414  void compute_lu_factorization ();
415 
426  void invert();
427 
442  std::vector<NumberType> eigenpairs_symmetric_by_index(const std::pair<unsigned int,unsigned int> &index_limits,
443  const bool compute_eigenvectors);
444 
453  std::vector<NumberType> eigenpairs_symmetric_by_value(const std::pair<NumberType,NumberType> &value_limits,
454  const bool compute_eigenvectors);
455 
470  std::vector<NumberType> eigenpairs_symmetric_by_index_MRRR(const std::pair<unsigned int,unsigned int> &index_limits,
471  const bool compute_eigenvectors);
472 
482  std::vector<NumberType> eigenpairs_symmetric_by_value_MRRR(const std::pair<NumberType,NumberType> &value_limits,
483  const bool compute_eigenvectors);
484 
508  std::vector<NumberType> compute_SVD(ScaLAPACKMatrix<NumberType> *U = nullptr,
509  ScaLAPACKMatrix<NumberType> *VT = nullptr);
510 
542  const bool transpose=false);
543 
562  unsigned int pseudoinverse(const NumberType ratio);
563 
577  NumberType reciprocal_condition_number(const NumberType a_norm) const;
578 
582  NumberType l1_norm() const;
583 
587  NumberType linfty_norm() const;
588 
592  NumberType frobenius_norm() const;
593 
597  size_type m() const;
598 
602  size_type n() const;
603 
607  unsigned int local_m() const;
608 
612  unsigned int local_n() const;
613 
617  unsigned int global_row(const unsigned int loc_row) const;
618 
622  unsigned int global_column(const unsigned int loc_column) const;
623 
627  NumberType local_el(const unsigned int loc_row, const unsigned int loc_column) const;
628 
632  NumberType &local_el(const unsigned int loc_row, const unsigned int loc_column);
633 
644  template <class InputVector>
645  void scale_columns(const InputVector &factors);
646 
657  template <class InputVector>
658  void scale_rows(const InputVector &factors);
659 
660 private:
661 
666  NumberType norm_symmetric(const char type) const;
667 
672  NumberType norm_general(const char type) const;
673 
684  std::vector<NumberType> eigenpairs_symmetric(const bool compute_eigenvectors,
685  const std::pair<unsigned int,unsigned int> &index_limits=
687  const std::pair<NumberType,NumberType> &value_limits=
688  std::make_pair(std::numeric_limits<NumberType>::quiet_NaN(),std::numeric_limits<NumberType>::quiet_NaN()));
689 
708  std::vector<NumberType> eigenpairs_symmetric_MRRR(const bool compute_eigenvectors,
709  const std::pair<unsigned int,unsigned int> &index_limits=
711  const std::pair<NumberType,NumberType> &value_limits=
712  std::make_pair(std::numeric_limits<NumberType>::quiet_NaN(),std::numeric_limits<NumberType>::quiet_NaN()));
713 
714  /*
715  * Stores the distributed matrix in @p filename
716  * using serial routines
717  */
718  void save_serial(const char *filename,
719  const std::pair<unsigned int,unsigned int> &chunk_size) const;
720 
721  /*
722  * Loads the distributed matrix from file @p filename
723  * using serial routines
724  */
725  void load_serial(const char *filename);
726 
727  /*
728  * Stores the distributed matrix in @p filename
729  * using parallel routines
730  */
731  void save_parallel(const char *filename,
732  const std::pair<unsigned int,unsigned int> &chunk_size) const;
733 
734  /*
735  * Loads the distributed matrix from file @p filename
736  * using parallel routines
737  */
738  void load_parallel(const char *filename);
739 
745 
751 
756  std::shared_ptr<const Utilities::MPI::ProcessGrid> grid;
757 
761  int n_rows;
762 
767 
772 
777 
782 
787 
791  int descriptor[9];
792 
796  mutable std::vector<NumberType> work;
797 
801  mutable std::vector<int> iwork;
802 
807  std::vector<int> ipiv;
808 
813  const char uplo;
814 
819  const int first_process_row;
820 
826 
831  const int submatrix_row;
832 
837  const int submatrix_column;
838 
843 };
844 
845 // ----------------------- inline functions ----------------------------
846 
847 #ifndef DOXYGEN
848 
849 template <typename NumberType>
850 inline
851 NumberType
852 ScaLAPACKMatrix<NumberType>::local_el(const unsigned int loc_row, const unsigned int loc_column) const
853 {
854  return (*this)(loc_row,loc_column);
855 }
856 
857 
858 
859 template <typename NumberType>
860 inline
861 NumberType &
862 ScaLAPACKMatrix<NumberType>::local_el(const unsigned int loc_row, const unsigned int loc_column)
863 {
864  return (*this)(loc_row,loc_column);
865 }
866 
867 
868 template <typename NumberType>
869 inline
870 unsigned int
872 {
873  return n_rows;
874 }
875 
876 
877 
878 template <typename NumberType>
879 inline
880 unsigned int
882 {
883  return n_columns;
884 }
885 
886 
887 
888 template <typename NumberType>
889 unsigned int
891 {
892  return n_local_rows;
893 }
894 
895 
896 
897 template <typename NumberType>
898 unsigned int
900 {
901  return n_local_columns;
902 }
903 
904 
905 #endif // DOXYGEN
906 
907 DEAL_II_NAMESPACE_CLOSE
908 
909 #endif // DEAL_II_WITH_SCALAPACK
910 
911 #endif
std::vector< NumberType > eigenpairs_symmetric_by_index_MRRR(const std::pair< unsigned int, unsigned int > &index_limits, const bool compute_eigenvectors)
Definition: scalapack.cc:1004
void copy_transposed(const ScaLAPACKMatrix< NumberType > &B)
Definition: scalapack.cc:499
ScaLAPACKMatrix(const size_type n_rows, const size_type n_columns, const std::shared_ptr< const Utilities::MPI::ProcessGrid > &process_grid, const size_type row_block_size=32, const size_type column_block_size=32, const LAPACKSupport::Property property=LAPACKSupport::Property::general)
Definition: scalapack.cc:73
std::vector< NumberType > eigenpairs_symmetric_by_value_MRRR(const std::pair< NumberType, NumberType > &value_limits, const bool compute_eigenvectors)
Definition: scalapack.cc:1024
void mmult(ScaLAPACKMatrix< NumberType > &C, const ScaLAPACKMatrix< NumberType > &B, const bool adding=false) const
Definition: scalapack.cc:634
static const unsigned int invalid_unsigned_int
Definition: types.h:173
unsigned int size_type
Definition: scalapack.h:112
LAPACKSupport::State state
Definition: scalapack.h:744
unsigned int global_row(const unsigned int loc_row) const
Definition: scalapack.cc:229
NumberType l1_norm() const
Definition: scalapack.cc:1392
void Tmmult(ScaLAPACKMatrix< NumberType > &C, const ScaLAPACKMatrix< NumberType > &B, const bool adding=false) const
Definition: scalapack.cc:647
ScaLAPACKMatrix< NumberType > & operator=(const FullMatrix< NumberType > &)
Definition: scalapack.cc:199
const char uplo
Definition: scalapack.h:813
unsigned int pseudoinverse(const NumberType ratio)
Definition: scalapack.cc:1308
void load(const char *filename)
Definition: scalapack.cc:1916
const TableIndices< N > & size() const
std::vector< NumberType > eigenpairs_symmetric(const bool compute_eigenvectors, const std::pair< unsigned int, unsigned int > &index_limits=std::make_pair(numbers::invalid_unsigned_int, numbers::invalid_unsigned_int), const std::pair< NumberType, NumberType > &value_limits=std::make_pair(std::numeric_limits< NumberType >::quiet_NaN(), std::numeric_limits< NumberType >::quiet_NaN()))
Definition: scalapack.cc:820
std::vector< NumberType > eigenpairs_symmetric_by_index(const std::pair< unsigned int, unsigned int > &index_limits, const bool compute_eigenvectors)
Definition: scalapack.cc:785
size_type n_rows() const
~ScaLAPACKMatrix()=default
void scale_columns(const InputVector &factors)
Definition: scalapack.cc:2303
std::shared_ptr< const Utilities::MPI::ProcessGrid > grid
Definition: scalapack.h:756
std::vector< int > iwork
Definition: scalapack.h:801
size_type n() const
const int submatrix_column
Definition: scalapack.h:837
void reinit(const size_type n_rows, const size_type n_columns, const std::shared_ptr< const Utilities::MPI::ProcessGrid > &process_grid, const size_type row_block_size=32, const size_type column_block_size=32, const LAPACKSupport::Property property=LAPACKSupport::Property::general)
Definition: scalapack.cc:104
void TmTmult(ScaLAPACKMatrix< NumberType > &C, const ScaLAPACKMatrix< NumberType > &B, const bool adding=false) const
Definition: scalapack.cc:673
void least_squares(ScaLAPACKMatrix< NumberType > &B, const bool transpose=false)
Definition: scalapack.cc:1254
void mTmult(ScaLAPACKMatrix< NumberType > &C, const ScaLAPACKMatrix< NumberType > &B, const bool adding=false) const
Definition: scalapack.cc:660
void scale_rows(const InputVector &factors)
Definition: scalapack.cc:2313
NumberType norm_symmetric(const char type) const
Definition: scalapack.cc:1466
void set_property(const LAPACKSupport::Property property)
Definition: scalapack.cc:172
std::vector< int > ipiv
Definition: scalapack.h:807
const int first_process_row
Definition: scalapack.h:819
LAPACKSupport::State get_state() const
Definition: scalapack.cc:190
unsigned int global_column(const unsigned int loc_column) const
Definition: scalapack.cc:241
void add(const ScaLAPACKMatrix< NumberType > &B, const NumberType a=0., const NumberType b=1., const bool transpose_B=false)
Definition: scalapack.cc:507
std::vector< NumberType > eigenpairs_symmetric_by_value(const std::pair< NumberType, NumberType > &value_limits, const bool compute_eigenvectors)
Definition: scalapack.cc:805
int column_block_size
Definition: scalapack.h:776
NumberType norm_general(const char type) const
Definition: scalapack.cc:1431
Threads::Mutex mutex
Definition: scalapack.h:842
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
const int submatrix_row
Definition: scalapack.h:831
void save(const char *filename, const std::pair< unsigned int, unsigned int > &chunk_size=std::make_pair(numbers::invalid_unsigned_int, numbers::invalid_unsigned_int)) const
Definition: scalapack.cc:1572
unsigned int local_m() const
std::vector< NumberType > eigenpairs_symmetric_MRRR(const bool compute_eigenvectors, const std::pair< unsigned int, unsigned int > &index_limits=std::make_pair(numbers::invalid_unsigned_int, numbers::invalid_unsigned_int), const std::pair< NumberType, NumberType > &value_limits=std::make_pair(std::numeric_limits< NumberType >::quiet_NaN(), std::numeric_limits< NumberType >::quiet_NaN()))
Definition: scalapack.cc:1039
size_type m() const
std::vector< NumberType > compute_SVD(ScaLAPACKMatrix< NumberType > *U=nullptr, ScaLAPACKMatrix< NumberType > *VT=nullptr)
Definition: scalapack.cc:1180
std::vector< NumberType > work
Definition: scalapack.h:796
void mult(const NumberType b, const ScaLAPACKMatrix< NumberType > &B, const NumberType c, ScaLAPACKMatrix< NumberType > &C, const bool transpose_A=false, const bool transpose_B=false) const
Definition: scalapack.cc:562
void copy_to(FullMatrix< NumberType > &matrix) const
Definition: scalapack.cc:253
unsigned int local_n() const
const int first_process_column
Definition: scalapack.h:825
LAPACKSupport::Property property
Definition: scalapack.h:750
void Tadd(const NumberType b, const ScaLAPACKMatrix< NumberType > &B)
Definition: scalapack.cc:553
NumberType local_el(const unsigned int loc_row, const unsigned int loc_column) const
NumberType reciprocal_condition_number(const NumberType a_norm) const
Definition: scalapack.cc:1356
NumberType linfty_norm() const
Definition: scalapack.cc:1405
int descriptor[9]
Definition: scalapack.h:791
void compute_cholesky_factorization()
Definition: scalapack.cc:686
NumberType frobenius_norm() const
Definition: scalapack.cc:1418
LAPACKSupport::Property get_property() const
Definition: scalapack.cc:181
void compute_lu_factorization()
Definition: scalapack.cc:708