forked from elemental/Elemental
-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Debugged recent modifications to tridiagonalization and reimplemented…
… correctness checks. --HG-- rename : src/lapack/Tridiag/ColReflector.cpp => src/lapack/Reflector/ColReflector.cpp rename : src/lapack/Tridiag/Reflector.cpp => src/lapack/Reflector/Reflector.cpp rename : src/lapack/Tridiag/RowReflector.cpp => src/lapack/Reflector/RowReflector.cpp
- Loading branch information
Showing
15 changed files
with
458 additions
and
259 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
File renamed without changes.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,205 @@ | ||
| /* | ||
| This file is part of Elemental, a library for distributed-memory dense | ||
| linear algebra. | ||
| Copyright (c) 1992-2008 The University of Tennessee | ||
| All rights reserved. | ||
| Copyright (c) 2009-2010 Jack Poulson <jack.poulson@gmail.com>. | ||
| All rights reserved. | ||
| This file is released under the terms of the license contained in the file | ||
| LICENSE-LAPACK-DERIV as it is a derivative work of the LAPACK routines | ||
| dlarfg.f and zlarfg.f. | ||
| */ | ||
| #include "elemental/lapack_internal.hpp" | ||
| using namespace std; | ||
| using namespace elemental; | ||
|
|
||
| template<typename R> | ||
| R | ||
| elemental::lapack::Reflector | ||
| ( Matrix<R>& chi, Matrix<R>& x ) | ||
| { | ||
| #ifndef RELEASE | ||
| PushCallStack("lapack::Reflector"); | ||
| #endif | ||
| if( x.Height() == 0 ) | ||
| { | ||
| chi(0,0) *= (R)-1; | ||
| #ifndef RELEASE | ||
| PopCallStack(); | ||
| #endif | ||
| return (R)2; | ||
| } | ||
|
|
||
| R norm = blas::Nrm2( x ); | ||
| R alpha = chi(0,0); | ||
|
|
||
| R beta; | ||
| if( alpha <= 0 ) | ||
| beta = wrappers::lapack::SafeNorm( alpha, norm ); | ||
| else | ||
| beta = -wrappers::lapack::SafeNorm( alpha, norm ); | ||
|
|
||
| R safeMin = numeric_limits<R>::min() / numeric_limits<R>::epsilon(); | ||
| int count = 0; | ||
| if( Abs( beta ) < safeMin ) | ||
| { | ||
| R invOfSafeMin = static_cast<R>(1) / safeMin; | ||
| do | ||
| { | ||
| ++count; | ||
| blas::Scal( invOfSafeMin, x ); | ||
| alpha *= invOfSafeMin; | ||
| beta *= invOfSafeMin; | ||
| } while( Abs( beta ) < safeMin ); | ||
|
|
||
| norm = blas::Nrm2( x ); | ||
| if( alpha <= 0 ) | ||
| beta = wrappers::lapack::SafeNorm( alpha, norm ); | ||
| else | ||
| beta = -wrappers::lapack::SafeNorm( alpha, norm ); | ||
| } | ||
|
|
||
| R tau = ( beta - alpha ) / beta; | ||
| blas::Scal( static_cast<R>(1)/(alpha-beta), x ); | ||
|
|
||
| for( int j=0; j<count; ++j ) | ||
| beta *= safeMin; | ||
| chi(0,0) = beta; | ||
| #ifndef RELEASE | ||
| PopCallStack(); | ||
| #endif | ||
| return tau; | ||
| } | ||
|
|
||
| #ifndef WITHOUT_COMPLEX | ||
| template<typename R> | ||
| complex<R> | ||
| elemental::lapack::Reflector | ||
| ( Matrix< complex<R> >& chi, Matrix< complex<R> >& x ) | ||
| { | ||
| #ifndef RELEASE | ||
| PushCallStack("lapack::Reflector"); | ||
| #endif | ||
| typedef complex<R> C; | ||
|
|
||
| R norm = blas::Nrm2( x ); | ||
| C alpha = chi(0,0); | ||
|
|
||
| if( norm == 0 && imag(alpha) == (R)0 ) | ||
| { | ||
| chi(0,0) *= (R)-1; | ||
| #ifndef RELEASE | ||
| PopCallStack(); | ||
| #endif | ||
| return (C)2; | ||
| } | ||
|
|
||
| R beta; | ||
| if( real(alpha) <= 0 ) | ||
| beta = wrappers::lapack::SafeNorm( real(alpha), imag(alpha), norm ); | ||
| else | ||
| beta = -wrappers::lapack::SafeNorm( real(alpha), imag(alpha), norm ); | ||
|
|
||
| R safeMin = numeric_limits<R>::min() / numeric_limits<R>::epsilon(); | ||
| int count = 0; | ||
| if( Abs( beta ) < safeMin ) | ||
| { | ||
| R invOfSafeMin = static_cast<R>(1) / safeMin; | ||
| do | ||
| { | ||
| ++count; | ||
| blas::Scal( (C)invOfSafeMin, x ); | ||
| alpha *= invOfSafeMin; | ||
| beta *= invOfSafeMin; | ||
| } while( Abs( beta ) < safeMin ); | ||
|
|
||
| norm = blas::Nrm2( x ); | ||
| if( real(alpha) <= 0 ) | ||
| beta = wrappers::lapack::SafeNorm | ||
| ( real(alpha), imag(alpha), norm ); | ||
| else | ||
| beta = -wrappers::lapack::SafeNorm | ||
| ( real(alpha), imag(alpha), norm ); | ||
| } | ||
|
|
||
| C tau = C( (beta-real(alpha))/beta, -imag(alpha)/beta ); | ||
| blas::Scal( static_cast<C>(1)/(alpha-beta), x ); | ||
|
|
||
| for( int j=0; j<count; ++j ) | ||
| beta *= safeMin; | ||
| chi(0,0) = beta; | ||
| #ifndef RELEASE | ||
| PopCallStack(); | ||
| #endif | ||
| return tau; | ||
| } | ||
| #endif // WITHOUT_COMPLEX | ||
|
|
||
| template<typename T> | ||
| T | ||
| elemental::lapack::Reflector | ||
| ( DistMatrix<T,MC,MR>& chi, DistMatrix<T,MC,MR>& x ) | ||
| { | ||
| #ifndef RELEASE | ||
| PushCallStack("lapack::Reflector"); | ||
| if( chi.GetGrid() != x.GetGrid() ) | ||
| throw logic_error( "chi and x must be distributed over the same grid" ); | ||
| if( chi.Height() != 1 || chi.Width() != 1 ) | ||
| throw logic_error( "chi must be a scalar." ); | ||
| if( x.Height() != 1 && x.Width() != 1 ) | ||
| throw logic_error( "x must be a vector." ); | ||
| #endif | ||
| if( max( x.Height(), x.Width() ) == 0 ) | ||
| return (T)0; | ||
|
|
||
| const Grid& g = x.GetGrid(); | ||
| T tau; | ||
| if( x.Width() == 1 ) | ||
| { | ||
| const bool thisIsMyColumn = ( g.MRRank() == x.RowAlignment() ); | ||
| if( thisIsMyColumn ) | ||
| tau = lapack::internal::ColReflector( chi, x ); | ||
| wrappers::mpi::Broadcast( &tau, 1, x.RowAlignment(), g.MRComm() ); | ||
| } | ||
| else | ||
| { | ||
| const bool thisIsMyRow = ( g.MCRank() == x.ColAlignment() ); | ||
| if( thisIsMyRow ) | ||
| tau = lapack::internal::RowReflector( chi, x ); | ||
| wrappers::mpi::Broadcast( &tau, 1, x.ColAlignment(), g.MCComm() ); | ||
| } | ||
| #ifndef RELEASE | ||
| PopCallStack(); | ||
| #endif | ||
| return tau; | ||
| } | ||
|
|
||
| template float elemental::lapack::Reflector | ||
| ( Matrix<float>& chi, Matrix<float>& x ); | ||
|
|
||
| template float elemental::lapack::Reflector | ||
| ( DistMatrix<float,MC,MR>& chi, DistMatrix<float,MC,MR>& x ); | ||
|
|
||
| template double elemental::lapack::Reflector | ||
| ( Matrix<double>& chi, Matrix<double>& x ); | ||
|
|
||
| template double elemental::lapack::Reflector | ||
| ( DistMatrix<double,MC,MR>& chi, DistMatrix<double,MC,MR>& x ); | ||
|
|
||
| #ifndef WITHOUT_COMPLEX | ||
| template scomplex elemental::lapack::Reflector | ||
| ( Matrix<scomplex>& chi, Matrix<scomplex>& x ); | ||
|
|
||
| template scomplex elemental::lapack::Reflector | ||
| ( DistMatrix<scomplex,MC,MR>& chi, DistMatrix<scomplex,MC,MR>& x ); | ||
|
|
||
| template dcomplex elemental::lapack::Reflector | ||
| ( Matrix<dcomplex>& chi, Matrix<dcomplex>& x ); | ||
|
|
||
| template dcomplex elemental::lapack::Reflector | ||
| ( DistMatrix<dcomplex,MC,MR>& chi, DistMatrix<dcomplex,MC,MR>& x ); | ||
| #endif | ||
|
|
File renamed without changes.
Oops, something went wrong.