cmatrix.go 5.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211
  1. // Copyright ©2013 The Gonum Authors. All rights reserved.
  2. // Use of this source code is governed by a BSD-style
  3. // license that can be found in the LICENSE file.
  4. package mat
  5. import (
  6. "math"
  7. "math/cmplx"
  8. "gonum.org/v1/gonum/floats"
  9. )
  10. // CMatrix is the basic matrix interface type for complex matrices.
  11. type CMatrix interface {
  12. // Dims returns the dimensions of a Matrix.
  13. Dims() (r, c int)
  14. // At returns the value of a matrix element at row i, column j.
  15. // It will panic if i or j are out of bounds for the matrix.
  16. At(i, j int) complex128
  17. // H returns the conjugate transpose of the Matrix. Whether H
  18. // returns a copy of the underlying data is implementation dependent.
  19. // This method may be implemented using the Conjugate type, which
  20. // provides an implicit matrix conjugate transpose.
  21. H() CMatrix
  22. }
  23. var (
  24. _ CMatrix = Conjugate{}
  25. _ Unconjugator = Conjugate{}
  26. )
  27. // Conjugate is a type for performing an implicit matrix conjugate transpose.
  28. // It implements the Matrix interface, returning values from the conjugate
  29. // transpose of the matrix within.
  30. type Conjugate struct {
  31. CMatrix CMatrix
  32. }
  33. // At returns the value of the element at row i and column j of the conjugate
  34. // transposed matrix, that is, row j and column i of the Matrix field.
  35. func (t Conjugate) At(i, j int) complex128 {
  36. z := t.CMatrix.At(j, i)
  37. return cmplx.Conj(z)
  38. }
  39. // Dims returns the dimensions of the transposed matrix. The number of rows returned
  40. // is the number of columns in the Matrix field, and the number of columns is
  41. // the number of rows in the Matrix field.
  42. func (t Conjugate) Dims() (r, c int) {
  43. c, r = t.CMatrix.Dims()
  44. return r, c
  45. }
  46. // H performs an implicit conjugate transpose by returning the Matrix field.
  47. func (t Conjugate) H() CMatrix {
  48. return t.CMatrix
  49. }
  50. // Unconjugate returns the Matrix field.
  51. func (t Conjugate) Unconjugate() CMatrix {
  52. return t.CMatrix
  53. }
  54. // Unconjugator is a type that can undo an implicit conjugate transpose.
  55. type Unconjugator interface {
  56. // Note: This interface is needed to unify all of the Conjugate types. In
  57. // the cmat128 methods, we need to test if the Matrix has been implicitly
  58. // transposed. If this is checked by testing for the specific Conjugate type
  59. // then the behavior will be different if the user uses H() or HTri() for a
  60. // triangular matrix.
  61. // Unconjugate returns the underlying Matrix stored for the implicit
  62. // conjugate transpose.
  63. Unconjugate() CMatrix
  64. }
  65. // useC returns a complex128 slice with l elements, using c if it
  66. // has the necessary capacity, otherwise creating a new slice.
  67. func useC(c []complex128, l int) []complex128 {
  68. if l <= cap(c) {
  69. return c[:l]
  70. }
  71. return make([]complex128, l)
  72. }
  73. // useZeroedC returns a complex128 slice with l elements, using c if it
  74. // has the necessary capacity, otherwise creating a new slice. The
  75. // elements of the returned slice are guaranteed to be zero.
  76. func useZeroedC(c []complex128, l int) []complex128 {
  77. if l <= cap(c) {
  78. c = c[:l]
  79. zeroC(c)
  80. return c
  81. }
  82. return make([]complex128, l)
  83. }
  84. // zeroC zeros the given slice's elements.
  85. func zeroC(c []complex128) {
  86. for i := range c {
  87. c[i] = 0
  88. }
  89. }
  90. // unconjugate unconjugates a matrix if applicable. If a is an Unconjugator, then
  91. // unconjugate returns the underlying matrix and true. If it is not, then it returns
  92. // the input matrix and false.
  93. func unconjugate(a CMatrix) (CMatrix, bool) {
  94. if ut, ok := a.(Unconjugator); ok {
  95. return ut.Unconjugate(), true
  96. }
  97. return a, false
  98. }
  99. // CEqual returns whether the matrices a and b have the same size
  100. // and are element-wise equal.
  101. func CEqual(a, b CMatrix) bool {
  102. ar, ac := a.Dims()
  103. br, bc := b.Dims()
  104. if ar != br || ac != bc {
  105. return false
  106. }
  107. // TODO(btracey): Add in fast-paths.
  108. for i := 0; i < ar; i++ {
  109. for j := 0; j < ac; j++ {
  110. if a.At(i, j) != b.At(i, j) {
  111. return false
  112. }
  113. }
  114. }
  115. return true
  116. }
  117. // CEqualApprox returns whether the matrices a and b have the same size and contain all equal
  118. // elements with tolerance for element-wise equality specified by epsilon. Matrices
  119. // with non-equal shapes are not equal.
  120. func CEqualApprox(a, b CMatrix, epsilon float64) bool {
  121. // TODO(btracey):
  122. ar, ac := a.Dims()
  123. br, bc := b.Dims()
  124. if ar != br || ac != bc {
  125. return false
  126. }
  127. for i := 0; i < ar; i++ {
  128. for j := 0; j < ac; j++ {
  129. if !cEqualWithinAbsOrRel(a.At(i, j), b.At(i, j), epsilon, epsilon) {
  130. return false
  131. }
  132. }
  133. }
  134. return true
  135. }
  136. // TODO(btracey): Move these into a cmplxs if/when we have one.
  137. func cEqualWithinAbsOrRel(a, b complex128, absTol, relTol float64) bool {
  138. if cEqualWithinAbs(a, b, absTol) {
  139. return true
  140. }
  141. return cEqualWithinRel(a, b, relTol)
  142. }
  143. // cEqualWithinAbs returns true if a and b have an absolute
  144. // difference of less than tol.
  145. func cEqualWithinAbs(a, b complex128, tol float64) bool {
  146. return a == b || cmplx.Abs(a-b) <= tol
  147. }
  148. const minNormalFloat64 = 2.2250738585072014e-308
  149. // cEqualWithinRel returns true if the difference between a and b
  150. // is not greater than tol times the greater value.
  151. func cEqualWithinRel(a, b complex128, tol float64) bool {
  152. if a == b {
  153. return true
  154. }
  155. if cmplx.IsNaN(a) || cmplx.IsNaN(b) {
  156. return false
  157. }
  158. // Cannot play the same trick as in floats because there are multiple
  159. // possible infinities.
  160. if cmplx.IsInf(a) {
  161. if !cmplx.IsInf(b) {
  162. return false
  163. }
  164. ra := real(a)
  165. if math.IsInf(ra, 0) {
  166. if ra == real(b) {
  167. return floats.EqualWithinRel(imag(a), imag(b), tol)
  168. }
  169. return false
  170. }
  171. if imag(a) == imag(b) {
  172. return floats.EqualWithinRel(ra, real(b), tol)
  173. }
  174. return false
  175. }
  176. if cmplx.IsInf(b) {
  177. return false
  178. }
  179. delta := cmplx.Abs(a - b)
  180. if delta <= minNormalFloat64 {
  181. return delta <= tol*minNormalFloat64
  182. }
  183. return delta/math.Max(cmplx.Abs(a), cmplx.Abs(b)) <= tol
  184. }