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)34inline 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)77inline 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)127inline 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