diff --git a/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/matricesCodegen.kt b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/matricesCodegen.kt
index 0dd880f7b..fd00a0255 100644
--- a/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/matricesCodegen.kt
+++ b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/matricesCodegen.kt
@@ -14,38 +14,53 @@ private fun KtPsiFactory.createMatrixClass(
     kotlinTypeName: String,
     kotlinTypeAlias: String = kotlinTypeName
 ) {
+    fun fn(pattern: String) = fn(pattern, cTypeName)
     val className = "Gsl${kotlinTypeAlias}Matrix"
     val structName = sn("gsl_matrixR", cTypeName)
 
     @Language("kotlin") val text = """internal class $className(
-    override val nativeHandle: CPointer<$structName>,
-    scope: DeferScope
+    override val rawNativeHandle: CPointer<$structName>,
+    scope: AutofreeScope,
 ) : GslMatrix<$kotlinTypeName, $structName>(scope) {
     override val rowNum: Int
-        get() = nativeHandleChecked().pointed.size1.toInt()
+        get() = nativeHandle.pointed.size1.toInt()
 
     override val colNum: Int
-        get() = nativeHandleChecked().pointed.size2.toInt()
+        get() = nativeHandle.pointed.size2.toInt()
+
+    override val rows: Buffer<Buffer<$kotlinTypeName>>
+        get() = VirtualBuffer(rowNum) { r ->
+            Gsl${kotlinTypeAlias}Vector(
+                ${fn("gsl_matrixRrow")}(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
+
+    override val columns: Buffer<Buffer<$kotlinTypeName>>
+        get() = VirtualBuffer(rowNum) { c ->
+            Gsl${kotlinTypeAlias}Vector(
+                ${fn("gsl_matrixRcolumn")}(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
 
     override operator fun get(i: Int, j: Int): $kotlinTypeName = 
-        ${fn("gsl_matrixRget", cTypeName)}(nativeHandleChecked(), i.toULong(), j.toULong())
+        ${fn("gsl_matrixRget")}(nativeHandle, i.toULong(), j.toULong())
 
     override operator fun set(i: Int, j: Int, value: ${kotlinTypeName}): Unit =
-        ${fn("gsl_matrixRset", cTypeName)}(nativeHandleChecked(), i.toULong(), j.toULong(), value)
+        ${fn("gsl_matrixRset")}(nativeHandle, i.toULong(), j.toULong(), value)
 
     override fun copy(): $className {
-        val new = requireNotNull(${fn("gsl_matrixRalloc", cTypeName)}(rowNum.toULong(), colNum.toULong()))
-        ${fn("gsl_matrixRmemcpy", cTypeName)}(new, nativeHandleChecked())
+        val new = requireNotNull(${fn("gsl_matrixRalloc")}(rowNum.toULong(), colNum.toULong()))
+        ${fn("gsl_matrixRmemcpy")}(new, nativeHandle)
         return $className(new, scope)
     }
 
-    override fun close(): Unit = ${fn("gsl_matrixRfree", cTypeName)}(nativeHandleChecked())
+    override fun close(): Unit = ${fn("gsl_matrixRfree")}(nativeHandle)
 
     override fun equals(other: Any?): Boolean {
         if (other is $className)
-            return ${
-        fn("gsl_matrixRequal", cTypeName)
-    }(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return ${fn("gsl_matrixRequal")}(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
@@ -64,7 +79,7 @@ fun matricesCodegen(outputFile: String, project: Project = createProject()) {
             f += createNewLine(2)
             f += createImportDirective(ImportPath.fromString("kotlinx.cinterop.*"))
             f += createNewLine(1)
-            f += createImportDirective(ImportPath.fromString("kscience.kmath.linear.*"))
+            f += createImportDirective(ImportPath.fromString("kscience.kmath.structures.*"))
             f += createNewLine(1)
             f += createImportDirective(ImportPath.fromString("org.gnu.gsl.*"))
             f += createNewLine(2)
diff --git a/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt
index 87acd258c..513d05727 100644
--- a/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt
+++ b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt
@@ -14,38 +14,36 @@ private fun KtPsiFactory.createVectorClass(
     kotlinTypeName: String,
     kotlinTypeAlias: String = kotlinTypeName
 ) {
+    fun fn(pattern: String) = fn(pattern, cTypeName)
     val className = "Gsl${kotlinTypeAlias}Vector"
     val structName = sn("gsl_vectorR", cTypeName)
 
     @Language("kotlin") val text =
-        """internal class $className(override val nativeHandle: CPointer<$structName>, scope: DeferScope) : 
+        """internal class $className(override val rawNativeHandle: CPointer<$structName>, scope: AutofreeScope) : 
     GslVector<$kotlinTypeName, $structName>(scope) {
     override val size: Int
-        get() = nativeHandleChecked().pointed.size.toInt()
+        get() = nativeHandle.pointed.size.toInt()
 
     override operator fun get(index: Int): $kotlinTypeName = 
-        ${fn("gsl_vectorRget", cTypeName)}(nativeHandleChecked(), index.toULong())
+        ${fn("gsl_vectorRget")}(nativeHandle, index.toULong())
 
     override operator fun set(index: Int, value: $kotlinTypeName): Unit = 
-        ${fn("gsl_vectorRset", cTypeName)}(nativeHandleChecked(), index.toULong(), value)
+        ${fn("gsl_vectorRset")}(nativeHandle, index.toULong(), value)
 
     override fun copy(): $className {
-        val new = requireNotNull(${fn("gsl_vectorRalloc", cTypeName)}(size.toULong()))
-        ${fn("gsl_vectorRmemcpy", cTypeName)}(new, nativeHandleChecked())
+        val new = requireNotNull(${fn("gsl_vectorRalloc")}(size.toULong()))
+        ${fn("gsl_vectorRmemcpy")}(new, nativeHandle)
         return ${className}(new, scope)
     }
 
     override fun equals(other: Any?): Boolean {
         if (other is $className)
-            return ${
-            fn("gsl_vectorRequal",
-                cTypeName)
-        }(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return ${fn("gsl_vectorRequal")}(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 
-    override fun close(): Unit = ${fn("gsl_vectorRfree", cTypeName)}(nativeHandleChecked())
+    override fun close(): Unit = ${fn("gsl_vectorRfree")}(nativeHandle)
 }"""
 
     f += createClass(text)
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
index 40c7cd5a4..ad1674c43 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
@@ -1,5 +1,6 @@
 package kscience.kmath.linear
 
+import kscience.kmath.misc.UnstableKMathAPI
 import kscience.kmath.operations.Ring
 import kscience.kmath.operations.SpaceOperations
 import kscience.kmath.operations.invoke
@@ -8,18 +9,22 @@ import kscience.kmath.structures.Buffer
 import kscience.kmath.structures.BufferFactory
 import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.asSequence
+import kotlin.reflect.KClass
 
 /**
- * Basic operations on matrices. Operates on [Matrix]
+ * Basic operations on matrices. Operates on [Matrix].
+ *
+ * @param T the type of items in the matrices.
+ * @param M the type of operated matrices.
  */
 public interface MatrixContext<T : Any, out M : Matrix<T>> : SpaceOperations<Matrix<T>> {
     /**
-     * Produce a matrix with this context and given dimensions
+     * Produces a matrix with this context and given dimensions.
      */
     public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): M
 
     /**
-     * Produce a point compatible with matrix space (and possibly optimized for it)
+     * Produces a point compatible with matrix space (and possibly optimized for it).
      */
     public fun point(size: Int, initializer: (Int) -> T): Point<T> = Buffer.boxing(size, initializer)
 
@@ -66,8 +71,19 @@ public interface MatrixContext<T : Any, out M : Matrix<T>> : SpaceOperations<Mat
      */
     public operator fun T.times(m: Matrix<T>): M = m * this
 
-    public companion object {
+    /**
+     * Gets a feature from the matrix. This function may return some additional features than
+     * [kscience.kmath.structures.NDStructure.getFeature].
+     *
+     * @param F the type of feature.
+     * @param m the matrix.
+     * @param type the [KClass] instance of [F].
+     * @return a feature object or `null` if it isn't present.
+     */
+    @UnstableKMathAPI
+    public fun <F : Any> getFeature(m: Matrix<T>, type: KClass<F>): F? = m.getFeature(type)
 
+    public companion object {
         /**
          * A structured matrix with custom buffer
          */
@@ -84,9 +100,31 @@ public interface MatrixContext<T : Any, out M : Matrix<T>> : SpaceOperations<Mat
     }
 }
 
+/**
+ * Gets a feature from the matrix. This function may return some additional features than
+ * [kscience.kmath.structures.NDStructure.getFeature].
+ *
+ * @param T the type of items in the matrices.
+ * @param M the type of operated matrices.
+ * @param F the type of feature.
+ * @receiver the [MatrixContext] of [T].
+ * @param m the matrix.
+ * @return a feature object or `null` if it isn't present.
+ */
+@UnstableKMathAPI
+public inline fun <T : Any, reified F : Any> MatrixContext<T, *>.getFeature(m: Matrix<T>): F? =
+    getFeature(m, F::class)
+
+/**
+ * Partial implementation of [MatrixContext] for matrices of [Ring].
+ *
+ * @param T the type of items in the matrices.
+ * @param R the type of ring of matrix elements.
+ * @param M the type of operated matrices.
+ */
 public interface GenericMatrixContext<T : Any, R : Ring<T>, out M : Matrix<T>> : MatrixContext<T, M> {
     /**
-     * The ring context for matrix elements
+     * The ring instance of [T].
      */
     public val elementContext: R
 
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt
index bfec6f871..ff19a4103 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt
@@ -241,11 +241,22 @@ public class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> {
     override fun copy(): MutableBuffer<T> = ArrayBuffer(array.copyOf())
 }
 
+
 /**
  * Returns an [ArrayBuffer] that wraps the original array.
  */
 public fun <T> Array<T>.asBuffer(): ArrayBuffer<T> = ArrayBuffer(this)
 
+/**
+ * Creates a new [ArrayBuffer] with the specified [size], where each element is calculated by calling the specified
+ * [init] function.
+ *
+ * The function [init] is called for each array element sequentially starting from the first one.
+ * It should return the value for an array element given its index.
+ */
+public inline fun <reified T> ArrayBuffer(size: Int, init: (Int) -> T): ArrayBuffer<T> =
+    Array(size) { i -> init(i) }.asBuffer()
+
 /**
  * Immutable wrapper for [MutableBuffer].
  *
diff --git a/kmath-gsl/src/nativeInterop/cinterop/libgsl.def b/kmath-gsl/src/nativeInterop/cinterop/libgsl.def
index bd3a861f0..c50dd41cd 100644
--- a/kmath-gsl/src/nativeInterop/cinterop/libgsl.def
+++ b/kmath-gsl/src/nativeInterop/cinterop/libgsl.def
@@ -1,5 +1,5 @@
 package=org.gnu.gsl
-headers=gsl/gsl_blas.h gsl/gsl_linalg.h
+headers=gsl/gsl_blas.h gsl/gsl_linalg.h gsl/gsl_permute_matrix.h
 compilerOpts.linux=-I/usr/include
 compilerOpts.osx=-I/usr/local -I/usr/local/include
 linkerOpts.linux=-L/usr/lib64 -L/usr/lib/x86_64-linux-gnu -lgsl
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslComplex.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslComplex.kt
index 0233ed8bc..6881bb0e4 100644
--- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslComplex.kt
+++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslComplex.kt
@@ -2,6 +2,8 @@ package kscience.kmath.gsl
 
 import kotlinx.cinterop.*
 import kscience.kmath.operations.Complex
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.VirtualBuffer
 import org.gnu.gsl.*
 
 internal fun CValue<gsl_complex>.toKMath(): Complex = useContents { Complex(dat[0], dat[1]) }
@@ -11,58 +13,74 @@ internal fun Complex.toGsl(): CValue<gsl_complex> = cValue {
     dat[1] = im
 }
 
-internal class GslComplexMatrix(override val nativeHandle: CPointer<gsl_matrix_complex>, scope: DeferScope) :
+internal class GslComplexMatrix(override val rawNativeHandle: CPointer<gsl_matrix_complex>, scope: AutofreeScope) :
     GslMatrix<Complex, gsl_matrix_complex>(scope) {
     override val rowNum: Int
-        get() = nativeHandleChecked().pointed.size1.toInt()
+        get() = nativeHandle.pointed.size1.toInt()
 
     override val colNum: Int
-        get() = nativeHandleChecked().pointed.size2.toInt()
+        get() = nativeHandle.pointed.size2.toInt()
+
+    override val rows: Buffer<Buffer<Complex>>
+        get() = VirtualBuffer(rowNum) { r ->
+            GslComplexVector(
+                gsl_matrix_complex_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
+
+    override val columns: Buffer<Buffer<Complex>>
+        get() = VirtualBuffer(rowNum) { c ->
+            GslComplexVector(
+                gsl_matrix_complex_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
 
     override operator fun get(i: Int, j: Int): Complex =
-        gsl_matrix_complex_get(nativeHandleChecked(), i.toULong(), j.toULong()).toKMath()
+        gsl_matrix_complex_get(nativeHandle, i.toULong(), j.toULong()).toKMath()
 
     override operator fun set(i: Int, j: Int, value: Complex): Unit =
-        gsl_matrix_complex_set(nativeHandleChecked(), i.toULong(), j.toULong(), value.toGsl())
+        gsl_matrix_complex_set(nativeHandle, i.toULong(), j.toULong(), value.toGsl())
 
     override fun copy(): GslComplexMatrix {
         val new = requireNotNull(gsl_matrix_complex_alloc(rowNum.toULong(), colNum.toULong()))
-        gsl_matrix_complex_memcpy(new, nativeHandleChecked())
+        gsl_matrix_complex_memcpy(new, nativeHandle)
         return GslComplexMatrix(new, scope)
     }
 
-    override fun close(): Unit = gsl_matrix_complex_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_matrix_complex_free(nativeHandle)
 
     override fun equals(other: Any?): Boolean {
         if (other is GslComplexMatrix)
-            return gsl_matrix_complex_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_matrix_complex_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 }
 
-internal class GslComplexVector(override val nativeHandle: CPointer<gsl_vector_complex>, scope: DeferScope) :
+internal class GslComplexVector(override val rawNativeHandle: CPointer<gsl_vector_complex>, scope: AutofreeScope) :
     GslVector<Complex, gsl_vector_complex>(scope) {
     override val size: Int
-        get() = nativeHandleChecked().pointed.size.toInt()
+        get() = nativeHandle.pointed.size.toInt()
 
-    override fun get(index: Int): Complex = gsl_vector_complex_get(nativeHandleChecked(), index.toULong()).toKMath()
+    override fun get(index: Int): Complex = gsl_vector_complex_get(nativeHandle, index.toULong()).toKMath()
 
     override fun set(index: Int, value: Complex): Unit =
-        gsl_vector_complex_set(nativeHandleChecked(), index.toULong(), value.toGsl())
+        gsl_vector_complex_set(nativeHandle, index.toULong(), value.toGsl())
 
     override fun copy(): GslComplexVector {
         val new = requireNotNull(gsl_vector_complex_alloc(size.toULong()))
-        gsl_vector_complex_memcpy(new, nativeHandleChecked())
+        gsl_vector_complex_memcpy(new, nativeHandle)
         return GslComplexVector(new, scope)
     }
 
     override fun equals(other: Any?): Boolean {
         if (other is GslComplexVector)
-            return gsl_vector_complex_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_vector_complex_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 
-    override fun close(): Unit = gsl_vector_complex_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_vector_complex_free(nativeHandle)
 }
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrix.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrix.kt
index 13283946e..c28236a67 100644
--- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrix.kt
+++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrix.kt
@@ -1,33 +1,40 @@
 package kscience.kmath.gsl
 
+import kotlinx.cinterop.AutofreeScope
 import kotlinx.cinterop.CStructVar
-import kotlinx.cinterop.DeferScope
 import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.NDStructure
+import kscience.kmath.structures.asSequence
 
 /**
  * Wraps gsl_matrix_* objects from GSL.
  */
-public abstract class GslMatrix<T : Any, H : CStructVar> internal constructor(scope: DeferScope) :
+public abstract class GslMatrix<T : Any, H : CStructVar> internal constructor(scope: AutofreeScope) :
     GslMemoryHolder<H>(scope), Matrix<T> {
     internal abstract operator fun set(i: Int, j: Int, value: T)
     internal abstract fun copy(): GslMatrix<T, H>
 
     public override fun equals(other: Any?): Boolean {
-        return NDStructure.equals(this, other as? NDStructure<*> ?: return false)
+        return NDStructure.contentEquals(this, other as? NDStructure<*> ?: return false)
     }
 
     public override fun hashCode(): Int {
         var ret = 7
-        val nRows = rowNum
-        val nCols = colNum
-        ret = ret * 31 + nRows
-        ret = ret * 31 + nCols
+        ret = ret * 31 + rowNum
+        ret = ret * 31 + colNum
 
-        for (row in 0 until nRows)
-            for (col in 0 until nCols)
+        for (row in 0 until rowNum)
+            for (col in 0 until colNum)
                 ret = ret * 31 + (11 * (row + 1) + 17 * (col + 1)) * this[row, col].hashCode()
 
         return ret
     }
+
+    public override fun toString(): String = if (rowNum <= 5 && colNum <= 5)
+        "Matrix(rowsNum = $rowNum, colNum = $colNum)\n" +
+                rows.asSequence().joinToString(prefix = "(", postfix = ")", separator = "\n ") { buffer ->
+                    buffer.asSequence().joinToString(separator = "\t") { it.toString() }
+                }
+    else
+        "Matrix(rowsNum = $rowNum, colNum = $colNum)"
 }
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt
index 95437c8c7..c32d852ba 100644
--- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt
+++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt
@@ -1,31 +1,33 @@
 package kscience.kmath.gsl
 
-import kotlinx.cinterop.CStructVar
-import kotlinx.cinterop.DeferScope
-import kotlinx.cinterop.memScoped
-import kotlinx.cinterop.pointed
-import kscience.kmath.linear.MatrixContext
-import kscience.kmath.linear.Point
+import kotlinx.cinterop.*
+import kscience.kmath.linear.*
+import kscience.kmath.misc.UnstableKMathAPI
 import kscience.kmath.operations.Complex
 import kscience.kmath.operations.ComplexField
 import kscience.kmath.operations.toComplex
 import kscience.kmath.structures.Matrix
+import kscience.kmath.structures.toList
 import org.gnu.gsl.*
+import kotlin.reflect.KClass
+import kotlin.reflect.cast
 
 internal inline fun <T : Any, H : CStructVar> GslMatrix<T, H>.fill(initializer: (Int, Int) -> T): GslMatrix<T, H> =
     apply {
-        (0 until rowNum).forEach { row -> (0 until colNum).forEach { col -> this[row, col] = initializer(row, col) } }
+        for (row in 0 until rowNum) {
+            for (col in 0 until colNum) this[row, col] = initializer(row, col)
+        }
     }
 
 internal inline fun <T : Any, H : CStructVar> GslVector<T, H>.fill(initializer: (Int) -> T): GslVector<T, H> =
-    apply { (0 until size).forEach { index -> this[index] = initializer(index) } }
+    apply {
+        for (index in 0 until size) this[index] = initializer(index)
+    }
 
 /**
  * Represents matrix context implementing where all the operations are delegated to GSL.
  */
-public abstract class GslMatrixContext<T : Any, H1 : CStructVar, H2 : CStructVar> internal constructor(
-    internal val scope: DeferScope,
-) : MatrixContext<T, GslMatrix<T, H1>> {
+public abstract class GslMatrixContext<T : Any, H1 : CStructVar, H2 : CStructVar> : MatrixContext<T, GslMatrix<T, H1>> {
     init {
         ensureHasGslErrorHandler()
     }
@@ -52,19 +54,22 @@ public abstract class GslMatrixContext<T : Any, H1 : CStructVar, H2 : CStructVar
     public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): GslMatrix<T, H1> =
         produceDirtyMatrix(rows, columns).fill(initializer)
 
-    public override fun point(size: Int, initializer: (Int) -> T): Point<T> =
+    public override fun point(size: Int, initializer: (Int) -> T): GslVector<T, H2> =
         produceDirtyVector(size).fill(initializer)
 }
 
 /**
  * Represents [Double] matrix context implementing where all the operations are delegated to GSL.
  */
-public class GslRealMatrixContext(scope: DeferScope) : GslMatrixContext<Double, gsl_matrix, gsl_vector>(scope) {
-    override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<Double, gsl_matrix> =
-        GslRealMatrix(nativeHandle = requireNotNull(gsl_matrix_alloc(rows.toULong(), columns.toULong())), scope = scope)
+public class GslRealMatrixContext(internal val scope: AutofreeScope) :
+    GslMatrixContext<Double, gsl_matrix, gsl_vector>() {
+    override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<Double, gsl_matrix> = GslRealMatrix(
+        rawNativeHandle = requireNotNull(gsl_matrix_alloc(rows.toULong(), columns.toULong())),
+        scope = scope,
+    )
 
     override fun produceDirtyVector(size: Int): GslVector<Double, gsl_vector> =
-        GslRealVector(nativeHandle = requireNotNull(gsl_vector_alloc(size.toULong())), scope = scope)
+        GslRealVector(rawNativeHandle = requireNotNull(gsl_vector_alloc(size.toULong())), scope = scope)
 
     public override fun Matrix<Double>.dot(other: Matrix<Double>): GslMatrix<Double, gsl_matrix> {
         val x = toGsl().nativeHandle
@@ -105,6 +110,58 @@ public class GslRealMatrixContext(scope: DeferScope) : GslMatrixContext<Double,
         gsl_matrix_sub(g1.nativeHandle, b.toGsl().nativeHandle)
         return g1
     }
+
+    @Suppress("IMPLICIT_CAST_TO_ANY")
+    @UnstableKMathAPI
+    public override fun <F : Any> getFeature(m: Matrix<Double>, type: KClass<F>): F? = when (type) {
+        LupDecompositionFeature::class -> object : LupDecompositionFeature<Double> {
+            private val lup by lazy {
+                val lu = m.toGsl().copy()
+                val n = lu.rowNum
+                val perm = gsl_permutation_alloc(n.toULong())
+
+                memScoped {
+                    val i = alloc<IntVar>()
+                    gsl_linalg_LU_decomp(lu.nativeHandle, perm, i.ptr)
+                }
+
+                val one = produce(n, n) { i, j -> if (i == j) 1.0 else 0.0 }
+                val p = produce(n, n) { _, _ -> 0.0 }
+
+                for (j in 0 until perm!!.pointed.size.toInt()) {
+                    println("7 $j")
+                    val k = gsl_permutation_get(perm, j.toULong()).toInt()
+                    val col = one.columns[k]
+                    println(p.toGsl())
+                    println(col.toList())
+                    gsl_matrix_set_col(p.nativeHandle, j.toULong(), col.toGsl().nativeHandle)
+                }
+
+                gsl_permutation_free(perm)
+                lu to p
+            }
+
+            override val l by lazy {
+                VirtualMatrix(lup.first.shape[0], lup.first.shape[1]) { i, j ->
+                    when {
+                        j < i -> lup.first[i, j]
+                        j == i -> 1.0
+                        else -> 0.0
+                    }
+                } + LFeature
+            }
+
+            override val u by lazy {
+                VirtualMatrix(
+                    lup.first.shape[0],
+                    lup.first.shape[1],
+                ) { i, j -> if (j >= i) lup.first[i, j] else 0.0 } + UFeature
+            }
+
+            override val p by lazy(lup::second)
+        }
+        else -> super.getFeature(m, type)
+    }?.let(type::cast)
 }
 
 /**
@@ -116,22 +173,22 @@ public fun <R> GslRealMatrixContext(block: GslRealMatrixContext.() -> R): R =
 /**
  * Represents [Float] matrix context implementing where all the operations are delegated to GSL.
  */
-public class GslFloatMatrixContext(scope: DeferScope) :
-    GslMatrixContext<Float, gsl_matrix_float, gsl_vector_float>(scope) {
+public class GslFloatMatrixContext(internal val scope: AutofreeScope) :
+    GslMatrixContext<Float, gsl_matrix_float, gsl_vector_float>() {
     override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<Float, gsl_matrix_float> = GslFloatMatrix(
-        nativeHandle = requireNotNull(gsl_matrix_float_alloc(rows.toULong(), columns.toULong())),
+        rawNativeHandle = requireNotNull(gsl_matrix_float_alloc(rows.toULong(), columns.toULong())),
         scope = scope,
     )
 
     override fun produceDirtyVector(size: Int): GslVector<Float, gsl_vector_float> =
-        GslFloatVector(nativeHandle = requireNotNull(value = gsl_vector_float_alloc(size.toULong())), scope = scope)
+        GslFloatVector(rawNativeHandle = requireNotNull(value = gsl_vector_float_alloc(size.toULong())), scope = scope)
 
     public override fun Matrix<Float>.dot(other: Matrix<Float>): GslMatrix<Float, gsl_matrix_float> {
         val x = toGsl().nativeHandle
         val a = other.toGsl().nativeHandle
         val result = requireNotNull(gsl_matrix_float_calloc(a.pointed.size1, a.pointed.size2))
         gsl_blas_sgemm(CblasNoTrans, CblasNoTrans, 1f, x, a, 1f, result)
-        return GslFloatMatrix(nativeHandle = result, scope = scope)
+        return GslFloatMatrix(rawNativeHandle = result, scope = scope)
     }
 
     public override fun Matrix<Float>.dot(vector: Point<Float>): GslVector<Float, gsl_vector_float> {
@@ -139,7 +196,7 @@ public class GslFloatMatrixContext(scope: DeferScope) :
         val a = vector.toGsl().nativeHandle
         val result = requireNotNull(gsl_vector_float_calloc(a.pointed.size))
         gsl_blas_sgemv(CblasNoTrans, 1f, x, a, 1f, result)
-        return GslFloatVector(nativeHandle = result, scope = scope)
+        return GslFloatVector(rawNativeHandle = result, scope = scope)
     }
 
     public override fun Matrix<Float>.times(value: Float): GslMatrix<Float, gsl_matrix_float> {
@@ -176,22 +233,22 @@ public fun <R> GslFloatMatrixContext(block: GslFloatMatrixContext.() -> R): R =
 /**
  * Represents [Complex] matrix context implementing where all the operations are delegated to GSL.
  */
-public class GslComplexMatrixContext(scope: DeferScope) :
-    GslMatrixContext<Complex, gsl_matrix_complex, gsl_vector_complex>(scope) {
+public class GslComplexMatrixContext(internal val scope: AutofreeScope) :
+    GslMatrixContext<Complex, gsl_matrix_complex, gsl_vector_complex>() {
     override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<Complex, gsl_matrix_complex> = GslComplexMatrix(
-        nativeHandle = requireNotNull(gsl_matrix_complex_alloc(rows.toULong(), columns.toULong())),
+        rawNativeHandle = requireNotNull(gsl_matrix_complex_alloc(rows.toULong(), columns.toULong())),
         scope = scope,
     )
 
     override fun produceDirtyVector(size: Int): GslVector<Complex, gsl_vector_complex> =
-        GslComplexVector(nativeHandle = requireNotNull(gsl_vector_complex_alloc(size.toULong())), scope = scope)
+        GslComplexVector(rawNativeHandle = requireNotNull(gsl_vector_complex_alloc(size.toULong())), scope = scope)
 
     public override fun Matrix<Complex>.dot(other: Matrix<Complex>): GslMatrix<Complex, gsl_matrix_complex> {
         val x = toGsl().nativeHandle
         val a = other.toGsl().nativeHandle
         val result = requireNotNull(gsl_matrix_complex_calloc(a.pointed.size1, a.pointed.size2))
         gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, ComplexField.one.toGsl(), x, a, ComplexField.one.toGsl(), result)
-        return GslComplexMatrix(nativeHandle = result, scope = scope)
+        return GslComplexMatrix(rawNativeHandle = result, scope = scope)
     }
 
     public override fun Matrix<Complex>.dot(vector: Point<Complex>): GslVector<Complex, gsl_vector_complex> {
@@ -225,6 +282,58 @@ public class GslComplexMatrixContext(scope: DeferScope) :
         gsl_matrix_complex_sub(g1.nativeHandle, b.toGsl().nativeHandle)
         return g1
     }
+
+    @Suppress("IMPLICIT_CAST_TO_ANY")
+    @UnstableKMathAPI
+    public override fun <F : Any> getFeature(m: Matrix<Complex>, type: KClass<F>): F? = when (type) {
+        LupDecompositionFeature::class -> object : LupDecompositionFeature<Complex> {
+            private val lup by lazy {
+                val lu = m.toGsl().copy()
+                val n = lu.rowNum
+                val perm = gsl_permutation_alloc(n.toULong())
+
+                memScoped {
+                    val i = alloc<IntVar>()
+                    gsl_linalg_complex_LU_decomp(lu.nativeHandle, perm, i.ptr)
+                }
+
+                val one = produce(n, n) { i, j -> if (i == j) 1.0.toComplex() else 0.0.toComplex() }
+                val p = produce(n, n) { _, _ -> 0.0.toComplex() }
+
+                for (j in 0 until perm!!.pointed.size.toInt()) {
+                    println("7 $j")
+                    val k = gsl_permutation_get(perm, j.toULong()).toInt()
+                    val col = one.columns[k]
+                    println(p.toGsl())
+                    println(col.toList())
+                    gsl_matrix_complex_set_col(p.nativeHandle, j.toULong(), col.toGsl().nativeHandle)
+                }
+
+                gsl_permutation_free(perm)
+                lu to p
+            }
+
+            override val l by lazy {
+                VirtualMatrix(lup.first.shape[0], lup.first.shape[1]) { i, j ->
+                    when {
+                        j < i -> lup.first[i, j]
+                        j == i -> 1.0.toComplex()
+                        else -> 0.0.toComplex()
+                    }
+                } + LFeature
+            }
+
+            override val u by lazy {
+                VirtualMatrix(
+                    lup.first.shape[0],
+                    lup.first.shape[1],
+                ) { i, j -> if (j >= i) lup.first[i, j] else 0.0.toComplex() } + UFeature
+            }
+
+            override val p by lazy(lup::second)
+        }
+        else -> super.getFeature(m, type)
+    }?.let(type::cast)
 }
 
 /**
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt
index d68a52a62..46d1ba60b 100644
--- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt
+++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt
@@ -1,5 +1,6 @@
 package kscience.kmath.gsl
 
+import kotlinx.cinterop.AutofreeScope
 import kotlinx.cinterop.CPointer
 import kotlinx.cinterop.CStructVar
 import kotlinx.cinterop.DeferScope
@@ -12,23 +13,26 @@ import kotlinx.cinterop.DeferScope
  *
  * @param scope the scope where this object is declared.
  */
-public abstract class GslMemoryHolder<H : CStructVar> internal constructor(internal val scope: DeferScope) {
-    internal abstract val nativeHandle: CPointer<H>
+public abstract class GslMemoryHolder<H : CStructVar> internal constructor(internal val scope: AutofreeScope) {
+    internal abstract val rawNativeHandle: CPointer<H>
     private var isClosed: Boolean = false
 
+    internal val nativeHandle: CPointer<H>
+        get() {
+            check(!isClosed) { "The use of GSL object that is closed." }
+            return rawNativeHandle
+        }
+
+    internal var shouldBeFreed = true
+
     init {
         ensureHasGslErrorHandler()
 
         scope.defer {
-            close()
+            if (shouldBeFreed) close()
             isClosed = true
         }
     }
 
-    internal fun nativeHandleChecked(): CPointer<H> {
-        check(!isClosed) { "The use of GSL object that is closed." }
-        return nativeHandle
-    }
-
     internal abstract fun close()
 }
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslVector.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslVector.kt
index 0033c47a1..8bc1240b4 100644
--- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslVector.kt
+++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslVector.kt
@@ -1,13 +1,13 @@
 package kscience.kmath.gsl
 
+import kotlinx.cinterop.AutofreeScope
 import kotlinx.cinterop.CStructVar
-import kotlinx.cinterop.DeferScope
 import kscience.kmath.linear.Point
 
 /**
  * Wraps gsl_vector_* objects from GSL.
  */
-public abstract class GslVector<T, H : CStructVar> internal constructor(scope: DeferScope) :
+public abstract class GslVector<T, H : CStructVar> internal constructor(scope: AutofreeScope) :
     GslMemoryHolder<H>(scope), Point<T> {
     internal abstract operator fun set(index: Int, value: T)
     internal abstract fun copy(): GslVector<T, H>
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Matrices.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Matrices.kt
index 07bc8af79..1a826d493 100644
--- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Matrices.kt
+++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Matrices.kt
@@ -1,260 +1,388 @@
 package kscience.kmath.gsl
 
 import kotlinx.cinterop.*
-import kscience.kmath.linear.*
+import kscience.kmath.structures.*
 import org.gnu.gsl.*
 
 internal class GslRealMatrix(
-    override val nativeHandle: CPointer<gsl_matrix>,
-    scope: DeferScope
+    override val rawNativeHandle: CPointer<gsl_matrix>,
+    scope: AutofreeScope,
 ) : GslMatrix<Double, gsl_matrix>(scope) {
     override val rowNum: Int
-        get() = nativeHandleChecked().pointed.size1.toInt()
+        get() = nativeHandle.pointed.size1.toInt()
 
     override val colNum: Int
-        get() = nativeHandleChecked().pointed.size2.toInt()
+        get() = nativeHandle.pointed.size2.toInt()
+
+    override val rows: Buffer<Buffer<Double>>
+        get() = VirtualBuffer(rowNum) { r ->
+            GslRealVector(
+                gsl_matrix_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
+
+    override val columns: Buffer<Buffer<Double>>
+        get() = VirtualBuffer(rowNum) { c ->
+            GslRealVector(
+                gsl_matrix_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
 
     override operator fun get(i: Int, j: Int): Double = 
-        gsl_matrix_get(nativeHandleChecked(), i.toULong(), j.toULong())
+        gsl_matrix_get(nativeHandle, i.toULong(), j.toULong())
 
     override operator fun set(i: Int, j: Int, value: Double): Unit =
-        gsl_matrix_set(nativeHandleChecked(), i.toULong(), j.toULong(), value)
+        gsl_matrix_set(nativeHandle, i.toULong(), j.toULong(), value)
 
     override fun copy(): GslRealMatrix {
         val new = requireNotNull(gsl_matrix_alloc(rowNum.toULong(), colNum.toULong()))
-        gsl_matrix_memcpy(new, nativeHandleChecked())
+        gsl_matrix_memcpy(new, nativeHandle)
         return GslRealMatrix(new, scope)
     }
 
-    override fun close(): Unit = gsl_matrix_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_matrix_free(nativeHandle)
 
     override fun equals(other: Any?): Boolean {
         if (other is GslRealMatrix)
-            return gsl_matrix_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_matrix_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 }
 
 internal class GslFloatMatrix(
-    override val nativeHandle: CPointer<gsl_matrix_float>,
-    scope: DeferScope
+    override val rawNativeHandle: CPointer<gsl_matrix_float>,
+    scope: AutofreeScope,
 ) : GslMatrix<Float, gsl_matrix_float>(scope) {
     override val rowNum: Int
-        get() = nativeHandleChecked().pointed.size1.toInt()
+        get() = nativeHandle.pointed.size1.toInt()
 
     override val colNum: Int
-        get() = nativeHandleChecked().pointed.size2.toInt()
+        get() = nativeHandle.pointed.size2.toInt()
+
+    override val rows: Buffer<Buffer<Float>>
+        get() = VirtualBuffer(rowNum) { r ->
+            GslFloatVector(
+                gsl_matrix_float_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
+
+    override val columns: Buffer<Buffer<Float>>
+        get() = VirtualBuffer(rowNum) { c ->
+            GslFloatVector(
+                gsl_matrix_float_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
 
     override operator fun get(i: Int, j: Int): Float = 
-        gsl_matrix_float_get(nativeHandleChecked(), i.toULong(), j.toULong())
+        gsl_matrix_float_get(nativeHandle, i.toULong(), j.toULong())
 
     override operator fun set(i: Int, j: Int, value: Float): Unit =
-        gsl_matrix_float_set(nativeHandleChecked(), i.toULong(), j.toULong(), value)
+        gsl_matrix_float_set(nativeHandle, i.toULong(), j.toULong(), value)
 
     override fun copy(): GslFloatMatrix {
         val new = requireNotNull(gsl_matrix_float_alloc(rowNum.toULong(), colNum.toULong()))
-        gsl_matrix_float_memcpy(new, nativeHandleChecked())
+        gsl_matrix_float_memcpy(new, nativeHandle)
         return GslFloatMatrix(new, scope)
     }
 
-    override fun close(): Unit = gsl_matrix_float_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_matrix_float_free(nativeHandle)
 
     override fun equals(other: Any?): Boolean {
         if (other is GslFloatMatrix)
-            return gsl_matrix_float_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_matrix_float_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 }
 
 internal class GslShortMatrix(
-    override val nativeHandle: CPointer<gsl_matrix_short>,
-    scope: DeferScope
+    override val rawNativeHandle: CPointer<gsl_matrix_short>,
+    scope: AutofreeScope,
 ) : GslMatrix<Short, gsl_matrix_short>(scope) {
     override val rowNum: Int
-        get() = nativeHandleChecked().pointed.size1.toInt()
+        get() = nativeHandle.pointed.size1.toInt()
 
     override val colNum: Int
-        get() = nativeHandleChecked().pointed.size2.toInt()
+        get() = nativeHandle.pointed.size2.toInt()
+
+    override val rows: Buffer<Buffer<Short>>
+        get() = VirtualBuffer(rowNum) { r ->
+            GslShortVector(
+                gsl_matrix_short_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
+
+    override val columns: Buffer<Buffer<Short>>
+        get() = VirtualBuffer(rowNum) { c ->
+            GslShortVector(
+                gsl_matrix_short_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
 
     override operator fun get(i: Int, j: Int): Short = 
-        gsl_matrix_short_get(nativeHandleChecked(), i.toULong(), j.toULong())
+        gsl_matrix_short_get(nativeHandle, i.toULong(), j.toULong())
 
     override operator fun set(i: Int, j: Int, value: Short): Unit =
-        gsl_matrix_short_set(nativeHandleChecked(), i.toULong(), j.toULong(), value)
+        gsl_matrix_short_set(nativeHandle, i.toULong(), j.toULong(), value)
 
     override fun copy(): GslShortMatrix {
         val new = requireNotNull(gsl_matrix_short_alloc(rowNum.toULong(), colNum.toULong()))
-        gsl_matrix_short_memcpy(new, nativeHandleChecked())
+        gsl_matrix_short_memcpy(new, nativeHandle)
         return GslShortMatrix(new, scope)
     }
 
-    override fun close(): Unit = gsl_matrix_short_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_matrix_short_free(nativeHandle)
 
     override fun equals(other: Any?): Boolean {
         if (other is GslShortMatrix)
-            return gsl_matrix_short_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_matrix_short_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 }
 
 internal class GslUShortMatrix(
-    override val nativeHandle: CPointer<gsl_matrix_ushort>,
-    scope: DeferScope
+    override val rawNativeHandle: CPointer<gsl_matrix_ushort>,
+    scope: AutofreeScope,
 ) : GslMatrix<UShort, gsl_matrix_ushort>(scope) {
     override val rowNum: Int
-        get() = nativeHandleChecked().pointed.size1.toInt()
+        get() = nativeHandle.pointed.size1.toInt()
 
     override val colNum: Int
-        get() = nativeHandleChecked().pointed.size2.toInt()
+        get() = nativeHandle.pointed.size2.toInt()
+
+    override val rows: Buffer<Buffer<UShort>>
+        get() = VirtualBuffer(rowNum) { r ->
+            GslUShortVector(
+                gsl_matrix_ushort_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
+
+    override val columns: Buffer<Buffer<UShort>>
+        get() = VirtualBuffer(rowNum) { c ->
+            GslUShortVector(
+                gsl_matrix_ushort_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
 
     override operator fun get(i: Int, j: Int): UShort = 
-        gsl_matrix_ushort_get(nativeHandleChecked(), i.toULong(), j.toULong())
+        gsl_matrix_ushort_get(nativeHandle, i.toULong(), j.toULong())
 
     override operator fun set(i: Int, j: Int, value: UShort): Unit =
-        gsl_matrix_ushort_set(nativeHandleChecked(), i.toULong(), j.toULong(), value)
+        gsl_matrix_ushort_set(nativeHandle, i.toULong(), j.toULong(), value)
 
     override fun copy(): GslUShortMatrix {
         val new = requireNotNull(gsl_matrix_ushort_alloc(rowNum.toULong(), colNum.toULong()))
-        gsl_matrix_ushort_memcpy(new, nativeHandleChecked())
+        gsl_matrix_ushort_memcpy(new, nativeHandle)
         return GslUShortMatrix(new, scope)
     }
 
-    override fun close(): Unit = gsl_matrix_ushort_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_matrix_ushort_free(nativeHandle)
 
     override fun equals(other: Any?): Boolean {
         if (other is GslUShortMatrix)
-            return gsl_matrix_ushort_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_matrix_ushort_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 }
 
 internal class GslLongMatrix(
-    override val nativeHandle: CPointer<gsl_matrix_long>,
-    scope: DeferScope
+    override val rawNativeHandle: CPointer<gsl_matrix_long>,
+    scope: AutofreeScope,
 ) : GslMatrix<Long, gsl_matrix_long>(scope) {
     override val rowNum: Int
-        get() = nativeHandleChecked().pointed.size1.toInt()
+        get() = nativeHandle.pointed.size1.toInt()
 
     override val colNum: Int
-        get() = nativeHandleChecked().pointed.size2.toInt()
+        get() = nativeHandle.pointed.size2.toInt()
+
+    override val rows: Buffer<Buffer<Long>>
+        get() = VirtualBuffer(rowNum) { r ->
+            GslLongVector(
+                gsl_matrix_long_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
+
+    override val columns: Buffer<Buffer<Long>>
+        get() = VirtualBuffer(rowNum) { c ->
+            GslLongVector(
+                gsl_matrix_long_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
 
     override operator fun get(i: Int, j: Int): Long = 
-        gsl_matrix_long_get(nativeHandleChecked(), i.toULong(), j.toULong())
+        gsl_matrix_long_get(nativeHandle, i.toULong(), j.toULong())
 
     override operator fun set(i: Int, j: Int, value: Long): Unit =
-        gsl_matrix_long_set(nativeHandleChecked(), i.toULong(), j.toULong(), value)
+        gsl_matrix_long_set(nativeHandle, i.toULong(), j.toULong(), value)
 
     override fun copy(): GslLongMatrix {
         val new = requireNotNull(gsl_matrix_long_alloc(rowNum.toULong(), colNum.toULong()))
-        gsl_matrix_long_memcpy(new, nativeHandleChecked())
+        gsl_matrix_long_memcpy(new, nativeHandle)
         return GslLongMatrix(new, scope)
     }
 
-    override fun close(): Unit = gsl_matrix_long_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_matrix_long_free(nativeHandle)
 
     override fun equals(other: Any?): Boolean {
         if (other is GslLongMatrix)
-            return gsl_matrix_long_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_matrix_long_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 }
 
 internal class GslULongMatrix(
-    override val nativeHandle: CPointer<gsl_matrix_ulong>,
-    scope: DeferScope
+    override val rawNativeHandle: CPointer<gsl_matrix_ulong>,
+    scope: AutofreeScope,
 ) : GslMatrix<ULong, gsl_matrix_ulong>(scope) {
     override val rowNum: Int
-        get() = nativeHandleChecked().pointed.size1.toInt()
+        get() = nativeHandle.pointed.size1.toInt()
 
     override val colNum: Int
-        get() = nativeHandleChecked().pointed.size2.toInt()
+        get() = nativeHandle.pointed.size2.toInt()
+
+    override val rows: Buffer<Buffer<ULong>>
+        get() = VirtualBuffer(rowNum) { r ->
+            GslULongVector(
+                gsl_matrix_ulong_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
+
+    override val columns: Buffer<Buffer<ULong>>
+        get() = VirtualBuffer(rowNum) { c ->
+            GslULongVector(
+                gsl_matrix_ulong_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
 
     override operator fun get(i: Int, j: Int): ULong = 
-        gsl_matrix_ulong_get(nativeHandleChecked(), i.toULong(), j.toULong())
+        gsl_matrix_ulong_get(nativeHandle, i.toULong(), j.toULong())
 
     override operator fun set(i: Int, j: Int, value: ULong): Unit =
-        gsl_matrix_ulong_set(nativeHandleChecked(), i.toULong(), j.toULong(), value)
+        gsl_matrix_ulong_set(nativeHandle, i.toULong(), j.toULong(), value)
 
     override fun copy(): GslULongMatrix {
         val new = requireNotNull(gsl_matrix_ulong_alloc(rowNum.toULong(), colNum.toULong()))
-        gsl_matrix_ulong_memcpy(new, nativeHandleChecked())
+        gsl_matrix_ulong_memcpy(new, nativeHandle)
         return GslULongMatrix(new, scope)
     }
 
-    override fun close(): Unit = gsl_matrix_ulong_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_matrix_ulong_free(nativeHandle)
 
     override fun equals(other: Any?): Boolean {
         if (other is GslULongMatrix)
-            return gsl_matrix_ulong_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_matrix_ulong_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 }
 
 internal class GslIntMatrix(
-    override val nativeHandle: CPointer<gsl_matrix_int>,
-    scope: DeferScope
+    override val rawNativeHandle: CPointer<gsl_matrix_int>,
+    scope: AutofreeScope,
 ) : GslMatrix<Int, gsl_matrix_int>(scope) {
     override val rowNum: Int
-        get() = nativeHandleChecked().pointed.size1.toInt()
+        get() = nativeHandle.pointed.size1.toInt()
 
     override val colNum: Int
-        get() = nativeHandleChecked().pointed.size2.toInt()
+        get() = nativeHandle.pointed.size2.toInt()
+
+    override val rows: Buffer<Buffer<Int>>
+        get() = VirtualBuffer(rowNum) { r ->
+            GslIntVector(
+                gsl_matrix_int_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
+
+    override val columns: Buffer<Buffer<Int>>
+        get() = VirtualBuffer(rowNum) { c ->
+            GslIntVector(
+                gsl_matrix_int_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
 
     override operator fun get(i: Int, j: Int): Int = 
-        gsl_matrix_int_get(nativeHandleChecked(), i.toULong(), j.toULong())
+        gsl_matrix_int_get(nativeHandle, i.toULong(), j.toULong())
 
     override operator fun set(i: Int, j: Int, value: Int): Unit =
-        gsl_matrix_int_set(nativeHandleChecked(), i.toULong(), j.toULong(), value)
+        gsl_matrix_int_set(nativeHandle, i.toULong(), j.toULong(), value)
 
     override fun copy(): GslIntMatrix {
         val new = requireNotNull(gsl_matrix_int_alloc(rowNum.toULong(), colNum.toULong()))
-        gsl_matrix_int_memcpy(new, nativeHandleChecked())
+        gsl_matrix_int_memcpy(new, nativeHandle)
         return GslIntMatrix(new, scope)
     }
 
-    override fun close(): Unit = gsl_matrix_int_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_matrix_int_free(nativeHandle)
 
     override fun equals(other: Any?): Boolean {
         if (other is GslIntMatrix)
-            return gsl_matrix_int_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_matrix_int_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 }
 
 internal class GslUIntMatrix(
-    override val nativeHandle: CPointer<gsl_matrix_uint>,
-    scope: DeferScope
+    override val rawNativeHandle: CPointer<gsl_matrix_uint>,
+    scope: AutofreeScope,
 ) : GslMatrix<UInt, gsl_matrix_uint>(scope) {
     override val rowNum: Int
-        get() = nativeHandleChecked().pointed.size1.toInt()
+        get() = nativeHandle.pointed.size1.toInt()
 
     override val colNum: Int
-        get() = nativeHandleChecked().pointed.size2.toInt()
+        get() = nativeHandle.pointed.size2.toInt()
+
+    override val rows: Buffer<Buffer<UInt>>
+        get() = VirtualBuffer(rowNum) { r ->
+            GslUIntVector(
+                gsl_matrix_uint_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
+
+    override val columns: Buffer<Buffer<UInt>>
+        get() = VirtualBuffer(rowNum) { c ->
+            GslUIntVector(
+                gsl_matrix_uint_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr,
+                scope,
+            ).apply { shouldBeFreed = false }
+        }
 
     override operator fun get(i: Int, j: Int): UInt = 
-        gsl_matrix_uint_get(nativeHandleChecked(), i.toULong(), j.toULong())
+        gsl_matrix_uint_get(nativeHandle, i.toULong(), j.toULong())
 
     override operator fun set(i: Int, j: Int, value: UInt): Unit =
-        gsl_matrix_uint_set(nativeHandleChecked(), i.toULong(), j.toULong(), value)
+        gsl_matrix_uint_set(nativeHandle, i.toULong(), j.toULong(), value)
 
     override fun copy(): GslUIntMatrix {
         val new = requireNotNull(gsl_matrix_uint_alloc(rowNum.toULong(), colNum.toULong()))
-        gsl_matrix_uint_memcpy(new, nativeHandleChecked())
+        gsl_matrix_uint_memcpy(new, nativeHandle)
         return GslUIntMatrix(new, scope)
     }
 
-    override fun close(): Unit = gsl_matrix_uint_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_matrix_uint_free(nativeHandle)
 
     override fun equals(other: Any?): Boolean {
         if (other is GslUIntMatrix)
-            return gsl_matrix_uint_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_matrix_uint_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt
index d7ed4718b..5b5236242 100644
--- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt
+++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt
@@ -3,219 +3,219 @@ package kscience.kmath.gsl
 import kotlinx.cinterop.*
 import org.gnu.gsl.*
 
-internal class GslRealVector(override val nativeHandle: CPointer<gsl_vector>, scope: DeferScope) : 
+internal class GslRealVector(override val rawNativeHandle: CPointer<gsl_vector>, scope: AutofreeScope) : 
     GslVector<Double, gsl_vector>(scope) {
     override val size: Int
-        get() = nativeHandleChecked().pointed.size.toInt()
+        get() = nativeHandle.pointed.size.toInt()
 
     override operator fun get(index: Int): Double = 
-        gsl_vector_get(nativeHandleChecked(), index.toULong())
+        gsl_vector_get(nativeHandle, index.toULong())
 
     override operator fun set(index: Int, value: Double): Unit = 
-        gsl_vector_set(nativeHandleChecked(), index.toULong(), value)
+        gsl_vector_set(nativeHandle, index.toULong(), value)
 
     override fun copy(): GslRealVector {
         val new = requireNotNull(gsl_vector_alloc(size.toULong()))
-        gsl_vector_memcpy(new, nativeHandleChecked())
+        gsl_vector_memcpy(new, nativeHandle)
         return GslRealVector(new, scope)
     }
 
     override fun equals(other: Any?): Boolean {
         if (other is GslRealVector)
-            return gsl_vector_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_vector_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 
-    override fun close(): Unit = gsl_vector_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_vector_free(nativeHandle)
 }
 
-internal class GslFloatVector(override val nativeHandle: CPointer<gsl_vector_float>, scope: DeferScope) : 
+internal class GslFloatVector(override val rawNativeHandle: CPointer<gsl_vector_float>, scope: AutofreeScope) : 
     GslVector<Float, gsl_vector_float>(scope) {
     override val size: Int
-        get() = nativeHandleChecked().pointed.size.toInt()
+        get() = nativeHandle.pointed.size.toInt()
 
     override operator fun get(index: Int): Float = 
-        gsl_vector_float_get(nativeHandleChecked(), index.toULong())
+        gsl_vector_float_get(nativeHandle, index.toULong())
 
     override operator fun set(index: Int, value: Float): Unit = 
-        gsl_vector_float_set(nativeHandleChecked(), index.toULong(), value)
+        gsl_vector_float_set(nativeHandle, index.toULong(), value)
 
     override fun copy(): GslFloatVector {
         val new = requireNotNull(gsl_vector_float_alloc(size.toULong()))
-        gsl_vector_float_memcpy(new, nativeHandleChecked())
+        gsl_vector_float_memcpy(new, nativeHandle)
         return GslFloatVector(new, scope)
     }
 
     override fun equals(other: Any?): Boolean {
         if (other is GslFloatVector)
-            return gsl_vector_float_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_vector_float_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 
-    override fun close(): Unit = gsl_vector_float_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_vector_float_free(nativeHandle)
 }
 
-internal class GslShortVector(override val nativeHandle: CPointer<gsl_vector_short>, scope: DeferScope) : 
+internal class GslShortVector(override val rawNativeHandle: CPointer<gsl_vector_short>, scope: AutofreeScope) : 
     GslVector<Short, gsl_vector_short>(scope) {
     override val size: Int
-        get() = nativeHandleChecked().pointed.size.toInt()
+        get() = nativeHandle.pointed.size.toInt()
 
     override operator fun get(index: Int): Short = 
-        gsl_vector_short_get(nativeHandleChecked(), index.toULong())
+        gsl_vector_short_get(nativeHandle, index.toULong())
 
     override operator fun set(index: Int, value: Short): Unit = 
-        gsl_vector_short_set(nativeHandleChecked(), index.toULong(), value)
+        gsl_vector_short_set(nativeHandle, index.toULong(), value)
 
     override fun copy(): GslShortVector {
         val new = requireNotNull(gsl_vector_short_alloc(size.toULong()))
-        gsl_vector_short_memcpy(new, nativeHandleChecked())
+        gsl_vector_short_memcpy(new, nativeHandle)
         return GslShortVector(new, scope)
     }
 
     override fun equals(other: Any?): Boolean {
         if (other is GslShortVector)
-            return gsl_vector_short_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_vector_short_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 
-    override fun close(): Unit = gsl_vector_short_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_vector_short_free(nativeHandle)
 }
 
-internal class GslUShortVector(override val nativeHandle: CPointer<gsl_vector_ushort>, scope: DeferScope) : 
+internal class GslUShortVector(override val rawNativeHandle: CPointer<gsl_vector_ushort>, scope: AutofreeScope) : 
     GslVector<UShort, gsl_vector_ushort>(scope) {
     override val size: Int
-        get() = nativeHandleChecked().pointed.size.toInt()
+        get() = nativeHandle.pointed.size.toInt()
 
     override operator fun get(index: Int): UShort = 
-        gsl_vector_ushort_get(nativeHandleChecked(), index.toULong())
+        gsl_vector_ushort_get(nativeHandle, index.toULong())
 
     override operator fun set(index: Int, value: UShort): Unit = 
-        gsl_vector_ushort_set(nativeHandleChecked(), index.toULong(), value)
+        gsl_vector_ushort_set(nativeHandle, index.toULong(), value)
 
     override fun copy(): GslUShortVector {
         val new = requireNotNull(gsl_vector_ushort_alloc(size.toULong()))
-        gsl_vector_ushort_memcpy(new, nativeHandleChecked())
+        gsl_vector_ushort_memcpy(new, nativeHandle)
         return GslUShortVector(new, scope)
     }
 
     override fun equals(other: Any?): Boolean {
         if (other is GslUShortVector)
-            return gsl_vector_ushort_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_vector_ushort_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 
-    override fun close(): Unit = gsl_vector_ushort_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_vector_ushort_free(nativeHandle)
 }
 
-internal class GslLongVector(override val nativeHandle: CPointer<gsl_vector_long>, scope: DeferScope) : 
+internal class GslLongVector(override val rawNativeHandle: CPointer<gsl_vector_long>, scope: AutofreeScope) : 
     GslVector<Long, gsl_vector_long>(scope) {
     override val size: Int
-        get() = nativeHandleChecked().pointed.size.toInt()
+        get() = nativeHandle.pointed.size.toInt()
 
     override operator fun get(index: Int): Long = 
-        gsl_vector_long_get(nativeHandleChecked(), index.toULong())
+        gsl_vector_long_get(nativeHandle, index.toULong())
 
     override operator fun set(index: Int, value: Long): Unit = 
-        gsl_vector_long_set(nativeHandleChecked(), index.toULong(), value)
+        gsl_vector_long_set(nativeHandle, index.toULong(), value)
 
     override fun copy(): GslLongVector {
         val new = requireNotNull(gsl_vector_long_alloc(size.toULong()))
-        gsl_vector_long_memcpy(new, nativeHandleChecked())
+        gsl_vector_long_memcpy(new, nativeHandle)
         return GslLongVector(new, scope)
     }
 
     override fun equals(other: Any?): Boolean {
         if (other is GslLongVector)
-            return gsl_vector_long_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_vector_long_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 
-    override fun close(): Unit = gsl_vector_long_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_vector_long_free(nativeHandle)
 }
 
-internal class GslULongVector(override val nativeHandle: CPointer<gsl_vector_ulong>, scope: DeferScope) : 
+internal class GslULongVector(override val rawNativeHandle: CPointer<gsl_vector_ulong>, scope: AutofreeScope) : 
     GslVector<ULong, gsl_vector_ulong>(scope) {
     override val size: Int
-        get() = nativeHandleChecked().pointed.size.toInt()
+        get() = nativeHandle.pointed.size.toInt()
 
     override operator fun get(index: Int): ULong = 
-        gsl_vector_ulong_get(nativeHandleChecked(), index.toULong())
+        gsl_vector_ulong_get(nativeHandle, index.toULong())
 
     override operator fun set(index: Int, value: ULong): Unit = 
-        gsl_vector_ulong_set(nativeHandleChecked(), index.toULong(), value)
+        gsl_vector_ulong_set(nativeHandle, index.toULong(), value)
 
     override fun copy(): GslULongVector {
         val new = requireNotNull(gsl_vector_ulong_alloc(size.toULong()))
-        gsl_vector_ulong_memcpy(new, nativeHandleChecked())
+        gsl_vector_ulong_memcpy(new, nativeHandle)
         return GslULongVector(new, scope)
     }
 
     override fun equals(other: Any?): Boolean {
         if (other is GslULongVector)
-            return gsl_vector_ulong_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_vector_ulong_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 
-    override fun close(): Unit = gsl_vector_ulong_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_vector_ulong_free(nativeHandle)
 }
 
-internal class GslIntVector(override val nativeHandle: CPointer<gsl_vector_int>, scope: DeferScope) : 
+internal class GslIntVector(override val rawNativeHandle: CPointer<gsl_vector_int>, scope: AutofreeScope) : 
     GslVector<Int, gsl_vector_int>(scope) {
     override val size: Int
-        get() = nativeHandleChecked().pointed.size.toInt()
+        get() = nativeHandle.pointed.size.toInt()
 
     override operator fun get(index: Int): Int = 
-        gsl_vector_int_get(nativeHandleChecked(), index.toULong())
+        gsl_vector_int_get(nativeHandle, index.toULong())
 
     override operator fun set(index: Int, value: Int): Unit = 
-        gsl_vector_int_set(nativeHandleChecked(), index.toULong(), value)
+        gsl_vector_int_set(nativeHandle, index.toULong(), value)
 
     override fun copy(): GslIntVector {
         val new = requireNotNull(gsl_vector_int_alloc(size.toULong()))
-        gsl_vector_int_memcpy(new, nativeHandleChecked())
+        gsl_vector_int_memcpy(new, nativeHandle)
         return GslIntVector(new, scope)
     }
 
     override fun equals(other: Any?): Boolean {
         if (other is GslIntVector)
-            return gsl_vector_int_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_vector_int_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 
-    override fun close(): Unit = gsl_vector_int_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_vector_int_free(nativeHandle)
 }
 
-internal class GslUIntVector(override val nativeHandle: CPointer<gsl_vector_uint>, scope: DeferScope) : 
+internal class GslUIntVector(override val rawNativeHandle: CPointer<gsl_vector_uint>, scope: AutofreeScope) : 
     GslVector<UInt, gsl_vector_uint>(scope) {
     override val size: Int
-        get() = nativeHandleChecked().pointed.size.toInt()
+        get() = nativeHandle.pointed.size.toInt()
 
     override operator fun get(index: Int): UInt = 
-        gsl_vector_uint_get(nativeHandleChecked(), index.toULong())
+        gsl_vector_uint_get(nativeHandle, index.toULong())
 
     override operator fun set(index: Int, value: UInt): Unit = 
-        gsl_vector_uint_set(nativeHandleChecked(), index.toULong(), value)
+        gsl_vector_uint_set(nativeHandle, index.toULong(), value)
 
     override fun copy(): GslUIntVector {
         val new = requireNotNull(gsl_vector_uint_alloc(size.toULong()))
-        gsl_vector_uint_memcpy(new, nativeHandleChecked())
+        gsl_vector_uint_memcpy(new, nativeHandle)
         return GslUIntVector(new, scope)
     }
 
     override fun equals(other: Any?): Boolean {
         if (other is GslUIntVector)
-            return gsl_vector_uint_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1
+            return gsl_vector_uint_equal(nativeHandle, other.nativeHandle) == 1
 
         return super.equals(other)
     }
 
-    override fun close(): Unit = gsl_vector_uint_free(nativeHandleChecked())
+    override fun close(): Unit = gsl_vector_uint_free(nativeHandle)
 }
 
diff --git a/kmath-gsl/src/nativeTest/kotlin/GslMatrixRealTest.kt b/kmath-gsl/src/nativeTest/kotlin/GslMatrixRealTest.kt
index b4b446462..c171f6d0a 100644
--- a/kmath-gsl/src/nativeTest/kotlin/GslMatrixRealTest.kt
+++ b/kmath-gsl/src/nativeTest/kotlin/GslMatrixRealTest.kt
@@ -3,6 +3,8 @@ package kscience.kmath.gsl
 import kscience.kmath.linear.RealMatrixContext
 import kscience.kmath.operations.invoke
 import kscience.kmath.structures.Matrix
+import kscience.kmath.structures.asIterable
+import kscience.kmath.structures.toList
 import kotlin.random.Random
 import kotlin.test.Test
 import kotlin.test.assertEquals
@@ -39,4 +41,22 @@ internal class GslMatrixRealTest {
         assertEquals(mat, mat3)
         assertEquals(mat2, mat3)
     }
+
+    @Test
+    fun rows() = GslRealMatrixContext {
+        val mat = produce(2, 2) { i, j -> i.toDouble() + j }
+
+        mat.rows.asIterable().zip(listOf(listOf(0.0, 1.0), listOf(1.0, 2.0))).forEach { (a, b) ->
+            assertEquals(a.toList(), b)
+        }
+    }
+
+    @Test
+    fun columns() = GslRealMatrixContext {
+        val mat = produce(2, 2) { i, j -> i.toDouble() + j }
+
+        mat.columns.asIterable().zip(listOf(listOf(0.0, 1.0), listOf(1.0, 2.0))).forEach { (a, b) ->
+            assertEquals(a.toList(), b)
+        }
+    }
 }