diff --git a/buildSrc/build.gradle.kts b/buildSrc/build.gradle.kts
new file mode 100644
index 000000000..7ca4df19d
--- /dev/null
+++ b/buildSrc/build.gradle.kts
@@ -0,0 +1,5 @@
+plugins {
+    `kotlin-dsl`
+}
+
+repositories.mavenCentral()
diff --git a/buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt b/buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt
new file mode 100644
index 000000000..e4b157973
--- /dev/null
+++ b/buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt
@@ -0,0 +1,415 @@
+/*
+ * Copyright 2018-2021 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+@file:Suppress("KDocUnresolvedReference")
+
+package space.kscience.kmath.ejml.codegen
+
+import org.intellij.lang.annotations.Language
+import java.io.File
+
+private fun Appendable.appendEjmlVector(type: String, ejmlMatrixType: String) {
+    @Language("kotlin") val text = """/**
+ * [EjmlVector] specialization for [$type].
+ */
+public class Ejml${type}Vector<out M : $ejmlMatrixType>(public override val origin: M) : EjmlVector<$type, M>(origin) {
+    init {
+        require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" }
+    }
+
+    public override operator fun get(index: Int): $type = origin[0, index]
+}"""
+    appendLine(text)
+    appendLine()
+}
+
+private fun Appendable.appendEjmlMatrix(type: String, ejmlMatrixType: String) {
+    val text = """/**
+ * [EjmlMatrix] specialization for [$type].
+ */
+public class Ejml${type}Matrix<out M : $ejmlMatrixType>(public override val origin: M) : EjmlMatrix<$type, M>(origin) {
+    public override operator fun get(i: Int, j: Int): $type = origin[i, j]
+}"""
+    appendLine(text)
+    appendLine()
+}
+
+private fun Appendable.appendEjmlLinearSpace(
+    type: String,
+    kmathAlgebra: String,
+    ejmlMatrixParentTypeMatrix: String,
+    ejmlMatrixType: String,
+    ejmlMatrixDenseType: String,
+    ops: String,
+    denseOps: String,
+    isDense: Boolean,
+) {
+    @Language("kotlin") val text = """/**
+ * [EjmlLinearSpace] implementation based on [CommonOps_$ops], [DecompositionFactory_${ops}] operations and
+ * [${ejmlMatrixType}] matrices.
+ */
+public object EjmlLinearSpace${ops} : EjmlLinearSpace<${type}, ${kmathAlgebra}, $ejmlMatrixType>() {
+    /**
+     * The [${kmathAlgebra}] reference.
+     */
+    public override val elementAlgebra: $kmathAlgebra get() = $kmathAlgebra
+
+    @Suppress("UNCHECKED_CAST")
+    public override fun Matrix<${type}>.toEjml(): Ejml${type}Matrix<${ejmlMatrixType}> = when {
+        this is Ejml${type}Matrix<*> && origin is $ejmlMatrixType -> this as Ejml${type}Matrix<${ejmlMatrixType}>
+        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
+    }
+
+    @Suppress("UNCHECKED_CAST")
+    public override fun Point<${type}>.toEjml(): Ejml${type}Vector<${ejmlMatrixType}> = when {
+        this is Ejml${type}Vector<*> && origin is $ejmlMatrixType -> this as Ejml${type}Vector<${ejmlMatrixType}>
+        else -> Ejml${type}Vector(${ejmlMatrixType}(size, 1).also {
+            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
+        })
+    }
+
+    public override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: ${kmathAlgebra}.(i: Int, j: Int) -> ${type},
+    ): Ejml${type}Matrix<${ejmlMatrixType}> = ${ejmlMatrixType}(rows, columns).also {
+        (0 until rows).forEach { row ->
+            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
+        }
+    }.wrapMatrix()
+
+    public override fun buildVector(
+        size: Int,
+        initializer: ${kmathAlgebra}.(Int) -> ${type},
+    ): Ejml${type}Vector<${ejmlMatrixType}> = Ejml${type}Vector(${ejmlMatrixType}(size, 1).also {
+        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
+    })
+
+    private fun <T : ${ejmlMatrixParentTypeMatrix}> T.wrapMatrix() = Ejml${type}Matrix(this)
+    private fun <T : ${ejmlMatrixParentTypeMatrix}> T.wrapVector() = Ejml${type}Vector(this)
+
+    public override fun Matrix<${type}>.unaryMinus(): Matrix<${type}> = this * elementAlgebra { -one }
+
+    public override fun Matrix<${type}>.dot(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
+        val out = ${ejmlMatrixType}(1, 1)
+        CommonOps_${ops}.mult(toEjml().origin, other.toEjml().origin, out)
+        return out.wrapMatrix()
+    }
+
+    public override fun Matrix<${type}>.dot(vector: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
+        val out = ${ejmlMatrixType}(1, 1)
+        CommonOps_${ops}.mult(toEjml().origin, vector.toEjml().origin, out)
+        return out.wrapVector()
+    }
+
+    public override operator fun Matrix<${type}>.minus(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
+        val out = ${ejmlMatrixType}(1, 1)
+
+        CommonOps_${ops}.add(
+            elementAlgebra.one, 
+            toEjml().origin, 
+            elementAlgebra { -one },
+            other.toEjml().origin, 
+            out,${
+        if (isDense) "" else
+            """
+            null, 
+            null,"""
+    }
+        )
+
+        return out.wrapMatrix()
+    }
+
+    public override operator fun Matrix<${type}>.times(value: ${type}): Ejml${type}Matrix<${ejmlMatrixType}> {
+        val res = ${ejmlMatrixType}(1, 1)
+        CommonOps_${ops}.scale(value, toEjml().origin, res)
+        return res.wrapMatrix()
+    }
+
+    public override fun Point<${type}>.unaryMinus(): Ejml${type}Vector<${ejmlMatrixType}> {
+        val res = ${ejmlMatrixType}(1, 1)
+        CommonOps_${ops}.changeSign(toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    public override fun Matrix<${type}>.plus(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
+        val out = ${ejmlMatrixType}(1, 1)
+        
+        CommonOps_${ops}.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin, 
+            out,${
+        if (isDense) "" else
+            """
+            null, 
+            null,"""
+    }
+        )
+
+        return out.wrapMatrix()
+    }
+
+    public override fun Point<${type}>.plus(other: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
+        val out = ${ejmlMatrixType}(1, 1)
+
+        CommonOps_${ops}.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin,
+            out,${
+        if (isDense) "" else
+            """
+            null, 
+            null,"""
+    }
+        )
+
+        return out.wrapVector()
+    }
+
+    public override fun Point<${type}>.minus(other: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
+        val out = ${ejmlMatrixType}(1, 1)
+
+        CommonOps_${ops}.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra { -one },
+            other.toEjml().origin,
+            out,${
+        if (isDense) "" else
+            """
+            null, 
+            null,"""
+    }
+        )
+
+        return out.wrapVector()
+    }
+
+    public override fun ${type}.times(m: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> = m * this
+
+    public override fun Point<${type}>.times(value: ${type}): Ejml${type}Vector<${ejmlMatrixType}> {
+        val res = ${ejmlMatrixType}(1, 1)
+        CommonOps_${ops}.scale(value, toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    public override fun ${type}.times(v: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> = v * this
+
+    @UnstableKMathAPI
+    public override fun <F : StructureFeature> getFeature(structure: Matrix<${type}>, type: KClass<out F>): F? {
+        structure.getFeature(type)?.let { return it }
+        val origin = structure.toEjml().origin
+
+        return when (type) {
+        ${
+        if (isDense)
+            """    InverseMatrixFeature::class -> object : InverseMatrixFeature<${type}> {
+                override val inverse: Matrix<${type}> by lazy {
+                    val res = origin.copy()
+                    CommonOps_${ops}.invert(res)
+                    res.wrapMatrix()
+                }
+            }
+
+            DeterminantFeature::class -> object : DeterminantFeature<${type}> {
+                override val determinant: $type by lazy { CommonOps_${ops}.det(origin) }
+            }
+
+            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<${type}> {
+                private val svd by lazy {
+                    DecompositionFactory_${ops}.svd(origin.numRows, origin.numCols, true, true, false)
+                        .apply { decompose(origin.copy()) }
+                }
+
+                override val u: Matrix<${type}> by lazy { svd.getU(null, false).wrapMatrix() }
+                override val s: Matrix<${type}> by lazy { svd.getW(null).wrapMatrix() }
+                override val v: Matrix<${type}> by lazy { svd.getV(null, false).wrapMatrix() }
+                override val singularValues: Point<${type}> by lazy { ${type}Buffer(svd.singularValues) }
+            }
+
+            QRDecompositionFeature::class -> object : QRDecompositionFeature<${type}> {
+                private val qr by lazy {
+                    DecompositionFactory_${ops}.qr().apply { decompose(origin.copy()) }
+                }
+
+                override val q: Matrix<${type}> by lazy {
+                    qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
+                }
+
+                override val r: Matrix<${type}> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
+            }
+
+            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<${type}> {
+                override val l: Matrix<${type}> by lazy {
+                    val cholesky =
+                        DecompositionFactory_${ops}.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
+
+                    cholesky.getT(null).wrapMatrix() + LFeature
+                }
+            }
+
+            LupDecompositionFeature::class -> object : LupDecompositionFeature<${type}> {
+                private val lup by lazy {
+                    DecompositionFactory_${ops}.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
+                }
+
+                override val l: Matrix<${type}> by lazy {
+                    lup.getLower(null).wrapMatrix() + LFeature
+                }
+
+                override val u: Matrix<${type}> by lazy {
+                    lup.getUpper(null).wrapMatrix() + UFeature
+                }
+
+                override val p: Matrix<${type}> by lazy { lup.getRowPivot(null).wrapMatrix() }
+            }""" else """    QRDecompositionFeature::class -> object : QRDecompositionFeature<$type> {
+                private val qr by lazy {
+                    DecompositionFactory_${ops}.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
+                }
+
+                override val q: Matrix<${type}> by lazy {
+                    qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
+                }
+
+                override val r: Matrix<${type}> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
+            }
+
+            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<${type}> {
+                override val l: Matrix<${type}> by lazy {
+                    val cholesky =
+                        DecompositionFactory_${ops}.cholesky().apply { decompose(origin.copy()) }
+
+                    (cholesky.getT(null) as ${ejmlMatrixParentTypeMatrix}).wrapMatrix() + LFeature
+                }
+            }
+
+            LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object :
+                LUDecompositionFeature<${type}>, DeterminantFeature<${type}>, InverseMatrixFeature<${type}> {
+                private val lu by lazy {
+                    DecompositionFactory_${ops}.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
+                }
+
+                override val l: Matrix<${type}> by lazy {
+                    lu.getLower(null).wrapMatrix() + LFeature
+                }
+
+                override val u: Matrix<${type}> by lazy {
+                    lu.getUpper(null).wrapMatrix() + UFeature
+                }
+
+                override val inverse: Matrix<${type}> by lazy {
+                    var a = origin
+                    val inverse = ${ejmlMatrixDenseType}(1, 1)
+                    val solver = LinearSolverFactory_${ops}.lu(FillReducing.NONE)
+                    if (solver.modifiesA()) a = a.copy()
+                    val i = CommonOps_${denseOps}.identity(a.numRows)
+                    solver.solve(i, inverse)
+                    inverse.wrapMatrix()
+                }
+
+                override val determinant: $type by lazy { elementAlgebra.number(lu.computeDeterminant().real) }
+            }"""
+    }
+
+            else -> null
+        }?.let(type::cast)
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p matrix.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<${type}>, b: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
+        val res = ${ejmlMatrixType}(1, 1)
+        CommonOps_${ops}.solve(${ejmlMatrixType}(a.toEjml().origin), ${ejmlMatrixType}(b.toEjml().origin), res)
+        return res.wrapMatrix()
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p vector.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<${type}>, b: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
+        val res = ${ejmlMatrixType}(1, 1)
+        CommonOps_${ops}.solve(${ejmlMatrixType}(a.toEjml().origin), ${ejmlMatrixType}(b.toEjml().origin), res)
+        return Ejml${type}Vector(res)
+    }
+}"""
+    appendLine(text)
+    appendLine()
+}
+
+
+/**
+ * Generates routine EJML classes.
+ */
+fun ejmlCodegen(outputFile: String): Unit = File(outputFile).run {
+    parentFile.mkdirs()
+
+    writer().use {
+        it.appendLine("/*")
+        it.appendLine(" * Copyright 2018-2021 KMath contributors.")
+        it.appendLine(" * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.")
+        it.appendLine(" */")
+        it.appendLine()
+        it.appendLine("/* This file is generated with buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt */")
+        it.appendLine()
+        it.appendLine("package space.kscience.kmath.ejml")
+        it.appendLine()
+        it.appendLine("import org.ejml.data.*")
+        it.appendLine("import space.kscience.kmath.linear.*")
+        it.appendLine("import space.kscience.kmath.operations.*")
+        it.appendLine("import space.kscience.kmath.structures.*")
+        it.appendLine("import space.kscience.kmath.misc.*")
+        it.appendLine("import kotlin.reflect.*")
+        it.appendLine("import org.ejml.dense.row.*")
+        it.appendLine("import org.ejml.dense.row.factory.*")
+        it.appendLine("import org.ejml.sparse.*")
+        it.appendLine("import org.ejml.sparse.csc.*")
+        it.appendLine("import org.ejml.sparse.csc.factory.*")
+        it.appendLine("import space.kscience.kmath.nd.*")
+        it.appendLine("import space.kscience.kmath.linear.Matrix")
+        it.appendLine()
+        it.appendEjmlVector("Double", "DMatrix")
+        it.appendEjmlVector("Float", "FMatrix")
+        it.appendEjmlMatrix("Double", "DMatrix")
+        it.appendEjmlMatrix("Float", "FMatrix")
+        it.appendEjmlLinearSpace("Double", "DoubleField", "DMatrix", "DMatrixRMaj", "DMatrixRMaj", "DDRM", "DDRM", true)
+        it.appendEjmlLinearSpace("Float", "FloatField", "FMatrix", "FMatrixRMaj", "FMatrixRMaj", "FDRM", "FDRM", true)
+
+        it.appendEjmlLinearSpace(
+            type = "Double",
+            kmathAlgebra = "DoubleField",
+            ejmlMatrixParentTypeMatrix = "DMatrix",
+            ejmlMatrixType = "DMatrixSparseCSC",
+            ejmlMatrixDenseType = "DMatrixRMaj",
+            ops = "DSCC",
+            denseOps = "DDRM",
+            isDense = false,
+        )
+
+        it.appendEjmlLinearSpace(
+            type = "Float",
+            kmathAlgebra = "FloatField",
+            ejmlMatrixParentTypeMatrix = "FMatrix",
+            ejmlMatrixType = "FMatrixSparseCSC",
+            ejmlMatrixDenseType = "FMatrixRMaj",
+            ops = "FSCC",
+            denseOps = "FDRM",
+            isDense = false,
+        )
+    }
+}
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt
index 4a0ca7dfe..37c93d249 100644
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt
@@ -75,6 +75,23 @@ public object LFeature : MatrixFeature
  */
 public object UFeature : MatrixFeature
 
+/**
+ * Matrices with this feature support LU factorization: *a = [l] &middot; [u]* where *a* is the owning matrix.
+ *
+ * @param T the type of matrices' items.
+ */
+public interface LUDecompositionFeature<T : Any> : MatrixFeature {
+    /**
+     * The lower triangular matrix in this decomposition. It may have [LFeature].
+     */
+    public val l: Matrix<T>
+
+    /**
+     * The upper triangular matrix in this decomposition. It may have [UFeature].
+     */
+    public val u: Matrix<T>
+}
+
 /**
  * Matrices with this feature support LU factorization with partial pivoting: *[p] &middot; a = [l] &middot; [u]* where
  * *a* is the owning matrix.
diff --git a/kmath-ejml/build.gradle.kts b/kmath-ejml/build.gradle.kts
index c8e2ecd8b..5107cfb68 100644
--- a/kmath-ejml/build.gradle.kts
+++ b/kmath-ejml/build.gradle.kts
@@ -1,3 +1,5 @@
+import space.kscience.kmath.ejml.codegen.ejmlCodegen
+
 plugins {
     kotlin("jvm")
     id("ru.mipt.npm.gradle.common")
@@ -5,6 +7,9 @@ plugins {
 
 dependencies {
     api("org.ejml:ejml-ddense:0.40")
+    api("org.ejml:ejml-fdense:0.40")
+    api("org.ejml:ejml-dsparse:0.40")
+    api("org.ejml:ejml-fsparse:0.40")
     api(project(":kmath-core"))
 }
 
@@ -14,19 +19,24 @@ readme {
 
     feature(
         id = "ejml-vector",
-        description = "Point implementations.",
         ref = "src/main/kotlin/space/kscience/kmath/ejml/EjmlVector.kt"
-    )
+    ) { "Point implementations." }
 
     feature(
         id = "ejml-matrix",
-        description = "Matrix implementation.",
         ref = "src/main/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt"
-    )
+    ) { "Matrix implementation." }
 
     feature(
         id = "ejml-linear-space",
-        description = "LinearSpace implementations.",
         ref = "src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt"
-    )
+    ) { "LinearSpace implementations." }
+}
+
+kotlin.sourceSets.main {
+    val codegen by tasks.creating {
+        ejmlCodegen(kotlin.srcDirs.first().absolutePath + "/space/kscience/kmath/ejml/_generated.kt")
+    }
+
+    kotlin.srcDirs(files().builtBy(codegen))
 }
diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt
index 71cae4829..f88e83369 100644
--- a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt
+++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt
@@ -5,19 +5,10 @@
 
 package space.kscience.kmath.ejml
 
-import org.ejml.data.DMatrix
-import org.ejml.data.DMatrixD1
-import org.ejml.data.DMatrixRMaj
-import org.ejml.dense.row.CommonOps_DDRM
-import org.ejml.dense.row.factory.DecompositionFactory_DDRM
-import space.kscience.kmath.linear.*
-import space.kscience.kmath.misc.UnstableKMathAPI
-import space.kscience.kmath.nd.StructureFeature
-import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.linear.LinearSpace
+import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.linear.Point
 import space.kscience.kmath.operations.Ring
-import space.kscience.kmath.structures.DoubleBuffer
-import kotlin.reflect.KClass
-import kotlin.reflect.cast
 
 /**
  * [LinearSpace] implementation specialized for a certain EJML type.
@@ -27,7 +18,7 @@ import kotlin.reflect.cast
  * @param M the EJML matrix type.
  * @author Iaroslav Postovalov
  */
-public abstract class EjmlLinearSpace<T : Any, out A : Ring<T>, M : org.ejml.data.Matrix> : LinearSpace<T, A> {
+public abstract class EjmlLinearSpace<T : Any, out A : Ring<T>, out M : org.ejml.data.Matrix> : LinearSpace<T, A> {
     /**
      * Converts this matrix to EJML one.
      */
@@ -46,209 +37,3 @@ public abstract class EjmlLinearSpace<T : Any, out A : Ring<T>, M : org.ejml.dat
 
     public abstract override fun buildVector(size: Int, initializer: A.(Int) -> T): EjmlVector<T, M>
 }
-
-/**
- * [EjmlLinearSpace] implementation based on [CommonOps_DDRM], [DecompositionFactory_DDRM] operations and
- * [DMatrixRMaj] matrices.
- *
- * @author Iaroslav Postovalov
- */
-public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, DoubleField, DMatrixRMaj>() {
-    /**
-     * The [DoubleField] reference.
-     */
-    public override val elementAlgebra: DoubleField get() = DoubleField
-
-    @Suppress("UNCHECKED_CAST")
-    public override fun Matrix<Double>.toEjml(): EjmlDoubleMatrix<DMatrixRMaj> = when {
-        this is EjmlDoubleMatrix<*> && origin is DMatrixRMaj -> this as EjmlDoubleMatrix<DMatrixRMaj>
-        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    @Suppress("UNCHECKED_CAST")
-    public override fun Point<Double>.toEjml(): EjmlDoubleVector<DMatrixRMaj> = when {
-        this is EjmlDoubleVector<*> && origin is DMatrixRMaj -> this as EjmlDoubleVector<DMatrixRMaj>
-        else -> EjmlDoubleVector(DMatrixRMaj(size, 1).also {
-            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
-        })
-    }
-
-    public override fun buildMatrix(
-        rows: Int,
-        columns: Int,
-        initializer: DoubleField.(i: Int, j: Int) -> Double,
-    ): EjmlDoubleMatrix<DMatrixRMaj> = EjmlDoubleMatrix(DMatrixRMaj(rows, columns).also {
-        (0 until rows).forEach { row ->
-            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
-        }
-    })
-
-    public override fun buildVector(
-        size: Int,
-        initializer: DoubleField.(Int) -> Double,
-    ): EjmlDoubleVector<DMatrixRMaj> = EjmlDoubleVector(DMatrixRMaj(size, 1).also {
-        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
-    })
-
-    private fun <T : DMatrix> T.wrapMatrix() = EjmlDoubleMatrix(this)
-    private fun <T : DMatrixD1> T.wrapVector() = EjmlDoubleVector(this)
-
-    public override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * (-1.0)
-
-    public override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.mult(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    public override fun Matrix<Double>.dot(vector: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.mult(toEjml().origin, vector.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    public override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.subtract(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    public override operator fun Matrix<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixRMaj> {
-        val res = this.toEjml().origin.copy()
-        CommonOps_DDRM.scale(value, res)
-        return res.wrapMatrix()
-    }
-
-    public override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixRMaj> {
-        val out = toEjml().origin.copy()
-        CommonOps_DDRM.changeSign(out)
-        return out.wrapVector()
-    }
-
-    public override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.add(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    public override fun Point<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.add(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    public override fun Point<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.subtract(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    public override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> = m * this
-
-    public override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixRMaj> {
-        val res = this.toEjml().origin.copy()
-        CommonOps_DDRM.scale(value, res)
-        return res.wrapVector()
-    }
-
-    public override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixRMaj> = v * this
-
-    @UnstableKMathAPI
-    public override fun <F : StructureFeature> getFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
-        // Return the feature if it is intrinsic to the structure
-        structure.getFeature(type)?.let { return it }
-
-        val origin = structure.toEjml().origin
-
-        return when (type) {
-            InverseMatrixFeature::class -> object : InverseMatrixFeature<Double> {
-                override val inverse: Matrix<Double> by lazy {
-                    val res = origin.copy()
-                    CommonOps_DDRM.invert(res)
-                    EjmlDoubleMatrix(res)
-                }
-            }
-
-            DeterminantFeature::class -> object : DeterminantFeature<Double> {
-                override val determinant: Double by lazy { CommonOps_DDRM.det(DMatrixRMaj(origin)) }
-            }
-
-            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Double> {
-                private val svd by lazy {
-                    DecompositionFactory_DDRM.svd(origin.numRows, origin.numCols, true, true, false)
-                        .apply { decompose(origin.copy()) }
-                }
-
-                override val u: Matrix<Double> by lazy { EjmlDoubleMatrix(svd.getU(null, false)) }
-                override val s: Matrix<Double> by lazy { EjmlDoubleMatrix(svd.getW(null)) }
-                override val v: Matrix<Double> by lazy { EjmlDoubleMatrix(svd.getV(null, false)) }
-                override val singularValues: Point<Double> by lazy { DoubleBuffer(svd.singularValues) }
-            }
-
-            QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
-                private val qr by lazy {
-                    DecompositionFactory_DDRM.qr().apply { decompose(origin.copy()) }
-                }
-
-                override val q: Matrix<Double> by lazy {
-                    EjmlDoubleMatrix(qr.getQ(null, false)) + OrthogonalFeature
-                }
-
-                override val r: Matrix<Double> by lazy { EjmlDoubleMatrix(qr.getR(null, false)) + UFeature }
-            }
-
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
-                override val l: Matrix<Double> by lazy {
-                    val cholesky =
-                        DecompositionFactory_DDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
-
-                    EjmlDoubleMatrix(cholesky.getT(null)) + LFeature
-                }
-            }
-
-            LupDecompositionFeature::class -> object : LupDecompositionFeature<Double> {
-                private val lup by lazy {
-                    DecompositionFactory_DDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
-                }
-
-                override val l: Matrix<Double> by lazy {
-                    EjmlDoubleMatrix(lup.getLower(null)) + LFeature
-                }
-
-                override val u: Matrix<Double> by lazy {
-                    EjmlDoubleMatrix(lup.getUpper(null)) + UFeature
-                }
-
-                override val p: Matrix<Double> by lazy { EjmlDoubleMatrix(lup.getRowPivot(null)) }
-            }
-
-            else -> null
-        }?.let(type::cast)
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
-     *
-     * @param a the base matrix.
-     * @param b n by p matrix.
-     * @return the solution for 'x' that is n by p.
-     */
-    public fun solve(a: Matrix<Double>, b: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        val res = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res)
-        return EjmlDoubleMatrix(res)
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
-     *
-     * @param a the base matrix.
-     * @param b n by p vector.
-     * @return the solution for 'x' that is n by p.
-     */
-    public fun solve(a: Matrix<Double>, b: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        val res = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res)
-        return EjmlDoubleVector(res)
-    }
-}
diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt
index 92c4d1cf0..fd43f295c 100644
--- a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt
+++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt
@@ -5,7 +5,6 @@
 
 package space.kscience.kmath.ejml
 
-import org.ejml.data.DMatrix
 import org.ejml.data.Matrix
 import space.kscience.kmath.nd.Structure2D
 
@@ -21,12 +20,3 @@ public abstract class EjmlMatrix<T, out M : Matrix>(public open val origin: M) :
     public override val rowNum: Int get() = origin.numRows
     public override val colNum: Int get() = origin.numCols
 }
-
-/**
- * [EjmlMatrix] specialization for [Double].
- *
- * @author Iaroslav Postovalov
- */
-public class EjmlDoubleMatrix<out M : DMatrix>(public override val origin: M) : EjmlMatrix<Double, M>(origin) {
-    public override operator fun get(i: Int, j: Int): Double = origin[i, j]
-}
diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlVector.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlVector.kt
index 81502d6d0..5d10d1fbb 100644
--- a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlVector.kt
+++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlVector.kt
@@ -5,7 +5,6 @@
 
 package space.kscience.kmath.ejml
 
-import org.ejml.data.DMatrixD1
 import org.ejml.data.Matrix
 import space.kscience.kmath.linear.Point
 
@@ -14,12 +13,12 @@ import space.kscience.kmath.linear.Point
  *
  * @param T the type of elements contained in the buffer.
  * @param M the type of EJML matrix.
- * @property origin The underlying matrix.
+ * @property origin The underlying matrix, must have only one row.
  * @author Iaroslav Postovalov
  */
 public abstract class EjmlVector<out T, out M : Matrix>(public open val origin: M) : Point<T> {
     public override val size: Int
-        get() = origin.numRows
+        get() = origin.numCols
 
     public override operator fun iterator(): Iterator<T> = object : Iterator<T> {
         private var cursor: Int = 0
@@ -34,12 +33,3 @@ public abstract class EjmlVector<out T, out M : Matrix>(public open val origin:
 
     public override fun toString(): String = "EjmlVector(origin=$origin)"
 }
-
-/**
- * [EjmlVector] specialization for [Double].
- *
- * @author Iaroslav Postovalov
- */
-public class EjmlDoubleVector<out M : DMatrixD1>(public override val origin: M) : EjmlVector<Double, M>(origin) {
-    public override operator fun get(index: Int): Double = origin[index]
-}
diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt
new file mode 100644
index 000000000..139c55697
--- /dev/null
+++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt
@@ -0,0 +1,995 @@
+/*
+ * Copyright 2018-2021 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+/* This file is generated with buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt */
+
+package space.kscience.kmath.ejml
+
+import org.ejml.data.*
+import org.ejml.dense.row.CommonOps_DDRM
+import org.ejml.dense.row.CommonOps_FDRM
+import org.ejml.dense.row.factory.DecompositionFactory_DDRM
+import org.ejml.dense.row.factory.DecompositionFactory_FDRM
+import org.ejml.sparse.FillReducing
+import org.ejml.sparse.csc.CommonOps_DSCC
+import org.ejml.sparse.csc.CommonOps_FSCC
+import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC
+import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC
+import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC
+import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC
+import space.kscience.kmath.linear.*
+import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.nd.StructureFeature
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.operations.FloatField
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.FloatBuffer
+import kotlin.reflect.KClass
+import kotlin.reflect.cast
+
+/**
+ * [EjmlVector] specialization for [Double].
+ */
+public class EjmlDoubleVector<out M : DMatrix>(public override val origin: M) : EjmlVector<Double, M>(origin) {
+    init {
+        require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" }
+    }
+
+    public override operator fun get(index: Int): Double = origin[0, index]
+}
+
+/**
+ * [EjmlVector] specialization for [Float].
+ */
+public class EjmlFloatVector<out M : FMatrix>(public override val origin: M) : EjmlVector<Float, M>(origin) {
+    init {
+        require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" }
+    }
+
+    public override operator fun get(index: Int): Float = origin[0, index]
+}
+
+/**
+ * [EjmlMatrix] specialization for [Double].
+ */
+public class EjmlDoubleMatrix<out M : DMatrix>(public override val origin: M) : EjmlMatrix<Double, M>(origin) {
+    public override operator fun get(i: Int, j: Int): Double = origin[i, j]
+}
+
+/**
+ * [EjmlMatrix] specialization for [Float].
+ */
+public class EjmlFloatMatrix<out M : FMatrix>(public override val origin: M) : EjmlMatrix<Float, M>(origin) {
+    public override operator fun get(i: Int, j: Int): Float = origin[i, j]
+}
+
+/**
+ * [EjmlLinearSpace] implementation based on [CommonOps_DDRM], [DecompositionFactory_DDRM] operations and
+ * [DMatrixRMaj] matrices.
+ */
+public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, DoubleField, DMatrixRMaj>() {
+    /**
+     * The [DoubleField] reference.
+     */
+    public override val elementAlgebra: DoubleField get() = DoubleField
+
+    @Suppress("UNCHECKED_CAST")
+    public override fun Matrix<Double>.toEjml(): EjmlDoubleMatrix<DMatrixRMaj> = when {
+        this is EjmlDoubleMatrix<*> && origin is DMatrixRMaj -> this as EjmlDoubleMatrix<DMatrixRMaj>
+        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
+    }
+
+    @Suppress("UNCHECKED_CAST")
+    public override fun Point<Double>.toEjml(): EjmlDoubleVector<DMatrixRMaj> = when {
+        this is EjmlDoubleVector<*> && origin is DMatrixRMaj -> this as EjmlDoubleVector<DMatrixRMaj>
+        else -> EjmlDoubleVector(DMatrixRMaj(size, 1).also {
+            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
+        })
+    }
+
+    public override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: DoubleField.(i: Int, j: Int) -> Double,
+    ): EjmlDoubleMatrix<DMatrixRMaj> = DMatrixRMaj(rows, columns).also {
+        (0 until rows).forEach { row ->
+            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
+        }
+    }.wrapMatrix()
+
+    public override fun buildVector(
+        size: Int,
+        initializer: DoubleField.(Int) -> Double,
+    ): EjmlDoubleVector<DMatrixRMaj> = EjmlDoubleVector(DMatrixRMaj(size, 1).also {
+        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
+    })
+
+    private fun <T : DMatrix> T.wrapMatrix() = EjmlDoubleMatrix(this)
+    private fun <T : DMatrix> T.wrapVector() = EjmlDoubleVector(this)
+
+    public override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
+
+    public override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.mult(toEjml().origin, other.toEjml().origin, out)
+        return out.wrapMatrix()
+    }
+
+    public override fun Matrix<Double>.dot(vector: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.mult(toEjml().origin, vector.toEjml().origin, out)
+        return out.wrapVector()
+    }
+
+    public override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+
+        CommonOps_DDRM.add(
+            elementAlgebra.one, 
+            toEjml().origin, 
+            elementAlgebra { -one },
+            other.toEjml().origin, 
+            out,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    public override operator fun Matrix<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixRMaj> {
+        val res = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.scale(value, toEjml().origin, res)
+        return res.wrapMatrix()
+    }
+
+    public override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixRMaj> {
+        val res = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.changeSign(toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    public override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+        
+        CommonOps_DDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin, 
+            out,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    public override fun Point<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+
+        CommonOps_DDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin,
+            out,
+        )
+
+        return out.wrapVector()
+    }
+
+    public override fun Point<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+
+        CommonOps_DDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra { -one },
+            other.toEjml().origin,
+            out,
+        )
+
+        return out.wrapVector()
+    }
+
+    public override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> = m * this
+
+    public override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixRMaj> {
+        val res = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.scale(value, toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    public override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixRMaj> = v * this
+
+    @UnstableKMathAPI
+    public override fun <F : StructureFeature> getFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
+        structure.getFeature(type)?.let { return it }
+        val origin = structure.toEjml().origin
+
+        return when (type) {
+            InverseMatrixFeature::class -> object : InverseMatrixFeature<Double> {
+                override val inverse: Matrix<Double> by lazy {
+                    val res = origin.copy()
+                    CommonOps_DDRM.invert(res)
+                    res.wrapMatrix()
+                }
+            }
+
+            DeterminantFeature::class -> object : DeterminantFeature<Double> {
+                override val determinant: Double by lazy { CommonOps_DDRM.det(origin) }
+            }
+
+            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Double> {
+                private val svd by lazy {
+                    DecompositionFactory_DDRM.svd(origin.numRows, origin.numCols, true, true, false)
+                        .apply { decompose(origin.copy()) }
+                }
+
+                override val u: Matrix<Double> by lazy { svd.getU(null, false).wrapMatrix() }
+                override val s: Matrix<Double> by lazy { svd.getW(null).wrapMatrix() }
+                override val v: Matrix<Double> by lazy { svd.getV(null, false).wrapMatrix() }
+                override val singularValues: Point<Double> by lazy { DoubleBuffer(svd.singularValues) }
+            }
+
+            QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
+                private val qr by lazy {
+                    DecompositionFactory_DDRM.qr().apply { decompose(origin.copy()) }
+                }
+
+                override val q: Matrix<Double> by lazy {
+                    qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
+                }
+
+                override val r: Matrix<Double> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
+            }
+
+            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
+                override val l: Matrix<Double> by lazy {
+                    val cholesky =
+                        DecompositionFactory_DDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
+
+                    cholesky.getT(null).wrapMatrix() + LFeature
+                }
+            }
+
+            LupDecompositionFeature::class -> object : LupDecompositionFeature<Double> {
+                private val lup by lazy {
+                    DecompositionFactory_DDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
+                }
+
+                override val l: Matrix<Double> by lazy {
+                    lup.getLower(null).wrapMatrix() + LFeature
+                }
+
+                override val u: Matrix<Double> by lazy {
+                    lup.getUpper(null).wrapMatrix() + UFeature
+                }
+
+                override val p: Matrix<Double> by lazy { lup.getRowPivot(null).wrapMatrix() }
+            }
+
+            else -> null
+        }?.let(type::cast)
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p matrix.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Double>, b: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+        val res = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res)
+        return res.wrapMatrix()
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p vector.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Double>, b: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+        val res = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res)
+        return EjmlDoubleVector(res)
+    }
+}
+
+/**
+ * [EjmlLinearSpace] implementation based on [CommonOps_FDRM], [DecompositionFactory_FDRM] operations and
+ * [FMatrixRMaj] matrices.
+ */
+public object EjmlLinearSpaceFDRM : EjmlLinearSpace<Float, FloatField, FMatrixRMaj>() {
+    /**
+     * The [FloatField] reference.
+     */
+    public override val elementAlgebra: FloatField get() = FloatField
+
+    @Suppress("UNCHECKED_CAST")
+    public override fun Matrix<Float>.toEjml(): EjmlFloatMatrix<FMatrixRMaj> = when {
+        this is EjmlFloatMatrix<*> && origin is FMatrixRMaj -> this as EjmlFloatMatrix<FMatrixRMaj>
+        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
+    }
+
+    @Suppress("UNCHECKED_CAST")
+    public override fun Point<Float>.toEjml(): EjmlFloatVector<FMatrixRMaj> = when {
+        this is EjmlFloatVector<*> && origin is FMatrixRMaj -> this as EjmlFloatVector<FMatrixRMaj>
+        else -> EjmlFloatVector(FMatrixRMaj(size, 1).also {
+            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
+        })
+    }
+
+    public override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: FloatField.(i: Int, j: Int) -> Float,
+    ): EjmlFloatMatrix<FMatrixRMaj> = FMatrixRMaj(rows, columns).also {
+        (0 until rows).forEach { row ->
+            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
+        }
+    }.wrapMatrix()
+
+    public override fun buildVector(
+        size: Int,
+        initializer: FloatField.(Int) -> Float,
+    ): EjmlFloatVector<FMatrixRMaj> = EjmlFloatVector(FMatrixRMaj(size, 1).also {
+        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
+    })
+
+    private fun <T : FMatrix> T.wrapMatrix() = EjmlFloatMatrix(this)
+    private fun <T : FMatrix> T.wrapVector() = EjmlFloatVector(this)
+
+    public override fun Matrix<Float>.unaryMinus(): Matrix<Float> = this * elementAlgebra { -one }
+
+    public override fun Matrix<Float>.dot(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.mult(toEjml().origin, other.toEjml().origin, out)
+        return out.wrapMatrix()
+    }
+
+    public override fun Matrix<Float>.dot(vector: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.mult(toEjml().origin, vector.toEjml().origin, out)
+        return out.wrapVector()
+    }
+
+    public override operator fun Matrix<Float>.minus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+
+        CommonOps_FDRM.add(
+            elementAlgebra.one, 
+            toEjml().origin, 
+            elementAlgebra { -one },
+            other.toEjml().origin, 
+            out,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    public override operator fun Matrix<Float>.times(value: Float): EjmlFloatMatrix<FMatrixRMaj> {
+        val res = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.scale(value, toEjml().origin, res)
+        return res.wrapMatrix()
+    }
+
+    public override fun Point<Float>.unaryMinus(): EjmlFloatVector<FMatrixRMaj> {
+        val res = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.changeSign(toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    public override fun Matrix<Float>.plus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+        
+        CommonOps_FDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin, 
+            out,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    public override fun Point<Float>.plus(other: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+
+        CommonOps_FDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin,
+            out,
+        )
+
+        return out.wrapVector()
+    }
+
+    public override fun Point<Float>.minus(other: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+
+        CommonOps_FDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra { -one },
+            other.toEjml().origin,
+            out,
+        )
+
+        return out.wrapVector()
+    }
+
+    public override fun Float.times(m: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> = m * this
+
+    public override fun Point<Float>.times(value: Float): EjmlFloatVector<FMatrixRMaj> {
+        val res = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.scale(value, toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    public override fun Float.times(v: Point<Float>): EjmlFloatVector<FMatrixRMaj> = v * this
+
+    @UnstableKMathAPI
+    public override fun <F : StructureFeature> getFeature(structure: Matrix<Float>, type: KClass<out F>): F? {
+        structure.getFeature(type)?.let { return it }
+        val origin = structure.toEjml().origin
+
+        return when (type) {
+            InverseMatrixFeature::class -> object : InverseMatrixFeature<Float> {
+                override val inverse: Matrix<Float> by lazy {
+                    val res = origin.copy()
+                    CommonOps_FDRM.invert(res)
+                    res.wrapMatrix()
+                }
+            }
+
+            DeterminantFeature::class -> object : DeterminantFeature<Float> {
+                override val determinant: Float by lazy { CommonOps_FDRM.det(origin) }
+            }
+
+            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Float> {
+                private val svd by lazy {
+                    DecompositionFactory_FDRM.svd(origin.numRows, origin.numCols, true, true, false)
+                        .apply { decompose(origin.copy()) }
+                }
+
+                override val u: Matrix<Float> by lazy { svd.getU(null, false).wrapMatrix() }
+                override val s: Matrix<Float> by lazy { svd.getW(null).wrapMatrix() }
+                override val v: Matrix<Float> by lazy { svd.getV(null, false).wrapMatrix() }
+                override val singularValues: Point<Float> by lazy { FloatBuffer(svd.singularValues) }
+            }
+
+            QRDecompositionFeature::class -> object : QRDecompositionFeature<Float> {
+                private val qr by lazy {
+                    DecompositionFactory_FDRM.qr().apply { decompose(origin.copy()) }
+                }
+
+                override val q: Matrix<Float> by lazy {
+                    qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
+                }
+
+                override val r: Matrix<Float> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
+            }
+
+            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Float> {
+                override val l: Matrix<Float> by lazy {
+                    val cholesky =
+                        DecompositionFactory_FDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
+
+                    cholesky.getT(null).wrapMatrix() + LFeature
+                }
+            }
+
+            LupDecompositionFeature::class -> object : LupDecompositionFeature<Float> {
+                private val lup by lazy {
+                    DecompositionFactory_FDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
+                }
+
+                override val l: Matrix<Float> by lazy {
+                    lup.getLower(null).wrapMatrix() + LFeature
+                }
+
+                override val u: Matrix<Float> by lazy {
+                    lup.getUpper(null).wrapMatrix() + UFeature
+                }
+
+                override val p: Matrix<Float> by lazy { lup.getRowPivot(null).wrapMatrix() }
+            }
+
+            else -> null
+        }?.let(type::cast)
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p matrix.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Float>, b: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
+        val res = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.solve(FMatrixRMaj(a.toEjml().origin), FMatrixRMaj(b.toEjml().origin), res)
+        return res.wrapMatrix()
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p vector.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Float>, b: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
+        val res = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.solve(FMatrixRMaj(a.toEjml().origin), FMatrixRMaj(b.toEjml().origin), res)
+        return EjmlFloatVector(res)
+    }
+}
+
+/**
+ * [EjmlLinearSpace] implementation based on [CommonOps_DSCC], [DecompositionFactory_DSCC] operations and
+ * [DMatrixSparseCSC] matrices.
+ */
+public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, DoubleField, DMatrixSparseCSC>() {
+    /**
+     * The [DoubleField] reference.
+     */
+    public override val elementAlgebra: DoubleField get() = DoubleField
+
+    @Suppress("UNCHECKED_CAST")
+    public override fun Matrix<Double>.toEjml(): EjmlDoubleMatrix<DMatrixSparseCSC> = when {
+        this is EjmlDoubleMatrix<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleMatrix<DMatrixSparseCSC>
+        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
+    }
+
+    @Suppress("UNCHECKED_CAST")
+    public override fun Point<Double>.toEjml(): EjmlDoubleVector<DMatrixSparseCSC> = when {
+        this is EjmlDoubleVector<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleVector<DMatrixSparseCSC>
+        else -> EjmlDoubleVector(DMatrixSparseCSC(size, 1).also {
+            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
+        })
+    }
+
+    public override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: DoubleField.(i: Int, j: Int) -> Double,
+    ): EjmlDoubleMatrix<DMatrixSparseCSC> = DMatrixSparseCSC(rows, columns).also {
+        (0 until rows).forEach { row ->
+            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
+        }
+    }.wrapMatrix()
+
+    public override fun buildVector(
+        size: Int,
+        initializer: DoubleField.(Int) -> Double,
+    ): EjmlDoubleVector<DMatrixSparseCSC> = EjmlDoubleVector(DMatrixSparseCSC(size, 1).also {
+        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
+    })
+
+    private fun <T : DMatrix> T.wrapMatrix() = EjmlDoubleMatrix(this)
+    private fun <T : DMatrix> T.wrapVector() = EjmlDoubleVector(this)
+
+    public override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
+
+    public override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.mult(toEjml().origin, other.toEjml().origin, out)
+        return out.wrapMatrix()
+    }
+
+    public override fun Matrix<Double>.dot(vector: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.mult(toEjml().origin, vector.toEjml().origin, out)
+        return out.wrapVector()
+    }
+
+    public override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+
+        CommonOps_DSCC.add(
+            elementAlgebra.one, 
+            toEjml().origin, 
+            elementAlgebra { -one },
+            other.toEjml().origin, 
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    public override operator fun Matrix<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixSparseCSC> {
+        val res = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.scale(value, toEjml().origin, res)
+        return res.wrapMatrix()
+    }
+
+    public override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixSparseCSC> {
+        val res = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.changeSign(toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    public override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+        
+        CommonOps_DSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin, 
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    public override fun Point<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+
+        CommonOps_DSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin,
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapVector()
+    }
+
+    public override fun Point<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+
+        CommonOps_DSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra { -one },
+            other.toEjml().origin,
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapVector()
+    }
+
+    public override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> = m * this
+
+    public override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixSparseCSC> {
+        val res = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.scale(value, toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    public override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> = v * this
+
+    @UnstableKMathAPI
+    public override fun <F : StructureFeature> getFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
+        structure.getFeature(type)?.let { return it }
+        val origin = structure.toEjml().origin
+
+        return when (type) {
+            QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
+                private val qr by lazy {
+                    DecompositionFactory_DSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
+                }
+
+                override val q: Matrix<Double> by lazy {
+                    qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
+                }
+
+                override val r: Matrix<Double> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
+            }
+
+            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
+                override val l: Matrix<Double> by lazy {
+                    val cholesky =
+                        DecompositionFactory_DSCC.cholesky().apply { decompose(origin.copy()) }
+
+                    (cholesky.getT(null) as DMatrix).wrapMatrix() + LFeature
+                }
+            }
+
+            LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object :
+                LUDecompositionFeature<Double>, DeterminantFeature<Double>, InverseMatrixFeature<Double> {
+                private val lu by lazy {
+                    DecompositionFactory_DSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
+                }
+
+                override val l: Matrix<Double> by lazy {
+                    lu.getLower(null).wrapMatrix() + LFeature
+                }
+
+                override val u: Matrix<Double> by lazy {
+                    lu.getUpper(null).wrapMatrix() + UFeature
+                }
+
+                override val inverse: Matrix<Double> by lazy {
+                    var a = origin
+                    val inverse = DMatrixRMaj(1, 1)
+                    val solver = LinearSolverFactory_DSCC.lu(FillReducing.NONE)
+                    if (solver.modifiesA()) a = a.copy()
+                    val i = CommonOps_DDRM.identity(a.numRows)
+                    solver.solve(i, inverse)
+                    inverse.wrapMatrix()
+                }
+
+                override val determinant: Double by lazy { elementAlgebra.number(lu.computeDeterminant().real) }
+            }
+
+            else -> null
+        }?.let(type::cast)
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p matrix.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Double>, b: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+        val res = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res)
+        return res.wrapMatrix()
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p vector.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Double>, b: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+        val res = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res)
+        return EjmlDoubleVector(res)
+    }
+}
+
+/**
+ * [EjmlLinearSpace] implementation based on [CommonOps_FSCC], [DecompositionFactory_FSCC] operations and
+ * [FMatrixSparseCSC] matrices.
+ */
+public object EjmlLinearSpaceFSCC : EjmlLinearSpace<Float, FloatField, FMatrixSparseCSC>() {
+    /**
+     * The [FloatField] reference.
+     */
+    public override val elementAlgebra: FloatField get() = FloatField
+
+    @Suppress("UNCHECKED_CAST")
+    public override fun Matrix<Float>.toEjml(): EjmlFloatMatrix<FMatrixSparseCSC> = when {
+        this is EjmlFloatMatrix<*> && origin is FMatrixSparseCSC -> this as EjmlFloatMatrix<FMatrixSparseCSC>
+        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
+    }
+
+    @Suppress("UNCHECKED_CAST")
+    public override fun Point<Float>.toEjml(): EjmlFloatVector<FMatrixSparseCSC> = when {
+        this is EjmlFloatVector<*> && origin is FMatrixSparseCSC -> this as EjmlFloatVector<FMatrixSparseCSC>
+        else -> EjmlFloatVector(FMatrixSparseCSC(size, 1).also {
+            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
+        })
+    }
+
+    public override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: FloatField.(i: Int, j: Int) -> Float,
+    ): EjmlFloatMatrix<FMatrixSparseCSC> = FMatrixSparseCSC(rows, columns).also {
+        (0 until rows).forEach { row ->
+            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
+        }
+    }.wrapMatrix()
+
+    public override fun buildVector(
+        size: Int,
+        initializer: FloatField.(Int) -> Float,
+    ): EjmlFloatVector<FMatrixSparseCSC> = EjmlFloatVector(FMatrixSparseCSC(size, 1).also {
+        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
+    })
+
+    private fun <T : FMatrix> T.wrapMatrix() = EjmlFloatMatrix(this)
+    private fun <T : FMatrix> T.wrapVector() = EjmlFloatVector(this)
+
+    public override fun Matrix<Float>.unaryMinus(): Matrix<Float> = this * elementAlgebra { -one }
+
+    public override fun Matrix<Float>.dot(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.mult(toEjml().origin, other.toEjml().origin, out)
+        return out.wrapMatrix()
+    }
+
+    public override fun Matrix<Float>.dot(vector: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.mult(toEjml().origin, vector.toEjml().origin, out)
+        return out.wrapVector()
+    }
+
+    public override operator fun Matrix<Float>.minus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+
+        CommonOps_FSCC.add(
+            elementAlgebra.one, 
+            toEjml().origin, 
+            elementAlgebra { -one },
+            other.toEjml().origin, 
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    public override operator fun Matrix<Float>.times(value: Float): EjmlFloatMatrix<FMatrixSparseCSC> {
+        val res = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.scale(value, toEjml().origin, res)
+        return res.wrapMatrix()
+    }
+
+    public override fun Point<Float>.unaryMinus(): EjmlFloatVector<FMatrixSparseCSC> {
+        val res = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.changeSign(toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    public override fun Matrix<Float>.plus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+        
+        CommonOps_FSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin, 
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    public override fun Point<Float>.plus(other: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+
+        CommonOps_FSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin,
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapVector()
+    }
+
+    public override fun Point<Float>.minus(other: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+
+        CommonOps_FSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra { -one },
+            other.toEjml().origin,
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapVector()
+    }
+
+    public override fun Float.times(m: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> = m * this
+
+    public override fun Point<Float>.times(value: Float): EjmlFloatVector<FMatrixSparseCSC> {
+        val res = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.scale(value, toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    public override fun Float.times(v: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> = v * this
+
+    @UnstableKMathAPI
+    public override fun <F : StructureFeature> getFeature(structure: Matrix<Float>, type: KClass<out F>): F? {
+        structure.getFeature(type)?.let { return it }
+        val origin = structure.toEjml().origin
+
+        return when (type) {
+            QRDecompositionFeature::class -> object : QRDecompositionFeature<Float> {
+                private val qr by lazy {
+                    DecompositionFactory_FSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
+                }
+
+                override val q: Matrix<Float> by lazy {
+                    qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
+                }
+
+                override val r: Matrix<Float> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
+            }
+
+            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Float> {
+                override val l: Matrix<Float> by lazy {
+                    val cholesky =
+                        DecompositionFactory_FSCC.cholesky().apply { decompose(origin.copy()) }
+
+                    (cholesky.getT(null) as FMatrix).wrapMatrix() + LFeature
+                }
+            }
+
+            LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object :
+                LUDecompositionFeature<Float>, DeterminantFeature<Float>, InverseMatrixFeature<Float> {
+                private val lu by lazy {
+                    DecompositionFactory_FSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
+                }
+
+                override val l: Matrix<Float> by lazy {
+                    lu.getLower(null).wrapMatrix() + LFeature
+                }
+
+                override val u: Matrix<Float> by lazy {
+                    lu.getUpper(null).wrapMatrix() + UFeature
+                }
+
+                override val inverse: Matrix<Float> by lazy {
+                    var a = origin
+                    val inverse = FMatrixRMaj(1, 1)
+                    val solver = LinearSolverFactory_FSCC.lu(FillReducing.NONE)
+                    if (solver.modifiesA()) a = a.copy()
+                    val i = CommonOps_FDRM.identity(a.numRows)
+                    solver.solve(i, inverse)
+                    inverse.wrapMatrix()
+                }
+
+                override val determinant: Float by lazy { elementAlgebra.number(lu.computeDeterminant().real) }
+            }
+
+            else -> null
+        }?.let(type::cast)
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p matrix.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Float>, b: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
+        val res = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.solve(FMatrixSparseCSC(a.toEjml().origin), FMatrixSparseCSC(b.toEjml().origin), res)
+        return res.wrapMatrix()
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p vector.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Float>, b: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
+        val res = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.solve(FMatrixSparseCSC(a.toEjml().origin), FMatrixSparseCSC(b.toEjml().origin), res)
+        return EjmlFloatVector(res)
+    }
+}
+
diff --git a/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlVectorTest.kt b/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlVectorTest.kt
index 9bf76033d..9592bfa6c 100644
--- a/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlVectorTest.kt
+++ b/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlVectorTest.kt
@@ -18,7 +18,7 @@ internal class EjmlVectorTest {
 
     private val randomMatrix: DMatrixRMaj
         get() {
-            val d = DMatrixRMaj(random.nextInt(2, 100), 1)
+            val d = DMatrixRMaj(1, random.nextInt(2, 100))
             RandomMatrices_DDRM.fillUniform(d, random.asJavaRandom())
             return d
         }
@@ -27,7 +27,7 @@ internal class EjmlVectorTest {
     fun size() {
         val m = randomMatrix
         val w = EjmlDoubleVector(m)
-        assertEquals(m.numRows, w.size)
+        assertEquals(m.numCols, w.size)
     }
 
     @Test
@@ -43,7 +43,7 @@ internal class EjmlVectorTest {
         val w = EjmlDoubleVector(m)
 
         assertEquals(
-            m.iterator(true, 0, 0, m.numRows - 1, 0).asSequence().toList(),
+            m.iterator(true, 0, 0, 0, m.numCols - 1).asSequence().toList(),
             w.iterator().asSequence().toList()
         )
     }