Skip to content

Commit

Permalink
cond: return value instead of output argument
Browse files Browse the repository at this point in the history
  • Loading branch information
@mgates3
mgates3 committed Nov 15, 2023
1 parent 86bb586 commit 17b05b6
Show file tree
Hide file tree
Showing 8 changed files with 130 additions and 119 deletions.
23 changes: 10 additions & 13 deletions include/slate/simplified_api.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -363,14 +363,13 @@ void lu_inverse_using_factor_out_of_place(

// gecondest
template <typename scalar_t>
void lu_condest_using_factor(
Norm norm,
blas::real_type<scalar_t> lu_rcondest_using_factor(
Norm in_norm,
Matrix<scalar_t>& A,
blas::real_type<scalar_t> Anorm,
blas::real_type<scalar_t> *rcond,
Options const& opts = Options())
{
gecondest(norm, A, Anorm, rcond, opts);
return gecondest( in_norm, A, Anorm, opts );
}

//-----------------------------------------
Expand Down Expand Up @@ -495,14 +494,13 @@ void chol_inverse_using_factor(

// pocondest
template <typename scalar_t>
void chol_condest_using_factor(
Norm norm,
blas::real_type<scalar_t> chol_rcondest_using_factor(
Norm in_norm,
HermitianMatrix<scalar_t>& A,
blas::real_type<scalar_t> Anorm,
blas::real_type<scalar_t> *rcond,
Options const& opts = Options())
{
pocondest(norm, A, Anorm, rcond, opts);
return pocondest( in_norm, A, Anorm, opts );
}

//-----------------------------------------
Expand Down Expand Up @@ -684,16 +682,15 @@ void lq_multiply_by_q(

// trcondest
template <typename scalar_t>
void triangular_condest(
Norm norm,
blas::real_type<scalar_t> triangular_rcondest(
Norm in_norm,
TriangularMatrix<scalar_t>& A,
blas::real_type<scalar_t> *rcond,
blas::real_type<scalar_t> Anorm,
Options const& opts = Options())
{
trcondest(norm, A, rcond, opts);
return trcondest( in_norm, A, Anorm, opts );
}


//------------------------------------------------------------------------------
// Symmetric/Hermitian Eigenvalues

Expand Down
40 changes: 25 additions & 15 deletions include/slate/slate.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1345,43 +1345,53 @@ void steqr2(
//-----------------------------------------
// gecondest()
template <typename scalar_t>
void gecondest(
Norm in_norm,
Matrix<scalar_t>& A,
blas::real_type<scalar_t> Anorm,
blas::real_type<scalar_t> *rcond,
Options const& opts = Options());
blas::real_type<scalar_t> gecondest(
Norm in_norm,
Matrix<scalar_t>& A,
blas::real_type<scalar_t> Anorm,
Options const& opts = Options());

template <typename scalar_t>
[[deprecated( "Pass Anorm by value instead. Will be removed 2024-11." )]]
[[deprecated( "Use rcond = gecondest(...) and pass Anorm by value instead. Will be removed 2024-11." )]]
void gecondest(
Norm in_norm,
Matrix<scalar_t>& A,
blas::real_type<scalar_t> *Anorm,
blas::real_type<scalar_t> *rcond,
Options const& opts = Options())
{
gecondest( in_norm, A, *Anorm, rcond, opts );
*rcond = gecondest( in_norm, A, *Anorm, opts );
}

//-----------------------------------------
// pocondest()
template <typename scalar_t>
void pocondest(
Norm in_norm,
HermitianMatrix<scalar_t>& A,
blas::real_type<scalar_t> Anorm,
blas::real_type<scalar_t> *rcond,
Options const& opts = Options());
blas::real_type<scalar_t> pocondest(
Norm in_norm,
HermitianMatrix<scalar_t>& A,
blas::real_type<scalar_t> Anorm,
Options const& opts = Options());

//-----------------------------------------
// trcondest()
template <typename scalar_t>
blas::real_type<scalar_t> trcondest(
Norm in_norm,
TriangularMatrix<scalar_t>& A,
blas::real_type<scalar_t> Anorm,
Options const& opts = Options());

template <typename scalar_t>
[[deprecated( "Use rcond = trcondest(...) and pass Anorm. Will be removed 2024-11." )]]
void trcondest(
Norm in_norm,
TriangularMatrix<scalar_t>& A,
blas::real_type<scalar_t> *rcond,
Options const& opts = Options());
Options const& opts = Options())
{
blas::real_type<scalar_t> Anorm = norm( in_norm, A, opts );
*rcond = trcondest( in_norm, A, Anorm, opts );
}

} // namespace slate

Expand Down
49 changes: 24 additions & 25 deletions src/gecondest.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,10 +38,6 @@ namespace slate {
/// If Norm::One, the 1-norm of the original matrix A.
/// If Norm::Inf, the infinity-norm of the original matrix A.
///
/// @param[in,out] rcond
/// The reciprocal of the condition number of the matrix A,
/// computed as stated above.
///
/// @param[in] opts
/// Additional options, as map of name = value pairs. Possible options:
/// - Option::Target:
Expand All @@ -51,15 +47,19 @@ namespace slate {
/// - HostBatch: batched BLAS on CPU host.
/// - Devices: batched BLAS on GPU device.
///
/// @return rcond
/// The reciprocal of the condition number of the matrix A,
/// computed as stated above.
/// rcond is used instead of cond to handle overflow better.
///
/// @ingroup cond_specialization
///
template <typename scalar_t>
void gecondest(
Norm in_norm,
Matrix<scalar_t>& A,
blas::real_type<scalar_t> Anorm,
blas::real_type<scalar_t> *rcond,
Options const& opts)
blas::real_type<scalar_t> gecondest(
Norm in_norm,
Matrix<scalar_t>& A,
blas::real_type<scalar_t> Anorm,
Options const& opts)
{
using blas::real;
using real_t = blas::real_type<scalar_t>;
Expand All @@ -74,17 +74,18 @@ void gecondest(
else {
slate_error("invalid norm.");
}
if (Anorm < 0 || std::isnan( Anorm )) {
slate_error( "invalid Anorm" );
}

int64_t m = A.m();

// Quick return
*rcond = 0.;
if (m <= 1) {
*rcond = 1.;
return;
return 1.;
}
else if (Anorm == 0.) {
return;
return 0.;
}

scalar_t alpha = 1.;
Expand Down Expand Up @@ -118,8 +119,7 @@ void gecondest(
MPI_Bcast( &isave[0], 4, MPI_INT64_T, X.tileRank(0, 0), A.mpiComm() );
MPI_Bcast( &kase, 1, MPI_INT, X.tileRank(0, 0), A.mpiComm() );

while (kase != 0)
{
while (kase != 0) {
if (kase == kase1) {
// Multiply by inv(L).
slate::trsm( Side::Left, alpha, L, X, opts );
Expand All @@ -144,42 +144,41 @@ void gecondest(

// Compute the estimate of the reciprocal condition number.
if (Ainvnorm != 0.0) {
*rcond = (1.0 / Ainvnorm) / Anorm;
return (1.0 / Ainvnorm) / Anorm;
}
else {
return 0.;
}
}

//------------------------------------------------------------------------------
// Explicit instantiations.
template
void gecondest<float>(
float gecondest<float>(
Norm in_norm,
Matrix<float>& A,
float Anorm,
float *rcond,
Options const& opts);

template
void gecondest<double>(
double gecondest<double>(
Norm in_norm,
Matrix<double>& A,
double Anorm,
double *rcond,
Options const& opts);

template
void gecondest< std::complex<float> >(
float gecondest< std::complex<float> >(
Norm in_norm,
Matrix< std::complex<float> >& A,
float Anorm,
float *rcond,
Options const& opts);

template
void gecondest< std::complex<double> >(
double gecondest< std::complex<double> >(
Norm in_norm,
Matrix< std::complex<double> >& A,
double Anorm,
double *rcond,
Options const& opts);

} // namespace slate
49 changes: 24 additions & 25 deletions src/pocondest.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,10 +39,6 @@ namespace slate {
/// If Norm::One, the 1-norm of the original matrix A.
/// If Norm::Inf, the infinity-norm of the original matrix A.
///
/// @param[in,out] rcond
/// The reciprocal of the condition number of the matrix A,
/// computed as stated above.
///
/// @param[in] opts
/// Additional options, as map of name = value pairs. Possible options:
/// - Option::Target:
Expand All @@ -52,15 +48,19 @@ namespace slate {
/// - HostBatch: batched BLAS on CPU host.
/// - Devices: batched BLAS on GPU device.
///
/// @return rcond
/// The reciprocal of the condition number of the matrix A,
/// computed as stated above.
/// rcond is used instead of cond to handle overflow better.
///
/// @ingroup cond_specialization
///
template <typename scalar_t>
void pocondest(
Norm in_norm,
HermitianMatrix<scalar_t>& A,
blas::real_type<scalar_t> Anorm,
blas::real_type<scalar_t> *rcond,
Options const& opts)
blas::real_type<scalar_t> pocondest(
Norm in_norm,
HermitianMatrix<scalar_t>& A,
blas::real_type<scalar_t> Anorm,
Options const& opts)
{
using blas::real;
using real_t = blas::real_type<scalar_t>;
Expand All @@ -69,17 +69,18 @@ void pocondest(
if (in_norm != Norm::One && in_norm != Norm::Inf) {
slate_error("invalid norm.");
}
if (Anorm < 0 || std::isnan( Anorm )) {
slate_error( "invalid Anorm" );
}

int64_t m = A.m();

// Quick return
*rcond = 0.;
if (m <= 1) {
*rcond = 1.;
return;
return 1.;
}
else if (Anorm == 0.) {
return;
return 0.;
}

real_t Ainvnorm = 0.0;
Expand Down Expand Up @@ -107,8 +108,7 @@ void pocondest(
MPI_Bcast( &isave[0], 4, MPI_INT64_T, X.tileRank(0, 0), A.mpiComm() );
MPI_Bcast( &kase, 1, MPI_INT, X.tileRank(0, 0), A.mpiComm() );

while (kase != 0)
{
while (kase != 0) {
// A is symmetric, so both cases are equivalent
potrs( A, X, opts );

Expand All @@ -119,42 +119,41 @@ void pocondest(

// Compute the estimate of the reciprocal condition number.
if (Ainvnorm != 0.0) {
*rcond = (1.0 / Ainvnorm) / Anorm;
return (1.0 / Ainvnorm) / Anorm;
}
else {
return 0.;
}
}

//------------------------------------------------------------------------------
// Explicit instantiations.
template
void pocondest<float>(
float pocondest<float>(
Norm in_norm,
HermitianMatrix<float>& A,
float Anorm,
float *rcond,
Options const& opts);

template
void pocondest<double>(
double pocondest<double>(
Norm in_norm,
HermitianMatrix<double>& A,
double Anorm,
double *rcond,
Options const& opts);

template
void pocondest< std::complex<float> >(
float pocondest< std::complex<float> >(
Norm in_norm,
HermitianMatrix< std::complex<float> >& A,
float Anorm,
float *rcond,
Options const& opts);

template
void pocondest< std::complex<double> >(
double pocondest< std::complex<double> >(
Norm in_norm,
HermitianMatrix< std::complex<double> >& A,
double Anorm,
double *rcond,
Options const& opts);

} // namespace slate
Loading

0 comments on commit 17b05b6

Please sign in to comment.