Conjugate transpose: Difference between revisions
Content added Content deleted
SqrtNegInf (talk | contribs) m (→{{header|Sidef}}: Fix link: Perl 6 --> Raku) |
(Added C++ solution) |
||
Line 318: | Line 318: | ||
Complex Matrix A is not normal</pre> |
Complex Matrix A is not normal</pre> |
||
=={{header|C++}}== |
|||
<lang cpp>#include <cassert> |
|||
#include <cmath> |
|||
#include <complex> |
|||
#include <iomanip> |
|||
#include <iostream> |
|||
#include <vector> |
|||
template <typename scalar_type> class complex_matrix { |
|||
public: |
|||
using element_type = std::complex<scalar_type>; |
|||
complex_matrix(size_t rows, size_t columns) |
|||
: rows_(rows), columns_(columns) { |
|||
elements_.resize(rows * columns); |
|||
} |
|||
complex_matrix(size_t rows, size_t columns, element_type value) |
|||
: rows_(rows), columns_(columns), elements_(rows * columns, value) {} |
|||
size_t rows() const { return rows_; } |
|||
size_t columns() const { return columns_; } |
|||
element_type* data() { return &elements_[0]; } |
|||
const element_type* data() const { return &elements_[0]; } |
|||
element_type* row_data(size_t row) { |
|||
assert(row < rows_); |
|||
return &elements_[row * columns_]; |
|||
} |
|||
const element_type* row_data(size_t row) const { |
|||
assert(row < rows_); |
|||
return &elements_[row * columns_]; |
|||
} |
|||
const element_type& at(size_t row, size_t column) const { |
|||
assert(row < rows_); |
|||
assert(column < columns_); |
|||
return elements_[index(row, column)]; |
|||
} |
|||
element_type& at(size_t row, size_t column) { |
|||
assert(row < rows_); |
|||
assert(column < columns_); |
|||
return elements_[index(row, column)]; |
|||
} |
|||
friend bool operator==(const complex_matrix& a, const complex_matrix& b) { |
|||
return a.rows_ == b.rows_ && a.columns_ == b.columns_ && |
|||
a.elements_ == b.elements_; |
|||
} |
|||
private: |
|||
size_t index(size_t row, size_t column) const { |
|||
return row * columns_ + column; |
|||
} |
|||
size_t rows_; |
|||
size_t columns_; |
|||
std::vector<element_type> elements_; |
|||
}; |
|||
template <typename scalar_type> |
|||
complex_matrix<scalar_type> product(const complex_matrix<scalar_type>& a, |
|||
const complex_matrix<scalar_type>& b) { |
|||
assert(a.columns() == b.rows()); |
|||
size_t arows = a.rows(); |
|||
size_t bcolumns = b.columns(); |
|||
size_t n = a.columns(); |
|||
complex_matrix<scalar_type> c(arows, bcolumns); |
|||
for (size_t i = 0; i < arows; ++i) { |
|||
const auto* aptr = a.row_data(i); |
|||
auto* cptr = c.row_data(i); |
|||
for (size_t j = 0; j < n; ++j) { |
|||
const auto* bptr = b.row_data(j); |
|||
for (size_t k = 0; k < bcolumns; ++k) |
|||
cptr[k] += aptr[j] * bptr[k]; |
|||
} |
|||
} |
|||
return c; |
|||
} |
|||
template <typename scalar_type> |
|||
complex_matrix<scalar_type> |
|||
conjugate_transpose(const complex_matrix<scalar_type>& a) { |
|||
size_t rows = a.rows(), columns = a.columns(); |
|||
complex_matrix<scalar_type> b(columns, rows); |
|||
for (size_t i = 0; i < columns; i++) { |
|||
auto* bptr = b.row_data(i); |
|||
for (size_t j = 0; j < rows; j++) { |
|||
*bptr++ = std::conj(a.at(j, i)); |
|||
} |
|||
} |
|||
return b; |
|||
} |
|||
template <typename scalar_type> |
|||
void print(std::ostream& out, const complex_matrix<scalar_type>& a) { |
|||
out << '['; |
|||
size_t rows = a.rows(), columns = a.columns(); |
|||
for (size_t row = 0; row < rows; ++row) { |
|||
if (row > 0) |
|||
out << ", "; |
|||
out << '['; |
|||
for (size_t column = 0; column < columns; ++column) { |
|||
if (column > 0) |
|||
out << ", "; |
|||
out << a.at(row, column); |
|||
} |
|||
out << ']'; |
|||
} |
|||
out << "]\n"; |
|||
} |
|||
template <typename scalar_type> |
|||
bool is_hermitian_matrix(const complex_matrix<scalar_type>& matrix) { |
|||
if (matrix.rows() != matrix.columns()) |
|||
return false; |
|||
return matrix == conjugate_transpose(matrix); |
|||
} |
|||
template <typename scalar_type> |
|||
bool is_normal_matrix(const complex_matrix<scalar_type>& matrix) { |
|||
if (matrix.rows() != matrix.columns()) |
|||
return false; |
|||
auto c = conjugate_transpose(matrix); |
|||
return product(c, matrix) == product(matrix, c); |
|||
} |
|||
bool is_equal(const std::complex<double>& a, double b) { |
|||
constexpr double e = 1e-15; |
|||
return std::abs(a.imag()) < e && std::abs(a.real() - b) < e; |
|||
} |
|||
template <typename scalar_type> |
|||
bool is_identity_matrix(const complex_matrix<scalar_type>& matrix) { |
|||
if (matrix.rows() != matrix.columns()) |
|||
return false; |
|||
size_t rows = matrix.rows(); |
|||
for (size_t i = 0; i < rows; ++i) { |
|||
for (size_t j = 0; j < rows; ++j) { |
|||
if (!is_equal(matrix.at(i, j), scalar_type(i == j ? 1 : 0))) |
|||
return false; |
|||
} |
|||
} |
|||
return true; |
|||
} |
|||
template <typename scalar_type> |
|||
bool is_unitary_matrix(const complex_matrix<scalar_type>& matrix) { |
|||
if (matrix.rows() != matrix.columns()) |
|||
return false; |
|||
auto c = conjugate_transpose(matrix); |
|||
auto p = product(c, matrix); |
|||
return is_identity_matrix(p) && p == product(matrix, c); |
|||
} |
|||
template <typename scalar_type> |
|||
void test(const complex_matrix<scalar_type>& matrix) { |
|||
std::cout << "Matrix: "; |
|||
print(std::cout, matrix); |
|||
std::cout << "Conjugate transpose: "; |
|||
print(std::cout, conjugate_transpose(matrix)); |
|||
std::cout << std::boolalpha; |
|||
std::cout << "Hermitian: " << is_hermitian_matrix(matrix) << '\n'; |
|||
std::cout << "Normal: " << is_normal_matrix(matrix) << '\n'; |
|||
std::cout << "Unitary: " << is_unitary_matrix(matrix) << '\n'; |
|||
} |
|||
int main() { |
|||
complex_matrix<double> matrix1(3, 3); |
|||
complex_matrix<double> matrix2(3, 3); |
|||
complex_matrix<double> matrix3(3, 3); |
|||
matrix1.at(0, 0) = {2, 0}; |
|||
matrix1.at(0, 1) = {2, 1}; |
|||
matrix1.at(0, 2) = {4, 0}; |
|||
matrix1.at(1, 0) = {2, -1}; |
|||
matrix1.at(1, 1) = {3, 0}; |
|||
matrix1.at(1, 2) = {0, 1}; |
|||
matrix1.at(2, 0) = {4, 0}; |
|||
matrix1.at(2, 1) = {0, -1}; |
|||
matrix1.at(2, 2) = {1, 0}; |
|||
double n = std::sqrt(0.5); |
|||
matrix2.at(0, 0) = {n, 0}; |
|||
matrix2.at(0, 1) = {n, 0}; |
|||
matrix2.at(0, 2) = {0, 0}; |
|||
matrix2.at(1, 0) = {0, -n}; |
|||
matrix2.at(1, 1) = {0, n}; |
|||
matrix2.at(1, 2) = {0, 0}; |
|||
matrix2.at(2, 0) = {0, 0}; |
|||
matrix2.at(2, 1) = {0, 0}; |
|||
matrix2.at(2, 2) = {0, 1}; |
|||
matrix3.at(0, 0) = {2, 2}; |
|||
matrix3.at(0, 1) = {3, 1}; |
|||
matrix3.at(0, 2) = {-3, 5}; |
|||
matrix3.at(1, 0) = {2, -1}; |
|||
matrix3.at(1, 1) = {4, 1}; |
|||
matrix3.at(1, 2) = {0, 0}; |
|||
matrix3.at(2, 0) = {7, -5}; |
|||
matrix3.at(2, 1) = {1, -4}; |
|||
matrix3.at(2, 2) = {1, 0}; |
|||
test(matrix1); |
|||
std::cout << '\n'; |
|||
test(matrix2); |
|||
std::cout << '\n'; |
|||
test(matrix3); |
|||
return 0; |
|||
}</lang> |
|||
{{out}} |
|||
<pre> |
|||
Matrix: [[(2,0), (2,1), (4,0)], [(2,-1), (3,0), (0,1)], [(4,0), (0,-1), (1,0)]] |
|||
Conjugate transpose: [[(2,-0), (2,1), (4,-0)], [(2,-1), (3,-0), (0,1)], [(4,-0), (0,-1), (1,-0)]] |
|||
Hermitian: true |
|||
Normal: true |
|||
Unitary: false |
|||
Matrix: [[(0.707107,0), (0.707107,0), (0,0)], [(0,-0.707107), (0,0.707107), (0,0)], [(0,0), (0,0), (0,1)]] |
|||
Conjugate transpose: [[(0.707107,-0), (0,0.707107), (0,-0)], [(0.707107,-0), (0,-0.707107), (0,-0)], [(0,-0), (0,-0), (0,-1)]] |
|||
Hermitian: false |
|||
Normal: true |
|||
Unitary: true |
|||
Matrix: [[(2,2), (3,1), (-3,5)], [(2,-1), (4,1), (0,0)], [(7,-5), (1,-4), (1,0)]] |
|||
Conjugate transpose: [[(2,-2), (2,1), (7,5)], [(3,-1), (4,-1), (1,4)], [(-3,-5), (0,-0), (1,-0)]] |
|||
Hermitian: false |
|||
Normal: false |
|||
Unitary: false |
|||
</pre> |
|||
=={{header|Common Lisp}}== |
=={{header|Common Lisp}}== |