mat: Add CDense type and basic methods (#846)

* mat: Add CDense type and basic methods

Updates #738.

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,
+	}
+}

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)
+				}
+			}
+		}
+	}
+}

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
+// At returns the value of the element at row i and column j of the conjugate
-// matrix, that is, row j and column i of the Matrix field.
+// 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

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 {

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 {