matrix_sparse_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_ITERATOR         STR_CONCAT(MATRIX_NAME, Iterator)
#define VECTOR_ITERATOR         STR_CONCAT(VECTOR_NAME, Iterator)
#define MATRIX_JOINT_ITERATOR   STR_CONCAT(MATRIX_NAME, JointIterator)
#define MATRIX_JOINT3_ITERATOR  STR_CONCAT(MATRIX_NAME, Joint3Iterator)

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

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
  tmp1       VECTOR_TYPE
  tmp2       VECTOR_TYPE
}

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

func NEW_MATRIX(rowIndices, colIndices []int, values []STORED_TYPE, rows, cols int) MATRIX_TYPE {
  m := NULL_MATRIX(rows, cols)
  if len(rowIndices) != len(colIndices) || len(colIndices) != len(values) {
    panic("number of row/col-indices does not match number of values")
  }
  for i := 0; i < len(colIndices); i++ {
    j1 := rowIndices[i]
    j2 := colIndices[i]
    if values[i] != 0.0 {
      m.At(j1, j2).SET_METHOD_NAME(values[i])
    }
  }
  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 it := matrix.ConstIterator(); it.Ok(); it.Next() {
    i, j := it.Index()
    r.AT(i,j).Set(it.GetConst())
  }
  return r
}

func (matrix MATRIX_TYPE) initTmp() {
  if matrix.tmp1 == nil || matrix.tmp1.Dim() < matrix.rows {
    matrix.tmp1 = NULL_VECTOR(matrix.rows)
  } else {
    matrix.tmp1 = matrix.tmp1.Slice(0, matrix.rows).(VECTOR_TYPE)
  }
  if matrix.tmp2 == nil || matrix.tmp2.Dim() < matrix.cols {
    matrix.tmp2 = NULL_VECTOR(matrix.cols)
  } else {
    matrix.tmp2 = matrix.tmp2.Slice(0, matrix.cols).(VECTOR_TYPE)
  }
}

/* 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,
    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))
  }
  return (matrix.rowOffset + i)*matrix.colMax + (matrix.colOffset + j)
}

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

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

func (matrix MATRIX_TYPE) AT(i, j int) SCALAR_TYPE {
  return matrix.values.AT(matrix.index(i, j))
}

func (matrix MATRIX_TYPE) ROW(i int) VECTOR_TYPE {
  var v VECTOR_TYPE
  v = NIL_VECTOR(matrix.cols)
  for j := 0; j < matrix.cols; j++ {
    if s := matrix.values.AT_(matrix.index(i, j)); !s.nullScalar() {
      v.AT(j).SET(s)
    }
  }
  return v
}

func (matrix MATRIX_TYPE) COL(j int) VECTOR_TYPE {
  var v VECTOR_TYPE
  v = NIL_VECTOR(matrix.rows)
  for i := 0; i < matrix.rows; i++ {
    if s := matrix.values.AT_(matrix.index(i, j)); !s.nullScalar() {
      v.AT(i).SET(s)
    }
  }
  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++ {
    if s := matrix.values.AT_(matrix.index(i, i)); !s.nullScalar() {
      v.At(i).Set(s)
    }
  }
  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 it := matrix.ConstIterator(); it.Ok(); it.Next() {
      i, j := it.Index()
      v.At(i*matrix.cols + j).Set(matrix.ConstAt(i, j))
    }
    return v
  } else {
    return 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 it := a.Iterator(); it.Ok(); it.Next() {
    i, j := it.Index()
    it.Get().Set(b.ConstAt(i, j))
  }
}

func (matrix MATRIX_TYPE) SetIdentity() {
  c := NewScalar(matrix.ElementType(), 1.0)
  for it := matrix.Iterator(); it.Ok(); it.Next() {
    i, j := it.Index()
    if i == j {
      it.Get().Set(c)
    } else {
      it.Get().Reset()
    }
  }
}

func (matrix MATRIX_TYPE) Reset() {
  for it := matrix.Iterator(); it.Ok(); it.Next() {
    it.Get().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.Swap(k1, k2)
}

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

func (matrix MATRIX_TYPE) Tip() {
  mn      := matrix.values.Dim()
  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.Swap(k, cycle)

      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.AT(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.ConstAt(matrix.index(i, j)).GetInt8()
}

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

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

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

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

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

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

func (matrix MATRIX_TYPE) ConstAt(i, j int) ConstScalar {
  return matrix.values.ConstAt(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 {
  var v VECTOR_TYPE
  i = matrix.index(i, 0)
  v = matrix.values.Slice(i, i + matrix.cols).(VECTOR_TYPE)
  return v
}

func (matrix MATRIX_TYPE) ConstCol(i int) ConstVector {
  return matrix.COL(i)
}

func (matrix MATRIX_TYPE) ConstDiag() ConstVector {
  return matrix.DIAG()
}

func (matrix MATRIX_TYPE) IsSymmetric(epsilon float64) bool {
  if n, m := matrix.Dims(); n != m {
    return false
  }
  for it := matrix.ConstIterator(); it.Ok(); it.Next() {
    i, j := it.Index()
    if !matrix.ConstAt(i,j).Equals(matrix.ConstAt(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 {
  m := &MATRIX_NAME{
    values    :  NULL_VECTOR(matrix.values.Dim()),
    rows      :  matrix.cols,
    cols      :  matrix.rows,
    rowOffset :  matrix.colOffset,
    rowMax    :  matrix.colMax,
    colOffset :  matrix.rowOffset,
    colMax    :  matrix.rowMax,
    tmp1      :  matrix.tmp2,
    tmp2      :  matrix.tmp1 }
  for k1, value := range matrix.values.values {
    // transform indices so that iterators operate correctly
    i1, j1 := matrix.ij(k1)
    k2     := m.index(j1, i1)
    m.values.values[k2] = value
    m.values.indexInsert(k2)
  }
  return m
}

func (matrix MATRIX_TYPE) ResetDerivatives() {
  for it := matrix.MagicIterator(); it.Ok(); it.Next() {
    it.GetMagic().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, v := range matrix.values.values {
    if err := v.SetVariable(i, matrix.values.Dim(), 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, "%d %d\n", m.rows, m.cols); err != nil {
    return err
  }
  for it := m.ITERATOR(); it.Ok(); it.Next() {
    i, j := it.Index()
    if _, err := fmt.Fprintf(w, "%d %d %v\n", i, j, it.GET()); err != nil {
      return err
    }
  }
  return nil
}

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

  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)
  }
  // scan header
  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 len(fields) != 2 {
      return fmt.Errorf("invalid sparse matrix format")
    }
    if v, err := strconv.ParseInt(fields[0], 10, 64); err != nil {
      return err
    } else {
      rows = int(v)
    }
    if v, err := strconv.ParseInt(fields[1], 10, 64); err != nil {
      return err
    } else {
      cols = int(v)
    }
    break
  }
  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 len(fields) != 3 {
      return fmt.Errorf("invalid sparse matrix format")
    }
    if v, err := strconv.ParseInt(fields[0], 10, 64); err != nil {
      return err
    } else {
      rowIndices = append(rowIndices, int(v))
    }
    if v, err := strconv.ParseInt(fields[1], 10, 64); err != nil {
      return err
    } else {
      colIndices = append(colIndices, int(v))
    }
    if v, err := strconv.ParseFloat(fields[2], 64); err != nil {
      return err
    } else {
      values = append(values, STORED_TYPE(v))
    }
  }
  *m = *NEW_MATRIX(rowIndices, colIndices, values, rows, cols)

  return nil
}

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

func (obj MATRIX_TYPE) MarshalJSON() ([]byte, error) {
  if obj.rowMax > obj.rows || obj.colMax > obj.cols {
    n, m := obj.Dims()
    tmp  := NULL_MATRIX(n, m)
    tmp.Set(obj)
    obj = tmp
  }
  k := []int{}
  v := []STORED_TYPE{}
  r := struct{Index []int; Value []STORED_TYPE; Rows int; Cols int}{}
  for it := obj.values.ConstIterator(); it.Ok(); it.Next() {
    k = append(k, it.Index())
    v = append(v, STORED_TYPE(it.GetConst().GET_METHOD_NAME()))
  }
  r.Index  = k
  r.Value  = v
  r.Rows   = obj.rows
  r.Cols   = obj.cols
  return json.MarshalIndent(r, "", "  ")
}

func (obj MATRIX_TYPE) UnmarshalJSON(data []byte) error {
  r := struct{Index []int; Value []STORED_TYPE; Rows int; Cols int}{}
  if err := json.Unmarshal(data, &r); err != nil {
    return err
  }
  if len(r.Index) != len(r.Value) {
    return fmt.Errorf("invalid sparse vector")
  }
  obj.values     = NEW_VECTOR(r.Index, r.Value, r.Rows*r.Cols)
  obj.rows       = r.Rows
  obj.rowMax     = r.Rows
  obj.rowOffset  = 0
  obj.cols       = r.Cols
  obj.colMax     = r.Cols
  obj.colOffset  = 0
  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.values.ITERATOR(), obj}
  return &r
}

func (obj MATRIX_TYPE) ITERATOR_FROM(i, j int) *MATRIX_ITERATOR {
  k := obj.index(i, j)
  r := MATRIX_ITERATOR{*obj.values.ITERATOR_FROM(k), obj}
  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
}

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

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

type MATRIX_ITERATOR struct {
  VECTOR_ITERATOR
  m MATRIX_TYPE
}

func (obj *MATRIX_ITERATOR) Index() (int, int) {
  return obj.m.ij(obj.VECTOR_ITERATOR.Index())
}

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

func (obj *MATRIX_ITERATOR) CloneConstIterator() MatrixConstIterator {
  return &MATRIX_ITERATOR{*obj.VECTOR_ITERATOR.Clone(), obj.m}
}

func (obj *MATRIX_ITERATOR) CloneMagicIterator() MatrixMagicIterator {
  return &MATRIX_ITERATOR{*obj.VECTOR_ITERATOR.Clone(), obj.m}
}

func (obj *MATRIX_ITERATOR) CloneIterator() MatrixIterator {
  return &MATRIX_ITERATOR{*obj.VECTOR_ITERATOR.Clone(), obj.m}
}

/* 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()
}

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

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

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

func (obj *MATRIX_JOINT3_ITERATOR) Ok() bool {
  return !(obj.s1 == nil || obj.s1.GET_METHOD_NAME() == 0.0) ||
         !(obj.s2 == nil || obj.s2.GET_METHOD_NAME() == 0.0) ||
         !(obj.s3 == nil || obj.s3.GET_METHOD_NAME() == 0.0)
}

func (obj *MATRIX_JOINT3_ITERATOR) Next() {
  ok1 := obj.it1.Ok()
  ok2 := obj.it2.Ok()
  ok3 := obj.it3.Ok()
  obj.s1 = nil
  obj.s2 = nil
  obj.s3 = 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  = i
      obj.j  = j
      obj.s1 = nil
      obj.s2 = obj.it2.GetConst()
    case obj.i == i && obj.j == j:
      obj.s2 = obj.it2.GetConst()
    }
  }
  if ok3 {
    i, j := obj.it3.Index()
    switch {
    case obj.i > i || (obj.i == i && obj.j > j) || (!ok1 && !ok2):
      obj.i  = i
      obj.j  = j
      obj.s1 = nil
      obj.s2 = nil
      obj.s3 = obj.it3.GetConst()
    case obj.i == i && obj.j == j:
      obj.s3 = obj.it3.GetConst()
    }
  }
  if obj.s1 != nil {
    obj.it1.Next()
  }
  if obj.s2 != nil {
    obj.it2.Next()
  } else {
    obj.s2 = CONST_SCALAR_TYPE(0.0)
  }
  if obj.s3 != nil {
    obj.it3.Next()
  } else {
    obj.s3 = CONST_SCALAR_TYPE(0.0)
  }
}

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

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