From 553c04a48c7266be3e9725fa34c1fe256c1bcd73 Mon Sep 17 00:00:00 2001 From: Mark Gates Date: Tue, 14 Nov 2023 21:10:09 -0700 Subject: [PATCH] examples: add cond --- examples/ex06_linear_system_lu.cc | 54 ++++++++++++++++++++++++- examples/ex07_linear_system_cholesky.cc | 46 ++++++++++++++++++++- 2 files changed, 97 insertions(+), 3 deletions(-) diff --git a/examples/ex06_linear_system_lu.cc b/examples/ex06_linear_system_lu.cc index 6bfb3f561..b2e559921 100644 --- a/examples/ex06_linear_system_lu.cc +++ b/examples/ex06_linear_system_lu.cc @@ -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 ); + } } //------------------------------------------------------------------------------ @@ -94,14 +96,19 @@ void test_lu_factor() int64_t n=1000, nrhs=100, nb=256; + //---------- begin factor1 slate::Matrix A( n, n, nb, grid_p, grid_q, MPI_COMM_WORLD ); slate::Matrix 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 ); @@ -109,6 +116,7 @@ void test_lu_factor() // traditional API slate::getrf( A, pivots ); // factor slate::getrs( A, pivots, B ); // solve + //---------- end factor2 } //------------------------------------------------------------------------------ @@ -140,6 +148,44 @@ void test_lu_inverse() //---------- end inverse2 } +//------------------------------------------------------------------------------ +template +void test_lu_cond() +{ + using real_t = blas::real_type; + + print_func( mpi_rank ); + + int64_t n=1000, nrhs=100, nb=256; + + //---------- begin cond1 + slate::Matrix 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 ) { @@ -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" ); @@ -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" ); @@ -190,6 +238,7 @@ int main( int argc, char** argv ) test_lu< std::complex >(); test_lu_factor< std::complex >(); test_lu_inverse< std::complex >(); + test_lu_cond< std::complex >(); } if (mpi_rank == 0) printf( "\n" ); @@ -199,6 +248,7 @@ int main( int argc, char** argv ) test_lu_factor< std::complex >(); test_lu_inverse< std::complex >(); test_lu_mixed< std::complex >(); + test_lu_cond< std::complex >(); } slate_mpi_call( diff --git a/examples/ex07_linear_system_cholesky.cc b/examples/ex07_linear_system_cholesky.cc index e8bad50e8..15c5dd076 100644 --- a/examples/ex07_linear_system_cholesky.cc +++ b/examples/ex07_linear_system_cholesky.cc @@ -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 ); + } } //------------------------------------------------------------------------------ @@ -138,6 +140,44 @@ void test_cholesky_inverse() //---------- end inverse2 } +//------------------------------------------------------------------------------ +template +void test_cholesky_cond() +{ + using real_t = blas::real_type; + + print_func( mpi_rank ); + + int64_t n=1000, nrhs=100, nb=256; + + //---------- begin cond1 + slate::HermitianMatrix + 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 ) { @@ -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" ); @@ -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" ); @@ -188,6 +230,7 @@ int main( int argc, char** argv ) test_cholesky< std::complex >(); test_cholesky_factor< std::complex >(); test_cholesky_inverse< std::complex >(); + test_cholesky_cond< std::complex >(); } if (mpi_rank == 0) printf( "\n" ); @@ -197,6 +240,7 @@ int main( int argc, char** argv ) test_cholesky_factor< std::complex >(); test_cholesky_inverse< std::complex >(); test_cholesky_mixed< std::complex >(); + test_cholesky_cond< std::complex >(); } slate_mpi_call(