19 #ifdef DEAL_II_WITH_SCALAPACK
23 # include <deal.II/base/mpi.templates.h>
25 # include <deal.II/lac/scalapack.templates.h>
27 # ifdef DEAL_II_WITH_HDF5
36 # ifdef DEAL_II_WITH_HDF5
38 template <
typename number>
50 return H5T_NATIVE_DOUBLE;
56 return H5T_NATIVE_FLOAT;
62 return H5T_NATIVE_INT;
68 return H5T_NATIVE_UINT;
74 return H5T_NATIVE_CHAR;
80 template <
typename NumberType>
84 const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
91 , first_process_column(0)
105 template <
typename NumberType>
108 const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
121 template <
typename NumberType>
123 const std::string & filename,
124 const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
130 , first_process_column(0)
132 , submatrix_column(1)
134 # ifndef DEAL_II_WITH_HDF5
142 "This function is only available when deal.II is configured with HDF5"));
156 hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
160 hid_t dataset = H5Dopen2(file,
"/matrix", H5P_DEFAULT);
164 hid_t filespace = H5Dget_space(dataset);
167 int rank = H5Sget_simple_extent_ndims(filespace);
170 status = H5Sget_simple_extent_dims(filespace, dims,
nullptr);
179 status = H5Sclose(filespace);
181 status = H5Dclose(dataset);
183 status = H5Fclose(file);
186 int ierr = MPI_Bcast(&
n_rows,
190 process_grid->mpi_communicator);
197 process_grid->mpi_communicator);
208 load(filename.c_str());
215 template <
typename NumberType>
220 const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
225 Assert(row_block_size_ > 0,
ExcMessage(
"Row block size has to be positive."));
226 Assert(column_block_size_ > 0,
227 ExcMessage(
"Column block size has to be positive."));
229 row_block_size_ <= n_rows_,
231 "Row block size can not be greater than the number of rows of the matrix"));
233 column_block_size_ <= n_columns_,
235 "Column block size can not be greater than the number of columns of the matrix"));
238 property = property_;
241 n_columns = n_columns_;
242 row_block_size = row_block_size_;
243 column_block_size = column_block_size_;
245 if (grid->mpi_process_is_active)
248 n_local_rows = numroc_(&n_rows,
250 &(grid->this_process_row),
252 &(grid->n_process_rows));
253 n_local_columns = numroc_(&n_columns,
255 &(grid->this_process_column),
256 &first_process_column,
257 &(grid->n_process_columns));
261 int lda =
std::max(1, n_local_rows);
264 descinit_(descriptor,
270 &first_process_column,
271 &(grid->blacs_context),
282 n_local_columns = -1;
289 template <
typename NumberType>
293 const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
302 template <
typename NumberType>
307 property = property_;
312 template <
typename NumberType>
321 template <
typename NumberType>
330 template <
typename NumberType>
343 if (grid->mpi_process_is_active)
345 for (
int i = 0; i < n_local_rows; ++i)
347 const int glob_i = global_row(i);
348 for (
int j = 0; j < n_local_columns; ++j)
350 const int glob_j = global_column(j);
351 local_el(i, j) =
matrix(glob_i, glob_j);
361 template <
typename NumberType>
364 const unsigned int rank)
366 if (n_rows * n_columns == 0)
374 ExcMessage(
"All processes have to call routine with identical rank"));
376 ExcMessage(
"All processes have to call routine with identical rank"));
391 MPI_Comm_group(this->grid->mpi_communicator, &group_A);
393 const std::vector<int> ranks(n, rank);
395 MPI_Group_incl(group_A, n, ranks.data(), &group_B);
396 MPI_Comm communicator_B;
399 const int ierr = MPI_Comm_create_group(this->grid->mpi_communicator,
404 int n_proc_rows_B = 1, n_proc_cols_B = 1;
405 int this_process_row_B = -1, this_process_column_B = -1;
406 int blacs_context_B = -1;
407 if (MPI_COMM_NULL != communicator_B)
410 blacs_context_B = Csys2blacs_handle(communicator_B);
411 const char *order =
"Col";
412 Cblacs_gridinit(&blacs_context_B, order, n_proc_rows_B, n_proc_cols_B);
413 Cblacs_gridinfo(blacs_context_B,
417 &this_process_column_B);
423 const bool mpi_process_is_active_B =
424 (this_process_row_B >= 0 && this_process_column_B >= 0);
427 std::vector<int> descriptor_B(9, -1);
428 const int first_process_row_B = 0, first_process_col_B = 0;
430 if (mpi_process_is_active_B)
433 int n_local_rows_B = numroc_(&n_rows,
436 &first_process_row_B,
438 int n_local_cols_B = numroc_(&n_columns,
440 &this_process_column_B,
441 &first_process_col_B,
445 (void)n_local_cols_B;
447 int lda =
std::max(1, n_local_rows_B);
449 descinit_(descriptor_B.data(),
454 &first_process_row_B,
455 &first_process_col_B,
461 if (this->grid->mpi_process_is_active)
464 NumberType *loc_vals_A =
465 this->
values.size() > 0 ? this->
values.data() :
nullptr;
466 const NumberType *loc_vals_B =
467 mpi_process_is_active_B ? &(B(0, 0)) :
nullptr;
481 &(this->grid->blacs_context));
483 if (mpi_process_is_active_B)
484 Cblacs_gridexit(blacs_context_B);
486 MPI_Group_free(&group_A);
487 MPI_Group_free(&group_B);
488 if (MPI_COMM_NULL != communicator_B)
496 template <
typename NumberType>
500 Assert(n_local_rows >= 0 && loc_row <
static_cast<unsigned int>(n_local_rows),
502 const int i = loc_row + 1;
505 &(grid->this_process_row),
507 &(grid->n_process_rows)) -
513 template <
typename NumberType>
517 Assert(n_local_columns >= 0 &&
518 loc_column <
static_cast<unsigned int>(n_local_columns),
520 const int j = loc_column + 1;
523 &(grid->this_process_column),
524 &first_process_column,
525 &(grid->n_process_columns)) -
531 template <
typename NumberType>
534 const unsigned int rank)
const
536 if (n_rows * n_columns == 0)
544 ExcMessage(
"All processes have to call routine with identical rank"));
546 ExcMessage(
"All processes have to call routine with identical rank"));
562 MPI_Comm_group(this->grid->mpi_communicator, &group_A);
564 const std::vector<int> ranks(n, rank);
566 MPI_Group_incl(group_A, n, ranks.data(), &group_B);
567 MPI_Comm communicator_B;
570 const int ierr = MPI_Comm_create_group(this->grid->mpi_communicator,
575 int n_proc_rows_B = 1, n_proc_cols_B = 1;
576 int this_process_row_B = -1, this_process_column_B = -1;
577 int blacs_context_B = -1;
578 if (MPI_COMM_NULL != communicator_B)
581 blacs_context_B = Csys2blacs_handle(communicator_B);
582 const char *order =
"Col";
583 Cblacs_gridinit(&blacs_context_B, order, n_proc_rows_B, n_proc_cols_B);
584 Cblacs_gridinfo(blacs_context_B,
588 &this_process_column_B);
594 const bool mpi_process_is_active_B =
595 (this_process_row_B >= 0 && this_process_column_B >= 0);
598 std::vector<int> descriptor_B(9, -1);
599 const int first_process_row_B = 0, first_process_col_B = 0;
601 if (mpi_process_is_active_B)
604 int n_local_rows_B = numroc_(&n_rows,
607 &first_process_row_B,
609 int n_local_cols_B = numroc_(&n_columns,
611 &this_process_column_B,
612 &first_process_col_B,
616 (void)n_local_cols_B;
618 int lda =
std::max(1, n_local_rows_B);
621 descinit_(descriptor_B.data(),
626 &first_process_row_B,
627 &first_process_col_B,
635 if (this->grid->mpi_process_is_active)
638 const NumberType *loc_vals_A =
639 this->
values.size() > 0 ? this->
values.data() :
nullptr;
640 NumberType *loc_vals_B = mpi_process_is_active_B ? &(B(0, 0)) :
nullptr;
652 &(this->grid->blacs_context));
654 if (mpi_process_is_active_B)
655 Cblacs_gridexit(blacs_context_B);
657 MPI_Group_free(&group_A);
658 MPI_Group_free(&group_B);
659 if (MPI_COMM_NULL != communicator_B)
665 template <
typename NumberType>
678 if (grid->mpi_process_is_active)
680 for (
int i = 0; i < n_local_rows; ++i)
682 const int glob_i = global_row(i);
683 for (
int j = 0; j < n_local_columns; ++j)
685 const int glob_j = global_column(j);
686 matrix(glob_i, glob_j) = local_el(i, j);
698 for (
unsigned int i = 0; i <
matrix.n(); ++i)
699 for (
unsigned int j = i + 1; j <
matrix.m(); ++j)
702 for (
unsigned int i = 0; i <
matrix.n(); ++i)
703 for (
unsigned int j = 0; j < i; ++j)
710 for (
unsigned int i = 0; i <
matrix.n(); ++i)
711 for (
unsigned int j = i + 1; j <
matrix.m(); ++j)
713 else if (uplo ==
'U')
714 for (
unsigned int i = 0; i <
matrix.n(); ++i)
715 for (
unsigned int j = 0; j < i; ++j)
722 template <
typename NumberType>
726 const std::pair<unsigned int, unsigned int> &offset_A,
727 const std::pair<unsigned int, unsigned int> &offset_B,
728 const std::pair<unsigned int, unsigned int> &submatrix_size)
const
731 if (submatrix_size.first == 0 || submatrix_size.second == 0)
744 int ierr, comparison;
745 ierr = MPI_Comm_compare(grid->mpi_communicator,
746 B.
grid->mpi_communicator,
749 Assert(comparison == MPI_IDENT,
750 ExcMessage(
"Matrix A and B must have a common MPI Communicator"));
758 int union_blacs_context = Csys2blacs_handle(this->grid->mpi_communicator);
759 const char *order =
"Col";
760 int union_n_process_rows =
762 int union_n_process_columns = 1;
763 Cblacs_gridinit(&union_blacs_context,
765 union_n_process_rows,
766 union_n_process_columns);
768 int n_grid_rows_A, n_grid_columns_A, my_row_A, my_column_A;
769 Cblacs_gridinfo(this->grid->blacs_context,
776 const bool in_context_A =
777 (my_row_A >= 0 && my_row_A < n_grid_rows_A) &&
778 (my_column_A >= 0 && my_column_A < n_grid_columns_A);
780 int n_grid_rows_B, n_grid_columns_B, my_row_B, my_column_B;
781 Cblacs_gridinfo(B.
grid->blacs_context,
788 const bool in_context_B =
789 (my_row_B >= 0 && my_row_B < n_grid_rows_B) &&
790 (my_column_B >= 0 && my_column_B < n_grid_columns_B);
792 const int n_rows_submatrix = submatrix_size.first;
793 const int n_columns_submatrix = submatrix_size.second;
796 int ia = offset_A.first + 1, ja = offset_A.second + 1;
797 int ib = offset_B.first + 1, jb = offset_B.second + 1;
799 std::array<int, 9> desc_A, desc_B;
801 const NumberType *loc_vals_A =
nullptr;
802 NumberType * loc_vals_B =
nullptr;
811 if (this->
values.size() != 0)
812 loc_vals_A = this->
values.data();
814 for (
unsigned int i = 0; i < desc_A.size(); ++i)
815 desc_A[i] = this->descriptor[i];
823 loc_vals_B = B.
values.data();
825 for (
unsigned int i = 0; i < desc_B.size(); ++i)
831 pgemr2d(&n_rows_submatrix,
832 &n_columns_submatrix,
841 &union_blacs_context);
846 Cblacs_gridexit(union_blacs_context);
851 template <
typename NumberType>
859 if (this->grid->mpi_process_is_active)
861 this->descriptor[0] == 1,
863 "Copying of ScaLAPACK matrices only implemented for dense matrices"));
864 if (dest.
grid->mpi_process_is_active)
868 "Copying of ScaLAPACK matrices only implemented for dense matrices"));
884 MPI_Group group_source, group_dest, group_union;
885 ierr = MPI_Comm_group(this->grid->mpi_communicator, &group_source);
887 ierr = MPI_Comm_group(dest.
grid->mpi_communicator, &group_dest);
889 ierr = MPI_Group_union(group_source, group_dest, &group_union);
891 MPI_Comm mpi_communicator_union;
906 ierr = MPI_Comm_create_group(MPI_COMM_WORLD,
909 &mpi_communicator_union);
917 int union_blacs_context = Csys2blacs_handle(mpi_communicator_union);
918 const char *order =
"Col";
919 int union_n_process_rows =
921 int union_n_process_columns = 1;
922 Cblacs_gridinit(&union_blacs_context,
924 union_n_process_rows,
925 union_n_process_columns);
927 const NumberType *loc_vals_source =
nullptr;
928 NumberType * loc_vals_dest =
nullptr;
930 if (this->grid->mpi_process_is_active && (this->values.size() > 0))
934 "source: process is active but local matrix empty"));
935 loc_vals_source = this->
values.data();
937 if (dest.
grid->mpi_process_is_active && (dest.
values.size() > 0))
942 "destination: process is active but local matrix empty"));
943 loc_vals_dest = dest.
values.data();
955 &union_blacs_context);
957 Cblacs_gridexit(union_blacs_context);
959 if (mpi_communicator_union != MPI_COMM_NULL)
961 ierr = MPI_Group_free(&group_source);
963 ierr = MPI_Group_free(&group_dest);
965 ierr = MPI_Group_free(&group_union);
970 if (this->grid->mpi_process_is_active)
979 template <
typename NumberType>
989 template <
typename NumberType>
992 const NumberType alpha,
993 const NumberType beta,
994 const bool transpose_B)
1016 ExcMessage(
"The matrices A and B need to have the same process grid"));
1018 if (this->grid->mpi_process_is_active)
1020 char trans_b = transpose_B ?
'T' :
'N';
1022 (this->
values.size() > 0) ? this->
values.data() :
nullptr;
1023 const NumberType *B_loc =
1045 template <
typename NumberType>
1050 add(B, 1, a,
false);
1055 template <
typename NumberType>
1065 template <
typename NumberType>
1071 const bool transpose_A,
1072 const bool transpose_B)
const
1075 ExcMessage(
"The matrices A and B need to have the same process grid"));
1077 ExcMessage(
"The matrices B and C need to have the same process grid"));
1081 if (!transpose_A && !transpose_B)
1085 Assert(this->n_rows ==
C.n_rows,
1089 Assert(this->row_block_size ==
C.row_block_size,
1096 else if (transpose_A && !transpose_B)
1100 Assert(this->n_columns ==
C.n_rows,
1104 Assert(this->column_block_size ==
C.row_block_size,
1111 else if (!transpose_A && transpose_B)
1115 Assert(this->n_rows ==
C.n_rows,
1119 Assert(this->row_block_size ==
C.row_block_size,
1131 Assert(this->n_columns ==
C.n_rows,
1135 Assert(this->column_block_size ==
C.row_block_size,
1143 if (this->grid->mpi_process_is_active)
1145 char trans_a = transpose_A ?
'T' :
'N';
1146 char trans_b = transpose_B ?
'T' :
'N';
1148 const NumberType *A_loc =
1149 (this->
values.size() > 0) ? this->
values.data() :
nullptr;
1150 const NumberType *B_loc =
1152 NumberType *C_loc = (
C.values.size() > 0) ?
C.values.data() :
nullptr;
1154 int n =
C.n_columns;
1155 int k = transpose_A ? this->n_rows : this->n_columns;
1164 &(this->submatrix_row),
1165 &(this->submatrix_column),
1174 &
C.submatrix_column,
1182 template <
typename NumberType>
1186 const bool adding)
const
1189 mult(1., B, 1.,
C,
false,
false);
1191 mult(1., B, 0,
C,
false,
false);
1196 template <
typename NumberType>
1200 const bool adding)
const
1203 mult(1., B, 1.,
C,
true,
false);
1205 mult(1., B, 0,
C,
true,
false);
1210 template <
typename NumberType>
1214 const bool adding)
const
1217 mult(1., B, 1.,
C,
false,
true);
1219 mult(1., B, 0,
C,
false,
true);
1224 template <
typename NumberType>
1228 const bool adding)
const
1231 mult(1., B, 1.,
C,
true,
true);
1233 mult(1., B, 0,
C,
true,
true);
1238 template <
typename NumberType>
1245 "Cholesky factorization can be applied to symmetric matrices only."));
1248 "Matrix has to be in Matrix state before calling this function."));
1250 if (grid->mpi_process_is_active)
1253 NumberType *A_loc = this->
values.data();
1271 template <
typename NumberType>
1277 "Matrix has to be in Matrix state before calling this function."));
1279 if (grid->mpi_process_is_active)
1282 NumberType *A_loc = this->
values.data();
1284 const int iarow = indxg2p_(&submatrix_row,
1286 &(grid->this_process_row),
1288 &(grid->n_process_rows));
1289 const int mp = numroc_(&n_rows,
1291 &(grid->this_process_row),
1293 &(grid->n_process_rows));
1294 ipiv.resize(mp + row_block_size);
1312 template <
typename NumberType>
1330 if (grid->mpi_process_is_active)
1332 const char uploTriangular =
1334 const char diag =
'N';
1336 NumberType *A_loc = this->
values.data();
1337 ptrtri(&uploTriangular,
1359 compute_cholesky_factorization();
1361 compute_lu_factorization();
1363 if (grid->mpi_process_is_active)
1366 NumberType *A_loc = this->
values.data();
1383 int lwork = -1, liwork = -1;
1401 lwork =
static_cast<int>(work[0]);
1404 iwork.resize(liwork);
1428 template <
typename NumberType>
1429 std::vector<NumberType>
1431 const std::pair<unsigned int, unsigned int> &index_limits,
1432 const bool compute_eigenvectors)
1438 std::pair<unsigned int, unsigned int> idx =
1439 std::make_pair(
std::min(index_limits.first, index_limits.second),
1440 std::max(index_limits.first, index_limits.second));
1443 if (idx.first == 0 && idx.second ==
static_cast<unsigned int>(n_rows - 1))
1444 return eigenpairs_symmetric(compute_eigenvectors);
1446 return eigenpairs_symmetric(compute_eigenvectors, idx);
1451 template <
typename NumberType>
1452 std::vector<NumberType>
1454 const std::pair<NumberType, NumberType> &value_limits,
1455 const bool compute_eigenvectors)
1457 Assert(!std::isnan(value_limits.first),
1459 Assert(!std::isnan(value_limits.second),
1462 std::pair<unsigned int, unsigned int> indices =
1466 return eigenpairs_symmetric(compute_eigenvectors, indices, value_limits);
1471 template <
typename NumberType>
1472 std::vector<NumberType>
1474 const bool compute_eigenvectors,
1475 const std::pair<unsigned int, unsigned int> &eigenvalue_idx,
1476 const std::pair<NumberType, NumberType> & eigenvalue_limits)
1480 "Matrix has to be in Matrix state before calling this function."));
1482 ExcMessage(
"Matrix has to be symmetric for this operation."));
1484 std::lock_guard<std::mutex> lock(mutex);
1486 const bool use_values = (std::isnan(eigenvalue_limits.first) ||
1487 std::isnan(eigenvalue_limits.second)) ?
1490 const bool use_indices =
1497 !(use_values && use_indices),
1499 "Prescribing both the index and value range for the eigenvalues is ambiguous"));
1503 std::unique_ptr<ScaLAPACKMatrix<NumberType>>
eigenvectors =
1504 compute_eigenvectors ?
1505 std::make_unique<ScaLAPACKMatrix<NumberType>>(n_rows,
1509 grid->n_process_rows, grid->n_process_columns, grid, 1, 1);
1516 std::vector<NumberType> ev(n_rows);
1518 if (grid->mpi_process_is_active)
1525 char jobz = compute_eigenvectors ?
'V' :
'N';
1528 bool all_eigenpairs =
true;
1529 NumberType vl = NumberType(), vu = NumberType();
1535 NumberType abstol = NumberType();
1542 NumberType orfac = 0;
1544 std::vector<int> ifail;
1551 std::vector<int> iclustr;
1557 std::vector<NumberType> gap(n_local_rows * n_local_columns);
1567 all_eigenpairs =
true;
1572 all_eigenpairs =
false;
1573 vl =
std::min(eigenvalue_limits.first, eigenvalue_limits.second);
1574 vu =
std::max(eigenvalue_limits.first, eigenvalue_limits.second);
1580 all_eigenpairs =
false;
1583 il =
std::min(eigenvalue_idx.first, eigenvalue_idx.second) + 1;
1584 iu =
std::max(eigenvalue_idx.first, eigenvalue_idx.second) + 1;
1586 NumberType *A_loc = this->
values.data();
1593 NumberType *eigenvectors_loc =
1594 (compute_eigenvectors ?
eigenvectors->values.data() :
nullptr);
1619 char cmach = compute_eigenvectors ?
'U' :
'S';
1620 plamch(&(this->grid->blacs_context), &cmach, abstol);
1622 ifail.resize(n_rows);
1623 iclustr.resize(2 * grid->n_process_rows * grid->n_process_columns);
1624 gap.resize(grid->n_process_rows * grid->n_process_columns);
1657 lwork =
static_cast<int>(work[0]);
1684 iwork.resize(liwork);
1721 if (compute_eigenvectors)
1722 this->
values.swap(eigenvectors->values);
1725 while (ev.size() >
static_cast<size_type>(m))
1731 grid->send_to_inactive(&m, 1);
1736 if (!grid->mpi_process_is_active)
1741 grid->send_to_inactive(ev.data(), ev.size());
1748 if (compute_eigenvectors)
1761 template <
typename NumberType>
1762 std::vector<NumberType>
1764 const std::pair<unsigned int, unsigned int> &index_limits,
1765 const bool compute_eigenvectors)
1771 const std::pair<unsigned int, unsigned int> idx =
1772 std::make_pair(
std::min(index_limits.first, index_limits.second),
1773 std::max(index_limits.first, index_limits.second));
1776 if (idx.first == 0 && idx.second ==
static_cast<unsigned int>(n_rows - 1))
1777 return eigenpairs_symmetric_MRRR(compute_eigenvectors);
1779 return eigenpairs_symmetric_MRRR(compute_eigenvectors, idx);
1784 template <
typename NumberType>
1785 std::vector<NumberType>
1787 const std::pair<NumberType, NumberType> &value_limits,
1788 const bool compute_eigenvectors)
1793 const std::pair<unsigned int, unsigned int> indices =
1797 return eigenpairs_symmetric_MRRR(compute_eigenvectors, indices, value_limits);
1802 template <
typename NumberType>
1803 std::vector<NumberType>
1805 const bool compute_eigenvectors,
1806 const std::pair<unsigned int, unsigned int> &eigenvalue_idx,
1807 const std::pair<NumberType, NumberType> & eigenvalue_limits)
1811 "Matrix has to be in Matrix state before calling this function."));
1813 ExcMessage(
"Matrix has to be symmetric for this operation."));
1815 std::lock_guard<std::mutex> lock(mutex);
1817 const bool use_values = (std::isnan(eigenvalue_limits.first) ||
1818 std::isnan(eigenvalue_limits.second)) ?
1821 const bool use_indices =
1828 !(use_values && use_indices),
1830 "Prescribing both the index and value range for the eigenvalues is ambiguous"));
1834 std::unique_ptr<ScaLAPACKMatrix<NumberType>>
eigenvectors =
1835 compute_eigenvectors ?
1836 std::make_unique<ScaLAPACKMatrix<NumberType>>(n_rows,
1840 grid->n_process_rows, grid->n_process_columns, grid, 1, 1);
1846 std::vector<NumberType> ev(n_rows);
1853 if (grid->mpi_process_is_active)
1860 char jobz = compute_eigenvectors ?
'V' :
'N';
1864 NumberType vl = NumberType(), vu = NumberType();
1879 vl =
std::min(eigenvalue_limits.first, eigenvalue_limits.second);
1880 vu =
std::max(eigenvalue_limits.first, eigenvalue_limits.second);
1888 il =
std::min(eigenvalue_idx.first, eigenvalue_idx.second) + 1;
1889 iu =
std::max(eigenvalue_idx.first, eigenvalue_idx.second) + 1;
1891 NumberType *A_loc = this->
values.data();
1899 NumberType *eigenvectors_loc =
1900 (compute_eigenvectors ?
eigenvectors->values.data() :
nullptr);
1931 lwork =
static_cast<int>(work[0]);
1934 iwork.resize(liwork);
1963 if (compute_eigenvectors)
1967 "psyevr failed to compute all eigenvectors for the selected eigenvalues"));
1972 if (compute_eigenvectors)
1973 this->
values.swap(eigenvectors->values);
1976 while (ev.size() >
static_cast<size_type>(m))
1982 grid->send_to_inactive(&m, 1);
1987 if (!grid->mpi_process_is_active)
1992 grid->send_to_inactive(ev.data(), ev.size());
1999 if (compute_eigenvectors)
2012 template <
typename NumberType>
2013 std::vector<NumberType>
2019 "Matrix has to be in Matrix state before calling this function."));
2020 Assert(row_block_size == column_block_size,
2023 const bool left_singluar_vectors = (
U !=
nullptr) ?
true :
false;
2024 const bool right_singluar_vectors = (VT !=
nullptr) ?
true :
false;
2026 if (left_singluar_vectors)
2031 Assert(row_block_size ==
U->row_block_size,
2033 Assert(column_block_size ==
U->column_block_size,
2035 Assert(grid->blacs_context ==
U->grid->blacs_context,
2038 if (right_singluar_vectors)
2048 Assert(grid->blacs_context == VT->
grid->blacs_context,
2050 VT->
grid->blacs_context));
2052 std::lock_guard<std::mutex> lock(mutex);
2054 std::vector<NumberType> sv(
std::min(n_rows, n_columns));
2056 if (grid->mpi_process_is_active)
2058 char jobu = left_singluar_vectors ?
'V' :
'N';
2059 char jobvt = right_singluar_vectors ?
'V' :
'N';
2060 NumberType *A_loc = this->
values.data();
2061 NumberType *U_loc = left_singluar_vectors ?
U->values.data() :
nullptr;
2062 NumberType *VT_loc = right_singluar_vectors ? VT->
values.data() :
nullptr;
2082 &
U->submatrix_column,
2093 lwork =
static_cast<int>(work[0]);
2107 &
U->submatrix_column,
2122 grid->send_to_inactive(sv.data(), sv.size());
2132 template <
typename NumberType>
2138 ExcMessage(
"The matrices A and B need to have the same process grid"));
2141 "Matrix has to be in Matrix state before calling this function."));
2144 "Matrix B has to be in Matrix state before calling this function."));
2157 Assert(row_block_size == column_block_size,
2159 "Use identical block sizes for rows and columns of matrix A"));
2162 "Use identical block sizes for rows and columns of matrix B"));
2165 "Use identical block-cyclic distribution for matrices A and B"));
2167 std::lock_guard<std::mutex> lock(mutex);
2169 if (grid->mpi_process_is_active)
2172 NumberType *A_loc = this->
values.data();
2173 NumberType *B_loc = B.
values.data();
2199 lwork =
static_cast<int>(work[0]);
2224 template <
typename NumberType>
2230 "Matrix has to be in Matrix state before calling this function."));
2231 Assert(row_block_size == column_block_size,
2233 "Use identical block sizes for rows and columns of matrix A"));
2235 ratio > 0. && ratio < 1.,
2237 "input parameter ratio has to be larger than zero and smaller than 1"));
2251 std::vector<NumberType> sv = this->compute_SVD(&
U, &VT);
2259 unsigned int n_sv = 1;
2260 std::vector<NumberType> inv_sigma;
2261 inv_sigma.push_back(1 / sv[0]);
2263 for (
unsigned int i = 1; i < sv.size(); ++i)
2264 if (sv[i] > sv[0] * ratio)
2267 inv_sigma.push_back(1 / sv[i]);
2289 std::make_pair(0, 0),
2290 std::make_pair(0, 0),
2291 std::make_pair(n_rows, n_sv));
2293 std::make_pair(0, 0),
2294 std::make_pair(0, 0),
2295 std::make_pair(n_sv, n_columns));
2304 VT_R.
mult(1, U_R, 0, *
this,
true,
true);
2311 template <
typename NumberType>
2314 const NumberType a_norm)
const
2318 "Matrix has to be in Cholesky state before calling this function."));
2319 std::lock_guard<std::mutex> lock(mutex);
2320 NumberType rcond = 0.;
2322 if (grid->mpi_process_is_active)
2324 int liwork = n_local_rows;
2325 iwork.resize(liwork);
2328 const NumberType *A_loc = this->
values.data();
2348 lwork =
static_cast<int>(
std::ceil(work[0]));
2367 grid->send_to_inactive(&rcond);
2373 template <
typename NumberType>
2377 const char type(
'O');
2380 return norm_symmetric(type);
2382 return norm_general(type);
2387 template <
typename NumberType>
2391 const char type(
'I');
2394 return norm_symmetric(type);
2396 return norm_general(type);
2401 template <
typename NumberType>
2405 const char type(
'F');
2408 return norm_symmetric(type);
2410 return norm_general(type);
2415 template <
typename NumberType>
2421 ExcMessage(
"norms can be called in matrix state only."));
2422 std::lock_guard<std::mutex> lock(mutex);
2423 NumberType res = 0.;
2425 if (grid->mpi_process_is_active)
2427 const int iarow = indxg2p_(&submatrix_row,
2429 &(grid->this_process_row),
2431 &(grid->n_process_rows));
2432 const int iacol = indxg2p_(&submatrix_column,
2434 &(grid->this_process_column),
2435 &first_process_column,
2436 &(grid->n_process_columns));
2437 const int mp0 = numroc_(&n_rows,
2439 &(grid->this_process_row),
2441 &(grid->n_process_rows));
2442 const int nq0 = numroc_(&n_columns,
2444 &(grid->this_process_column),
2446 &(grid->n_process_columns));
2452 if (type ==
'O' || type ==
'1')
2454 else if (type ==
'I')
2458 const NumberType *A_loc = this->
values.begin();
2468 grid->send_to_inactive(&res);
2474 template <
typename NumberType>
2480 ExcMessage(
"norms can be called in matrix state only."));
2482 ExcMessage(
"Matrix has to be symmetric for this operation."));
2483 std::lock_guard<std::mutex> lock(mutex);
2484 NumberType res = 0.;
2486 if (grid->mpi_process_is_active)
2491 ilcm_(&(grid->n_process_rows), &(grid->n_process_columns));
2492 const int v2 = lcm / (grid->n_process_rows);
2494 const int IAROW = indxg2p_(&submatrix_row,
2496 &(grid->this_process_row),
2498 &(grid->n_process_rows));
2499 const int IACOL = indxg2p_(&submatrix_column,
2501 &(grid->this_process_column),
2502 &first_process_column,
2503 &(grid->n_process_columns));
2504 const int Np0 = numroc_(&n_columns ,
2506 &(grid->this_process_row),
2508 &(grid->n_process_rows));
2509 const int Nq0 = numroc_(&n_columns ,
2511 &(grid->this_process_column),
2513 &(grid->n_process_columns));
2515 const int v1 = iceil_(&Np0, &row_block_size);
2516 const int ldw = (n_local_rows == n_local_columns) ?
2518 row_block_size * iceil_(&
v1, &v2);
2521 (type ==
'M' || type ==
'F' || type ==
'E') ? 0 : 2 * Nq0 + Np0 + ldw;
2523 const NumberType *A_loc = this->
values.begin();
2533 grid->send_to_inactive(&res);
2539 # ifdef DEAL_II_WITH_HDF5
2545 create_HDF5_state_enum_id(hid_t &state_enum_id)
2551 herr_t status = H5Tenum_insert(state_enum_id,
"cholesky", &val);
2554 status = H5Tenum_insert(state_enum_id,
"eigenvalues", &val);
2557 status = H5Tenum_insert(state_enum_id,
"inverse_matrix", &val);
2560 status = H5Tenum_insert(state_enum_id,
"inverse_svd", &val);
2563 status = H5Tenum_insert(state_enum_id,
"lu", &val);
2566 status = H5Tenum_insert(state_enum_id,
"matrix", &val);
2569 status = H5Tenum_insert(state_enum_id,
"svd", &val);
2572 status = H5Tenum_insert(state_enum_id,
"unusable", &val);
2577 create_HDF5_property_enum_id(hid_t &property_enum_id)
2582 herr_t status = H5Tenum_insert(property_enum_id,
"diagonal", &prop);
2585 status = H5Tenum_insert(property_enum_id,
"general", &prop);
2588 status = H5Tenum_insert(property_enum_id,
"hessenberg", &prop);
2591 status = H5Tenum_insert(property_enum_id,
"lower_triangular", &prop);
2594 status = H5Tenum_insert(property_enum_id,
"symmetric", &prop);
2597 status = H5Tenum_insert(property_enum_id,
"upper_triangular", &prop);
2606 template <
typename NumberType>
2609 const std::string & filename,
2610 const std::pair<unsigned int, unsigned int> &
chunk_size)
const
2612 # ifndef DEAL_II_WITH_HDF5
2618 std::pair<unsigned int, unsigned int> chunks_size_ =
chunk_size;
2624 chunks_size_.first = n_rows;
2625 chunks_size_.second = 1;
2627 Assert(chunks_size_.first > 0,
2628 ExcMessage(
"The row chunk size must be larger than 0."));
2630 Assert(chunks_size_.second > 0,
2631 ExcMessage(
"The column chunk size must be larger than 0."));
2634 # ifdef H5_HAVE_PARALLEL
2636 save_parallel(filename, chunks_size_);
2640 save_serial(filename, chunks_size_);
2648 template <
typename NumberType>
2651 const std::string & filename,
2652 const std::pair<unsigned int, unsigned int> &
chunk_size)
const
2654 # ifndef DEAL_II_WITH_HDF5
2669 const auto column_grid =
2670 std::make_shared<Utilities::MPI::ProcessGrid>(this->grid->mpi_communicator,
2674 const int MB = n_rows, NB = n_columns;
2680 if (tmp.
grid->mpi_process_is_active)
2686 H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
2689 hsize_t chunk_dims[2];
2694 hid_t data_property = H5Pcreate(H5P_DATASET_CREATE);
2695 status = H5Pset_chunk(data_property, 2, chunk_dims);
2702 dims[0] = n_columns;
2704 hid_t dataspace_id = H5Screate_simple(2, dims,
nullptr);
2708 hid_t dataset_id = H5Dcreate2(file_id,
2718 dataset_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, tmp.
values.data());
2723 hid_t state_enum_id, property_enum_id;
2724 internal::create_HDF5_state_enum_id(state_enum_id);
2725 internal::create_HDF5_property_enum_id(property_enum_id);
2728 hsize_t dims_state[1];
2730 hid_t state_enum_dataspace = H5Screate_simple(1, dims_state,
nullptr);
2732 hid_t state_enum_dataset = H5Dcreate2(file_id,
2735 state_enum_dataspace,
2740 status = H5Dwrite(state_enum_dataset,
2749 hsize_t dims_property[1];
2750 dims_property[0] = 1;
2751 hid_t property_enum_dataspace =
2752 H5Screate_simple(1, dims_property,
nullptr);
2754 hid_t property_enum_dataset = H5Dcreate2(file_id,
2757 property_enum_dataspace,
2762 status = H5Dwrite(property_enum_dataset,
2771 status = H5Dclose(dataset_id);
2773 status = H5Dclose(state_enum_dataset);
2775 status = H5Dclose(property_enum_dataset);
2779 status = H5Sclose(dataspace_id);
2781 status = H5Sclose(state_enum_dataspace);
2783 status = H5Sclose(property_enum_dataspace);
2787 status = H5Tclose(state_enum_id);
2789 status = H5Tclose(property_enum_id);
2793 status = H5Pclose(data_property);
2797 status = H5Fclose(file_id);
2805 template <
typename NumberType>
2808 const std::string & filename,
2809 const std::pair<unsigned int, unsigned int> &
chunk_size)
const
2811 # ifndef DEAL_II_WITH_HDF5
2819 MPI_Info info = MPI_INFO_NULL;
2827 const auto column_grid =
2828 std::make_shared<Utilities::MPI::ProcessGrid>(this->grid->mpi_communicator,
2832 const int MB = n_rows;
2853 NumberType *data = (tmp.
values.size() > 0) ? tmp.
values.data() :
nullptr;
2860 hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
2861 status = H5Pset_fapl_mpio(plist_id, tmp.
grid->mpi_communicator, info);
2866 H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
2867 status = H5Pclose(plist_id);
2876 hid_t filespace = H5Screate_simple(2, dims,
nullptr);
2879 hsize_t chunk_dims[2];
2883 plist_id = H5Pcreate(H5P_DATASET_CREATE);
2884 H5Pset_chunk(plist_id, 2, chunk_dims);
2886 hid_t dset_id = H5Dcreate2(
2887 file_id,
"/matrix", type_id, filespace, H5P_DEFAULT, plist_id, H5P_DEFAULT);
2889 status = H5Sclose(filespace);
2892 status = H5Pclose(plist_id);
2901 proc_n_local_rows.data(),
2904 tmp.
grid->mpi_communicator);
2909 proc_n_local_columns.data(),
2912 tmp.
grid->mpi_communicator);
2915 const unsigned int my_rank(
2924 hid_t memspace = H5Screate_simple(2, count,
nullptr);
2926 hsize_t offset[2] = {0};
2927 for (
unsigned int i = 0; i < my_rank; ++i)
2928 offset[0] += proc_n_local_columns[i];
2931 filespace = H5Dget_space(dset_id);
2932 status = H5Sselect_hyperslab(
2933 filespace, H5S_SELECT_SET, offset,
nullptr, count,
nullptr);
2937 plist_id = H5Pcreate(H5P_DATASET_XFER);
2938 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT);
2942 if (tmp.
values.size() > 0)
2944 status = H5Dwrite(dset_id, type_id, memspace, filespace, plist_id, data);
2948 status = H5Dclose(dset_id);
2950 status = H5Sclose(filespace);
2952 status = H5Sclose(memspace);
2954 status = H5Pclose(plist_id);
2956 status = H5Fclose(file_id);
2961 ierr = MPI_Barrier(tmp.
grid->mpi_communicator);
2965 if (tmp.
grid->this_mpi_process == 0)
2968 hid_t file_id_reopen =
2969 H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
2973 hid_t state_enum_id, property_enum_id;
2974 internal::create_HDF5_state_enum_id(state_enum_id);
2975 internal::create_HDF5_property_enum_id(property_enum_id);
2978 hsize_t dims_state[1];
2980 hid_t state_enum_dataspace = H5Screate_simple(1, dims_state,
nullptr);
2982 hid_t state_enum_dataset = H5Dcreate2(file_id_reopen,
2985 state_enum_dataspace,
2990 status = H5Dwrite(state_enum_dataset,
2999 hsize_t dims_property[1];
3000 dims_property[0] = 1;
3001 hid_t property_enum_dataspace =
3002 H5Screate_simple(1, dims_property,
nullptr);
3004 hid_t property_enum_dataset = H5Dcreate2(file_id_reopen,
3007 property_enum_dataspace,
3012 status = H5Dwrite(property_enum_dataset,
3020 status = H5Dclose(state_enum_dataset);
3022 status = H5Dclose(property_enum_dataset);
3024 status = H5Sclose(state_enum_dataspace);
3026 status = H5Sclose(property_enum_dataspace);
3028 status = H5Tclose(state_enum_id);
3030 status = H5Tclose(property_enum_id);
3032 status = H5Fclose(file_id_reopen);
3041 template <
typename NumberType>
3045 # ifndef DEAL_II_WITH_HDF5
3049 # ifdef H5_HAVE_PARALLEL
3051 load_parallel(filename);
3055 load_serial(filename);
3062 template <
typename NumberType>
3066 # ifndef DEAL_II_WITH_HDF5
3077 const auto one_grid =
3078 std::make_shared<Utilities::MPI::ProcessGrid>(this->grid->mpi_communicator,
3082 const int MB = n_rows, NB = n_columns;
3086 int property_int = -1;
3090 if (tmp.
grid->mpi_process_is_active)
3095 hid_t file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
3098 hid_t dataset_id = H5Dopen2(file_id,
"/matrix", H5P_DEFAULT);
3104 hid_t datatype = H5Dget_type(dataset_id);
3105 H5T_class_t t_class_in = H5Tget_class(datatype);
3108 t_class_in == t_class,
3110 "The data type of the matrix to be read does not match the archive"));
3113 hid_t dataspace_id = H5Dget_space(dataset_id);
3115 const int ndims = H5Sget_simple_extent_ndims(dataspace_id);
3119 H5Sget_simple_extent_dims(dataspace_id, dims,
nullptr);
3121 static_cast<int>(dims[0]) == n_columns,
3123 "The number of columns of the matrix does not match the content of the archive"));
3125 static_cast<int>(dims[1]) == n_rows,
3127 "The number of rows of the matrix does not match the content of the archive"));
3130 status = H5Dread(dataset_id,
3140 hid_t state_enum_id, property_enum_id;
3141 internal::create_HDF5_state_enum_id(state_enum_id);
3142 internal::create_HDF5_property_enum_id(property_enum_id);
3145 hid_t dataset_state_id = H5Dopen2(file_id,
"/state", H5P_DEFAULT);
3146 hid_t datatype_state = H5Dget_type(dataset_state_id);
3147 H5T_class_t t_class_state = H5Tget_class(datatype_state);
3150 hid_t dataset_property_id = H5Dopen2(file_id,
"/property", H5P_DEFAULT);
3151 hid_t datatype_property = H5Dget_type(dataset_property_id);
3152 H5T_class_t t_class_property = H5Tget_class(datatype_property);
3156 hid_t dataspace_state = H5Dget_space(dataset_state_id);
3157 hid_t dataspace_property = H5Dget_space(dataset_property_id);
3159 const int ndims_state = H5Sget_simple_extent_ndims(dataspace_state);
3161 const int ndims_property = H5Sget_simple_extent_ndims(dataspace_property);
3164 hsize_t dims_state[1];
3165 H5Sget_simple_extent_dims(dataspace_state, dims_state,
nullptr);
3167 hsize_t dims_property[1];
3168 H5Sget_simple_extent_dims(dataspace_property, dims_property,
nullptr);
3172 status = H5Dread(dataset_state_id,
3182 state_int =
static_cast<int>(tmp.
state);
3184 status = H5Dread(dataset_property_id,
3194 property_int =
static_cast<int>(tmp.
property);
3197 status = H5Sclose(dataspace_id);
3199 status = H5Sclose(dataspace_state);
3201 status = H5Sclose(dataspace_property);
3205 status = H5Tclose(datatype);
3207 status = H5Tclose(state_enum_id);
3209 status = H5Tclose(property_enum_id);
3213 status = H5Dclose(dataset_state_id);
3215 status = H5Dclose(dataset_id);
3217 status = H5Dclose(dataset_property_id);
3221 status = H5Fclose(file_id);
3225 tmp.
grid->send_to_inactive(&state_int, 1);
3228 tmp.
grid->send_to_inactive(&property_int, 1);
3240 template <
typename NumberType>
3244 # ifndef DEAL_II_WITH_HDF5
3248 # ifndef H5_HAVE_PARALLEL
3254 MPI_Info info = MPI_INFO_NULL;
3261 const auto column_grid =
3262 std::make_shared<Utilities::MPI::ProcessGrid>(this->grid->mpi_communicator,
3266 const int MB = n_rows;
3275 NumberType *data = (tmp.
values.size() > 0) ? tmp.
values.data() :
nullptr;
3280 hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
3281 status = H5Pset_fapl_mpio(plist_id, tmp.
grid->mpi_communicator, info);
3286 hid_t file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, plist_id);
3287 status = H5Pclose(plist_id);
3291 hid_t dataset_id = H5Dopen2(file_id,
"/matrix", H5P_DEFAULT);
3298 hid_t datatype_inp = H5Dget_type(dataset_id);
3299 H5T_class_t t_class_inp = H5Tget_class(datatype_inp);
3300 H5T_class_t t_class = H5Tget_class(datatype);
3302 t_class_inp == t_class,
3304 "The data type of the matrix to be read does not match the archive"));
3308 hid_t dataspace_id = H5Dget_space(dataset_id);
3310 const int ndims = H5Sget_simple_extent_ndims(dataspace_id);
3314 status = H5Sget_simple_extent_dims(dataspace_id, dims,
nullptr);
3317 static_cast<int>(dims[0]) == n_columns,
3319 "The number of columns of the matrix does not match the content of the archive"));
3321 static_cast<int>(dims[1]) == n_rows,
3323 "The number of rows of the matrix does not match the content of the archive"));
3331 proc_n_local_rows.data(),
3334 tmp.
grid->mpi_communicator);
3339 proc_n_local_columns.data(),
3342 tmp.
grid->mpi_communicator);
3345 const unsigned int my_rank(
3355 hsize_t offset[2] = {0};
3356 for (
unsigned int i = 0; i < my_rank; ++i)
3357 offset[0] += proc_n_local_columns[i];
3360 status = H5Sselect_hyperslab(
3361 dataspace_id, H5S_SELECT_SET, offset,
nullptr, count,
nullptr);
3365 hid_t memspace = H5Screate_simple(2, count,
nullptr);
3369 H5Dread(dataset_id, datatype, memspace, dataspace_id, H5P_DEFAULT, data);
3373 hid_t state_enum_id, property_enum_id;
3374 internal::create_HDF5_state_enum_id(state_enum_id);
3375 internal::create_HDF5_property_enum_id(property_enum_id);
3378 hid_t dataset_state_id = H5Dopen2(file_id,
"/state", H5P_DEFAULT);
3379 hid_t datatype_state = H5Dget_type(dataset_state_id);
3380 H5T_class_t t_class_state = H5Tget_class(datatype_state);
3383 hid_t dataset_property_id = H5Dopen2(file_id,
"/property", H5P_DEFAULT);
3384 hid_t datatype_property = H5Dget_type(dataset_property_id);
3385 H5T_class_t t_class_property = H5Tget_class(datatype_property);
3389 hid_t dataspace_state = H5Dget_space(dataset_state_id);
3390 hid_t dataspace_property = H5Dget_space(dataset_property_id);
3392 const int ndims_state = H5Sget_simple_extent_ndims(dataspace_state);
3394 const int ndims_property = H5Sget_simple_extent_ndims(dataspace_property);
3397 hsize_t dims_state[1];
3398 H5Sget_simple_extent_dims(dataspace_state, dims_state,
nullptr);
3400 hsize_t dims_property[1];
3401 H5Sget_simple_extent_dims(dataspace_property, dims_property,
nullptr);
3406 dataset_state_id, state_enum_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, &tmp.
state);
3409 status = H5Dread(dataset_property_id,
3418 status = H5Sclose(memspace);
3420 status = H5Dclose(dataset_id);
3422 status = H5Dclose(dataset_state_id);
3424 status = H5Dclose(dataset_property_id);
3426 status = H5Sclose(dataspace_id);
3428 status = H5Sclose(dataspace_state);
3430 status = H5Sclose(dataspace_property);
3434 status = H5Tclose(state_enum_id);
3436 status = H5Tclose(property_enum_id);
3438 status = H5Fclose(file_id);
3454 template <
typename NumberType>
3462 for (
unsigned int i = 0; i <
matrix.local_n(); ++i)
3464 const NumberType s = factors[
matrix.global_column(i)];
3466 for (
unsigned int j = 0; j <
matrix.local_m(); ++j)
3467 matrix.local_el(j, i) *= s;
3471 template <
typename NumberType>
3479 for (
unsigned int i = 0; i <
matrix.local_m(); ++i)
3481 const NumberType s = factors[
matrix.global_row(i)];
3483 for (
unsigned int j = 0; j <
matrix.local_n(); ++j)
3484 matrix.local_el(i, j) *= s;
3493 template <
typename NumberType>
3494 template <
class InputVector>
3498 if (this->grid->mpi_process_is_active)
3504 template <
typename NumberType>
3505 template <
class InputVector>
3509 if (this->grid->mpi_process_is_active)
3516 # include "scalapack.inst"
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
std::vector< NumberType > compute_SVD(ScaLAPACKMatrix< NumberType > *U=nullptr, ScaLAPACKMatrix< NumberType > *VT=nullptr)
std::vector< NumberType > eigenpairs_symmetric_by_value(const std::pair< NumberType, NumberType > &value_limits, const bool compute_eigenvectors)
NumberType frobenius_norm() const
unsigned int pseudoinverse(const NumberType ratio)
std::vector< NumberType > eigenpairs_symmetric_by_value_MRRR(const std::pair< NumberType, NumberType > &value_limits, const bool compute_eigenvectors)
void copy_from(const LAPACKFullMatrix< NumberType > &matrix, const unsigned int rank)
void save_parallel(const std::string &filename, const std::pair< unsigned int, unsigned int > &chunk_size) const
void least_squares(ScaLAPACKMatrix< NumberType > &B, const bool transpose=false)
ScaLAPACKMatrix< NumberType > & operator=(const FullMatrix< NumberType > &)
void Tadd(const NumberType b, const ScaLAPACKMatrix< NumberType > &B)
void mTmult(ScaLAPACKMatrix< NumberType > &C, const ScaLAPACKMatrix< NumberType > &B, const bool adding=false) const
void add(const ScaLAPACKMatrix< NumberType > &B, const NumberType a=0., const NumberType b=1., const bool transpose_B=false)
LAPACKSupport::State get_state() const
LAPACKSupport::Property get_property() const
void Tmmult(ScaLAPACKMatrix< NumberType > &C, const ScaLAPACKMatrix< NumberType > &B, const bool adding=false) 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()))
void scale_rows(const InputVector &factors)
std::shared_ptr< const Utilities::MPI::ProcessGrid > grid
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)
void load(const std::string &filename)
const int submatrix_column
void mmult(ScaLAPACKMatrix< NumberType > &C, const ScaLAPACKMatrix< NumberType > &B, const bool adding=false) 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()))
void save_serial(const std::string &filename, const std::pair< unsigned int, unsigned int > &chunk_size) const
NumberType norm_general(const char type) const
void save(const std::string &filename, const std::pair< unsigned int, unsigned int > &chunk_size=std::make_pair(numbers::invalid_unsigned_int, numbers::invalid_unsigned_int)) const
void load_parallel(const std::string &filename)
NumberType l1_norm() const
void compute_lu_factorization()
NumberType norm_symmetric(const char type) const
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
LAPACKSupport::Property property
void set_property(const LAPACKSupport::Property property)
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)
void load_serial(const std::string &filename)
std::vector< NumberType > eigenpairs_symmetric_by_index_MRRR(const std::pair< unsigned int, unsigned int > &index_limits, const bool compute_eigenvectors)
NumberType reciprocal_condition_number(const NumberType a_norm) const
LAPACKSupport::State state
void TmTmult(ScaLAPACKMatrix< NumberType > &C, const ScaLAPACKMatrix< NumberType > &B, const bool adding=false) const
std::vector< NumberType > eigenpairs_symmetric_by_index(const std::pair< unsigned int, unsigned int > &index_limits, const bool compute_eigenvectors)
unsigned int global_column(const unsigned int loc_column) const
void copy_to(FullMatrix< NumberType > &matrix) const
unsigned int global_row(const unsigned int loc_row) const
void compute_cholesky_factorization()
void copy_transposed(const ScaLAPACKMatrix< NumberType > &B)
NumberType linfty_norm() const
void scale_columns(const InputVector &factors)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
AlignedVector< T > values
typename AlignedVector< T >::size_type size_type
void reinit(const size_type size1, const size_type size2, const bool omit_default_initialization=false)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNeedsHDF5()
static ::ExceptionBase & ExcErrorCode(std::string arg1, types::blas_int arg2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Expression ceil(const Expression &x)
@ cholesky
Contents is a Cholesky decomposition.
@ lu
Contents is an LU decomposition.
@ matrix
Contents is actually a matrix.
@ unusable
Contents is something useless.
@ inverse_matrix
Contents is the inverse of a matrix.
@ svd
Matrix contains singular value decomposition,.
@ inverse_svd
Matrix is the inverse of a singular value decomposition.
@ eigenvalues
Eigenvalue vector is filled.
@ symmetric
Matrix is symmetric.
@ hessenberg
Matrix is in upper Hessenberg form.
@ diagonal
Matrix is diagonal.
@ upper_triangular
Matrix is upper triangular.
@ lower_triangular
Matrix is lower triangular.
@ general
No special properties.
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
void free_communicator(MPI_Comm &mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
T min(const T &t, const MPI_Comm &mpi_communicator)
T sum(const T &t, const MPI_Comm &mpi_communicator)
const MPI_Datatype mpi_type_id_for_type
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
T max(const T &t, const MPI_Comm &mpi_communicator)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
static const unsigned int invalid_unsigned_int
hid_t hdf5_type_id(const number *)