void compute_cholesky_factorization ();
/**
- * Invert the matrix by first computing a Cholesky factorization and then
- * building the actual inverse using <code>pXpotri</code>. The inverse is stored
- * in this object.
+ * Compute the LU factorization of the matrix using ScaLAPACK
+ * function <code>pXgetrf</code> and partial pivoting with row interchanges.
+ * The result of the factorization is stored in this object.
+ */
+ void compute_lu_factorization ();
+
+ /**
+ * Invert the matrix by first computing a Cholesky for symmetric matrices
+ * or a LU factorization for general matrices and then
+ * building the actual inverse using <code>pXpotri</code> or <code>pXgetri</code>.
+ *
+ * If a Cholesky or LU factorization has been applied previously,
+ * <code>pXpotri</code> or <code>pXgetri</code> are called directly.
+ *
+ * The inverse is stored in this object.
*/
void invert();
*/
mutable std::vector<int> iwork;
+ /**
+ * Integer array holding pivoting information required
+ * by ScaLAPACK's matrix factorization routines.
+ */
+ std::vector<int> ipiv;
+
/**
* A character to define where elements are stored in case
* ScaLAPACK operations support this.
const int *DESCA,
int *INFO);
+ /**
+ * Computes an LU factorization of a general distributed matrix sub( A )
+ * using partial pivoting with row interchanges.
+ *
+ * http://www.netlib.org/scalapack/explore-html/df/dfe/pdgetrf_8f_source.html
+ * https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lgetrf.htm
+ */
+ void pdgetrf_(const int *m,
+ const int *n,
+ double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *ipiv,
+ int *INFO);
+ void psgetrf_(const int *m,
+ const int *n,
+ float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *ipiv,
+ int *INFO);
/**
* Compute the inverse of a real symmetric positive definite
const int *DESCA,
int *INFO);
+ /**
+ * PDGETRI computes the inverse of a distributed matrix using the LU
+ * factorization computed by PDGETRF. This method inverts U and then
+ * computes the inverse of sub( A ) = A(IA:IA+N-1,JA:JA+N-1) denoted
+ * InvA by solving the system InvA*L = inv(U) for InvA.
+ *
+ * http://www.netlib.org/scalapack/explore-html/d3/df3/pdgetri_8f_source.html
+ * https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lgetri.htm
+ */
+ void pdgetri_(const int *N,
+ double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ const int *ipiv,
+ double *work,
+ int *lwork,
+ int *iwork,
+ int *liwork,
+ int *info);
+ void psgetri_(const int *N,
+ float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ const int *ipiv,
+ float *work,
+ int *lwork,
+ int *iwork,
+ int *liwork,
+ int *info);
+
/**
* Estimate the reciprocal of the condition number (in the
* l1-norm) of a real symmetric positive definite distributed matrix
}
+template <typename number>
+inline void pgetrf(const int * /*m*/,
+ const int * /*n*/,
+ number * /*A*/,
+ const int * /*IA*/,
+ const int * /*JA*/,
+ const int * /*DESCA*/,
+ int * /*ipiv*/,
+ int * /*INFO*/)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void pgetrf(const int *m,
+ const int *n,
+ double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *ipiv,
+ int *INFO)
+{
+ pdgetrf_(m,n,A,IA,JA,DESCA,ipiv,INFO);
+}
+
+inline void pgetrf(const int *m,
+ const int *n,
+ float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *ipiv,
+ int *INFO)
+{
+ psgetrf_(m,n,A,IA,JA,DESCA,ipiv,INFO);
+}
+
+
template <typename number>
inline void ppotri(const char * /*UPLO*/,
const int * /*N*/,
}
+template <typename number>
+inline void pgetri(const int * /*N*/,
+ number * /*A*/,
+ const int * /*IA*/,
+ const int * /*JA*/,
+ const int * /*DESCA*/,
+ const int * /*ipiv*/,
+ number * /*work*/,
+ int * /*lwork*/,
+ int * /*iwork*/,
+ int * /*liwork*/,
+ int * /*info*/)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void pgetri(const int *N,
+ double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ const int *ipiv,
+ double *work,
+ int *lwork,
+ int *iwork,
+ int *liwork,
+ int *info)
+{
+ pdgetri_(N,A,IA,JA,DESCA,ipiv,work,lwork,iwork,liwork,info);
+}
+
+inline void pgetri(const int *N,
+ float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ const int *ipiv,
+ float *work,
+ int *lwork,
+ int *iwork,
+ int *liwork,
+ int *info)
+{
+ psgetri_(N,A,IA,JA,DESCA,ipiv,work,lwork,iwork,liwork,info);
+}
+
+
template <typename number>
inline void ppocon(const char * /*uplo*/,
const int * /*N*/,
template <typename NumberType>
void ScaLAPACKMatrix<NumberType>::compute_cholesky_factorization()
{
- Assert (n_columns == n_rows,
- ExcMessage("Cholesky factorization can be applied to SPD matrices only."));
+ Assert (n_columns==n_rows && property == LAPACKSupport::Property::symmetric,
+ ExcMessage("Cholesky factorization can be applied to symmetric matrices only."));
+ Assert (state == LAPACKSupport::matrix,
+ ExcMessage("Matrix has to be in Matrix state before calling this function."));
if (grid->mpi_process_is_active)
{
ppotrf(&uplo,&n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor,&info);
AssertThrow (info==0, LAPACKSupport::ExcErrorCode("ppotrf", info));
}
- property = (uplo=='L' ? LAPACKSupport::lower_triangular : LAPACKSupport::upper_triangular);
state = LAPACKSupport::cholesky;
+ property = (uplo=='L' ? LAPACKSupport::lower_triangular : LAPACKSupport::upper_triangular);
}
template <typename NumberType>
-void ScaLAPACKMatrix<NumberType>::invert()
+void ScaLAPACKMatrix<NumberType>::compute_lu_factorization()
{
- if (state == LAPACKSupport::matrix)
- compute_cholesky_factorization();
+ Assert (state == LAPACKSupport::matrix,
+ ExcMessage("Matrix has to be in Matrix state before calling this function."));
if (grid->mpi_process_is_active)
{
int info = 0;
NumberType *A_loc = &this->values[0];
- ppotri (&uplo,&n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor,&info);
- AssertThrow (info==0, LAPACKSupport::ExcErrorCode("ppotri", info));
+
+ const int iarow = indxg2p_(&submatrix_row, &row_block_size, &(grid->this_process_row), &first_process_row, &(grid->n_process_rows));
+ const int mp = numroc_(&n_rows, &row_block_size, &(grid->this_process_row), &iarow, &(grid->n_process_rows));
+ ipiv.resize(mp+row_block_size);
+
+ pgetrf(&n_rows,&n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor,ipiv.data(),&info);
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgetrf", info));
+ }
+ state = LAPACKSupport::State::lu;
+ property = LAPACKSupport::Property::general;
+}
+
+
+
+template <typename NumberType>
+void ScaLAPACKMatrix<NumberType>::invert()
+{
+ // Check whether matrix is symmetric and save flag.
+ // If a Cholesky factorization has been applied previously,
+ // the original matrix was symmetric.
+ const bool is_symmetric = (property == LAPACKSupport::symmetric || state == LAPACKSupport::State::cholesky);
+
+ // Matrix is neither in Cholesky nor LU state.
+ // Compute the required factorizations based on the property of the matrix.
+ if ( !(state==LAPACKSupport::State::lu || state==LAPACKSupport::State::cholesky) )
+ {
+ if (is_symmetric)
+ compute_cholesky_factorization();
+ else
+ compute_lu_factorization();
+ }
+ if (grid->mpi_process_is_active)
+ {
+ int info = 0;
+ NumberType *A_loc = &this->values[0];
+
+ if (is_symmetric)
+ {
+ ppotri(&uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, &info);
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("ppotri", info));
+ }
+ else
+ {
+ int lwork=-1, liwork=-1;
+ work.resize(1);
+ iwork.resize(1);
+
+ pgetri(&n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor,
+ ipiv.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
+
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgetri", info));
+ lwork=work[0];
+ liwork=iwork[0];
+ work.resize(lwork);
+ iwork.resize(liwork);
+
+ pgetri(&n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor,
+ ipiv.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
+
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgetri", info));
+ }
}
- state = LAPACKSupport::inverse_matrix;
+ state = LAPACKSupport::State::inverse_matrix;
}
--- /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 inverse via LU decomposition in ScaLAPACK vs FullMatrix for general matrices
+ *
+ * relevant tests:
+ * full_matrix/full_matrix_invert_01
+ * full_matrix/full_matrix_56
+ * full_matrix/full_matrix_53
+ * full_matrix/full_matrix_42
+ * full_matrix/full_matrix_41
+ */
+
+#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)
+{
+ 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_in(size), inverse(size), full_out(size), diff(size), prod1(size), prod2(size), one(size);
+ create_spd (full_in);
+
+ std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size,block_size);
+ ScaLAPACKMatrix<NumberType> scalapack_matrix (size, grid, block_size,
+ LAPACKSupport::Property::general);
+
+ pcout << size << " " << block_size << " " << grid->get_process_grid_rows() << " " << grid->get_process_grid_columns() << std::endl;
+
+ // add a skew-symmetric matrix to the s.p.d matrix
+ // the positive definiteness is inherited from the symmetric part
+ for (unsigned int i=1; i<size; ++i)
+ for (unsigned int j=i+1; j<size; ++j)
+ {
+ full_in(i,j) *= 1.3;
+ full_in(j,i) *= -1.3;
+ }
+ one = 0.;
+ for (unsigned int i=0; i<size; ++i)
+ one(i,i) = 1.;
+ // invert via Lapack
+ inverse.invert(full_in);
+ inverse.mmult(prod1, full_in);
+ prod1.add(-1., one);
+ const NumberType lapack_error = prod1.linfty_norm();
+
+ // estimated condition number from 1-norm:
+ const NumberType k = full_in.l1_norm() * inverse.l1_norm();
+ const NumberType tol = k * 1000 * std::numeric_limits<NumberType>::epsilon();
+
+ // invert via ScaLAPACK
+ scalapack_matrix = full_in;
+ scalapack_matrix.invert();
+ scalapack_matrix.copy_to(full_out);
+ full_out.mmult(prod2, full_in);
+ prod2.add(-1.,one);
+ const NumberType error = prod2.linfty_norm();
+
+ if ( error > tol && this_mpi_process == 0)
+ {
+ diff = 0;
+ diff.add( 1., inverse);
+ diff.add(-1., full_out);
+
+ std::cout << "Norm of the error " << error
+ << " is more than the threshold " << tol
+ << " . Norm of the A^{-1}*A using Lapack is " << lapack_error << std::endl
+ << "===== Expected to have:" << std::endl;
+ inverse.print_formatted (std::cout);
+ std::cout << "===== But got:" << std::endl;
+ full_out.print_formatted(std::cout);
+ std::cout << "===== Difference:" << std::endl;
+ diff.print_formatted(std::cout);
+ AssertThrow (false, dealii::ExcInternalError());
+ }
+ else
+ pcout << "Ok" << std::endl;
+}
+
+
+
+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}};
+
+ for (const auto &s : sizes)
+ for (const auto &b : blocks)
+ if (b <= s)
+ test<float>(s,b);
+
+ for (const auto &s : sizes)
+ for (const auto &b : blocks)
+ if (b <= s)
+ test<double>(s,b);
+}
--- /dev/null
+32 32 1 1
+Ok
+64 32 1 1
+Ok
+64 64 1 1
+Ok
+120 32 1 1
+Ok
+120 64 1 1
+Ok
+320 32 1 1
+Ok
+320 64 1 1
+Ok
+640 32 1 1
+Ok
+640 64 1 1
+Ok
+32 32 1 1
+Ok
+64 32 1 1
+Ok
+64 64 1 1
+Ok
+120 32 1 1
+Ok
+120 64 1 1
+Ok
+320 32 1 1
+Ok
+320 64 1 1
+Ok
+640 32 1 1
+Ok
+640 64 1 1
+Ok
--- /dev/null
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 3 3
+Ok
+120 64 2 2
+Ok
+320 32 3 3
+Ok
+320 64 3 3
+Ok
+640 32 3 3
+Ok
+640 64 3 3
+Ok
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 3 3
+Ok
+120 64 2 2
+Ok
+320 32 3 3
+Ok
+320 64 3 3
+Ok
+640 32 3 3
+Ok
+640 64 3 3
+Ok
--- /dev/null
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 2 2
+Ok
+120 64 2 2
+Ok
+320 32 2 2
+Ok
+320 64 2 2
+Ok
+640 32 2 2
+Ok
+640 64 2 2
+Ok
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 2 2
+Ok
+120 64 2 2
+Ok
+320 32 2 2
+Ok
+320 64 2 2
+Ok
+640 32 2 2
+Ok
+640 64 2 2
+Ok
projects / dealii.git / commitdiff