From 1c1db70354c3ba6dd08bdd2304a62f95ec61e964 Mon Sep 17 00:00:00 2001 From: Benjamin Brands Date: Mon, 15 Jan 2018 14:29:07 +0100 Subject: [PATCH] added least squares routine and tests --- include/deal.II/lac/scalapack.h | 40 +++++- include/deal.II/lac/scalapack.templates.h | 95 +++++++++++++ source/lac/scalapack.cc | 48 +++++++ tests/scalapack/scalapack_09.cc | 137 +++++++++++++++++++ tests/scalapack/scalapack_09.mpirun=1.output | 4 + tests/scalapack/scalapack_09.mpirun=2.output | 4 + tests/scalapack/scalapack_09.mpirun=4.output | 4 + 7 files changed, 331 insertions(+), 1 deletion(-) create mode 100644 tests/scalapack/scalapack_09.cc create mode 100644 tests/scalapack/scalapack_09.mpirun=1.output create mode 100644 tests/scalapack/scalapack_09.mpirun=2.output create mode 100644 tests/scalapack/scalapack_09.mpirun=4.output diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h index c58936b101..405b1248fd 100644 --- a/include/deal.II/lac/scalapack.h +++ b/include/deal.II/lac/scalapack.h @@ -215,12 +215,50 @@ public: * have to be constructed with the same process grid and block cyclic distribution. * If right singular vectors are required matrices A and VT * have to be constructed with the same process grid and block cyclic distribution. - */ + */ std::vector compute_SVD(ScaLAPACKMatrix &U, ScaLAPACKMatrix &VT, const bool left_singluar_vectors=false, const bool right_singluar_vectors=false); + + /** + * Function solves overdetermined or underdetermined real linear + * systems involving an M-by-N matrix A, or its transpose, using a QR or LQ factorization of A. + * + * It is assumed that A has full rank: \f$rank(A) = \min(M,N)\f$. + * Upon exit the columns of B contain the solutions and + * the following options are supported: + * - 1. If transpose==false and \f$M \geq N\f$: least squares solution of overdetermined system + * \f$\min \Vert B - A X\Vert\f$. + * + * Upon exit the rows 0 to N-1 contain the least square solution vectors. The residual sum of squares + * for each column is given by the sum of squares of elements N to M-1 in that column + * + * - 2. If transpose==false and \f$M < N\f$: find minimum norm solutions of underdetermined systems + * \f$A X = B\f$. + * + * Upon exit the columns of B contain the minimum norm solution vectors + * + * - 3. If transpose==true and and \f$M \geq N\f$: find minimum norm solutions of underdetermined system + * \f$ A^\top X = B\f$ + * + * Upon exit the columns of B contain the minimum norm solution vectors + * + * - 4. If transpose==true and \f$M < N\f$: least squares solution of overdetermined system + * \f$\min \Vert B - A^\top X\Vert\f$. + * + * Upon exit the rows 0 to M-1 contain the least square solution vectors. The residual sum of squares + * for each column is given by the sum of squares of elements M to N-1 in that column + * . + * If transpose==false B is M x NRHS matrix, otherwise it is NxNRHS. + * The matrices A and B must have an identical block cyclic distribution for rows and columns + */ + void least_squares(ScaLAPACKMatrix &B, + const bool transpose=false); + + + /** * Estimate the the condition number of a SPD matrix in the $l_1$-norm. * The matrix has to be in the Cholesky state (see compute_cholesky_factorization()). diff --git a/include/deal.II/lac/scalapack.templates.h b/include/deal.II/lac/scalapack.templates.h index 9ba8796d8c..173ff298a5 100644 --- a/include/deal.II/lac/scalapack.templates.h +++ b/include/deal.II/lac/scalapack.templates.h @@ -647,6 +647,42 @@ extern "C" float *work, int *lwork, int *info); + + /* + * P_GELS solves overdetermined or underdetermined real linear + * systems involving an M-by-N matrix A, or its transpose, + * using a QR or LQ factorization of A. It is assumed that A has full rank. + */ + void pdgels_(const char *trans, + const int *m, + const int *n, + const int *nrhs, + double *A, + const int *ia, + const int *ja, + const int *desca, + double *B, + const int *ib, + const int *jb, + const int *descb, + double *work, + int *lwork, + int *info); + void psgels_(const char *trans, + const int *m, + const int *n, + const int *nrhs, + float *A, + const int *ia, + const int *ja, + const int *desca, + float *B, + const int *ib, + const int *jb, + const int *descb, + float *work, + int *lwork, + int *info); } @@ -1394,6 +1430,65 @@ inline void pgesvd(const char *jobu, psgesvd_(jobu,jobvt,m,n,A,ia,ja,desca,S,U,iu,ju,descu,VT,ivt,jvt,descvt,work,lwork,info); } + +template +inline void pgels(const char *trans, + const int *m, + const int *n, + const int *nrhs, + number *A, + const int *ia, + const int *ja, + const int *desca, + number *B, + const int *ib, + const int *jb, + const int *descb, + number *work, + int *lwork, + int *info) +{ + Assert (false, dealii::ExcNotImplemented()); +} + +inline void pgels(const char *trans, + const int *m, + const int *n, + const int *nrhs, + double *A, + const int *ia, + const int *ja, + const int *desca, + double *B, + const int *ib, + const int *jb, + const int *descb, + double *work, + int *lwork, + int *info) +{ + pdgels_(trans,m,n,nrhs,A,ia,ja,desca,B,ib,jb,descb,work,lwork,info); +} + +inline void pgels(const char *trans, + const int *m, + const int *n, + const int *nrhs, + float *A, + const int *ia, + const int *ja, + const int *desca, + float *B, + const int *ib, + const int *jb, + const int *descb, + float *work, + int *lwork, + int *info) +{ + psgels_(trans,m,n,nrhs,A,ia,ja,desca,B,ib,jb,descb,work,lwork,info); +} + #endif // DEAL_II_WITH_SCALAPACK #endif // dealii_scalapack_templates_h diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc index ea9d7dfa4c..832199c2bb 100644 --- a/source/lac/scalapack.cc +++ b/source/lac/scalapack.cc @@ -571,6 +571,54 @@ std::vector ScaLAPACKMatrix::compute_SVD(ScaLAPACKMatrix +template +void ScaLAPACKMatrix::least_squares(ScaLAPACKMatrix &B, + const bool transpose) +{ + Assert(grid==B.grid,ExcMessage("The matrices A and B need to have the same process grid")); + Assert (state == LAPACKSupport::matrix, + ExcMessage("Matrix has to be in Matrix state before calling this function.")); + Assert (B.state == LAPACKSupport::matrix, + ExcMessage("Matrix B has to be in Matrix state before calling this function.")); + + transpose ? + (Assert(n_columns==B.n_rows,ExcDimensionMismatch(n_columns,B.n_rows))) : + (Assert(n_rows==B.n_rows,ExcDimensionMismatch(n_rows,B.n_rows))); + + //see https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lgels.htm + Assert(row_block_size==column_block_size,ExcMessage("Use identical block sizes for rows and columns of matrix A")); + Assert(B.row_block_size==B.column_block_size,ExcMessage("Use identical block sizes for rows and columns of matrix B")); + Assert(row_block_size==B.row_block_size,ExcMessage("Use identical block-cyclic distribution for matrices A and B")); + + Threads::Mutex::ScopedLock lock (mutex); + + if (grid->mpi_process_is_active) + { + char trans = transpose ? 'T' : 'N'; + NumberType *A_loc = & this->values[0]; + NumberType *B_loc = & B.values[0]; + int info = 0; + /* + * by setting lwork to -1 a workspace query for optimal length of work is performed + */ + int lwork=-1; + work.resize(1); + + pgels(&trans,&n_rows,&n_columns,&B.n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor, + B_loc,&B.submatrix_row,&B.submatrix_column,B.descriptor,&work[0],&lwork,&info); + AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgels", info)); + + lwork=work[0]; + work.resize(lwork); + + pgels(&trans,&n_rows,&n_columns,&B.n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor, + B_loc,&B.submatrix_row,&B.submatrix_column,B.descriptor,&work[0],&lwork,&info); + AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgels", info)); + } +} + + + template NumberType ScaLAPACKMatrix::reciprocal_condition_number(const NumberType a_norm) const { diff --git a/tests/scalapack/scalapack_09.cc b/tests/scalapack/scalapack_09.cc new file mode 100644 index 0000000000..c7b3ac7b11 --- /dev/null +++ b/tests/scalapack/scalapack_09.cc @@ -0,0 +1,137 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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" + +// test least_squares(ScaLAPACKMatrix&,const bool) + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + + +template +void test(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)); + + std::shared_ptr grid_2d = std::make_shared(mpi_communicator,4,3,block_size,block_size); + + //examples from https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lgels.htm + FullMatrix full_A_I(4,3), full_B_I(4,5), full_X_I(4,5); + + //FIXME: Add more tests for different cases!!! + pcout << "Solving least squares problem ||B - A*X||" << std::endl; + + full_A_I(0,0)=1.; + full_A_I(0,1)=-2.; + full_A_I(0,2)=-1.; + full_A_I(1,0)=2.; + full_A_I(1,1)=0.; + full_A_I(1,2)=1.; + full_A_I(2,0)=2.; + full_A_I(2,1)=-4.; + full_A_I(2,2)=2.; + full_A_I(3,0)=4.; + full_A_I(3,1)=0.; + full_A_I(3,2)=0.; + + full_B_I(0,0)=-1.; + full_B_I(0,1)=-2.; + full_B_I(0,2)=-7.; + full_B_I(0,3)=0.; + full_B_I(0,4)=-5.; + full_B_I(1,0)=1.; + full_B_I(1,1)=3.; + full_B_I(1,2)=4.; + full_B_I(1,3)=3.; + full_B_I(1,4)=5.; + full_B_I(2,0)=1.; + full_B_I(2,1)=0.; + full_B_I(2,2)=4.; + full_B_I(2,3)=2.; + full_B_I(2,4)=2.; + full_B_I(3,0)=-2.; + full_B_I(3,1)=4.; + full_B_I(3,2)=4.; + full_B_I(3,3)=0.; + full_B_I(3,4)=4.; + + full_X_I(0,0)=-0.4; + full_X_I(0,1)=1.; + full_X_I(0,2)=0.8; + full_X_I(0,3)=0.2; + full_X_I(0,4)=1.; + full_X_I(1,0)=0.; + full_X_I(1,1)=1.; + full_X_I(1,2)=1.5; + full_X_I(1,3)=0.; + full_X_I(1,4)=1.5; + full_X_I(2,0)=1.; + full_X_I(2,1)=1.; + full_X_I(2,2)=4.; + full_X_I(2,3)=1.; + full_X_I(2,4)=3.; + full_X_I(3,0)=-1.; + full_X_I(3,1)=0.; + full_X_I(3,2)=2.; + full_X_I(3,3)=-2.; + full_X_I(3,4)=0.; + + //compute eigenpairs of s.p.d matrix + ScaLAPACKMatrix scalapack_A (4,3,grid_2d,block_size,block_size); + ScaLAPACKMatrix scalapack_B (4,5,grid_2d,block_size,block_size); + scalapack_A.set_property(LAPACKSupport::Property::general); + scalapack_A = full_A_I; + scalapack_B = full_B_I; + scalapack_A.least_squares(scalapack_B,false); + FullMatrix result(4,5); + scalapack_B.copy_to(result); + + result.add(-1,full_X_I); + AssertThrow(result.frobenius_norm() blocks = {{1,2}}; + const double tol_double = 1e-10; + const float tol_float = 1e-5; + + for (const auto &b : blocks) + test(b,tol_double); + + for (const auto &b : blocks) + test(b,tol_float); +} diff --git a/tests/scalapack/scalapack_09.mpirun=1.output b/tests/scalapack/scalapack_09.mpirun=1.output new file mode 100644 index 0000000000..28791e6d80 --- /dev/null +++ b/tests/scalapack/scalapack_09.mpirun=1.output @@ -0,0 +1,4 @@ +Solving least squares problem ||B - A*X|| +Solving least squares problem ||B - A*X|| +Solving least squares problem ||B - A*X|| +Solving least squares problem ||B - A*X|| diff --git a/tests/scalapack/scalapack_09.mpirun=2.output b/tests/scalapack/scalapack_09.mpirun=2.output new file mode 100644 index 0000000000..28791e6d80 --- /dev/null +++ b/tests/scalapack/scalapack_09.mpirun=2.output @@ -0,0 +1,4 @@ +Solving least squares problem ||B - A*X|| +Solving least squares problem ||B - A*X|| +Solving least squares problem ||B - A*X|| +Solving least squares problem ||B - A*X|| diff --git a/tests/scalapack/scalapack_09.mpirun=4.output b/tests/scalapack/scalapack_09.mpirun=4.output new file mode 100644 index 0000000000..28791e6d80 --- /dev/null +++ b/tests/scalapack/scalapack_09.mpirun=4.output @@ -0,0 +1,4 @@ +Solving least squares problem ||B - A*X|| +Solving least squares problem ||B - A*X|| +Solving least squares problem ||B - A*X|| +Solving least squares problem ||B - A*X|| -- 2.39.5