diff --git a/CMakeLists.txt b/CMakeLists.txt index 8a6d8951a8..22df7a284b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -265,29 +265,17 @@ if(PROJECT_SOURCE_DIR STREQUAL PROJECT_BINARY_DIR) message(FATAL_ERROR "In-source build attempted; please clean the CMake cache and then switch to an out-of-source build, e.g.,\nrm CMakeCache.txt && rm -Rf CMakeFiles/\nmkdir build/ && cd build/ && cmake ..") endif() -# Get the Git revision -include(GetGitRevisionDescription) -get_git_head_revision(GIT_REFSPEC GIT_SHA1) - -# Ensure that the build type is set to either Release or Debug -if(CMAKE_BUILD_TYPE STREQUAL "Release") - # This option is okay as-is - set(CMAKE_BUILD_TYPE Release) -elseif(CMAKE_BUILD_TYPE STREQUAL "Debug") - # This option is okay as-is - set(CMAKE_BUILD_TYPE Debug) -elseif(CMAKE_BUILD_TYPE STREQUAL "RelWithDebInfo") - message(WARNING "RelWithDebInfo not supported; switching to Release") - set(CMAKE_BUILD_TYPE Release) -elseif(CMAKE_BUILD_TYPE STREQUAL "MinSizeRel") - message(WARNING "MinSizeRel not supported; switching to Release") - set(CMAKE_BUILD_TYPE Release) -else() - message(WARNING "Build mode not specified, defaulting to Release build.") - set(CMAKE_BUILD_TYPE Release) +if (NOT CMAKE_BUILD_TYPE) + set( CMAKE_BUILD_TYPE Release) +endif() + +set(ACCEPTED_BUILD_TYPES Release Debug RelWithDebInfo MinSizeRel) +list(FIND ACCEPTED_BUILD_TYPES ${CMAKE_BUILD_TYPE} IS_BUILD_TYPE_ACCEPTED) +if(${IS_BUILD_TYPE_ACCEPTED} EQUAL -1) + message(FATAL_ERROR "CMAKE_BUILD_TYPE of ${CMAKE_BUILD_TYPE} is not accepted.") endif() -if(CMAKE_BUILD_TYPE STREQUAL "Release") +if(CMAKE_BUILD_TYPE STREQUAL "Release" OR CMAKE_BUILD_TYPE STREQUAL "MinSizeRel") set(EL_RELEASE TRUE) else() set(EL_RELEASE FALSE) diff --git a/include/El/lapack_like/spectral.hpp b/include/El/lapack_like/spectral.hpp index 89f3c2afbb..1d1951b9a4 100644 --- a/include/El/lapack_like/spectral.hpp +++ b/include/El/lapack_like/spectral.hpp @@ -2012,6 +2012,7 @@ DistMatrix HessenbergSpectralCloud #include #include #include +//#include #include #include diff --git a/include/El/lapack_like/spectral/TSVD.hpp b/include/El/lapack_like/spectral/TSVD.hpp new file mode 100644 index 0000000000..15bd5fa9cc --- /dev/null +++ b/include/El/lapack_like/spectral/TSVD.hpp @@ -0,0 +1,368 @@ +/*** +* +* TSVD Based off of implementation by Andreas Noack. See: http://github.com/andreasnoack/TSVD.jl +*/ +#ifndef EL_LAPACKLIKE_SPECTRAL_TSVD_HPP +#define EL_LAPACKLIKE_SPECTRAL_TSVD_HPP +#include +#include + +namespace El { + +template +Real CopySign(const Real& x, const Real& y) { return y >= Real(0) ? Abs(x) : -Abs(x); } + +template +struct BidiagonalMatrix { + typedef Matrix > RealMatrix; + + BidiagonalMatrix(El::Int N) : mainDiag(N, 1), superDiag(N - 1, 1) {} + + RealMatrix mainDiag; + RealMatrix superDiag; +}; //BidiagonalMatrix + + +struct BidiagInfo { + Int numVecsReorth; +}; //struct BidiagInfo + +template +struct BidiagCtrl { + bool reorthIn = false; + Int minBlockSize=50; + Base reorthogTol = Pow(El::limits::Epsilon>(), El::Base(.7)); + Base convValTol = Pow(El::limits::Epsilon>(), El::Base(.7)); + Base convVecTol = Pow(El::limits::Epsilon>(), El::Base(.5)); + bool partialProject = false; +}; //end struct BidiagCtrl + +template +struct ReorthogonalizationData { + ReorthogonalizationData(El::Int N, F mu, F nu, F tau_): tau(tau_) { + muList.reserve(N); + nuList.reserve(N); + muList.push_back(mu); + muList.push_back(1); + nuList.push_back(nu); + nuList.push_back(1); + } + std::vector> muList, nuList; + F tau; +}; + +//TODO: Make batch with Gemv +template +Int Reorthogonalize(DistMatrix & Q, + Int j, + std::vector>& termList, + const BidiagCtrl& ctrl, + DistMatrix & x, + Matrix >& diagList) { + typedef Base Real; + const Real eps = limits::Epsilon(); + Int numVecsReorth = 0; + termList.emplace_back(Real(1)); + for (Int i = 0; i < j; ++i) { + auto qi = Q(ALL, IR(i)); + Axpy(-Dot(qi, x), qi, x); + termList[i] = eps; + numVecsReorth++; + } + Real alpha = Nrm2(x); + diagList.Set(j, 0, alpha); + Scale(Real(1) / alpha, x); + return numVecsReorth; +} + +template +bool omegaRecurrenceV(const Int j, + const El::Base& tau, + const BidiagonalMatrix & B, + const std::vector>& muList, + const El::Base& alpha, + const El::Base& beta, + const BidiagCtrl& ctrl, + std::vector>& nuList) { + auto& superDiag = B.superDiag; + auto& mainDiag = B.mainDiag; + bool foundInaccurate = false; //run omega recurrence + for (Int i = 0; i < j; ++i) { + auto nu = superDiag.Get(i, 0) * muList[i + 1] + + mainDiag.Get(i, 0) * muList[i] - beta * nuList[i]; + nu = ( nu + CopySign(tau, nu) ) / alpha; + foundInaccurate |= ( Abs(nu) > ctrl.reorthogTol ); + nuList[i] = nu; + } + return foundInaccurate; +} + + +template +bool omegaRecurrenceU(const Int j, + const El::Base& tau, + const BidiagonalMatrix & B, + const std::vector>& nuList, + const El::Base& alpha, + const El::Base& beta, + const BidiagCtrl& ctrl, + std::vector>& muList) { + auto& superDiag = B.superDiag; + auto& mainDiag = B.mainDiag; + bool foundInaccurate = false; //run omega recurrence + for (Int i = 0; i <= j; ++i) { + auto mu = mainDiag.Get(i, 0) * nuList[i] - alpha * muList[i]; + if ( i > 0 ) { + mu += superDiag.Get(i - 1, 0) * nuList[i - 1]; + } + mu = ( mu + CopySign(tau, mu) ) / beta; + foundInaccurate |= ( Abs(mu) > ctrl.reorthogTol ); + muList[i] = mu; + }; + return foundInaccurate; +} + + +template +BidiagInfo +GolubKahan( + El::Int iter, + El::Int k, + const ForwardOperator& A, const AdjointOperator& AAdj, + BidiagonalMatrix & B, El::DistMatrix& U, El::DistMatrix& V, + ReorthogonalizationData& reorthData, const El::BidiagCtrl& ctrl) { + + typedef El::Base Real; + typedef El::DistMatrix DM; + const Real eps = El::limits::Epsilon(); + + bool reorthB = ctrl.reorthIn; + auto& nuList = reorthData.nuList; + auto& muList = reorthData.muList; + El::Int numVecsReorth = 0; + //std::vector> maxMuList, maxNuList; + //maxMuList.reserve(steps); + //maxNuList.reserve(steps); + + DM v_j = V(ALL, IR(iter)); + DM u_j = U(ALL, IR(iter + 1)); + DM v_prev(v_j); + DM u_prev(u_j); + Real beta_j = 0; + Real alpha = B.superDiag.Get(iter, 0); + //auto absComp = [](const Real& a, const Real& b) { return Abs(a) < Abs(b); }; + for (Int j = iter + 1; j < k; ++j) { + + //The u step + // apply the operator + u_prev = u_j; + u_j = A * v_prev; + Axpy(-alpha, u_prev, u_j); + beta_j = Nrm2(u_j); + Scale(F(1.0)/beta_j, u_j); + + //reorthData.tau = Max(reorthData.tau, eps * ( alpha + beta )); + ////compute omega recurrence + //auto foundInaccurateU = omegaRecurrenceU(j, reorthData.tau, B, nuList, alpha, beta, ctrl, muList); + //maxMuList.emplace_back(*std::max_element(muList.begin(), muList.end(), absComp)); + //muList.emplace_back(1); + //if ( reorthB || foundInaccurateU ) { + numVecsReorth += Reorthogonalize(U, j, muList, ctrl, u_j, B.superDiag); + U(ALL, IR(j)) = u_j; + //} + + v_prev = v_j; + v_j = AAdj * u_j; + Axpy(-beta_j, v_prev, v_j); + Real alpha = Nrm2(v_j); + Scale(F(1.0)/alpha, v_j); + + //reorthData.tau = Max(reorthData.tau, eps * ( alpha + beta )); + //bool foundInaccurate; + //Real maxElt; + //auto foundInaccurateV = omegaRecurrenceV(j, reorthData.tau, B, muList, alpha, beta, ctrl, nuList); + //maxNuList.emplace_back(*std::max_element(nuList.begin(), nuList.end(), absComp)); + ////Reorthogonalize if necessary. + //if ( reorthB || foundInaccurateV ) { + numVecsReorth += Reorthogonalize(V, j, nuList, ctrl, v_j, B.mainDiag); + V(ALL, IR(j)) = v_j; + //} + } + BidiagInfo info; + info.numVecsReorth = numVecsReorth; + return info; +} + +template +El::Matrix BidiagonalSVD(const BidiagonalMatrix& B, + Int n, Int numColsVT, Int numRowsU, + El::DistMatrix& U, El::DistMatrix VT) { + auto s = B.mainDiag; + auto superDiagCopy = B.superDiag; + lapack::BidiagSVDQRAlg('U', n, numColsVT, numRowsU, + s.Buffer(), superDiagCopy.Buffer(), + VT.Buffer(), VT.LDim(), + U.Buffer(), U.LDim()); + return s; +} + +template +bool hasConverged(Int nVals, + Int newIndex, + const Matrix >& s, + const Matrix >& sOld, + const DistMatrix & VTHat, + const DistMatrix & UHat, + const El::Base& beta, + const BidiagCtrl& ctrl) { + for (Int j = 0; j < nVals; ++j) { + if ( Abs(s.Get(j, 0) - sOld.Get(j, 0)) > ctrl.convValTol ) { + return false; + } + if ( beta * Abs(VTHat.Get(j, newIndex)) > ctrl.convVecTol ) { + return false; + } + if ( beta * Abs(UHat.Get(newIndex, j)) > ctrl.convVecTol ) { + return false; + } + } + return true; +} + +template< typename F, + typename ForwardOperator, + typename AdjointOperator> +ReorthogonalizationData InitialIteration(const ForwardOperator& A, const AdjointOperator& AAdj, + El::AbstractDistMatrix& initialVec, + DistMatrix& UTilde, DistMatrix& VTilde, + BidiagonalMatrix& B){ + typedef El::Base Real; + + const Real eps = El::limits::Epsilon(); + + Real initNorm = Nrm2(initialVec); + Scale(Real(1)/initNorm, initialVec); + UTilde(El::ALL, 0) = initialVec; + DistMatrix x = UTilde(El::ALL, IR(0)); + VTilde(El::ALL, 0) = AAdj * x; + + auto v = VTilde(El::ALL, 0); + + Real alpha = El::Nrm2(v); //record the norm, for reorth. + El::Scale(Real(1)/alpha, v); //make v a unit vector + + UTilde(El::ALL, 1) = A * v; + + auto u0 = UTilde(El::ALL, IR(0)); + auto u1 = UTilde(El::ALL, IR(1)); + + El::Axpy(-alpha, u0, u1); + Real beta = El::Nrm2(u1); + El::Scale(Real(1)/beta,u1); + auto tau = eps*(alpha + beta); + B.mainDiag.Set(0,0, alpha); + B.superDiag.Set(0,0, beta); + //1) Compute nu, a tolerance for reorthoganlization + Real nu = 1 + tau/alpha; //think of nu = 1+\epsilon1 + Real mu = tau/beta; + //first argument is maxIter which at the time of typing this was VTilde.Height and UTilde.Height(). + //If this changes revisit. + ReorthogonalizationData reorthData(VTilde.Height(), nu, mu, tau); + + return reorthData; +} +/** +* TSVD +*/ +template< typename F, + typename ForwardOperator, + typename AdjointOperator> +inline El::Int +TSVD( + const ForwardOperator& A, + const AdjointOperator& AAdj, + El::Int nVals, + El::AbstractDistMatrix& U, + El::AbstractDistMatrix& S, + El::AbstractDistMatrix& V, + El::AbstractDistMatrix& initialVec, //should be of length m + El::Int maxIter=El::Int(1000), + const El::BidiagCtrl& ctrl=El::BidiagCtrl()) //requires InfiniteCharacteristic && ForwardOperator::ScalarType is F + { + typedef El::DistMatrix DM; + typedef El::Matrix> RealMatrix; + typedef El::Base Real; + + const Real eps = El::limits::Epsilon(); + const El::Grid& g = initialVec.Grid(); + const El::Int m = A.Height(); + const El::Int n = A.Width(); + + maxIter = El::Min(maxIter,El::Min(m,n)); + + DM VTilde( n, maxIter, g); + DM UTilde( m, maxIter, g); + + BidiagonalMatrix B(maxIter+1); + auto reorthData = InitialIteration(A, AAdj, initialVec, UTilde, VTilde, B); + El::Display(UTilde, "UTilde"); + El::Display(VTilde, "VTilde"); + + RealMatrix sOld( maxIter+1, 1); + + El::Int blockSize = 1;//El::Max(nVals, ctrl.minBlockSize); + for(Int i = 0; i < maxIter; i+=blockSize){ + Int numSteps = Min( blockSize, maxIter-i); + Int k = i+numSteps; + //Write U_kA'V_k = B_k where k is i+numStep and B_k is diagonal. + GolubKahan( i, k, A, AAdj, B, UTilde, VTilde, reorthData, ctrl); + El::Display(UTilde(ALL, Range(0,k)), "UTilde"); + El::Display(VTilde(ALL, Range(0,k)), "VTilde"); + + //C++17 this should be auto [UHat, S, VThat] = BidiagonalSVD(B,n,m); + DM VTHat(k, nVals, g); + DM UHat(nVals, k, g); + auto s = BidiagonalSVD(B,k,nVals,nVals,UHat,VTHat); + Display(UHat, "UHat"); + Display(VTHat, "VTHat"); + Real beta = B.superDiag.Get(k-1,0); + if( hasConverged(nVals, k, s, sOld,VTHat, UHat, beta, ctrl)){ + std::cout << "Apparently, we have converged" << std::endl; + El::Gemm(El::NORMAL, El::NORMAL, Real(1), UTilde, UHat, Real(0), U); + El::Gemm(El::NORMAL, El::ADJOINT, Real(1), VTilde, VTHat, Real(0), V); + DistMatrix S_STAR_STAR( S.Grid()); + S_STAR_STAR.LockedAttach( S.Height(), S.Width(), S.Grid(), S.ColAlign(), S.RowAlign(), s.Buffer(), s.LDim(), S.Root()); + Copy(S_STAR_STAR, S); + return k; + } + sOld = s; + reorthData.tau = eps*s.Get(0,0); + } + return Int(-1); +} + +/** +* TSVD don't provide initial guess. +*/ +template< typename F, + typename ForwardOperator, + typename AdjointOperator> +inline Int +TSVD( const ForwardOperator& A, + const AdjointOperator& AAdj, + Int nVals, + AbstractDistMatrix& U, + AbstractDistMatrix& S, + AbstractDistMatrix& V, + Int maxIter=Int(1000), + const BidiagCtrl& ctrl=BidiagCtrl()) { + DistMatrix initialVec(U.Grid()); + F center = rand(); + Uniform( initialVec, A.Height(), 1, center); + return TSVD(A, AAdj, nVals, U,S,V, initialVec, maxIter, ctrl); +} + +} //end namespace El +#endif //ifndef EL_LAPACKLIKE_SPECTRAL_TSVD diff --git a/tests/lapack_like/TSVD.cpp b/tests/lapack_like/TSVD.cpp new file mode 100644 index 0000000000..83c29466be --- /dev/null +++ b/tests/lapack_like/TSVD.cpp @@ -0,0 +1,128 @@ +/* + Copyright (c) 2009-2015, Jack Poulson + All rights reserved. + + This file is part of Elemental and is under the BSD 2-Clause License, + which can be found in the LICENSE file in the root directory, or at + http://opensource.org/licenses/BSD-2-Clause +*/ +#include "El.hpp" +#include +using namespace El; + +template< typename Matrix> +class AOp { + + public: + AOp( const Matrix& A_, Orientation o_ = El::NORMAL) : A( A_), o(o_) {} + + template< typename V> + V operator*( V& v) const { + V w; + Gemv(o, 1.0, A, v, w); + return w; + } + + Int Height() const { + if( o == El::NORMAL){ return A.Height(); } + return A.Width(); + } + + Int Width() const { + if( o == El::NORMAL){ return A.Width(); } + return A.Height(); + } + + const Matrix& A; + const Orientation o = El::NORMAL; +}; //Matrix Operator + + +template +void TestTSVD( const DistMatrix& A, Int k=3){ + typedef Base Real; + const auto& g = A.Grid(); + int m = A.Height(); + int n = A.Width(); + DistMatrix Uk( g), Sk( g), Vk( g); + DistMatrix U( g), V( g); + DistMatrix, VR, STAR> S( g); + mpi::Barrier( g.Comm() ); + if( g.Rank() == 0 ) + Output(" Starting SVD factorization..."); + SVD( A, U, S, V); + mpi::Barrier( g.Comm() ); + std::cout << "Singular Values: " << std::flush; + Display(S(Range(0,k), 0), "Singular Values"); + std::cout << std::endl; + if( g.Rank() == 0 ) + Output(" Starting TSVD factorization..."); + const double startTime = mpi::Time(); + typedef AOp> Operator; + Operator Aop( A); + Operator AAdjop( A, El::ADJOINT); + TSVD( Aop, AAdjop, k, Uk, Sk, Vk, 5); + const double runTime = mpi::Time() - startTime; + Output("TSVD Time = ",runTime); + Real eps = limits::Epsilon(); + for( Int i = 0; i < k; ++i){ + std::cout << i << ": "; + const Real& A = Sk.Get(i,0); + const Real& B = S.Get(i,0); + Output( "Sk = ", A, "SVD: S_k", B); + auto eigenvalue_error = Abs(A - B); + if( eigenvalue_error > eps){ + RuntimeError("Convergence Issue in TSVD. Eigenvalues are not close enough."); + } + F left_svec_error = F(0); + F right_svec_error = F(0); + for (int j = 0; j <= m; ++j){ + auto t = Abs(Uk.Get( j, i)) - Abs(U.Get( j, i)); + left_svec_error += t*t; + auto s = Abs(Vk.Get( j, i)) - Abs(V.Get( j, i)); + right_svec_error += s*s; + } + left_svec_error = Sqrt(left_svec_error); + right_svec_error = Sqrt(right_svec_error); + if( left_svec_error > eps || right_svec_error > eps){ + RuntimeError("Convergence Issue in TSVD. Eigenvector Entries do not match."); + } + } +} + +int +main( int argc, char* argv[] ) +{ + Environment env( argc, argv ); + mpi::Comm comm = mpi::COMM_WORLD; + const Int commRank = mpi::Rank( comm ); + + try + { + const bool colMajor = Input("--colMajor","column-major ordering?",true); + const Int m = Input("--height","height of matrix",3); + const Int n = Input("--width","width of matrix",3); + const Int nb = Input("--nb","algorithmic blocksize",96); + const bool testCorrectness = Input + ("--correctness","test correctness?",true); + const bool print = Input("--print","print matrices?",false); + ProcessInput(); + PrintInputReport(); + + const GridOrder order = ( colMajor ? COLUMN_MAJOR : ROW_MAJOR ); + const Grid g( comm, order ); + SetBlocksize( nb ); + ComplainIfDebug(); + if( commRank == 0 ) + Output("Will test TSQR"); + + if( commRank == 0 ) + Output("Testing with doubles:"); + DistMatrix A; + Laplacian(A, m, n); + TestTSVD( A, 3); + } + catch( exception& e ) { ReportException(e); } + + return 0; +} diff --git a/travis/install-mpi.sh b/travis/install-mpi.sh index 0d3513c141..31532abd2a 100644 --- a/travis/install-mpi.sh +++ b/travis/install-mpi.sh @@ -2,7 +2,7 @@ set -e case $1 in mpich) set -x; - sudo apt-get install -q mpich libmpich-dev;; + sudo apt-get install -q --force-yes mpich libmpich-dev;; openmpi) set -x; sudo apt-get install openmpi-bin openmpi-dev;; *)