* 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<NumberType> compute_SVD(ScaLAPACKMatrix<NumberType> &U,
ScaLAPACKMatrix<NumberType> &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<NumberType> &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()).
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);
}
psgesvd_(jobu,jobvt,m,n,A,ia,ja,desca,S,U,iu,ju,descu,VT,ivt,jvt,descvt,work,lwork,info);
}
+
+template <typename number>
+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
+template <typename NumberType>
+void ScaLAPACKMatrix<NumberType>::least_squares(ScaLAPACKMatrix<NumberType> &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 <typename NumberType>
NumberType ScaLAPACKMatrix<NumberType>::reciprocal_condition_number(const NumberType a_norm) const
{
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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<NumberType>&,const bool)
+
+#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/base/process_grid.h>
+
+#include <deal.II/lac/scalapack.h>
+
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include <memory>
+
+
+template <typename NumberType>
+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<Utilities::MPI::ProcessGrid> grid_2d = std::make_shared<Utilities::MPI::ProcessGrid>(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<NumberType> 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<NumberType> scalapack_A (4,3,grid_2d,block_size,block_size);
+ ScaLAPACKMatrix<NumberType> 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<NumberType> result(4,5);
+ scalapack_B.copy_to(result);
+
+ result.add(-1,full_X_I);
+ AssertThrow(result.frobenius_norm()<tol,ExcMessage("solution deviates from reference"));
+}
+
+
+
+int main (int argc,char **argv)
+{
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
+
+ const std::vector<unsigned int> blocks = {{1,2}};
+ const double tol_double = 1e-10;
+ const float tol_float = 1e-5;
+
+ for (const auto &b : blocks)
+ test<double>(b,tol_double);
+
+ for (const auto &b : blocks)
+ test<float>(b,tol_float);
+}
--- /dev/null
+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||
--- /dev/null
+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||
--- /dev/null
+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||
projects / dealii.git / commitdiff