matrix_dense_real_template.in

/* -*- mode: go; -*-
 *
 * Copyright (C) 2015-2020 Philipp Benner
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "macros.h"

#define MATRIX_CONST_ITERATOR  STR_CONCAT(MATRIX_NAME, ConstIterator)
#define MATRIX_JOINT_ITERATOR  STR_CONCAT(MATRIX_NAME, JointIterator)
#define MATRIX_ITERATOR        STR_CONCAT(MATRIX_NAME, Iterator)
#define VECTOR_CONST_ITERATOR  STR_CONCAT(VECTOR_NAME, ConstIterator)
#define VECTOR_ITERATOR        STR_CONCAT(VECTOR_NAME, Iterator)

#ifdef STORE_PTR
# define STORED_SCALAR_TYPE *SCALAR_NAME
#else
# define STORED_SCALAR_TYPE  SCALAR_NAME
#endif

/* -------------------------------------------------------------------------- */

package autodiff

/* -------------------------------------------------------------------------- */

import "fmt"
import "bytes"
import "bufio"
import "compress/gzip"
import "encoding/json"
import "io"
import "os"
import "strconv"
import "strings"
import "unsafe"

/* matrix type declaration
 * -------------------------------------------------------------------------- */

type MATRIX_NAME struct {
  values     VECTOR_TYPE
  rows       int
  cols       int
  rowOffset  int
  rowMax     int
  colOffset  int
  colMax     int
  transposed bool
  tmp1       VECTOR_TYPE
  tmp2       VECTOR_TYPE
}

/* constructors
 * -------------------------------------------------------------------------- */

func NEW_MATRIX(values []STORED_TYPE, rows, cols int) MATRIX_TYPE {
  m := NIL_MATRIX(rows, cols)
  v := m.values
  if len(values) == 1 {
    for i := 0; i < rows*cols; i++ {
#ifdef STORE_PTR
      v[i] =  NEW_SCALAR(values[0])
#else
      v[i] = *NEW_SCALAR(values[0])
#endif
    }
  } else if len(values) == rows*cols {
    for i := 0; i < rows*cols; i++ {
#ifdef STORE_PTR
      v[i] =  NEW_SCALAR(values[i])
#else
      v[i] = *NEW_SCALAR(values[i])
#endif
    }
  } else {
    panic("NewMatrix(): Matrix dimension does not fit input values!")
  }
  m.initTmp()

  return m
}

func NULL_MATRIX(rows, cols int) MATRIX_TYPE {
  m := MATRIX_NAME{}
  m.values    = NULL_VECTOR(rows*cols)
  m.rows      = rows
  m.cols      = cols
  m.rowOffset = 0
  m.rowMax    = rows
  m.colOffset = 0
  m.colMax    = cols
  m.initTmp()
  return &m
}

func NIL_MATRIX(rows, cols int) MATRIX_TYPE {
  m := MATRIX_NAME{}
  m.values    = NIL_VECTOR(rows*cols)
  m.rows      = rows
  m.cols      = cols
  m.rowOffset = 0
  m.rowMax    = rows
  m.colOffset = 0
  m.colMax    = cols
  return &m
}

func AS_MATRIX(matrix ConstMatrix) MATRIX_TYPE {
  switch matrix_ := matrix.(type) {
  case MATRIX_TYPE:
    return matrix_.Clone()
  }
  n, m := matrix.Dims()
  r := NULL_MATRIX(n, m)
  for i := 0; i < n; i++ {
    for j := 0; j < m; j++ {
      r.AT(i,j).Set(matrix.ConstAt(i,j))
    }
  }
  return r
}

func (matrix MATRIX_TYPE) initTmp() {
  if len(matrix.tmp1) < matrix.rows {
    matrix.tmp1 = NULL_VECTOR(matrix.rows)
  } else {
    matrix.tmp1 = matrix.tmp1[0:matrix.rows]
  }
  if len(matrix.tmp2) < matrix.cols {
    matrix.tmp2 = NULL_VECTOR(matrix.cols)
  } else {
    matrix.tmp2 = matrix.tmp2[0:matrix.cols]
  }
}

/* cloning
 * -------------------------------------------------------------------------- */

// Clone matrix including data.
func (matrix MATRIX_TYPE) Clone() MATRIX_TYPE {
  return &MATRIX_NAME{
    values    : matrix.values.Clone(),
    rows      : matrix.rows,
    cols      : matrix.cols,
    transposed: matrix.transposed,
    rowOffset : matrix.rowOffset,
    rowMax    : matrix.rowMax,
    colOffset : matrix.colOffset,
    colMax    : matrix.colMax,
    tmp1      : matrix.tmp1.Clone(),
    tmp2      : matrix.tmp2.Clone() }
}

/* indexing
 * -------------------------------------------------------------------------- */

func (matrix MATRIX_TYPE) index(i, j int) int {
  if i < 0 || j < 0 || i >= matrix.rows || j >= matrix.cols {
    panic(fmt.Errorf("index (%d,%d) out of bounds for matrix of dimension %dx%d", i, j, matrix.rows, matrix.cols))
  }
  if matrix.transposed {
    return (matrix.colOffset + j)*matrix.rowMax + (matrix.rowOffset + i)
  } else {
    return (matrix.rowOffset + i)*matrix.colMax + (matrix.colOffset + j)
  }
}

func (matrix MATRIX_TYPE) ij(k int) (int, int) {
  if matrix.transposed {
    i := (k%matrix.colMax) - matrix.colOffset
    j := (k/matrix.colMax) - matrix.rowOffset
    return i, j
  } else {
    i := (k/matrix.rowMax) - matrix.rowOffset
    j := (k%matrix.rowMax) - matrix.colOffset
    return i, j
  }
}

/* native matrix methods
 * -------------------------------------------------------------------------- */

func (matrix MATRIX_TYPE) AT(i, j int) SCALAR_TYPE {
#ifdef STORE_PTR
  return matrix.values[matrix.index(i, j)]
#else
  return &matrix.values[matrix.index(i, j)]
#endif
}

func (matrix MATRIX_TYPE) ROW(i int) VECTOR_TYPE {
  v := NIL_VECTOR(matrix.cols)
  for j := 0; j < matrix.cols; j++ {
    v[j] = matrix.values[matrix.index(i, j)].Clone()
  }
  return v
}

func (matrix MATRIX_TYPE) COL(j int) VECTOR_TYPE {
  v := NIL_VECTOR(matrix.rows)
  for i := 0; i < matrix.rows; i++ {
    v[i] = matrix.values[matrix.index(i, j)].Clone()
  }
  return v
}

func (matrix MATRIX_TYPE) DIAG() VECTOR_TYPE {
  n, m := matrix.Dims()
  if n != m {
    panic("Diag(): not a square matrix!")
  }
  v := NIL_VECTOR(n)
  for i := 0; i < n; i++ {
    v[i] = matrix.values[matrix.index(i, i)].Clone()
  }
  return v
}

func (matrix MATRIX_TYPE) SLICE(rfrom, rto, cfrom, cto int) MATRIX_TYPE {
  m := *matrix
  m.rowOffset += rfrom
  m.rows       = rto - rfrom
  m.colOffset += cfrom
  m.cols       = cto - cfrom
  // crop tmp vectors
  m.initTmp()
  return &m
}

func (matrix MATRIX_TYPE) STR_CONCAT(As, VECTOR_NAME)() VECTOR_TYPE {
  if matrix.cols < matrix.colMax - matrix.colOffset ||
    (matrix.rows < matrix.rowMax - matrix.rowOffset) {
    n, m := matrix.Dims()
    v := NIL_VECTOR(n*m)
    for i := 0; i < n; i++ {
      for j := 0; j < m; j++ {
#ifdef STORE_PTR
        v[i*matrix.cols + j] =  matrix.AT(i, j)
#else
        v[i*matrix.cols + j] = *matrix.AT(i, j)
#endif
      }
    }
    return v
  } else {
    return VECTOR_TYPE(matrix.values)
  }
}

/* matrix interface
 * -------------------------------------------------------------------------- */

func (matrix MATRIX_TYPE) CloneMatrix() Matrix {
  return matrix.Clone()
}

func (matrix MATRIX_TYPE) At(i, j int) Scalar {
  return matrix.AT(i, j)
}

func (a MATRIX_TYPE) Set(b ConstMatrix) {
  n1, m1 := a.Dims()
  n2, m2 := b.Dims()
  if n1 != n2 || m1 != m2 {
    panic("Copy(): Matrix dimension does not match!")
  }
  for i := 0; i < n1; i++ {
    for j := 0; j < m1; j++ {
      a.At(i, j).Set(b.ConstAt(i, j))
    }
  }
}

func (matrix MATRIX_TYPE) SetIdentity() {
  n, m := matrix.Dims()
  c := NewScalar(matrix.ElementType(), 1.0)
  for i := 0; i < n; i++ {
    for j := 0; j < m; j++ {
      if i == j {
        matrix.At(i, j).Set(c)
      } else {
        matrix.At(i, j).Reset()
      }
    }
  }
}

func (matrix MATRIX_TYPE) Reset() {
  for i := 0; i < len(matrix.values); i++ {
    matrix.values[i].Reset()
  }
}

func (matrix MATRIX_TYPE) Row(i int) Vector {
  return matrix.ROW(i)
}

func (matrix MATRIX_TYPE) Col(j int) Vector {
  return matrix.COL(j)
}

func (matrix MATRIX_TYPE) Diag() Vector {
  return matrix.DIAG()
}

func (matrix MATRIX_TYPE) Slice(rfrom, rto, cfrom, cto int) Matrix {
  return matrix.SLICE(rfrom, rto, cfrom, cto)
}

func (matrix MATRIX_TYPE) Swap(i1, j1, i2, j2 int) {
  k1 := matrix.index(i1, j1)
  k2 := matrix.index(i2, j2)
  matrix.values[k1], matrix.values[k2] = matrix.values[k2], matrix.values[k1]
}

func (matrix MATRIX_TYPE) T() Matrix {
  return matrix.MagicT()
}

func (matrix MATRIX_TYPE) Tip() {
  mn      := len(matrix.values)
  visited := make([]bool, mn)
  k       := 0
  for cycle := 1; cycle < mn; cycle++ {
    if visited[cycle] {
      continue
    }
    k = cycle

    for {
      if k != mn-1 {
        k = matrix.rows*k % (mn-1)
      }
      visited[k] = true
      // swap
      matrix.values[k], matrix.values[cycle] = matrix.values[cycle], matrix.values[k]

      if k == cycle {
        break
      }
    }
  }
  matrix.rows,      matrix.cols      = matrix.cols,      matrix.rows
  matrix.rowOffset, matrix.colOffset = matrix.colOffset, matrix.rowOffset
  matrix.rowMax,    matrix.colMax    = matrix.colMax,    matrix.rowMax
  matrix.tmp1,      matrix.tmp2      = matrix.tmp2,      matrix.tmp1
}


func (matrix MATRIX_TYPE) AsVector() Vector {
  return matrix.STR_CONCAT(As, VECTOR_NAME)()
}

func (matrix MATRIX_TYPE) storageLocation() uintptr {
  return uintptr(unsafe.Pointer(&matrix.values[0]))
}

/* const interface
 * -------------------------------------------------------------------------- */

func (matrix MATRIX_TYPE) CloneConstMatrix() ConstMatrix {
  return matrix.Clone()
}

func (matrix MATRIX_TYPE) Dims() (int, int) {
  if matrix == nil {
    return 0, 0
  } else {
    return matrix.rows, matrix.cols
  }
}

func (matrix MATRIX_TYPE) Int8At(i, j int) int8 {
  return matrix.values[matrix.index(i, j)].GetInt8()
}

func (matrix MATRIX_TYPE) Int16At(i, j int) int16 {
  return matrix.values[matrix.index(i, j)].GetInt16()
}

func (matrix MATRIX_TYPE) Int32At(i, j int) int32 {
  return matrix.values[matrix.index(i, j)].GetInt32()
}

func (matrix MATRIX_TYPE) Int64At(i, j int) int64 {
  return matrix.values[matrix.index(i, j)].GetInt64()
}

func (matrix MATRIX_TYPE) IntAt(i, j int) int {
  return matrix.values[matrix.index(i, j)].GetInt()
}

func (matrix MATRIX_TYPE) Float32At(i, j int) float32 {
  return matrix.values[matrix.index(i, j)].GetFloat32()
}

func (matrix MATRIX_TYPE) Float64At(i, j int) float64 {
  return matrix.values[matrix.index(i, j)].GetFloat64()
}

func (matrix MATRIX_TYPE) ConstAt(i, j int) ConstScalar {
  return matrix.values[matrix.index(i, j)]
}

func (matrix MATRIX_TYPE) ConstSlice(rfrom, rto, cfrom, cto int) ConstMatrix {
  return matrix.Slice(rfrom, rto, cfrom, cto)
}

func (matrix MATRIX_TYPE) ConstRow(i int) ConstVector {
  // no cloning required...
  var v VECTOR_TYPE
  if matrix.transposed {
    v = NIL_VECTOR(matrix.cols)
    for j := 0; j < matrix.cols; j++ {
      v[j] = matrix.values[matrix.index(i, j)]
    }
  } else {
    i = matrix.index(i, 0)
    v = matrix.values[i:i + matrix.cols]
  }
  return v
}

func (matrix MATRIX_TYPE) ConstCol(j int) ConstVector {
  // no cloning required...
  var v VECTOR_TYPE
  if matrix.transposed {
    j = matrix.index(0, j)
    v = matrix.values[j:j + matrix.rows]
  } else {
    v = NIL_VECTOR(matrix.rows)
    for i := 0; i < matrix.rows; i++ {
      v[i] = matrix.values[matrix.index(i, j)]
    }
  }
  return v
}

func (matrix MATRIX_TYPE) ConstDiag() ConstVector {
  // no cloning required...
  n, m := matrix.Dims()
  if n != m {
    panic("Diag(): not a square matrix!")
  }
  v := NIL_VECTOR(n)
  for i := 0; i < n; i++ {
    v[i] = matrix.values[matrix.index(i, i)]
  }
  return v
}

func (matrix MATRIX_TYPE) IsSymmetric(epsilon float64) bool {
  n, m := matrix.Dims()
  if n != m {
    return false
  }
  for i := 0; i < n; i++ {
    for j := i+1; j < m; j++ {
      if !matrix.At(i,j).Equals(matrix.At(j,i), 1e-12) {
        return false
      }
    }
  }
  return true
}

func (matrix MATRIX_TYPE) AsConstVector() ConstVector {
  return matrix.STR_CONCAT(As, VECTOR_NAME)()
}

/* magic interface
 * -------------------------------------------------------------------------- */

func (matrix MATRIX_TYPE) CloneMagicMatrix() MagicMatrix {
  return matrix.Clone()
}

func (matrix MATRIX_TYPE) MagicAt(i, j int) MagicScalar {
  return matrix.AT(i, j)
}

func (matrix MATRIX_TYPE) MagicSlice(rfrom, rto, cfrom, cto int) MagicMatrix {
  return matrix.SLICE(rfrom, rto, cfrom, cto)
}

func (matrix MATRIX_TYPE) MagicT() MagicMatrix {
  return &MATRIX_NAME{
    values    :  matrix.values,
    rows      :  matrix.cols,
    cols      :  matrix.rows,
    transposed: !matrix.transposed,
    rowOffset :  matrix.colOffset,
    rowMax    :  matrix.colMax,
    colOffset :  matrix.rowOffset,
    colMax    :  matrix.rowMax,
    tmp1      :  matrix.tmp2,
    tmp2      :  matrix.tmp1 }
}

func (matrix MATRIX_TYPE) ResetDerivatives() {
  for i := 0; i < len(matrix.values); i++ {
    matrix.values[i].ResetDerivatives()
  }
}

func (matrix MATRIX_TYPE) AsMagicVector() MagicVector {
  return matrix.STR_CONCAT(As, VECTOR_NAME)()
}

/* implement MagicScalarContainer
 * -------------------------------------------------------------------------- */

func (matrix MATRIX_TYPE) Map(f func(Scalar)) {
  n, m := matrix.Dims()
  for i := 0; i < n; i++ {
    for j := 0; j < m; j++ {
      f(matrix.At(i, j))
    }
  }
}

func (matrix MATRIX_TYPE) MapSet(f func(ConstScalar) Scalar) {
  n, m := matrix.Dims()
  for i := 0; i < n; i++ {
    for j := 0; j < m; j++ {
      matrix.At(i,j).Set(f(matrix.ConstAt(i, j)))
    }
  }
}

func (matrix MATRIX_TYPE) Reduce(f func(Scalar, ConstScalar) Scalar, r Scalar) Scalar {
  n, m := matrix.Dims()
  for i := 0; i < n; i++ {
    for j := 0; j < m; j++ {
      r = f(r, matrix.ConstAt(i, j))
    }
  }
  return r
}

func (matrix MATRIX_TYPE) ElementType() ScalarType {
  return SCALAR_REFLECT_TYPE
}

// Treat all elements as variables for automatic differentiation. This method should only be called on a single vector or matrix. If multiple matrices should be treated as variables, then a single matrix must be allocated first and sliced after calling this method.
func (matrix MATRIX_TYPE) Variables(order int) error {
  for i, _ := range matrix.values {
    if err := matrix.values[i].SetVariable(i, len(matrix.values), order); err != nil {
      return err
    }
  }
  return nil
}

/* permutations
 * -------------------------------------------------------------------------- */

func (matrix MATRIX_TYPE) SwapRows(i, j int) error {
  n, m := matrix.Dims()
  if n != m {
    return fmt.Errorf("SymmetricPermutation(): matrix is not a square matrix")
  }
  for k := 0; k < m; k++ {
    matrix.Swap(i, k, j, k)
  }
  return nil
}

func (matrix MATRIX_TYPE) SwapColumns(i, j int) error {
  n, m := matrix.Dims()
  if n != m {
    return fmt.Errorf("SymmetricPermutation(): matrix is not a square matrix")
  }
  for k := 0; k < n; k++ {
    matrix.Swap(k, i, k, j)
  }
  return nil
}

func (matrix MATRIX_TYPE) PermuteRows(pi []int) error {
  n, m := matrix.Dims()
  if n != m {
    return fmt.Errorf("SymmetricPermutation(): matrix is not a square matrix")
  }
  // permute matrix
  for i := 0; i < n; i++ {
    if pi[i] < 0 || pi[i] > n {
      return fmt.Errorf("SymmetricPermutation(): invalid permutation")
    }
    if i != pi[i] && pi[i] > i {
      matrix.SwapRows(i, pi[i])
    }
  }
  return nil
}

func (matrix MATRIX_TYPE) PermuteColumns(pi []int) error {
  n, m := matrix.Dims()
  if n != m {
    return fmt.Errorf("SymmetricPermutation(): matrix is not a square matrix")
  }
  // permute matrix
  for i := 0; i < m; i++ {
    if pi[i] < 0 || pi[i] > n {
      return fmt.Errorf("SymmetricPermutation(): invalid permutation")
    }
    if i != pi[i] && pi[i] > i {
      matrix.SwapColumns(i, pi[i])
    }
  }
  return nil
}

func (matrix MATRIX_TYPE) SymmetricPermutation(pi []int) error {
  n, m := matrix.Dims()
  if n != m {
    return fmt.Errorf("SymmetricPermutation(): matrix is not a square matrix")
  }
  for i := 0; i < n; i++ {
    if pi[i] < 0 || pi[i] > n {
      return fmt.Errorf("SymmetricPermutation(): invalid permutation")
    }
    if pi[i] > i {
      // permute rows
      matrix.SwapRows(i, pi[i])
      // permute colums
      matrix.SwapColumns(i, pi[i])
    }
  }
  return nil
}

/* type conversion
 * -------------------------------------------------------------------------- */

func (m MATRIX_TYPE) String() string {
  var buffer bytes.Buffer

  buffer.WriteString("[")
  for i := 0; i < m.rows; i++ {
    if i != 0 {
      buffer.WriteString(",\n ")
    }
    buffer.WriteString("[")
    for j := 0; j < m.cols; j++ {
      if j != 0 {
        buffer.WriteString(", ")
      }
      buffer.WriteString(m.ConstAt(i,j).String())
    }
    buffer.WriteString("]")
  }
  buffer.WriteString("]")

  return buffer.String()
}

func (a MATRIX_TYPE) Table() string {
  var buffer bytes.Buffer

  n, m := a.Dims()

  for i := 0; i < n; i++ {
    if i != 0 {
      buffer.WriteString("\n")
    }
    for j := 0; j < m; j++ {
      if j != 0 {
        buffer.WriteString(" ")
      }
      buffer.WriteString(a.ConstAt(i,j).String())
    }
  }

  return buffer.String()
}

func (m MATRIX_TYPE) Export(filename string) error {
  f, err := os.Create(filename)
  if err != nil {
    return err
  }
  defer f.Close()

  w := bufio.NewWriter(f)
  defer w.Flush()

  if _, err := fmt.Fprintf(w, "%s\n", m.Table()); err != nil {
    return err
  }
  return nil
}

func (m MATRIX_TYPE) Import(filename string) error {
  values := []STORED_TYPE{}
  rows   := 0
  cols   := 0

  var reader *bufio.Reader
  // open file
  f, err := os.Open(filename)
  if err != nil {
    return err
  }
  defer f.Close()
  isgzip, err := isGzip(filename)
  if err != nil {
    return err
  }
  // check if file is gzipped
  if isgzip {
    g, err := gzip.NewReader(f)
    if err != nil {
      return err
    }
    defer g.Close()
    reader = bufio.NewReader(g)
  } else {
    reader = bufio.NewReader(f)
  }

  for i_ := 1;; i_++ {
    l, err := bufioReadLine(reader)
    if err == io.EOF {
      break
    }
    if err != nil {
      return err
    }
    if len(l) == 0 {
      continue
    }
    fields := strings.Fields(l)
    if cols == 0 {
      cols = len(fields)
    }
    if cols != len(fields) {
      return fmt.Errorf("invalid table")
    }
    for i := 0; i < len(fields); i++ {
      value, err := strconv.ParseFloat(fields[i], 64)
      if err != nil {
        return fmt.Errorf("invalid table")
      }
      values = append(values, STORED_TYPE(value))
    }
    rows++
  }
  *m = *NEW_MATRIX(values, rows, cols)

  return nil
}

/* json
 * -------------------------------------------------------------------------- */

func (obj MATRIX_TYPE) MarshalJSON() ([]byte, error) {
  if obj.transposed || obj.rowMax > obj.rows || obj.colMax > obj.cols {
    n, m := obj.Dims()
    tmp  := NULL_MATRIX(n, m)
    tmp.Set(obj)
    obj = tmp
  }
  r := struct{Values []STORED_SCALAR_TYPE; Rows int; Cols int}{}
  r.Values = obj.values
  r.Rows   = obj.rows
  r.Cols   = obj.cols
  return json.MarshalIndent(r, "", "  ")
}

func (obj MATRIX_TYPE) UnmarshalJSON(data []byte) error {
  r := struct{Values []STORED_SCALAR_TYPE; Rows int; Cols int}{}
  if err := json.Unmarshal(data, &r); err != nil {
    return err
  }
  obj.values = NIL_VECTOR(len(r.Values))
  for i := 0; i < len(r.Values); i++ {
    obj.values[i] = r.Values[i]
  }
  obj.rows       = r.Rows
  obj.rowMax     = r.Rows
  obj.rowOffset  = 0
  obj.cols       = r.Cols
  obj.colMax     = r.Cols
  obj.colOffset  = 0
  obj.transposed = false
  obj.initTmp()
  return nil
}


/* iterator methods
 * -------------------------------------------------------------------------- */

func (obj MATRIX_TYPE) ConstIterator() MatrixConstIterator {
  return obj.ITERATOR()
}

func (obj MATRIX_TYPE) ConstIteratorFrom(i, j int) MatrixConstIterator {
  return obj.ITERATOR_FROM(i, j)
}

func (obj MATRIX_TYPE) MagicIterator() MatrixMagicIterator {
  return obj.ITERATOR()
}

func (obj MATRIX_TYPE) MagicIteratorFrom(i, j int) MatrixMagicIterator {
  return obj.ITERATOR_FROM(i, j)
}

func (obj MATRIX_TYPE) Iterator() MatrixIterator {
  return obj.ITERATOR()
}

func (obj MATRIX_TYPE) IteratorFrom(i, j int) MatrixIterator {
  return obj.ITERATOR_FROM(i, j)
}

func (obj MATRIX_TYPE) JointIterator(b ConstMatrix) MatrixJointIterator {
  return obj.JOINT_ITERATOR(b)
}

func (obj MATRIX_TYPE) ITERATOR() *MATRIX_ITERATOR {
  r := MATRIX_ITERATOR{obj, 0, -1}
  r.Next()
  return &r
}

func (obj MATRIX_TYPE) ITERATOR_FROM(i, j int) *MATRIX_ITERATOR {
  r := MATRIX_ITERATOR{obj, i, j-1}
  r.Next()
  return &r
}

func (obj MATRIX_TYPE) JOINT_ITERATOR(b ConstMatrix) *MATRIX_JOINT_ITERATOR {
  r := MATRIX_JOINT_ITERATOR{obj.ITERATOR(), b.ConstIterator(), -1, -1, nil, nil}
  r.Next()
  return &r
}

/* iterator
 * -------------------------------------------------------------------------- */

type MATRIX_ITERATOR struct {
  m MATRIX_TYPE
  i, j int
}

func (obj *MATRIX_ITERATOR) Get() Scalar {
  return obj.GET()
}

func (obj *MATRIX_ITERATOR) GetConst() ConstScalar {
  return obj.GET()
}

func (obj *MATRIX_ITERATOR) GetMagic() MagicScalar {
  return obj.GET()
}

func (obj *MATRIX_ITERATOR) GET() SCALAR_TYPE {
  return obj.m.AT(obj.i, obj.j)
}

func (obj *MATRIX_ITERATOR) Ok() bool {
  return obj.i < obj.m.rowMax && obj.j < obj.m.colMax
}

func (obj *MATRIX_ITERATOR) next() {
  if obj.j == obj.m.cols-1 {
    obj.i = obj.i + 1
    obj.j = 0
  } else {
    obj.j = obj.j + 1
  }
}

func (obj *MATRIX_ITERATOR) Next() {
  obj.next()
  for obj.Ok() && obj.GET().nullScalar() {
    obj.next()
  }
}

func (obj *MATRIX_ITERATOR) Index() (int, int) {
  return obj.i, obj.j
}

func (obj *MATRIX_ITERATOR) Clone() *MATRIX_ITERATOR {
  return &MATRIX_ITERATOR{obj.m, obj.i, obj.j}
}

func (obj *MATRIX_ITERATOR) CloneIterator() MatrixIterator {
  return &MATRIX_ITERATOR{obj.m, obj.i, obj.j}
}

func (obj *MATRIX_ITERATOR) CloneConstIterator() MatrixConstIterator {
  return &MATRIX_ITERATOR{obj.m, obj.i, obj.j}
}

func (obj *MATRIX_ITERATOR) CloneMagicIterator() MatrixMagicIterator {
  return &MATRIX_ITERATOR{obj.m, obj.i, obj.j}
}

/* joint iterator
 * -------------------------------------------------------------------------- */

type MATRIX_JOINT_ITERATOR struct {
  it1 *MATRIX_ITERATOR
  it2  MatrixConstIterator
  i, j int
  s1   SCALAR_TYPE
  s2   ConstScalar
}

func (obj *MATRIX_JOINT_ITERATOR) Index() (int, int) {
  return obj.i, obj.j
}

func (obj *MATRIX_JOINT_ITERATOR) Ok() bool {
  return !(obj.s1 == nil || obj.s1.GET_METHOD_NAME() == STORED_TYPE(0)) ||
         !(obj.s2 == nil || obj.s2.GET_METHOD_NAME() == STORED_TYPE(0))
}

func (obj *MATRIX_JOINT_ITERATOR) Next() {
  ok1 := obj.it1.Ok()
  ok2 := obj.it2.Ok()
  obj.s1 = nil
  obj.s2 = nil
  if ok1 {
    obj.i, obj.j = obj.it1.Index()
    obj.s1       = obj.it1.GET()
  }
  if ok2 {
    i, j := obj.it2.Index()
    switch {
    case obj.i > i || (obj.i == i && obj.j > j) || !ok1:
      obj.i, obj.j = i, j
      obj.s1       = nil
      obj.s2       = obj.it2.GetConst()
    case obj.i == i && obj.j == j:
      obj.s2       = obj.it2.GetConst()
    }
  }
  if obj.s1 != nil {
    obj.it1.Next()
  }
  if obj.s2 != nil {
    obj.it2.Next()
  } else {
    obj.s2 = CONST_SCALAR_TYPE(0.0)
  }
}

func (obj *MATRIX_JOINT_ITERATOR) Get() (Scalar, ConstScalar) {
  if obj.s1 == nil {
    return nil, obj.s2
  } else {
    return obj.s1, obj.s2
  }
}

func (obj *MATRIX_JOINT_ITERATOR) GetConst() (ConstScalar, ConstScalar) {
  if obj.s1 == nil {
    return nil, obj.s2
  } else {
    return obj.s1, obj.s2
  }
}

func (obj *MATRIX_JOINT_ITERATOR) GET() (SCALAR_TYPE, ConstScalar) {
  return obj.s1, obj.s2
}

func (obj *MATRIX_JOINT_ITERATOR) Clone() *MATRIX_JOINT_ITERATOR {
  r := MATRIX_JOINT_ITERATOR{}
  r.it1 = obj.it1.Clone()
  r.it2 = obj.it2.CloneConstIterator()
  r.i   = obj.i
  r.j   = obj.j
  r.s1  = obj.s1
  r.s2  = obj.s2
  return &r
}

func (obj *MATRIX_JOINT_ITERATOR) CloneJointIterator() MatrixJointIterator {
  return obj.Clone()
}

func (obj *MATRIX_JOINT_ITERATOR) CloneConstJointIterator() MatrixConstJointIterator {
  return obj.Clone()
}