Example 04: Matrix Norms

This example demonstrates how to compute various matrix norms.

Key Concepts

  1. Norm Types: * Norm::One: Maximum column sum. * Norm::Inf: Maximum row sum. * Norm::Max: Maximum absolute element. * Norm::Fro: Frobenius norm (square root of sum of squares).

  2. Polymorphism: The slate::norm function works for all matrix types (General, Symmetric, Hermitian, Triangular, Trapezoid).

C++ Example

General Matrix Norms (Lines 27-41)

slate::Matrix<scalar_type> A( m, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
// ... insert tiles and fill data ...

real_type A_norm_one = slate::norm( slate::Norm::One, A );
real_type A_norm_inf = slate::norm( slate::Norm::Inf, A );
// ...

Computes norms for a standard general matrix A.

  • Norm::One (1-norm): Maximum absolute column sum. \(\max_j \sum_i |a_{ij}|\).

  • Norm::Inf (Infinity-norm): Maximum absolute row sum. \(\max_i \sum_j |a_{ij}|\).

  • Norm::Max (Max-norm): Maximum absolute value of any element. \(\max_{i,j} |a_{ij}|\).

  • Norm::Fro (Frobenius-norm): Square root of the sum of squares of all elements. \(\sqrt{\sum_{i,j} |a_{ij}|^2}\).

Symmetric Matrix Norms (Lines 59-73)

slate::SymmetricMatrix<scalar_type> S( ... );
// ...
real_type S_norm_one = slate::norm( slate::Norm::One, S );

slate::norm automatically handles symmetric matrices. It understands that elements in the unreferenced triangle are equal to their symmetric counterparts. For a symmetric matrix, the 1-norm and Infinity-norm are equal.

Hermitian Matrix Norms (Lines 88-96)

Similar to symmetric matrices, but for Hermitian matrices (where symmetric elements are complex conjugates). The logic for norms accounts for the magnitudes properly.

Triangular and Trapezoid Norms (Lines 111-143)

slate::TriangularMatrix<scalar_type> T( ... );
real_type T_norm_one = slate::norm( slate::Norm::One, T );

For triangular (and trapezoid) matrices, slate::norm respects the structure. Elements in the empty part of the matrix (e.g., upper triangle for a Lower triangular matrix) are treated as zero and do not contribute to the norm.

  1// ex04_norm.cc
  2// matrix norms
  3
  4/// !!!   Lines between `//---------- begin label`          !!!
  5/// !!!             and `//---------- end label`            !!!
  6/// !!!   are included in the SLATE Users' Guide.           !!!
  7
  8#include <slate/slate.hh>
  9
 10#include "util.hh"
 11
 12int mpi_size = 0;
 13int mpi_rank = 0;
 14int grid_p = 0;
 15int grid_q = 0;
 16
 17//------------------------------------------------------------------------------
 18template <typename scalar_type>
 19void test_general_norm()
 20{
 21    using real_type = blas::real_type<scalar_type>;
 22
 23    print_func( mpi_rank );
 24
 25    int64_t m=2000, n=1000, nb=256;
 26
 27    //---------- begin genorm1
 28    slate::Matrix<scalar_type> A( m, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
 29    // ...
 30
 31    //---------- end genorm1
 32
 33    A.insertLocalTiles();
 34    random_matrix( A );
 35
 36    //---------- begin genorm2
 37    real_type A_norm_one = slate::norm( slate::Norm::One, A );
 38    real_type A_norm_inf = slate::norm( slate::Norm::Inf, A );
 39    real_type A_norm_max = slate::norm( slate::Norm::Max, A );
 40    real_type A_norm_fro = slate::norm( slate::Norm::Fro, A );
 41    //---------- end genorm2
 42    printf( "rank %d: norms: one %12.6f, inf %12.6f, max %12.6f, fro %12.6f\n",
 43            mpi_rank, A_norm_one, A_norm_inf, A_norm_max, A_norm_fro );
 44}
 45
 46//------------------------------------------------------------------------------
 47template <typename scalar_type>
 48void test_symmetric_norm()
 49{
 50    using real_type = blas::real_type<scalar_type>;
 51
 52    print_func( mpi_rank );
 53
 54    int64_t n=1000, nb=256;
 55
 56    //---------- begin synorm1
 57
 58    // norm() is overloaded for all matrix types: Symmetric, Triangular, etc.
 59    slate::SymmetricMatrix<scalar_type>
 60        S( slate::Uplo::Lower, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
 61    // ...
 62
 63    //---------- end synorm1
 64
 65    S.insertLocalTiles();
 66    random_matrix( S );
 67
 68    //---------- begin synorm2
 69    real_type S_norm_one = slate::norm( slate::Norm::One, S );
 70    real_type S_norm_inf = slate::norm( slate::Norm::Inf, S );
 71    real_type S_norm_max = slate::norm( slate::Norm::Max, S );
 72    real_type S_norm_fro = slate::norm( slate::Norm::Fro, S );
 73    //---------- end synorm2
 74    printf( "rank %d: norms: one %12.6f, inf %12.6f, max %12.6f, fro %12.6f\n",
 75            mpi_rank, S_norm_one, S_norm_inf, S_norm_max, S_norm_fro );
 76}
 77
 78//------------------------------------------------------------------------------
 79template <typename scalar_type>
 80void test_hermitian_norm()
 81{
 82    using real_type = blas::real_type<scalar_type>;
 83
 84    print_func( mpi_rank );
 85
 86    int64_t n=1000, nb=256;
 87
 88    slate::HermitianMatrix<scalar_type>
 89        H( slate::Uplo::Lower, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
 90    H.insertLocalTiles();
 91    random_matrix( H );
 92
 93    real_type H_norm_one = slate::norm( slate::Norm::One, H );
 94    real_type H_norm_inf = slate::norm( slate::Norm::Inf, H );
 95    real_type H_norm_max = slate::norm( slate::Norm::Max, H );
 96    real_type H_norm_fro = slate::norm( slate::Norm::Fro, H );
 97    printf( "rank %d: norms: one %12.6f, inf %12.6f, max %12.6f, fro %12.6f\n",
 98            mpi_rank, H_norm_one, H_norm_inf, H_norm_max, H_norm_fro );
 99}
100
101//------------------------------------------------------------------------------
102template <typename scalar_type>
103void test_triangular_norm()
104{
105    using real_type = blas::real_type<scalar_type>;
106
107    print_func( mpi_rank );
108
109    int64_t n=1000, nb=256;
110
111    slate::TriangularMatrix<scalar_type>
112        T( slate::Uplo::Lower, slate::Diag::NonUnit, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
113    T.insertLocalTiles();
114    random_matrix( T );
115
116    real_type T_norm_one = slate::norm( slate::Norm::One, T );
117    real_type T_norm_inf = slate::norm( slate::Norm::Inf, T );
118    real_type T_norm_max = slate::norm( slate::Norm::Max, T );
119    real_type T_norm_fro = slate::norm( slate::Norm::Fro, T );
120    printf( "rank %d: norms: one %12.6f, inf %12.6f, max %12.6f, fro %12.6f\n",
121            mpi_rank, T_norm_one, T_norm_inf, T_norm_max, T_norm_fro );
122}
123
124//------------------------------------------------------------------------------
125template <typename scalar_type>
126void test_trapezoid_norm()
127{
128    using real_type = blas::real_type<scalar_type>;
129
130    print_func( mpi_rank );
131
132    int64_t m=2000, n=1000, nb=256;
133
134    slate::TrapezoidMatrix<scalar_type>
135        Tz( slate::Uplo::Lower, slate::Diag::NonUnit, m, n, nb,
136            grid_p, grid_q, MPI_COMM_WORLD );
137    Tz.insertLocalTiles();
138    random_matrix( Tz );
139
140    real_type Tz_norm_one = slate::norm( slate::Norm::One, Tz );
141    real_type Tz_norm_inf = slate::norm( slate::Norm::Inf, Tz );
142    real_type Tz_norm_max = slate::norm( slate::Norm::Max, Tz );
143    real_type Tz_norm_fro = slate::norm( slate::Norm::Fro, Tz );
144    printf( "rank %d: norms: one %12.6f, inf %12.6f, max %12.6f, fro %12.6f\n",
145            mpi_rank, Tz_norm_one, Tz_norm_inf, Tz_norm_max, Tz_norm_fro );
146}
147
148//------------------------------------------------------------------------------
149int main( int argc, char** argv )
150{
151    try {
152        // Parse command line to set types for s, d, c, z precisions.
153        bool types[ 4 ];
154        parse_args( argc, argv, types );
155
156        int provided = 0;
157        slate_mpi_call(
158            MPI_Init_thread( &argc, &argv, MPI_THREAD_MULTIPLE, &provided ) );
159        assert( provided == MPI_THREAD_MULTIPLE );
160
161        slate_mpi_call(
162            MPI_Comm_size( MPI_COMM_WORLD, &mpi_size ) );
163
164        slate_mpi_call(
165            MPI_Comm_rank( MPI_COMM_WORLD, &mpi_rank ) );
166
167        // Determine p-by-q grid for this MPI size.
168        grid_size( mpi_size, &grid_p, &grid_q );
169        if (mpi_rank == 0) {
170            printf( "mpi_size %d, grid_p %d, grid_q %d\n",
171                    mpi_size, grid_p, grid_q );
172        }
173
174        // so random_matrix is different on different ranks.
175        srand( 100 * mpi_rank );
176
177        if (types[ 0 ]) {
178            test_general_norm   < float >();
179            test_symmetric_norm < float >();
180            test_hermitian_norm < float >();
181            test_triangular_norm< float >();
182            test_trapezoid_norm < float >();
183        }
184        if (mpi_rank == 0)
185            printf( "\n" );
186
187        if (types[ 1 ]) {
188            test_general_norm   < double >();
189            test_symmetric_norm < double >();
190            test_hermitian_norm < double >();
191            test_triangular_norm< double >();
192            test_trapezoid_norm < double >();
193        }
194        if (mpi_rank == 0)
195            printf( "\n" );
196
197        if (types[ 2 ]) {
198            test_general_norm   < std::complex<float> >();
199            test_symmetric_norm < std::complex<float> >();
200            test_hermitian_norm < std::complex<float> >();
201            test_triangular_norm< std::complex<float> >();
202            test_trapezoid_norm < std::complex<float> >();
203        }
204        if (mpi_rank == 0)
205            printf( "\n" );
206
207        if (types[ 3 ]) {
208            test_general_norm   < std::complex<double> >();
209            test_symmetric_norm < std::complex<double> >();
210            test_hermitian_norm < std::complex<double> >();
211            test_triangular_norm< std::complex<double> >();
212            test_trapezoid_norm < std::complex<double> >();
213        }
214
215        slate_mpi_call(
216            MPI_Finalize() );
217    }
218    catch (std::exception const& ex) {
219        fprintf( stderr, "%s", ex.what() );
220        return 1;
221    }
222    return 0;
223}