1 // Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
2 // other Axom Project Developers. See the top-level LICENSE file for details.
3 //
4 // SPDX-License-Identifier: (BSD-3-Clause)
5 
6 #ifndef AXOM_MATRIX_NORMS_HPP_
7 #define AXOM_MATRIX_NORMS_HPP_
8 
9 #include "axom/core/numerics/Matrix.hpp"      // for numerics::Matrix
10 #include "axom/core/utilities/Utilities.hpp"  // for utilities::abs()
11 
12 namespace axom
13 {
14 namespace numerics
15 {
16 namespace internal
17 {
18 /*!
19  * \brief Computes the p1-norm of a given \f$ M \times N \f$ matrix
20  *
21  * \param [in] A the matrix whose norm is computed
22  * \return norm the computed p1-norm of the matrix
23  *
24  * \note The p1-norm is computed by taking the maximum absolute column sum,
25  *  given by:
26  *  \f[
27  *     ||A||_1 = \max\limits_{1 \le j \le N}( \sum_{i=1}^M | a_{ij} | )
28  *  \f]
29  *
30  * \pre A.getNumRows() >= 2
31  * \pre A.getNumCols() >= 2
32  */
33 template <typename T>
matrix_p1_norm(const Matrix<T> & A)34 inline T matrix_p1_norm(const Matrix<T>& A)
35 {
36   const int numRows = A.getNumRows();
37   const int numCols = A.getNumColumns();
38 
39   if(numRows < 2 || numCols < 2)
40   {
41     /* short-circuit and return -1.0 to indicate error */
42     return -1.0;
43   }
44 
45   T norm = -1.0;
46   for(IndexType j = 0; j < numCols; ++j)
47   {
48     const T* col = A.getColumn(j);
49     T col_sum = 0.0;
50     for(IndexType i = 0; i < numRows; ++i)
51     {
52       col_sum += utilities::abs(col[i]);
53     }  // END for all rows
54 
55     norm = (col_sum > norm) ? col_sum : norm;
56   }  // END for all columns
57 
58   return norm;
59 }
60 
61 /*!
62  * \brief Computes the infinity-norm of a given \f$ M \times N \f$ matrix
63  *
64  * \param [in] A the matrix whose norm is computed
65  * \return norm the computed infinity-norm of the matrix.
66  *
67  * \note The infinity-norm is computed by taking the maximum absolute row sum,
68  *  given by:
69  *  \f[
70  *    ||A||_\infty = \max\limits_{1 \le j \le M}( \sum_{i=1}^N | a_{ij} | )
71  *  \f]
72  *
73  * \pre A.getNumRows() >= 2
74  * \pre A.getNumCols() >= 2
75  */
76 template <typename T>
matrix_infty_norm(const Matrix<T> & A)77 inline T matrix_infty_norm(const Matrix<T>& A)
78 {
79   const int numRows = A.getNumRows();
80   const int numCols = A.getNumColumns();
81 
82   if(numRows < 2 || numCols < 2)
83   {
84     /* short-circuit and return -1.0 to indicate error */
85     return -1.0;
86   }
87 
88   T norm = -1.0;
89   for(IndexType i = 0; i < numRows; ++i)
90   {
91     // The lines bracketed by rowsum_start and rowsum_end (prepended with
92     // an underscore) are used in the Sphinx documentation of Matrix.
93     // _rowsum_start
94     IndexType p = 0;
95     IndexType N = 0;
96     const T* row = A.getRow(i, p, N);
97 
98     T row_sum = 0.0;
99     for(IndexType j = 0; j < N; j += p)
100     {
101       row_sum += utilities::abs(row[j]);
102     }  // END for all columns
103     // _rowsum_end
104 
105     norm = (row_sum > norm) ? row_sum : norm;
106   }  // END for all rows
107 
108   return norm;
109 }
110 
111 /*!
112  * \brief Computes the frobenius norm of a given \f$ M \times N \f$ matrix
113  *
114  * \param [in] A the matrix whose norm is computed
115  * \return norm the computed frobenius norm of the matrix.
116  *
117  * \note The frobenius norm of an \f$ M \times N \f$ matrix, \f$ A \f$ is
118  *  defined as follows:
119  *  \f[
120  *    ||A||_F = \sqrt{ \sum_{i=1}^M \sum_{j=1}^N ({a_{ij}})^2  }
121  *  \f]
122  *
123  * \pre A.getNumRows() >= 2
124  * \pre A.getNumCols() >= 2
125  */
126 template <typename T>
matrix_frobenious_norm(const Matrix<T> & A)127 inline T matrix_frobenious_norm(const Matrix<T>& A)
128 {
129   AXOM_STATIC_ASSERT_MSG(
130     std::is_floating_point<T>::value,
131     "T is required to be a floating type for computing the frobenious norm");
132 
133   const int numRows = A.getNumRows();
134   const int numCols = A.getNumColumns();
135 
136   if(numRows < 2 || numCols < 2)
137   {
138     /* short-circuit and return -1.0 to indicate error */
139     return -1.0;
140   }
141 
142   T norm = 0.0;
143   for(IndexType i = 0; i < numRows; ++i)
144   {
145     for(IndexType j = 0; j < numCols; ++j)
146     {
147       const T abs_a_ij = utilities::abs(A(i, j));
148       norm += abs_a_ij * abs_a_ij;
149     }  // END for all columns
150   }    // END for all rows
151 
152   norm = sqrt(norm);
153   return norm;
154 }
155 
156 } /* end namespace internal */
157 } /* end namespace numerics */
158 } /* end namespace axom */
159 
160 #endif /* AXOM_MATRIX_NORMS_HPP_ */
161