diff --git a/mat/cdense.go b/mat/cdense.go new file mode 100644 index 00000000..a6835a45 --- /dev/null +++ b/mat/cdense.go @@ -0,0 +1,59 @@ +// Copyright ©2019 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package mat + +import "gonum.org/v1/gonum/blas/cblas128" + +// Dense is a dense matrix representation with complex data. +type CDense struct { + mat cblas128.General + + capRows, capCols int +} + +// Dims returns the number of rows and columns in the matrix. +func (m *CDense) Dims() (r, c int) { + return m.mat.Rows, m.mat.Cols +} + +// H performs an implicit conjugate transpose by returning the receiver inside a +// Conjugate. +func (m *CDense) H() CMatrix { + return Conjugate{m} +} + +// NewCDense creates a new complex Dense matrix with r rows and c columns. +// If data == nil, a new slice is allocated for the backing slice. +// If len(data) == r*c, data is used as the backing slice, and changes to the +// elements of the returned CDense will be reflected in data. +// If neither of these is true, NewCDense will panic. +// NewCDense will panic if either r or c is zero. +// +// The data must be arranged in row-major order, i.e. the (i*c + j)-th +// element in the data slice is the {i, j}-th element in the matrix. +func NewCDense(r, c int, data []complex128) *CDense { + if r <= 0 || c <= 0 { + if r == 0 || c == 0 { + panic(ErrZeroLength) + } + panic("mat: negative dimension") + } + if data != nil && r*c != len(data) { + panic(ErrShape) + } + if data == nil { + data = make([]complex128, r*c) + } + return &CDense{ + mat: cblas128.General{ + Rows: r, + Cols: c, + Stride: c, + Data: data, + }, + capRows: r, + capCols: c, + } +} diff --git a/mat/cdense_test.go b/mat/cdense_test.go new file mode 100644 index 00000000..53585064 --- /dev/null +++ b/mat/cdense_test.go @@ -0,0 +1,48 @@ +// Copyright ©2019 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package mat + +import "testing" + +func TestCDenseNewAtSet(t *testing.T) { + for cas, test := range []struct { + a []complex128 + rows, cols int + }{ + { + a: []complex128{0, 0, 0, + 0, 0, 0, + 0, 0, 0}, + rows: 3, + cols: 3, + }, + } { + aCopy := make([]complex128, len(test.a)) + copy(aCopy, test.a) + mZero := NewCDense(test.rows, test.cols, nil) + rows, cols := mZero.Dims() + if rows != test.rows { + t.Errorf("unexpected number of rows for test %d: got: %d want: %d", cas, rows, test.rows) + } + if cols != test.cols { + t.Errorf("unexpected number of cols for test %d: got: %d want: %d", cas, cols, test.cols) + } + m := NewCDense(test.rows, test.cols, aCopy) + for i := 0; i < test.rows; i++ { + for j := 0; j < test.cols; j++ { + v := m.At(i, j) + idx := i*test.rows + j + if v != test.a[idx] { + t.Errorf("unexpected get value for test %d at i=%d, j=%d: got: %v, want: %v", cas, i, j, v, test.a[idx]) + } + add := complex(float64(i+1), float64(j+1)) + m.Set(i, j, v+add) + if m.At(i, j) != test.a[idx]+add { + t.Errorf("unexpected set value for test %d at i=%d, j=%d: got: %v, want: %v", cas, i, j, v, test.a[idx]+add) + } + } + } + } +} diff --git a/mat/cmatrix.go b/mat/cmatrix.go index fa8e135b..bff38c20 100644 --- a/mat/cmatrix.go +++ b/mat/cmatrix.go @@ -4,6 +4,8 @@ package mat +import "math/cmplx" + // CMatrix is the basic matrix interface type for complex matrices. type CMatrix interface { // Dims returns the dimensions of a Matrix. @@ -32,11 +34,11 @@ type Conjugate struct { CMatrix CMatrix } -// At returns the value of the element at row i and column j of the transposed -// matrix, that is, row j and column i of the Matrix field. +// At returns the value of the element at row i and column j of the conjugate +// transposed matrix, that is, row j and column i of the Matrix field. func (t Conjugate) At(i, j int) complex128 { z := t.CMatrix.At(j, i) - return complex(real(z), -imag(z)) + return cmplx.Conj(z) } // Dims returns the dimensions of the transposed matrix. The number of rows returned diff --git a/mat/index_bound_checks.go b/mat/index_bound_checks.go index af593b37..59815a67 100644 --- a/mat/index_bound_checks.go +++ b/mat/index_bound_checks.go @@ -38,6 +38,36 @@ func (m *Dense) set(i, j int, v float64) { m.mat.Data[i*m.mat.Stride+j] = v } +// At returns the element at row i, column j. +func (m *CDense) At(i, j int) complex128 { + return m.at(i, j) +} + +func (m *CDense) at(i, j int) complex128 { + if uint(i) >= uint(m.mat.Rows) { + panic(ErrRowAccess) + } + if uint(j) >= uint(m.mat.Cols) { + panic(ErrColAccess) + } + return m.mat.Data[i*m.mat.Stride+j] +} + +// Set sets the element at row i, column j to the value v. +func (m *CDense) Set(i, j int, v complex128) { + m.set(i, j, v) +} + +func (m *CDense) set(i, j int, v complex128) { + if uint(i) >= uint(m.mat.Rows) { + panic(ErrRowAccess) + } + if uint(j) >= uint(m.mat.Cols) { + panic(ErrColAccess) + } + m.mat.Data[i*m.mat.Stride+j] = v +} + // At returns the element at row i. // It panics if i is out of bounds or if j is not zero. func (v *VecDense) At(i, j int) float64 { diff --git a/mat/index_no_bound_checks.go b/mat/index_no_bound_checks.go index b2606768..051f8437 100644 --- a/mat/index_no_bound_checks.go +++ b/mat/index_no_bound_checks.go @@ -38,6 +38,36 @@ func (m *Dense) set(i, j int, v float64) { m.mat.Data[i*m.mat.Stride+j] = v } +// At returns the element at row i, column j. +func (m *CDense) At(i, j int) complex128 { + if uint(i) >= uint(m.mat.Rows) { + panic(ErrRowAccess) + } + if uint(j) >= uint(m.mat.Cols) { + panic(ErrColAccess) + } + return m.at(i, j) +} + +func (m *CDense) at(i, j int) complex128 { + return m.mat.Data[i*m.mat.Stride+j] +} + +// Set sets the element at row i, column j to the value v. +func (m *CDense) Set(i, j int, v complex128) { + if uint(i) >= uint(m.mat.Rows) { + panic(ErrRowAccess) + } + if uint(j) >= uint(m.mat.Cols) { + panic(ErrColAccess) + } + m.set(i, j, v) +} + +func (m *CDense) set(i, j int, v complex128) { + m.mat.Data[i*m.mat.Stride+j] = v +} + // At returns the element at row i. // It panics if i is out of bounds or if j is not zero. func (v *VecDense) At(i, j int) float64 {