diff --git a/.travis.yml b/.travis.yml
index 3db93d4..eebfa6c 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,5 +1,4 @@
 sudo: false
 language: ruby
 rvm:
-  - 2.5.0
-before_install: gem install bundler -v 1.15.4
+  - ruby-head
diff --git a/lib/matrix.rb b/lib/matrix.rb
index 5cc9697..7f338bb 100644
--- a/lib/matrix.rb
+++ b/lib/matrix.rb
@@ -145,8 +145,8 @@ def Matrix.identity(n)
     scalar(n, 1)
   end
   class << Matrix
-    alias unit identity
-    alias I identity
+    alias_method :unit, :identity
+    alias_method :I, :identity
   end
 
   #
@@ -289,10 +289,9 @@ def initialize(rows, column_count = rows[0].size)
     @column_count = column_count
   end
 
-  def new_matrix(rows, column_count = rows[0].size) # :nodoc:
+  private def new_matrix(rows, column_count = rows[0].size) # :nodoc:
     self.class.send(:new, rows, column_count) # bypass privacy of Matrix.new
   end
-  private :new_matrix
 
   #
   # Returns element (+i+,+j+) of the matrix.  That is: row +i+, column +j+.
@@ -303,12 +302,107 @@ def [](i, j)
   alias element []
   alias component []
 
+  #
+  # :call-seq:
+  #   matrix[range, range] = matrix/element
+  #   matrix[range, integer] = vector/column_matrix/element
+  #   matrix[integer, range] = vector/row_matrix/element
+  #   matrix[integer, integer] = element
+  #
+  # Set element or elements of matrix.
   def []=(i, j, v)
-    @rows[i][j] = v
+    raise FrozenError, "can't modify frozen Matrix" if frozen?
+    rows = check_range(i, :row) or row = check_int(i, :row)
+    columns = check_range(j, :column) or column = check_int(j, :column)
+    if rows && columns
+      set_row_and_col_range(rows, columns, v)
+    elsif rows
+      set_row_range(rows, column, v)
+    elsif columns
+      set_col_range(row, columns, v)
+    else
+      set_value(row, column, v)
+    end
   end
   alias set_element []=
   alias set_component []=
-  private :[]=, :set_element, :set_component
+  private :set_element, :set_component
+
+  # Returns range or nil
+  private def check_range(val, direction)
+    return unless val.is_a?(Range)
+    count = direction == :row ? row_count : column_count
+    CoercionHelper.check_range(val, count, direction)
+  end
+
+  private def check_int(val, direction)
+    count = direction == :row ? row_count : column_count
+    CoercionHelper.check_int(val, count, direction)
+  end
+
+  private def set_value(row, col, value)
+    raise ErrDimensionMismatch, "Expected a a value, got a #{value.class}" if value.respond_to?(:to_matrix)
+
+    @rows[row][col] = value
+  end
+
+  private def set_row_and_col_range(row_range, col_range, value)
+    if value.is_a?(Matrix)
+      if row_range.size != value.row_count || col_range.size != value.column_count
+        raise ErrDimensionMismatch, [
+          'Expected a Matrix of dimensions',
+          "#{row_range.size}x#{col_range.size}",
+          'got',
+          "#{value.row_count}x#{value.column_count}",
+        ].join(' ')
+      end
+      source = value.instance_variable_get :@rows
+      row_range.each_with_index do |row, i|
+        @rows[row][col_range] = source[i]
+      end
+    elsif value.is_a?(Vector)
+      raise ErrDimensionMismatch, 'Expected a Matrix or a value, got a Vector'
+    else
+      value_to_set = Array.new(col_range.size, value)
+      row_range.each do |i|
+        @rows[i][col_range] = value_to_set
+      end
+    end
+  end
+
+  private def set_row_range(row_range, col, value)
+    if value.is_a?(Vector)
+      Matrix.Raise ErrDimensionMismatch unless row_range.size == value.size
+      set_column_vector(row_range, col, value)
+    elsif value.is_a?(Matrix)
+      Matrix.Raise ErrDimensionMismatch unless value.column_count == 1
+      value = value.column(0)
+      Matrix.Raise ErrDimensionMismatch unless row_range.size == value.size
+      set_column_vector(row_range, col, value)
+    else
+      @rows[row_range].each{|e| e[col] = value }
+    end
+  end
+
+  private def set_column_vector(row_range, col, value)
+    value.each_with_index do |e, index|
+      r = row_range.begin + index
+      @rows[r][col] = e
+    end
+  end
+
+  private def set_col_range(row, col_range, value)
+    value = if value.is_a?(Vector)
+      value.to_a
+    elsif value.is_a?(Matrix)
+      Matrix.Raise ErrDimensionMismatch unless value.row_count == 1
+      value.row(0).to_a
+    else
+      Array.new(col_range.size, value)
+    end
+    Matrix.Raise ErrDimensionMismatch unless col_range.size == value.size
+    @rows[row][col_range] = value
+  end
 
   #
   # Returns the number of rows.
@@ -361,16 +455,48 @@ def column(j) # :yield: e
   #
   # Returns a matrix that is the result of iteration of the given block over all
   # elements of the matrix.
+  # Elements can be restricted by passing an argument:
+  # * :all (default): yields all elements
+  # * :diagonal: yields only elements on the diagonal
+  # * :off_diagonal: yields all elements except on the diagonal
+  # * :lower: yields only elements on or below the diagonal
+  # * :strict_lower: yields only elements below the diagonal
+  # * :strict_upper: yields only elements above the diagonal
+  # * :upper: yields only elements on or above the diagonal
   #   Matrix[ [1,2], [3,4] ].collect { |e| e**2 }
   #     => 1  4
   #        9 16
   #
-  def collect(&block) # :yield: e
-    return to_enum(:collect) unless block_given?
-    rows = @rows.collect{|row| row.collect(&block)}
-    new_matrix rows, column_count
+  def collect(which = :all, &block) # :yield: e
+    return to_enum(:collect, which) unless block_given?
+    dup.collect!(which, &block)
+  end
+  alias_method :map, :collect
+
+  #
+  # Invokes the given block for each element of matrix, replacing the element with the value
+  # returned by the block.
+  # Elements can be restricted by passing an argument:
+  # * :all (default): yields all elements
+  # * :diagonal: yields only elements on the diagonal
+  # * :off_diagonal: yields all elements except on the diagonal
+  # * :lower: yields only elements on or below the diagonal
+  # * :strict_lower: yields only elements below the diagonal
+  # * :strict_upper: yields only elements above the diagonal
+  # * :upper: yields only elements on or above the diagonal
+  #
+  def collect!(which = :all)
+    return to_enum(:collect!, which) unless block_given?
+    raise FrozenError, "can't modify frozen Matrix" if frozen?
+    each_with_index(which){ |e, row_index, col_index| @rows[row_index][col_index] = yield e }
+  end
+
+  alias map! collect!
+
+  def freeze
+    @rows.freeze
+    super
   end
-  alias map collect
 
   #
   # Yields all elements of the matrix, starting with those of the first row,
@@ -812,6 +938,7 @@ def antisymmetric?
     end
     true
   end
+  alias_method :skew_symmetric?, :antisymmetric?
 
   #
   # Returns +true+ if this is a unitary matrix
@@ -865,12 +992,11 @@ def eql?(other)
   end
 
   #
-  # Returns a clone of the matrix, so that the contents of each do not reference
-  # identical objects.
-  # There should be no good reason to do this since Matrices are immutable.
+  # Called for dup & clone.
   #
-  def clone
-    new_matrix @rows.map(&:dup), column_count
+  private def initialize_copy(m)
+    super
+    @rows = @rows.map(&:dup) unless frozen?
   end
 
   #
@@ -1012,9 +1138,9 @@ def inverse
     Matrix.Raise ErrDimensionMismatch unless square?
     self.class.I(row_count).send(:inverse_from, self)
   end
-  alias inv inverse
+  alias_method :inv, :inverse
 
-  def inverse_from(src) # :nodoc:
+  private def inverse_from(src) # :nodoc:
     last = row_count - 1
     a = src.to_a
 
@@ -1057,7 +1183,6 @@ def inverse_from(src) # :nodoc:
     end
     self
   end
-  private :inverse_from
 
   #
   # Matrix exponentiation.
@@ -1164,7 +1289,7 @@ def determinant
   # with smaller bignums (if any), while a matrix of Float will usually have
   # intermediate results with better precision.
   #
-  def determinant_bareiss
+  private def determinant_bareiss
     size = row_count
     last = size - 1
     a = to_a
@@ -1190,7 +1315,6 @@ def determinant_bareiss
     end
     sign * pivot
   end
-  private :determinant_bareiss
 
   #
   # deprecated; use Matrix#determinant
@@ -1199,7 +1323,7 @@ def determinant_e
     warn "Matrix#determinant_e is deprecated; use #determinant", uplevel: 1
     determinant
   end
-  alias det_e determinant_e
+  alias_method :det_e, :determinant_e
 
   #
   # Returns a new matrix resulting by stacking horizontally
@@ -1276,7 +1400,7 @@ def trace
       tr + @rows[i][i]
     end
   end
-  alias tr trace
+  alias_method :tr, :trace
 
   #
   # Returns the transpose of the matrix.
@@ -1292,7 +1416,7 @@ def transpose
     return self.class.empty(column_count, 0) if row_count.zero?
     new_matrix @rows.transpose, row_count
   end
-  alias t transpose
+  alias_method :t, :transpose
 
   #
   # Returns a new matrix resulting by stacking vertically
@@ -1321,7 +1445,7 @@ def vstack(*matrices)
   def eigensystem
     EigenvalueDecomposition.new(self)
   end
-  alias eigen eigensystem
+  alias_method :eigen, :eigensystem
 
   #
   # Returns the LUP decomposition of the matrix; see +LUPDecomposition+.
@@ -1336,7 +1460,7 @@ def eigensystem
   def lup
     LUPDecomposition.new(self)
   end
-  alias lup_decomposition lup
+  alias_method :lup_decomposition, :lup
 
   #--
   # COMPLEX ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
@@ -1354,7 +1478,7 @@ def lup
   def conjugate
     collect(&:conjugate)
   end
-  alias conj conjugate
+  alias_method :conj, :conjugate
 
   #
   # Returns the imaginary part of the matrix.
@@ -1368,7 +1492,7 @@ def conjugate
   def imaginary
     collect(&:imaginary)
   end
-  alias imag imaginary
+  alias_method :imag, :imaginary
 
   #
   # Returns the real part of the matrix.
@@ -1392,7 +1516,7 @@ def real
   def rect
     [real, imag]
   end
-  alias rectangular rect
+  alias_method :rectangular, :rect
 
   #--
   # CONVERTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
@@ -1505,7 +1629,7 @@ module ConversionHelper # :nodoc:
     # Converts the obj to an Array. If copy is set to true
     # a copy of obj will be made if necessary.
     #
-    def convert_to_array(obj, copy = false) # :nodoc:
+    private def convert_to_array(obj, copy = false) # :nodoc:
       case obj
       when Array
         copy ? obj.dup : obj
@@ -1521,7 +1645,6 @@ def convert_to_array(obj, copy = false) # :nodoc:
         converted
       end
     end
-    private :convert_to_array
   end
 
   extend ConversionHelper
@@ -1531,14 +1654,13 @@ module CoercionHelper # :nodoc:
     # Applies the operator +oper+ with argument +obj+
     # through coercion of +obj+
     #
-    def apply_through_coercion(obj, oper)
+    private def apply_through_coercion(obj, oper)
       coercion = obj.coerce(self)
       raise TypeError unless coercion.is_a?(Array) && coercion.length == 2
       coercion[0].public_send(oper, coercion[1])
     rescue
       raise TypeError, "#{obj.inspect} can't be coerced into #{self.class}"
     end
-    private :apply_through_coercion
 
     #
     # Helper method to coerce a value into a specific class.
@@ -1566,6 +1688,26 @@ def self.coerce_to_int(obj)
     def self.coerce_to_matrix(obj)
       coerce_to(obj, Matrix, :to_matrix)
     end
+
+    # Returns `nil` for non Ranges
+    # Checks range validity, return canonical range with 0 <= begin <= end < count
+    def self.check_range(val, count, kind)
+      canonical = (val.begin + (val.begin < 0 ? count : 0))..
+                  (val.end ? val.end + (val.end < 0 ? count : 0) - (val.exclude_end? ? 1 : 0)
+                           : count - 1)
+      unless 0 <= canonical.begin && canonical.begin <= canonical.end && canonical.end < count
+        raise IndexError, "given range #{val} is outside of #{kind} dimensions: 0...#{count}"
+      end
+      canonical
+    end
+
+    def self.check_int(val, count, kind)
+      val = CoercionHelper.coerce_to_int(val)
+      if val >= count || val < -count
+        raise IndexError, "given #{kind} #{val} is outside of #{-count}...#{count}"
+      end
+      val
+    end
   end
 
   include CoercionHelper
@@ -1660,6 +1802,9 @@ def **(other)
 # To access elements:
 # * #[](i)
 #
+# To set elements:
+# * #[]=(i, v)
+#
 # To enumerate the elements:
 # * #each2(v)
 # * #collect2(v)
@@ -1682,8 +1827,10 @@ def **(other)
 # * #inner_product(v), dot(v)
 # * #cross_product(v), cross(v)
 # * #collect
+# * #collect!
 # * #magnitude
 # * #map
+# * #map!
 # * #map2(v)
 # * #norm
 # * #normalize
@@ -1762,7 +1909,11 @@ def initialize(array)
   # ACCESSING
 
   #
-  # Returns element number +i+ (starting at zero) of the vector.
+  # :call-seq:
+  #   vector[range]
+  #   vector[integer]
+  #
+  # Returns element or elements of the vector.
   #
   def [](i)
     @elements[i]
@@ -1770,12 +1921,44 @@ def [](i)
   alias element []
   alias component []
 
+  #
+  # :call-seq:
+  #   vector[range] = new_vector
+  #   vector[range] = row_matrix
+  #   vector[range] = new_element
+  #   vector[integer] = new_element
+  #
+  # Set element or elements of vector.
+  #
   def []=(i, v)
-    @elements[i]= v
+    raise FrozenError, "can't modify frozen Vector" if frozen?
+    if i.is_a?(Range)
+      range = Matrix::CoercionHelper.check_range(i, size, :vector)
+      set_range(range, v)
+    else
+      index = Matrix::CoercionHelper.check_int(i, size, :index)
+      set_value(index, v)
+    end
   end
   alias set_element []=
   alias set_component []=
-  private :[]=, :set_element, :set_component
+  private :set_element, :set_component
+
+  private def set_value(index, value)
+    @elements[index] = value
+  end
+
+  private def set_range(range, value)
+    if value.is_a?(Vector)
+      raise ArgumentError, "vector to be set has wrong size" unless range.size == value.size
+      @elements[range] = value.elements
+    elsif value.is_a?(Matrix)
+      Matrix.Raise ErrDimensionMismatch unless value.row_count == 1
+      @elements[range] = value.row(0).elements
+    else
+      @elements[range] = Array.new(range.size, value)
+    end
+  end
 
   # Returns a vector with entries rounded to the given precision
   # (see Float#round)
@@ -1872,6 +2055,20 @@ def zero?
     all?(&:zero?)
   end
 
+  def freeze
+    @elements.freeze
+    super
+  end
+
+  #
+  # Called for dup & clone.
+  #
+  private def initialize_copy(v)
+    super
+    @elements = @elements.dup unless frozen?
+  end
+
+
   #--
   # COMPARING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
   #++
@@ -1889,13 +2086,6 @@ def eql?(other)
     @elements.eql? other.elements
   end
 
-  #
-  # Returns a copy of the vector.
-  #
-  def clone
-    self.class.elements(@elements)
-  end
-
   #
   # Returns a hash-code for the vector.
   #
@@ -2044,7 +2234,18 @@ def collect(&block) # :yield: e
     els = @elements.collect(&block)
     self.class.elements(els, false)
   end
-  alias map collect
+  alias_method :map, :collect
+
+  #
+  # Like Array#collect!
+  #
+  def collect!(&block)
+    return to_enum(:collect!) unless block_given?
+    raise FrozenError, "can't modify frozen Vector" if frozen?
+    @elements.collect!(&block)
+    self
+  end
+  alias map! collect!
 
   #
   # Returns the modulus (Pythagorean distance) of the vector.
@@ -2053,8 +2254,8 @@ def collect(&block) # :yield: e
   def magnitude
     Math.sqrt(@elements.inject(0) {|v, e| v + e.abs2})
   end
-  alias r magnitude
-  alias norm magnitude
+  alias_method :r, :magnitude
+  alias_method :norm, :magnitude
 
   #
   # Like Vector#collect2, but returns a Vector instead of an Array.
@@ -2080,7 +2281,7 @@ def normalize
   end
 
   #
-  # Returns an angle with another vector. Result is within the [0...Math::PI].
+  # Returns an angle with another vector. Result is within the [0..Math::PI].
   #   Vector[1,0].angle_with(Vector[0,1])
   #   # => Math::PI / 2
   #
@@ -2089,8 +2290,12 @@ def angle_with(v)
     Vector.Raise ErrDimensionMismatch if size != v.size
     prod = magnitude * v.magnitude
     raise ZeroVectorError, "Can't get angle of zero vector" if prod == 0
-
-    Math.acos( inner_product(v) / prod )
+    dot = inner_product(v)
+    if dot.abs >= prod
+      dot.positive? ? 0 : Math::PI
+    else
+      Math.acos(dot / prod)
+    end
   end
 
   #--
diff --git a/lib/matrix/eigenvalue_decomposition.rb b/lib/matrix/eigenvalue_decomposition.rb
index 919db9e..bf66376 100644
--- a/lib/matrix/eigenvalue_decomposition.rb
+++ b/lib/matrix/eigenvalue_decomposition.rb
@@ -43,7 +43,7 @@ def initialize(a)
     def eigenvector_matrix
       Matrix.send(:new, build_eigenvectors.transpose)
     end
-    alias v eigenvector_matrix
+    alias_method :v, :eigenvector_matrix
 
     # Returns the inverse of the eigenvector matrix +V+
     #
@@ -52,7 +52,7 @@ def eigenvector_matrix_inv
       r = r.transpose.inverse unless @symmetric
       r
     end
-    alias v_inv eigenvector_matrix_inv
+    alias_method :v_inv, :eigenvector_matrix_inv
 
     # Returns the eigenvalues in an array
     #
@@ -73,7 +73,7 @@ def eigenvectors
     def eigenvalue_matrix
       Matrix.diagonal(*eigenvalues)
     end
-    alias d eigenvalue_matrix
+    alias_method :d, :eigenvalue_matrix
 
     # Returns [eigenvector_matrix, eigenvalue_matrix, eigenvector_matrix_inv]
     #
@@ -82,8 +82,8 @@ def to_ary
     end
     alias_method :to_a, :to_ary
 
-  private
-    def build_eigenvectors
+
+    private def build_eigenvectors
       # JAMA stores complex eigenvectors in a strange way
       # See http://web.archive.org/web/20111016032731/http://cio.nist.gov/esd/emaildir/lists/jama/msg01021.html
       @e.each_with_index.map do |imag, i|
@@ -96,9 +96,10 @@ def build_eigenvectors
         end
       end
     end
+
     # Complex scalar division.
 
-    def cdiv(xr, xi, yr, yi)
+    private def cdiv(xr, xi, yr, yi)
       if (yr.abs > yi.abs)
         r = yi/yr
         d = yr + r*yi
@@ -113,7 +114,7 @@ def cdiv(xr, xi, yr, yi)
 
     # Symmetric Householder reduction to tridiagonal form.
 
-    def tridiagonalize
+    private def tridiagonalize
 
       #  This is derived from the Algol procedures tred2 by
       #  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
@@ -231,7 +232,7 @@ def tridiagonalize
 
     # Symmetric tridiagonal QL algorithm.
 
-    def diagonalize
+    private def diagonalize
       #  This is derived from the Algol procedures tql2, by
       #  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
       #  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
@@ -350,7 +351,7 @@ def diagonalize
 
     # Nonsymmetric reduction to Hessenberg form.
 
-    def reduce_to_hessenberg
+    private def reduce_to_hessenberg
       #  This is derived from the Algol procedures orthes and ortran,
       #  by Martin and Wilkinson, Handbook for Auto. Comp.,
       #  Vol.ii-Linear Algebra, and the corresponding
@@ -440,11 +441,9 @@ def reduce_to_hessenberg
       end
     end
 
-
-
     # Nonsymmetric reduction from Hessenberg to real Schur form.
 
-    def hessenberg_to_real_schur
+    private def hessenberg_to_real_schur
 
       #  This is derived from the Algol procedure hqr2,
       #  by Martin and Wilkinson, Handbook for Auto. Comp.,
diff --git a/test/matrix/test_matrix.rb b/test/matrix/test_matrix.rb
index d40e0c0..3ecb9d2 100644
--- a/test/matrix/test_matrix.rb
+++ b/test/matrix/test_matrix.rb
@@ -283,7 +283,18 @@ def test_column
   end
 
   def test_collect
-    assert_equal(Matrix[[1, 4, 9], [16, 25, 36]], @m1.collect {|x| x ** 2 })
+    m1 = Matrix.zero(2,2)
+    m2 = Matrix.build(3,4){|row, col| 1}
+
+    assert_equal(Matrix[[5, 5, 5, 5], [5, 5, 5, 5], [5, 5, 5, 5]], m2.collect{|e| e * 5})
+    assert_equal(Matrix[[7, 0],[0, 7]], m1.collect(:diagonal){|e| e + 7})
+    assert_equal(Matrix[[0, 5],[5, 0]], m1.collect(:off_diagonal){|e| e + 5})
+    assert_equal(Matrix[[8, 1, 1, 1], [8, 8, 1, 1], [8, 8, 8, 1]], m2.collect(:lower){|e| e + 7})
+    assert_equal(Matrix[[1, 1, 1, 1], [-11, 1, 1, 1], [-11, -11, 1, 1]], m2.collect(:strict_lower){|e| e - 12})
+    assert_equal(Matrix[[1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]], m2.collect(:strict_upper){|e| e ** 2})
+    assert_equal(Matrix[[-1, -1, -1, -1], [1, -1, -1, -1], [1, 1, -1, -1]], m2.collect(:upper){|e| -e})
+    assert_raise(ArgumentError) {m1.collect(:test){|e| e + 7}}
+    assert_not_equal(m2, m2.collect {|e| e * 2 })
   end
 
   def test_minor
@@ -623,6 +634,134 @@ def test_combine
     assert_equal Matrix.empty(0,3), Matrix.combine(Matrix.empty(0,3), Matrix.empty(0,3)) { raise }
   end
 
+  def test_set_element
+    src = Matrix[
+      [1, 2, 3, 4],
+      [5, 6, 7, 8],
+      [9, 10, 11, 12],
+    ]
+    rows = {
+      range:   [1..2, 1...3, 1..-1, -2..2, 1.., 1..., -2.., -2...],
+      int:     [2, -1],
+      invalid: [-4, 4, -4..2, 2..-4, 0...0, 2..0, -4..],
+    }
+    columns = {
+      range:   [2..3, 2...4, 2..-1, -2..3, 2.., 2..., -2..., -2..],
+      int:     [3, -1],
+      invalid: [-5, 5, -5..2, 2..-5, 0...0, -5..],
+    }
+    values = {
+      element: 42,
+      matrix:  Matrix[[20, 21], [22, 23]],
+      vector:  Vector[30, 31],
+      row:     Matrix[[60, 61]],
+      column:  Matrix[[50], [51]],
+      mismatched_matrix: Matrix.identity(3),
+      mismatched_vector: Vector[0, 1, 2, 3],
+    }
+    solutions = {
+      [:int, :int] => {
+        element: Matrix[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 42]],
+      },
+      [:range , :int] => {
+        element: Matrix[[1, 2, 3, 4], [5, 6, 7, 42], [9, 10, 11, 42]],
+        column:  Matrix[[1, 2, 3, 4], [5, 6, 7, 50], [9, 10, 11, 51]],
+        vector:  Matrix[[1, 2, 3, 4], [5, 6, 7, 30], [9, 10, 11, 31]],
+      },
+      [:int, :range] => {
+        element: Matrix[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 42, 42]],
+        row:     Matrix[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 60, 61]],
+        vector:  Matrix[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 30, 31]],
+      },
+      [:range , :range] => {
+        element: Matrix[[1, 2, 3, 4], [5, 6, 42, 42], [9, 10, 42, 42]],
+        matrix:  Matrix[[1, 2, 3, 4], [5, 6, 20, 21], [9, 10, 22, 23]],
+      },
+    }
+    solutions.default = Hash.new(IndexError)
+
+    rows.each do |row_style, row_arguments|
+      row_arguments.each do |row_argument|
+        columns.each do |column_style, column_arguments|
+          column_arguments.each do |column_argument|
+            values.each do |value_type, value|
+              expected = solutions[[row_style, column_style]][value_type] || Matrix::ErrDimensionMismatch
+
+              result = src.clone
+              begin
+                result[row_argument, column_argument] = value
+                assert_equal expected, result,
+                  "m[#{row_argument.inspect}][#{column_argument.inspect}] = #{value.inspect} failed"
+              rescue Exception => e
+                raise if e.class != expected
+              end
+            end
+          end
+        end
+      end
+    end
+  end
+
+  def test_map!
+    m1 = Matrix.zero(2,2)
+    m2 = Matrix.build(3,4){|row, col| 1}
+    m3 = Matrix.zero(3,5).freeze
+    m4 = Matrix.empty.freeze
+
+    assert_equal Matrix[[5, 5, 5, 5], [5, 5, 5, 5], [5, 5, 5, 5]], m2.map!{|e| e * 5}
+    assert_equal Matrix[[7, 0],[0, 7]], m1.map!(:diagonal){|e| e + 7}
+    assert_equal Matrix[[7, 5],[5, 7]], m1.map!(:off_diagonal){|e| e + 5}
+    assert_equal Matrix[[12, 5, 5, 5], [12, 12, 5, 5], [12, 12, 12, 5]], m2.map!(:lower){|e| e + 7}
+    assert_equal Matrix[[12, 5, 5, 5], [0, 12, 5, 5], [0, 0, 12, 5]], m2.map!(:strict_lower){|e| e - 12}
+    assert_equal Matrix[[12, 25, 25, 25], [0, 12, 25, 25], [0, 0, 12, 25]], m2.map!(:strict_upper){|e| e ** 2}
+    assert_equal Matrix[[-12, -25, -25, -25], [0, -12, -25, -25], [0, 0, -12, -25]], m2.map!(:upper){|e| -e}
+    assert_equal m1, m1.map!{|e| e ** 2 }
+    assert_equal m2, m2.map!(:lower){ |e| e - 3 }
+    assert_raise(ArgumentError) {m1.map!(:test){|e| e + 7}}
+    assert_raise(FrozenError) { m3.map!{|e| e * 2} }
+    assert_raise(FrozenError) { m4.map!{} }
+  end
+
+  def test_freeze
+    m = Matrix[[1, 2, 3],[4, 5, 6]]
+    f = m.freeze
+    assert_equal true, f.frozen?
+    assert m.equal?(f)
+    assert m.equal?(f.freeze)
+    assert_raise(FrozenError){ m[0, 1] = 56 }
+    assert_equal m.dup, m
+  end
+
+  def test_clone
+    a = Matrix[[4]]
+    def a.foo
+      42
+    end
+
+    m = a.clone
+    m[0, 0] = 2
+    assert_equal a, m * 2
+    assert_equal 42, m.foo
+
+    a.freeze
+    m = a.clone
+    assert m.frozen?
+    assert_equal 42, m.foo
+  end
+
+  def test_dup
+    a = Matrix[[4]]
+    def a.foo
+      42
+    end
+    a.freeze
+
+    m = a.dup
+    m[0, 0] = 2
+    assert_equal a, m * 2
+    assert !m.respond_to?(:foo)
+  end
+
   def test_eigenvalues_and_eigenvectors_symmetric
     m = Matrix[
       [8, 1],
diff --git a/test/matrix/test_vector.rb b/test/matrix/test_vector.rb
index 5278512..8b771ef 100644
--- a/test/matrix/test_vector.rb
+++ b/test/matrix/test_vector.rb
@@ -27,6 +27,108 @@ def test_basis
     assert_raise(ArgumentError) { Vector.basis(index: 3) }
   end
 
+  def test_get_element
+    assert_equal(@v1[0..], [1, 2, 3])
+    assert_equal(@v1[0..1], [1, 2])
+    assert_equal(@v1[2], 3)
+    assert_equal(@v1[4], nil)
+  end
+
+  def test_set_element
+
+    assert_block do
+      v = Vector[5, 6, 7, 8, 9]
+      v[1..2] = Vector[1, 2]
+      v == Vector[5, 1, 2, 8, 9]
+    end
+
+    assert_block do
+      v = Vector[6, 7, 8]
+      v[1..2] = Matrix[[1, 3]]
+      v == Vector[6, 1, 3]
+    end
+
+    assert_block do
+      v = Vector[1, 2, 3, 4, 5, 6]
+      v[0..2] = 8
+      v == Vector[8, 8, 8, 4, 5, 6]
+    end
+
+    assert_block do
+      v = Vector[1, 3, 4, 5]
+      v[2] = 5
+      v == Vector[1, 3, 5, 5]
+    end
+
+    assert_block do
+      v = Vector[2, 3, 5]
+      v[-2] = 13
+      v == Vector[2, 13, 5]
+    end
+
+    assert_block do
+      v = Vector[4, 8, 9, 11, 30]
+      v[1..-2] = Vector[1, 2, 3]
+      v == Vector[4, 1, 2, 3, 30]
+    end
+
+    assert_raise(IndexError) {Vector[1, 3, 4, 5][5..6] = 17}
+    assert_raise(IndexError) {Vector[1, 3, 4, 5][6] = 17}
+    assert_raise(Matrix::ErrDimensionMismatch) {Vector[1, 3, 4, 5][0..2] = Matrix[[1], [2], [3]]}
+    assert_raise(ArgumentError) {Vector[1, 2, 3, 4, 5, 6][0..2] = Vector[1, 2, 3, 4, 5, 6]}
+    assert_raise(FrozenError) { Vector[7, 8, 9].freeze[0..1] = 5}
+  end
+
+  def test_map!
+    v1 = Vector[1, 2, 3]
+    v2 = Vector[1, 3, 5].freeze
+    v3 = Vector[].freeze
+    assert_equal Vector[1, 4, 9], v1.map!{|e| e ** 2}
+    assert_equal v1, v1.map!{|e| e - 8}
+    assert_raise(FrozenError) { v2.map!{|e| e + 2 }}
+    assert_raise(FrozenError){ v3.map!{} }
+  end
+
+  def test_freeze
+    v = Vector[1,2,3]
+    f = v.freeze
+    assert_equal true, f.frozen?
+    assert v.equal?(f)
+    assert v.equal?(f.freeze)
+    assert_raise(FrozenError){ v[1] = 56 }
+    assert_equal v.dup, v
+  end
+
+  def test_clone
+    a = Vector[4]
+    def a.foo
+      42
+    end
+
+    v = a.clone
+    v[0] = 2
+    assert_equal a, v * 2
+    assert_equal 42, v.foo
+
+    a.freeze
+    v = a.clone
+    assert v.frozen?
+    assert_equal 42, v.foo
+  end
+
+  def test_dup
+    a = Vector[4]
+    def a.foo
+      42
+    end
+    a.freeze
+
+    v = a.dup
+    v[0] = 2
+    assert_equal a, v * 2
+    assert !v.respond_to?(:foo)
+  end
+
   def test_identity
     assert_same @v1, @v1
     assert_not_same @v1, @v2
@@ -223,6 +325,8 @@ def test_angle_with
     assert_in_epsilon(Math::PI/2, Vector[1, 0].angle_with(Vector[0, -1]))
     assert_in_epsilon(Math::PI/4, Vector[2, 2].angle_with(Vector[0, 1]))
     assert_in_delta(0.0, Vector[1, 1].angle_with(Vector[1, 1]), 0.00001)
+    assert_equal(Vector[6, 6].angle_with(Vector[7, 7]), 0.0)
+    assert_equal(Vector[6, 6].angle_with(Vector[-7, -7]), Math::PI)
 
     assert_raise(Vector::ZeroVectorError) { Vector[1, 1].angle_with(Vector[0, 0]) }
     assert_raise(Vector::ZeroVectorError) { Vector[0, 0].angle_with(Vector[1, 1]) }