+++ /dev/null
-// Copyright ©2015 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 "fmt"
-
-// Product calculates the product of the given factors and places the result in
-// the receiver. The order of multiplication operations is optimized to minimize
-// the number of floating point operations on the basis that all matrix
-// multiplications are general.
-func (m *Dense) Product(factors ...Matrix) {
- // The operation order optimisation is the naive O(n^3) dynamic
- // programming approach and does not take into consideration
- // finer-grained optimisations that might be available.
- //
- // TODO(kortschak) Consider using the O(nlogn) or O(mlogn)
- // algorithms that are available. e.g.
- //
- // e.g. http://www.jofcis.com/publishedpapers/2014_10_10_4299_4306.pdf
- //
- // In the case that this is replaced, retain this code in
- // tests to compare against.
-
- r, c := m.Dims()
- switch len(factors) {
- case 0:
- if r != 0 || c != 0 {
- panic(ErrShape)
- }
- return
- case 1:
- m.reuseAs(factors[0].Dims())
- m.Copy(factors[0])
- return
- case 2:
- // Don't do work that we know the answer to.
- m.Mul(factors[0], factors[1])
- return
- }
-
- p := newMultiplier(m, factors)
- p.optimize()
- result := p.multiply()
- m.reuseAs(result.Dims())
- m.Copy(result)
- putWorkspace(result)
-}
-
-// debugProductWalk enables debugging output for Product.
-const debugProductWalk = false
-
-// multiplier performs operation order optimisation and tree traversal.
-type multiplier struct {
- // factors is the ordered set of
- // factors to multiply.
- factors []Matrix
- // dims is the chain of factor
- // dimensions.
- dims []int
-
- // table contains the dynamic
- // programming costs and subchain
- // division indices.
- table table
-}
-
-func newMultiplier(m *Dense, factors []Matrix) *multiplier {
- // Check size early, but don't yet
- // allocate data for m.
- r, c := m.Dims()
- fr, fc := factors[0].Dims() // newMultiplier is only called with len(factors) > 2.
- if !m.IsZero() {
- if fr != r {
- panic(ErrShape)
- }
- if _, lc := factors[len(factors)-1].Dims(); lc != c {
- panic(ErrShape)
- }
- }
-
- dims := make([]int, len(factors)+1)
- dims[0] = r
- dims[len(dims)-1] = c
- pc := fc
- for i, f := range factors[1:] {
- cr, cc := f.Dims()
- dims[i+1] = cr
- if pc != cr {
- panic(ErrShape)
- }
- pc = cc
- }
-
- return &multiplier{
- factors: factors,
- dims: dims,
- table: newTable(len(factors)),
- }
-}
-
-// optimize determines an optimal matrix multiply operation order.
-func (p *multiplier) optimize() {
- if debugProductWalk {
- fmt.Printf("chain dims: %v\n", p.dims)
- }
- const maxInt = int(^uint(0) >> 1)
- for f := 1; f < len(p.factors); f++ {
- for i := 0; i < len(p.factors)-f; i++ {
- j := i + f
- p.table.set(i, j, entry{cost: maxInt})
- for k := i; k < j; k++ {
- cost := p.table.at(i, k).cost + p.table.at(k+1, j).cost + p.dims[i]*p.dims[k+1]*p.dims[j+1]
- if cost < p.table.at(i, j).cost {
- p.table.set(i, j, entry{cost: cost, k: k})
- }
- }
- }
- }
-}
-
-// multiply walks the optimal operation tree found by optimize,
-// leaving the final result in the stack. It returns the
-// product, which may be copied but should be returned to
-// the workspace pool.
-func (p *multiplier) multiply() *Dense {
- result, _ := p.multiplySubchain(0, len(p.factors)-1)
- if debugProductWalk {
- r, c := result.Dims()
- fmt.Printf("\tpop result (%d×%d) cost=%d\n", r, c, p.table.at(0, len(p.factors)-1).cost)
- }
- return result.(*Dense)
-}
-
-func (p *multiplier) multiplySubchain(i, j int) (m Matrix, intermediate bool) {
- if i == j {
- return p.factors[i], false
- }
-
- a, aTmp := p.multiplySubchain(i, p.table.at(i, j).k)
- b, bTmp := p.multiplySubchain(p.table.at(i, j).k+1, j)
-
- ar, ac := a.Dims()
- br, bc := b.Dims()
- if ac != br {
- // Panic with a string since this
- // is not a user-facing panic.
- panic(ErrShape.Error())
- }
-
- if debugProductWalk {
- fmt.Printf("\tpush f[%d] (%d×%d)%s * f[%d] (%d×%d)%s\n",
- i, ar, ac, result(aTmp), j, br, bc, result(bTmp))
- }
-
- r := getWorkspace(ar, bc, false)
- r.Mul(a, b)
- if aTmp {
- putWorkspace(a.(*Dense))
- }
- if bTmp {
- putWorkspace(b.(*Dense))
- }
- return r, true
-}
-
-type entry struct {
- k int // is the chain subdivision index.
- cost int // cost is the cost of the operation.
-}
-
-// table is a row major n×n dynamic programming table.
-type table struct {
- n int
- entries []entry
-}
-
-func newTable(n int) table {
- return table{n: n, entries: make([]entry, n*n)}
-}
-
-func (t table) at(i, j int) entry { return t.entries[i*t.n+j] }
-func (t table) set(i, j int, e entry) { t.entries[i*t.n+j] = e }
-
-type result bool
-
-func (r result) String() string {
- if r {
- return " (popped result)"
- }
- return ""
-}
OSDN Git Service
