Skip to content

Commit

Permalink
examples: add cond
Browse files Browse the repository at this point in the history
  • Loading branch information
@mgates3
mgates3 committed Nov 15, 2023
1 parent 17b05b6 commit 553c04a
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 3 deletions.
54 changes: 52 additions & 2 deletions examples/ex06_linear_system_lu.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,9 @@ void test_lu_mixed()
slate::gesv_mixed_gmres( A, pivots, B1, X1, iters ); // only one RHS
//---------- end mixed2

printf( "rank %d: iters %d\n", mpi_rank, iters );
if (mpi_rank == 0) {
printf( "rank %d: iters %d\n", mpi_rank, iters );
}
}

//------------------------------------------------------------------------------
Expand All @@ -94,21 +96,27 @@ void test_lu_factor()

int64_t n=1000, nrhs=100, nb=256;

//---------- begin factor1
slate::Matrix<scalar_type> A( n, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
slate::Matrix<scalar_type> B( n, nrhs, nb, grid_p, grid_q, MPI_COMM_WORLD );
slate::Pivots pivots;
// ...
//---------- end factor1

A.insertLocalTiles();
B.insertLocalTiles();
random_matrix( A );
random_matrix( B );
slate::Pivots pivots;

//---------- begin factor2
// simplified API
slate::lu_factor( A, pivots );
slate::lu_solve_using_factor( A, pivots, B );

// traditional API
slate::getrf( A, pivots ); // factor
slate::getrs( A, pivots, B ); // solve
//---------- end factor2
}

//------------------------------------------------------------------------------
Expand Down Expand Up @@ -140,6 +148,44 @@ void test_lu_inverse()
//---------- end inverse2
}

//------------------------------------------------------------------------------
template <typename scalar_type>
void test_lu_cond()
{
using real_t = blas::real_type<scalar_type>;

print_func( mpi_rank );

int64_t n=1000, nrhs=100, nb=256;

//---------- begin cond1
slate::Matrix<scalar_type> A( n, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
slate::Pivots pivots;
// ...
//---------- end cond1

A.insertLocalTiles();
random_matrix( A );

//---------- begin cond2

// Compute A_norm before factoring.
real_t A_norm = slate::norm( slate::Norm::One, A );

// Factor using lu_factor or lu_solve.
slate::lu_factor( A, pivots );

// reciprocal condition number, 1 / (||A|| * ||A^{-1}||)
real_t A_rcond = slate::lu_rcondest_using_factor( slate::Norm::One, A, A_norm );
real_t A_cond = 1. / A_rcond;
//---------- end cond2

if (mpi_rank == 0) {
printf( "rank %d: norm %.2e, rcond %.2e, cond %.2e\n",
mpi_rank, A_norm, A_rcond, 1 / A_rcond );
}
}

//------------------------------------------------------------------------------
int main( int argc, char** argv )
{
Expand Down Expand Up @@ -173,6 +219,7 @@ int main( int argc, char** argv )
test_lu< float >();
test_lu_factor< float >();
test_lu_inverse< float >();
test_lu_cond< float >();
}
if (mpi_rank == 0)
printf( "\n" );
Expand All @@ -182,6 +229,7 @@ int main( int argc, char** argv )
test_lu_factor< double >();
test_lu_inverse< double >();
test_lu_mixed< double >();
test_lu_cond< double >();
}
if (mpi_rank == 0)
printf( "\n" );
Expand All @@ -190,6 +238,7 @@ int main( int argc, char** argv )
test_lu< std::complex<float> >();
test_lu_factor< std::complex<float> >();
test_lu_inverse< std::complex<float> >();
test_lu_cond< std::complex<float> >();
}
if (mpi_rank == 0)
printf( "\n" );
Expand All @@ -199,6 +248,7 @@ int main( int argc, char** argv )
test_lu_factor< std::complex<double> >();
test_lu_inverse< std::complex<double> >();
test_lu_mixed< std::complex<double> >();
test_lu_cond< std::complex<double> >();
}

slate_mpi_call(
Expand Down
46 changes: 45 additions & 1 deletion examples/ex07_linear_system_cholesky.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,9 @@ void test_cholesky_mixed()
slate::posv_mixed_gmres( A, B1, X1, iters ); // only one RHS
//---------- end mixed2

printf( "rank %d: iters %d\n", mpi_rank, iters );
if (mpi_rank == 0) {
printf( "rank %d: iters %d\n", mpi_rank, iters );
}
}

//------------------------------------------------------------------------------
Expand Down Expand Up @@ -138,6 +140,44 @@ void test_cholesky_inverse()
//---------- end inverse2
}

//------------------------------------------------------------------------------
template <typename scalar_type>
void test_cholesky_cond()
{
using real_t = blas::real_type<scalar_type>;

print_func( mpi_rank );

int64_t n=1000, nrhs=100, nb=256;

//---------- begin cond1
slate::HermitianMatrix<scalar_type>
A( slate::Uplo::Lower, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
// ...
//---------- end cond1

A.insertLocalTiles();
random_matrix_diag_dominant( A );

//---------- begin cond2

// Compute A_norm before factoring.
real_t A_norm = slate::norm( slate::Norm::One, A );

// Factor using chol_factor or chol_solve.
slate::chol_factor( A );

// reciprocal condition number, 1 / (||A|| * ||A^{-1}||)
real_t A_rcond = slate::chol_rcondest_using_factor( slate::Norm::One, A, A_norm );
real_t A_cond = 1. / A_rcond;
//---------- end cond2

if (mpi_rank == 0) {
printf( "rank %d: norm %.2e, rcond %.2e, cond %.2e\n",
mpi_rank, A_norm, A_rcond, 1 / A_rcond );
}
}

//------------------------------------------------------------------------------
int main( int argc, char** argv )
{
Expand Down Expand Up @@ -171,6 +211,7 @@ int main( int argc, char** argv )
test_cholesky< float >();
test_cholesky_factor< float >();
test_cholesky_inverse< float >();
test_cholesky_cond< float >();
}
if (mpi_rank == 0)
printf( "\n" );
Expand All @@ -180,6 +221,7 @@ int main( int argc, char** argv )
test_cholesky_factor< double >();
test_cholesky_inverse< double >();
test_cholesky_mixed< double >();
test_cholesky_cond< double >();
}
if (mpi_rank == 0)
printf( "\n" );
Expand All @@ -188,6 +230,7 @@ int main( int argc, char** argv )
test_cholesky< std::complex<float> >();
test_cholesky_factor< std::complex<float> >();
test_cholesky_inverse< std::complex<float> >();
test_cholesky_cond< std::complex<float> >();
}
if (mpi_rank == 0)
printf( "\n" );
Expand All @@ -197,6 +240,7 @@ int main( int argc, char** argv )
test_cholesky_factor< std::complex<double> >();
test_cholesky_inverse< std::complex<double> >();
test_cholesky_mixed< std::complex<double> >();
test_cholesky_cond< std::complex<double> >();
}

slate_mpi_call(
Expand Down

0 comments on commit 553c04a

Please sign in to comment.