From 0c14bcea33e22ed59df9a423cb4ce62dfb5150aa Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Tue, 5 Apr 2022 11:14:49 +0200 Subject: [PATCH 1/2] initial version --- include/alp/algorithms/hessenberg_qr.cpp | 163 ++++++++++++++++++++ include/alp/algorithms/svd.cpp | 184 +++++++++++++++++++++++ 2 files changed, 347 insertions(+) create mode 100644 include/alp/algorithms/hessenberg_qr.cpp create mode 100644 include/alp/algorithms/svd.cpp diff --git a/include/alp/algorithms/hessenberg_qr.cpp b/include/alp/algorithms/hessenberg_qr.cpp new file mode 100644 index 0000000..422027f --- /dev/null +++ b/include/alp/algorithms/hessenberg_qr.cpp @@ -0,0 +1,163 @@ + +/* + * Copyright 2021 Huawei Technologies Co., Ltd. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include + +#include + +using namespace alp; + +/** + * @brief Computes Hessenberg QR decomposition \f$H = QR\f$ + * where \a H is tridiagonal symmetric matrix, + * where \a Q ortogonal matrix + * where \a R upper triangular matrix. + * + * @tparam D Data element type + * @tparam Ring Type of the semiring used in the computation + * @tparam Minus Type minus operator used in the computation + * @tparam Divide Type of divide operator used in the computation + * @param[in] H input tridiagonal symmetric matrix H + * @param[out] R output upper triangular matrix (R) + * @param[out] Q output orthogonal matrix + * @return RC SUCCESS if the execution was correct + * + */ +template< + typename D = double, + class Ring = Semiring< operators::add< D >, operators::mul< D >, identities::zero, identities::one >, + class Minus = operators::subtract< D >, + class Divide = operators::divide< D > > +RC hessenberg_qr( + Matrix< D, structures::Hessenberg, Dense > &H, + Matrix< D, structures::Orthogonal, Dense > &Q, + Matrix< D, structures::UpperTriangular, Dense > &R, + const Ring & ring = Ring(), + const Minus & minus = Minus(), + const Divide & divide = Divide() ) { + + const Scalar< D > zero( ring.template getZero< D >() ); + + // n=len(H) + const size_t n = nrows( H ); + + // Q=identity(n) + rc = set( Q , structures::constant::I( n ) ); + + // for i in range(0,n-1): + for( size_t i = 0; i < n - 1; ++i ) { + // c,s=zrotg(*(H[i:i+2,i])) + zrotg(H[i,i],H[i+1,i],c,s,) ; //?? what to do with this + //D=[[c,s],[-conjugate(s),c]] + + // R=identity(n).astype(complex) + Matrix< D, structures::Symmetric, Dense > R( n ); + rc = set( Qk, structures::constant::I( n ) ); + + // R[i:i+2,i:i+2]=D + auto Rview = get_view( R, utils::range( i , i + 2 ), utils::range( i , i + 2 ) ); + rc = set( Rview, D ); + + // H=R.dot(H) + Matrix< D, structures::Square, Dense > Tmp( n ); + rc = set( Tmp, zero ); + rc = mxm( Tmp, R, H, ring ); + rc = set( H, Tmp ); + + // Q=Q.dot(conjugate(R.T)) + rc = set( Tmp, zero ); + auto Rt_view = get_view< view::transpose > + ( R ); + rc = mxm( Tmp, Q, Rt_view, ring ); + rc = set( Q, Tmp ); + + } + // return(H,Q) + + + RC rc = SUCCESS; + + + return rc; +} + +void alp_program( const size_t & unit, alp::RC & rc ) { + alp::Semiring< alp::operators::add< double >, alp::operators::mul< double >, alp::identities::zero, alp::identities::one > ring; + + // dimensions of rectangular matrix H + size_t N = 10 * unit; + + // dimensions of views over H + size_t n = 2 * unit; + + alp::Matrix< double, structures::Hessenberg > H( N ); + alp::Matrix< double, structures::Orthogonal > Q( N ); + alp::Matrix< D, structures::UpperTriangular, Dense > R( N ); + + rc = hessenberg_qr( H, Q, R, ring ); + +} + +int main( int argc, char ** argv ) { + // defaults + bool printUsage = false; + size_t in = 100; + + // error checking + if( argc > 2 ) { + printUsage = true; + } + if( argc == 2 ) { + size_t read; + std::istringstream ss( argv[ 1 ] ); + if( ! ( ss >> read ) ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else if( ! ss.eof() ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else if( read % 2 != 0 ) { + std::cerr << "Given value for n is odd\n"; + printUsage = true; + } else { + // all OK + in = read; + } + } + if( printUsage ) { + std::cerr << "Usage: " << argv[ 0 ] << " [n]\n"; + std::cerr << " -n (optional, default is 100): an even integer, the " + "test size.\n"; + return 1; + } + + std::cout << "This is functional test " << argv[ 0 ] << "\n"; + alp::Launcher< AUTOMATIC > launcher; + alp::RC out; + if( launcher.exec( &alp_program, in, out, true ) != SUCCESS ) { + std::cerr << "Launching test FAILED\n"; + return 255; + } + if( out != SUCCESS ) { + std::cerr << "Test FAILED (" << alp::toString( out ) << ")" << std::endl; + } else { + std::cout << "Test OK" << std::endl; + } + return 0; +} diff --git a/include/alp/algorithms/svd.cpp b/include/alp/algorithms/svd.cpp new file mode 100644 index 0000000..59617e1 --- /dev/null +++ b/include/alp/algorithms/svd.cpp @@ -0,0 +1,184 @@ + +/* + * Copyright 2021 Huawei Technologies Co., Ltd. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include + +#include + +using namespace alp; + +/** + * @brief Computes Householder tridiagonalization \f$H = USVh^T\f$ + * where \a H is real rectangular n x m matrix, where n>=m + * where \a U is real rectangular orthogonal n x m matrix + * where \a S is positive diagoanl m x m matrix, with + * singular values s_1 >= s_2 >= >= s_m on diagonal + * where \a Vh is real square orthogonal m x m matrix . + * + * @tparam D Data element type + * @tparam Ring Type of the semiring used in the computation + * @tparam Minus Type minus operator used in the computation + * @tparam Divide Type of divide operator used in the computation + * @param[out] U output real rectangular orthogonal n x m matrix + * @param[out] S output is positive diagoanl m x m matrix + * @param[out] Vh output real square orthogonal m x m matrix + * @param[in] H input real rectangular n x m matrix, where n>=m + * @param[in] ring A semiring for operations + * @return RC SUCCESS if the execution was correct + * + */ +template< + typename D = double, + class Ring = Semiring< operators::add< D >, operators::mul< D >, identities::zero, identities::one >, + class Minus = operators::subtract< D >, + class Divide = operators::divide< D > > +RC svd( + Matrix< D, structures::Orthogonal, Dense > &U, + Matrix< D, structures::Diagonal, Dense > &S, + Matrix< D, structures::SquareOrthogonal, Dense > &Vh, + const Matrix< D, structures::General, Dense > &H, + size_t imax = 100, + const Ring & ring = Ring(), + const Minus & minus = Minus(), + const Divide & divide = Divide() ) { + + + // n,m=shape(H) + const size_t n = nrows( H ); + const size_t m = ncolumns( H ); + + + + // Atmp=zeros((n+m,n+m)) + // Atmp[-n:,:m]=H + // A=A_r+A_r.T + Matrix< D, structures::Symmetric, Dense > Atmp( (n+m,n+m) ); + ... + + // Q,Atmp=house_trid(Atmp) + rc = zhetrd( Q, T, H, ring ); + + // Atmp=Atmp+identity((n+m)) + + // for i in range(imax): + // R,QQ=Hessenberg_QR(Atmp) + // Atmp=((R)).dot(QQ) + // Q=Q.dot(QQ) + for( size_t k = 0; k < imax; ++k ) { + hessenberg_qr( H, Q, ring ) + + } + + + // this strange sort should give permuataions (ii): s1,s2,s3...sm,-s1,-s2,-s3...-sm,0,0,0 + // e=real(diag(Atmp))-1 + // ####### + // ii=argsort(abs(e))[::-1] + // jj=abs(e[ii])>1.e-10 + // kk=argsort(e[ii][jj])[::-1] + // kk[len(kk)//2:]=kk[len(kk)//2:][::-1] + // ii[:len(kk)]=ii[kk] + + + // use premuations ii to get submatrices + // ####### + // e=e[ii] + // v=Q[:,ii] + // ####### + // Vh=sqrt(2)*v[:m,:m].T + // U=sqrt(2)*v[-n:,:n] + // ####### + // return(U[:, :m],e[:m],Vh) + + + + RC rc = SUCCESS; + + + return rc; +} + +void alp_program( const size_t & unit, alp::RC & rc ) { + alp::Semiring< alp::operators::add< double >, alp::operators::mul< double >, alp::identities::zero, alp::identities::one > ring; + + // dimensions of rectangular matrix H + size_t N = 10 * unit; + size_t M = 10 * unit; + + // dimensions of views over H + size_t n = 2 * unit; + size_t m = 2 * unit; + + alp::Matrix< double, structures::General > H( N, M ); + + alp::Matrix< D, structures::Orthogonal, Dense > U ( N, M ); + alp::Matrix< D, structures::Diagonal, Dense > S ( M ); + alp::Matrix< D, structures::SquareOrthogonal, Dense > Vh( M, M); + + rc = svd( U, S, Vh, H, (size_t)100, ring ); + +} + +int main( int argc, char ** argv ) { + // defaults + bool printUsage = false; + size_t in = 100; + + // error checking + if( argc > 2 ) { + printUsage = true; + } + if( argc == 2 ) { + size_t read; + std::istringstream ss( argv[ 1 ] ); + if( ! ( ss >> read ) ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else if( ! ss.eof() ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else if( read % 2 != 0 ) { + std::cerr << "Given value for n is odd\n"; + printUsage = true; + } else { + // all OK + in = read; + } + } + if( printUsage ) { + std::cerr << "Usage: " << argv[ 0 ] << " [n]\n"; + std::cerr << " -n (optional, default is 100): an even integer, the " + "test size.\n"; + return 1; + } + + std::cout << "This is functional test " << argv[ 0 ] << "\n"; + alp::Launcher< AUTOMATIC > launcher; + alp::RC out; + if( launcher.exec( &alp_program, in, out, true ) != SUCCESS ) { + std::cerr << "Launching test FAILED\n"; + return 255; + } + if( out != SUCCESS ) { + std::cerr << "Test FAILED (" << alp::toString( out ) << ")" << std::endl; + } else { + std::cout << "Test OK" << std::endl; + } + return 0; +} -- Gitee From 49f7d44b7357f3d8e37fb3c00314c96d9d7bdf43 Mon Sep 17 00:00:00 2001 From: Daniele Giuseppe Spampinato Date: Tue, 5 Apr 2022 12:15:17 +0200 Subject: [PATCH 2/2] Minor fixes --- include/alp/algorithms/hessenberg_qr.cpp | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/include/alp/algorithms/hessenberg_qr.cpp b/include/alp/algorithms/hessenberg_qr.cpp index 422027f..7abc668 100644 --- a/include/alp/algorithms/hessenberg_qr.cpp +++ b/include/alp/algorithms/hessenberg_qr.cpp @@ -59,31 +59,33 @@ RC hessenberg_qr( // Q=identity(n) rc = set( Q , structures::constant::I( n ) ); - + + Matrix< D, structures::General, Dense > Tmp( n ); + // for i in range(0,n-1): for( size_t i = 0; i < n - 1; ++i ) { // c,s=zrotg(*(H[i:i+2,i])) zrotg(H[i,i],H[i+1,i],c,s,) ; //?? what to do with this //D=[[c,s],[-conjugate(s),c]] + // WE could express the below as: auto G = structures::constant::Givens(a, b); // R=identity(n).astype(complex) Matrix< D, structures::Symmetric, Dense > R( n ); - rc = set( Qk, structures::constant::I( n ) ); + rc = set( R, structures::constant::I( n ) ); // R[i:i+2,i:i+2]=D auto Rview = get_view( R, utils::range( i , i + 2 ), utils::range( i , i + 2 ) ); rc = set( Rview, D ); - - // H=R.dot(H) - Matrix< D, structures::Square, Dense > Tmp( n ); + // closes Givens comment + + // H=R.dot(H) rc = set( Tmp, zero ); rc = mxm( Tmp, R, H, ring ); rc = set( H, Tmp ); // Q=Q.dot(conjugate(R.T)) rc = set( Tmp, zero ); - auto Rt_view = get_view< view::transpose > - ( R ); + auto Rt_view = get_view< view::transpose >( R ); rc = mxm( Tmp, Q, Rt_view, ring ); rc = set( Q, Tmp ); -- Gitee