private:
+ /**
+ * Calculate the norm of a distributed symmetric dense matrix using ScaLAPACK's
+ * internal function.
+ */
+ NumberType norm_symmetric(const char type) const;
+
/**
* Calculate the norm of a distributed dense matrix using ScaLAPACK's
* internal function.
*/
- NumberType norm(const char type) const;
+ NumberType norm_general(const char type) const;
/**
* Computing selected eigenvalues and, optionally, the eigenvectors.
* or the element of largest absolute value of a distributed matrix
*/
double pdlange_(char const *norm,
- int const &m,
- int const &n,
- double *A,
- int const &ia,
- int const &ja,
- int *desca,
+ const int *m,
+ const int *n,
+ const double *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
double *work);
- float pslange_(char const *norm,
- int const &m,
- int const &n,
- float *A,
- int const &ia,
- int const &ja,
- int *desca,
+ float pslange_(const char *norm,
+ const int *m,
+ const int *n,
+ const float *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
float *work);
/**
template <typename number>
-inline number plange(char const * /*norm*/,
- int const &/*m*/,
- int const &/*n*/,
- number * /*A*/,
- int const &/*ia*/,
- int const &/*ja*/,
- int * /*desca*/,
+inline number plange(const char * /*norm*/,
+ const int * /*m*/,
+ const int * /*n*/,
+ const number * /*A*/,
+ const int * /*ia*/,
+ const int * /*ja*/,
+ const int * /*desca*/,
number * /*work*/)
{
Assert (false, dealii::ExcNotImplemented());
}
-inline double plange(char const *norm,
- int const &m,
- int const &n,
- double *A,
- int const &ia,
- int const &ja,
- int *desca,
+inline double plange(const char *norm,
+ const int *m,
+ const int *n,
+ const double *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
double *work)
{
return pdlange_(norm, m, n, A, ia, ja,desca, work);
}
-inline float plange(char const *norm,
- int const &m,
- int const &n,
- float *A,
- int const &ia,
- int const &ja,
- int *desca,
+inline float plange(const char *norm,
+ const int *m,
+ const int *n,
+ const float *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
float *work)
{
return pslange_(norm, m, n, A, ia, ja,desca, work);
NumberType ScaLAPACKMatrix<NumberType>::l1_norm() const
{
const char type('O');
- return norm(type);
+
+ if (property == LAPACKSupport::symmetric)
+ return norm_symmetric(type);
+ else
+ return norm_general(type);
}
NumberType ScaLAPACKMatrix<NumberType>::linfty_norm() const
{
const char type('I');
- return norm(type);
+
+ if (property == LAPACKSupport::symmetric)
+ return norm_symmetric(type);
+ else
+ return norm_general(type);
}
NumberType ScaLAPACKMatrix<NumberType>::frobenius_norm() const
{
const char type('F');
- return norm(type);
+
+ if (property == LAPACKSupport::symmetric)
+ return norm_symmetric(type);
+ else
+ return norm_general(type);
+}
+
+
+
+template <typename NumberType>
+NumberType ScaLAPACKMatrix<NumberType>::norm_general(const char type) const
+{
+ Assert (state == LAPACKSupport::matrix ||
+ state == LAPACKSupport::inverse_matrix,
+ ExcMessage("norms can be called in matrix state only."));
+ Threads::Mutex::ScopedLock lock (mutex);
+ NumberType res = 0.;
+
+ if (grid->mpi_process_is_active)
+ {
+ const int iarow = indxg2p_(&submatrix_row, &row_block_size, &(grid->this_process_row), &first_process_row, &(grid->n_process_rows));
+ const int iacol = indxg2p_(&submatrix_column, &column_block_size, &(grid->this_process_column), &first_process_column, &(grid->n_process_columns));
+ const int mp0 = numroc_(&n_rows, &row_block_size, &(grid->this_process_row), &iarow, &(grid->n_process_rows));
+ const int nq0 = numroc_(&n_columns, &column_block_size, &(grid->this_process_column), &iacol, &(grid->n_process_columns));
+
+ // type='M': compute largest absolute value
+ // type='F' || type='E': compute Frobenius norm
+ // type='0' || type='1': compute infinity norm
+ int lwork=0; // for type == 'M' || type == 'F' || type == 'E'
+ if (type=='O' || type=='1')
+ lwork = nq0;
+ else if (type=='I')
+ lwork = mp0;
+
+ work.resize(lwork);
+ const NumberType *A_loc = this->values.begin();
+ res = plange(&type, &n_rows, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, work.data());
+ }
+ grid->send_to_inactive(&res);
+ return res;
}
template <typename NumberType>
-NumberType ScaLAPACKMatrix<NumberType>::norm(const char type) const
+NumberType ScaLAPACKMatrix<NumberType>::norm_symmetric(const char type) const
{
Assert (state == LAPACKSupport::matrix ||
state == LAPACKSupport::inverse_matrix,
ExcMessage("norms can be called in matrix state only."));
+ Assert (property == LAPACKSupport::symmetric,
+ ExcMessage("Matrix has to be symmetric for this operation."));
Threads::Mutex::ScopedLock lock (mutex);
NumberType res = 0.;
0 :
2*Nq0+Np0+ldw;
work.resize(lwork);
- const NumberType *A_loc = &this->values[0];
- res = plansy(&type, &uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, &work[0]);
+ const NumberType *A_loc = this->values.begin();
+ res = plansy(&type, &uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, work.data());
}
grid->send_to_inactive(&res);
return res;
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, 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.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include "../tests.h"
+#include "../lapack/create_matrix.h"
+
+// test norms of ScaLAPACK vs FullMatrix for general matrices
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/multithread_info.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/lac/scalapack.h>
+
+#include <fstream>
+#include <iostream>
+
+template <typename NumberType>
+void test(const unsigned int size, const unsigned int block_size, const NumberType tol)
+{
+ MPI_Comm mpi_communicator(MPI_COMM_WORLD);
+ const unsigned int n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator));
+ const unsigned int this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator));
+
+ ConditionalOStream pcout (std::cout, (this_mpi_process ==0));
+
+ // Create SPD matrices of requested size:
+ FullMatrix<NumberType> full(size);
+
+ std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size,block_size);
+ ScaLAPACKMatrix<NumberType> scalapack (size, size, grid, block_size, block_size,
+ LAPACKSupport::Property::general);
+
+ pcout << size << " " << block_size << " " << grid->get_process_grid_rows() << " " << grid->get_process_grid_columns() << std::endl;
+
+ create_random (full);
+ scalapack = full;
+
+ const NumberType l1 = full.l1_norm();
+ const NumberType linfty = full.linfty_norm();
+ const NumberType frobenius = full.frobenius_norm();
+
+ // local result on this core:
+ const NumberType s_l1 = scalapack.l1_norm();
+ const NumberType s_linfty = scalapack.linfty_norm();
+ const NumberType s_frobenius = scalapack.frobenius_norm();
+
+ // make sure we have the same result on all cores, do average:
+ const NumberType as_l1 = dealii::Utilities::MPI::sum(s_l1, mpi_communicator)/n_mpi_processes;
+ const NumberType as_linfty = dealii::Utilities::MPI::sum(s_linfty, mpi_communicator)/n_mpi_processes;
+ const NumberType as_frobenius = dealii::Utilities::MPI::sum(s_frobenius, mpi_communicator)/n_mpi_processes;
+
+ pcout << l1 << " " << s_l1 << " " << as_l1 << std::endl
+ << linfty << " " << s_linfty << " " << as_linfty << std::endl
+ << frobenius << " " << s_frobenius << " " << as_frobenius << std::endl;
+
+ AssertThrow (std::abs(l1 -as_l1) < tol*std::abs(l1), dealii::ExcInternalError());
+ AssertThrow (std::abs(linfty -as_linfty) < tol*std::abs(linfty), dealii::ExcInternalError());
+ AssertThrow (std::abs(frobenius -as_frobenius) < tol*std::abs(frobenius), dealii::ExcInternalError());
+}
+
+
+
+int main (int argc,char **argv)
+{
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
+
+ const std::vector<unsigned int> sizes = {{32,64,120,320,640}};
+ const std::vector<unsigned int> blocks = {{32,64}};
+
+ const double tol_double = 1e-10;
+ const float tol_float = 1e-5;
+
+ /*for (const auto &s : sizes)
+ for (const auto &b : blocks)
+ if (b <= s)
+ test<float>(s,b,tol_float);*/
+
+ for (const auto &s : sizes)
+ for (const auto &b : blocks)
+ if (b <= s)
+ test<double>(s,b,tol_double);
+}
--- /dev/null
+32 32 1 1
+19.2741 19.2741 19.2741
+20.6796 20.6796 20.6796
+18.5785 18.5785 18.5785
+64 32 1 1
+38.5552 38.5552 38.5552
+38.2739 38.2739 38.2739
+36.8389 36.8389 36.8389
+64 64 1 1
+37.8671 37.8671 37.8671
+37.8657 37.8657 37.8657
+36.7171 36.7171 36.7171
+120 32 1 1
+69.0663 69.0663 69.0663
+68.0208 68.0208 68.0208
+69.1394 69.1394 69.1394
+120 64 1 1
+68.2078 68.2078 68.2078
+66.8542 66.8542 66.8542
+69.0112 69.0112 69.0112
+320 32 1 1
+174.387 174.387 174.387
+172.758 172.758 172.758
+184.876 184.876 184.876
+320 64 1 1
+182.948 182.948 182.948
+174.945 174.945 174.945
+184.521 184.521 184.521
+640 32 1 1
+338.935 338.935 338.935
+342.898 342.898 342.898
+369.65 369.65 369.65
+640 64 1 1
+342.421 342.421 342.421
+341.617 341.617 341.617
+369.497 369.497 369.497
--- /dev/null
+32 32 1 1
+19.2741 19.2741 19.2741
+20.6796 20.6796 20.6796
+18.5785 18.5785 18.5785
+64 32 2 2
+38.5552 38.5552 38.5552
+38.2739 38.2739 38.2739
+36.8389 36.8389 36.8389
+64 64 1 1
+37.8671 37.8671 37.8671
+37.8657 37.8657 37.8657
+36.7171 36.7171 36.7171
+120 32 3 3
+69.0663 69.0663 69.0663
+68.0208 68.0208 68.0208
+69.1394 69.1394 69.1394
+120 64 2 2
+68.2078 68.2078 68.2078
+66.8542 66.8542 66.8542
+69.0112 69.0112 69.0112
+320 32 3 3
+174.387 174.387 174.387
+172.758 172.758 172.758
+184.876 184.876 184.876
+320 64 3 3
+182.948 182.948 182.948
+174.945 174.945 174.945
+184.521 184.521 184.521
+640 32 3 3
+338.935 338.935 338.935
+342.898 342.898 342.898
+369.65 369.65 369.65
+640 64 3 3
+342.421 342.421 342.421
+341.617 341.617 341.617
+369.497 369.497 369.497
--- /dev/null
+32 32 1 1
+19.2741 19.2741 19.2741
+20.6796 20.6796 20.6796
+18.5785 18.5785 18.5785
+64 32 2 2
+38.5552 38.5552 38.5552
+38.2739 38.2739 38.2739
+36.8389 36.8389 36.8389
+64 64 1 1
+37.8671 37.8671 37.8671
+37.8657 37.8657 37.8657
+36.7171 36.7171 36.7171
+120 32 2 2
+69.0663 69.0663 69.0663
+68.0208 68.0208 68.0208
+69.1394 69.1394 69.1394
+120 64 2 2
+68.2078 68.2078 68.2078
+66.8542 66.8542 66.8542
+69.0112 69.0112 69.0112
+320 32 2 2
+174.387 174.387 174.387
+172.758 172.758 172.758
+184.876 184.876 184.876
+320 64 2 2
+182.948 182.948 182.948
+174.945 174.945 174.945
+184.521 184.521 184.521
+640 32 2 2
+338.935 338.935 338.935
+342.898 342.898 342.898
+369.65 369.65 369.65
+640 64 2 2
+342.421 342.421 342.421
+341.617 341.617 341.617
+369.497 369.497 369.497
projects / dealii.git / commitdiff