diff --git a/doc/DaphneDSL/Builtins.md b/doc/DaphneDSL/Builtins.md
index c9214a1c3..110f15822 100644
--- a/doc/DaphneDSL/Builtins.md
+++ b/doc/DaphneDSL/Builtins.md
@@ -330,6 +330,23 @@ The following built-in functions all follow the same scheme:
 | `cumMin` | cumulative minimum |
 | `cumMax` | cumulative maximum |
 
+## Map
+
+Standard element-wise mapping, as well as row- and column-wise mapping is supported.
+
+- **`map`**`(arg:matrix, func:str)`
+
+    Element-wise mapping over a *(n x m)* matrix `arg` using a user-defined function `func` written in DaphneDSL. Applies the given UDF to each element of the given matrix.
+
+- **`map`**`(arg:matrix, func:str, axis:si64[, udfReturnsScalar:bool])`
+
+    Row- or column-wise mapping over a *(n x m)* matrix `arg` using a user-defined function `func` written in DaphneDSL. Applies the given UDF to each row/column of the given matrix. If the input of the UDF is a row matrix, the output can be a row or a scalar; if the input is a column matrix, the output can be a column or a scalar.
+    
+    - `axis` == 0: Map an entire row of the input matrix to an entire row of the output matrix; the result is a *(n x ?)* matrix
+    - `axis` == 1: Map an entire column of the input matrix to an entire column of the output matrix; the result is a *(? x m)* matrix
+    - `udfReturnsScalar` == false (optional): The given UDF `func` returns a matrix (default for row-/column-wise map), must match the UDFs output type; the result is as previously described
+    - `udfReturnsScalar` == true: The given UDF `func` returns a scalar; the result is a *(n x 1)* (column) or matrix *(1 x m)* (row) matrix (depending on `axis`)
+
 ## Reorganization
 
 - **`reshape`**`(arg:matrix, numRows:size, numCols:size)`
@@ -701,4 +718,4 @@ These must be provided in a separate [`.meta`-file](/doc/FileMetaDataFormat.md).
 - **`remove`**`(lst:list, idx:size)`
 
     Removes the element at position `idx` (counting starts at zero) from the given list `lst`.
-    Returns (1) the result as a new list (the argument list stays unchanged), and (2) the removed element.
\ No newline at end of file
+    Returns (1) the result as a new list (the argument list stays unchanged), and (2) the removed element.
diff --git a/src/compiler/lowering/LowerToLLVMPass.cpp b/src/compiler/lowering/LowerToLLVMPass.cpp
index a30dc4199..262cadc59 100644
--- a/src/compiler/lowering/LowerToLLVMPass.cpp
+++ b/src/compiler/lowering/LowerToLLVMPass.cpp
@@ -508,11 +508,17 @@ class MapOpLowering : public OpConversionPattern<daphne::MapOp> {
         // Pointer to UDF
         callee << "__void";
 
+        // Axis
+        callee << "__int64_t";
+
+        // udfReturnsScalar
+        callee << "__bool";
+
         // get pointer to UDF
         LLVM::LLVMFuncOp udfFuncOp = module.lookupSymbol<LLVM::LLVMFuncOp>(op.getFunc());
         auto udfFnPtr = rewriter.create<LLVM::AddressOfOp>(loc, udfFuncOp);
 
-        std::vector<Value> kernelOperands{op.getArg(), udfFnPtr};
+        std::vector<Value> kernelOperands{op.getArg(), udfFnPtr, op.getAxis(), op.getUdfReturnsScalar()};
 
         auto kernel = rewriter.create<daphne::CallKernelOp>(loc, callee.str(), kernelOperands, op->getResultTypes());
         rewriter.replaceOp(op, kernel.getResults());
diff --git a/src/compiler/lowering/SpecializeGenericFunctionsPass.cpp b/src/compiler/lowering/SpecializeGenericFunctionsPass.cpp
index e6b98bc8a..a085cb4e1 100644
--- a/src/compiler/lowering/SpecializeGenericFunctionsPass.cpp
+++ b/src/compiler/lowering/SpecializeGenericFunctionsPass.cpp
@@ -369,8 +369,19 @@ class SpecializeGenericFunctionsPass : public PassWrapper<SpecializeGenericFunct
                 // mapped on
                 mlir::Type opTy = mapOp.getArg().getType();
                 auto inpMatrixTy = opTy.dyn_cast<daphne::MatrixType>();
-                func::FuncOp specializedFunc =
-                    createOrReuseSpecialization(inpMatrixTy.getElementType(), {}, calledFunction, mapOp.getLoc());
+                int64_t axis = CompilerUtils::constantOrThrow<int64_t>(mapOp.getAxis(), "map axis must be a constant.");
+                func::FuncOp specializedFunc;
+                // Set function input type based on given axis
+                if (axis == 0) { // row-wise map
+                    specializedFunc =
+                        createOrReuseSpecialization(inpMatrixTy.withShape(1, -1), {}, calledFunction, mapOp.getLoc());
+                } else if (axis == 1) { // column-wise map
+                    specializedFunc =
+                        createOrReuseSpecialization(inpMatrixTy.withShape(-1, 1), {}, calledFunction, mapOp.getLoc());
+                } else { // element-wise
+                    specializedFunc =
+                        createOrReuseSpecialization(inpMatrixTy.getElementType(), {}, calledFunction, mapOp.getLoc());
+                }
                 mapOp.setFuncAttr(specializedFunc.getSymNameAttr());
 
                 // We only allow functions that return exactly one result for
@@ -385,16 +396,45 @@ class SpecializeGenericFunctionsPass : public PassWrapper<SpecializeGenericFunct
 
                 // Get current mapOp result matrix type and fix it if needed.
                 // If we fixed something we rerun inference of the whole
-                // function
+                // function.
                 daphne::MatrixType resMatrixTy = mapOp.getType().dyn_cast<daphne::MatrixType>();
                 mlir::Type funcResTy = specializedFunc.getFunctionType().getResult(0);
+                bool madeChanges = false;
+
+                auto udfReturnsScalar = CompilerUtils::constantOrThrow<bool>(
+                    mapOp.getUdfReturnsScalar(), "map parameter udfReturnsScalar must be a bool.");
+
+                // If the specialized function returns a scalar, the previously
+                // unknown dimension is set to one, and if the specialized function
+                // returns a matrix, this dimension is still unknown.
+                if (axis == 0 || axis == 1) {
+                    if (dyn_cast_or_null<daphne::MatrixType>(funcResTy) && !udfReturnsScalar) { // Matrix -> Matrix
+                        // Set function result type to the matrix's element
+                        // type for further processing
+                        funcResTy = dyn_cast<daphne::MatrixType>(funcResTy).getElementType();
+                    } else if (!dyn_cast_or_null<daphne::MatrixType>(funcResTy) &&
+                               udfReturnsScalar) { // Matrix -> Scalar
+                        if (axis == 0)
+                            resMatrixTy = resMatrixTy.withShape(inpMatrixTy.getNumRows(), 1);
+                        else if (axis == 1)
+                            resMatrixTy = resMatrixTy.withShape(1, inpMatrixTy.getNumCols());
+                        madeChanges = true;
+                    } else { // udfReturnsScalar does not match funcResTy
+                        throw ErrorHandler::compilerError(
+                            mapOp.getOperation(), "SpecializeGenericFunctionsPass",
+                            "map parameter udfReturnsScalar does not match the output type of the provided function.");
+                    }
+                }
+
+                // The matrix that results from the mapOp has the same
+                // element-type returned by the specialized function
+                if (resMatrixTy.getElementType() != funcResTy) {
+                    resMatrixTy = resMatrixTy.withElementType(funcResTy);
+                    madeChanges = true;
+                }
 
-                // The matrix that results from the mapOp has the same dimension
-                // as the input matrix and the element-type returned by the
-                // specialized function
-                if (resMatrixTy.getNumCols() != inpMatrixTy.getNumCols() ||
-                    resMatrixTy.getNumRows() != inpMatrixTy.getNumRows() || resMatrixTy.getElementType() != funcResTy) {
-                    mapOp.getResult().setType(inpMatrixTy.withElementType(funcResTy));
+                if (madeChanges) {
+                    mapOp.getResult().setType(resMatrixTy);
                     inferTypesInFunction(function);
                 }
 
diff --git a/src/ir/daphneir/DaphneInferShapeOpInterface.cpp b/src/ir/daphneir/DaphneInferShapeOpInterface.cpp
index c12aac089..1018bfba7 100644
--- a/src/ir/daphneir/DaphneInferShapeOpInterface.cpp
+++ b/src/ir/daphneir/DaphneInferShapeOpInterface.cpp
@@ -568,6 +568,30 @@ std::vector<std::pair<ssize_t, ssize_t>> daphne::RecodeOp::inferShape() {
     return {{resNumRows, resNumCols}, {dictNumRows, dictNumCols}};
 }
 
+std::vector<std::pair<ssize_t, ssize_t>> daphne::MapOp::inferShape() {
+    mlir::Type opTy = getArg().getType();
+    auto inpMatrixTy = opTy.dyn_cast<daphne::MatrixType>();
+
+    // For element-wise mapOp, the result matrix has the same
+    // dimension as the input matrix
+    ssize_t resNumRows = inpMatrixTy.getNumRows();
+    ssize_t resNumCols = inpMatrixTy.getNumCols();
+
+    int64_t axis = CompilerUtils::constantOrThrow<int64_t>(getAxis(), "map axis must be a constant.");
+
+    // For row- and column-wise mapOp, the result matrix does not
+    // have the same dimensions as the input matrix.
+    // During a row-wise mapOp the number of rows stays the same
+    // as the input matrix, and during a column-wise mapOp the
+    // number of columns stays the same.
+    if (axis == 0)
+        resNumCols = -1;
+    else if (axis == 1)
+        resNumRows = -1;
+
+    return {{resNumRows, resNumCols}};
+}
+
 // ****************************************************************************
 // Shape inference trait implementations
 // ****************************************************************************
diff --git a/src/ir/daphneir/DaphneOps.td b/src/ir/daphneir/DaphneOps.td
index c8c4958d2..4cc9f3549 100644
--- a/src/ir/daphneir/DaphneOps.td
+++ b/src/ir/daphneir/DaphneOps.td
@@ -1918,9 +1918,9 @@ def Daphne_DistributedPipelineOp : Daphne_Op<"distributedPipeline", [AttrSizedOp
 // ****************************************************************************
 // Higher-order operations
 // ****************************************************************************
-def Daphne_MapOp : Daphne_Op<"map", [ShapeFromArg]> {
+def Daphne_MapOp : Daphne_Op<"map", [DeclareOpInterfaceMethods<InferShapeOpInterface>]> {
     let summary = "Applies a user defined function to elements of a matrix.";
-    let arguments = (ins MatrixOrU:$arg, SymbolNameAttr:$func);
+    let arguments = (ins MatrixOrU:$arg, SymbolNameAttr:$func, SI64:$axis, BoolScalar:$udfReturnsScalar);
     let results = (outs MatrixOrU:$res);
 }
 
diff --git a/src/parser/daphnedsl/DaphneDSLBuiltins.cpp b/src/parser/daphnedsl/DaphneDSLBuiltins.cpp
index a2577995a..38c0670d6 100644
--- a/src/parser/daphnedsl/DaphneDSLBuiltins.cpp
+++ b/src/parser/daphnedsl/DaphneDSLBuiltins.cpp
@@ -1299,14 +1299,32 @@ antlrcpp::Any DaphneDSLBuiltins::build(mlir::Location loc, const std::string &fu
     // ****************************************************************************
 
     if (func == "map") {
-        checkNumArgsExact(loc, func, numArgs, 2);
+        checkNumArgsBetween(loc, func, numArgs, 2, 4);
+
         mlir::Value source = args[0];
 
         auto co = args[1].getDefiningOp<mlir::daphne::ConstantOp>();
         mlir::Attribute attr = co.getValue();
 
-        return static_cast<mlir::Value>(
-            builder.create<MapOp>(loc, source.getType(), source, attr.dyn_cast<mlir::StringAttr>()));
+        // Default values, if not given
+        mlir::Value axis = builder.create<ConstantOp>(loc, int64_t(-1));
+        mlir::Value udfReturnsScalar = builder.create<ConstantOp>(loc, false);
+
+        if (numArgs >= 3) { // axis is given
+            int64_t axisInt =
+                CompilerUtils::constantOrThrow<int64_t>(args[2], "third argument of map must be a constant");
+            if (axisInt == 0 || axisInt == 1)
+                axis = args[2];
+            else
+                throw ErrorHandler::compilerError(loc, "DSLBuiltins", "invalid axis for aggregation.");
+        }
+
+        if (numArgs == 4) { // udfReturnsScalar is given
+            udfReturnsScalar = args[3];
+        }
+
+        return static_cast<mlir::Value>(builder.create<MapOp>(
+            loc, source.getType(), source, attr.dyn_cast<mlir::StringAttr>(), axis, udfReturnsScalar));
     }
 
     // ****************************************************************************
diff --git a/src/parser/daphnedsl/DaphneDSLVisitor.cpp b/src/parser/daphnedsl/DaphneDSLVisitor.cpp
index c3dc2e628..cbadc1ae2 100644
--- a/src/parser/daphnedsl/DaphneDSLVisitor.cpp
+++ b/src/parser/daphnedsl/DaphneDSLVisitor.cpp
@@ -892,9 +892,9 @@ antlrcpp::Any DaphneDSLVisitor::handleMapOpCall(DaphneDSLGrammarParser::CallExpr
         throw ErrorHandler::compilerError(loc, "DSLVisitor",
                                           "called 'handleMapOpCall' for function " + func + " instead of 'map'");
 
-    if (ctx->expr().size() != 2) {
+    if (ctx->expr().size() < 2 || ctx->expr().size() > 4) {
         throw ErrorHandler::compilerError(loc, "DSLVisitor",
-                                          "built-in function 'map' expects exactly 2 argument(s), but got " +
+                                          "built-in function 'map' expects 2-4 argument(s), but got " +
                                               std::to_string(ctx->expr().size()));
     }
 
@@ -918,6 +918,16 @@ antlrcpp::Any DaphneDSLVisitor::handleMapOpCall(DaphneDSLGrammarParser::CallExpr
     args.push_back(
         static_cast<mlir::Value>(builder.create<mlir::daphne::ConstantOp>(loc, maybeUDF->getSymName().str())));
 
+    if (ctx->expr().size() >= 3) {
+        auto axis = valueOrErrorOnVisit(ctx->expr(2));
+        args.push_back(axis);
+
+        if (ctx->expr().size() == 4) {
+            auto udfReturnsScalar = utils.castBoolIf(valueOrErrorOnVisit(ctx->expr(3)));
+            args.push_back(udfReturnsScalar);
+        }
+    }
+
     // Create DaphneIR operation for the built-in function.
     return builtins.build(loc, func, args);
 }
diff --git a/src/runtime/local/kernels/Map.h b/src/runtime/local/kernels/Map.h
index 8f2382353..5f9f29c6c 100644
--- a/src/runtime/local/kernels/Map.h
+++ b/src/runtime/local/kernels/Map.h
@@ -29,18 +29,17 @@
 // ****************************************************************************
 
 template <class DTRes, class DTArg> struct Map {
-    // We could have a more specialized function pointer here i.e.
-    // (DTRes::VT)(*func)(DTArg::VT). The problem is that this is currently not
-    // supported by kernels.json.
-    static void apply(DTRes *&res, const DTArg *arg, void *func, DCTX(ctx)) = delete;
+    static void apply(DTRes *&res, const DTArg *arg, void *func, const int64_t axis, const bool udfReturnsScalar,
+                      DCTX(ctx)) = delete;
 };
 
 // ****************************************************************************
 // Convenience function
 // ****************************************************************************
 
-template <class DTRes, class DTArg> void map(DTRes *&res, const DTArg *arg, void *func, DCTX(ctx)) {
-    Map<DTRes, DTArg>::apply(res, arg, func, ctx);
+template <class DTRes, class DTArg>
+void map(DTRes *&res, const DTArg *arg, void *func, const int64_t axis, const bool udfReturnsScalar, DCTX(ctx)) {
+    Map<DTRes, DTArg>::apply(res, arg, func, axis, udfReturnsScalar, ctx);
 }
 
 // ****************************************************************************
@@ -52,23 +51,95 @@ template <class DTRes, class DTArg> void map(DTRes *&res, const DTArg *arg, void
 // ----------------------------------------------------------------------------
 
 template <typename VTRes, typename VTArg> struct Map<DenseMatrix<VTRes>, DenseMatrix<VTArg>> {
-    static void apply(DenseMatrix<VTRes> *&res, const DenseMatrix<VTArg> *arg, void *func, DCTX(ctx)) {
+    static void apply(DenseMatrix<VTRes> *&res, const DenseMatrix<VTArg> *arg, void *func, const int64_t axis,
+                      const bool udfReturnsScalar, DCTX(ctx)) {
         const size_t numRows = arg->getNumRows();
         const size_t numCols = arg->getNumCols();
 
-        if (res == nullptr)
-            res = DataObjectFactory::create<DenseMatrix<VTRes>>(numRows, numCols, false);
-
-        auto udf = reinterpret_cast<VTRes (*)(VTArg)>(func);
+        VTRes *valuesRes = nullptr;
 
-        const VTArg *valuesArg = arg->getValues();
-        VTRes *valuesRes = res->getValues();
+        if (axis != 0 && axis != 1) { // element-wise
+            res = DataObjectFactory::create<DenseMatrix<VTRes>>(numRows, numCols, false);
+            valuesRes = res->getValues();
+        }
 
-        for (size_t r = 0; r < numRows; r++) {
-            for (size_t c = 0; c < numCols; c++)
-                valuesRes[c] = udf(valuesArg[c]);
-            valuesArg += arg->getRowSkip();
-            valuesRes += res->getRowSkip();
+        auto udfElem = reinterpret_cast<VTRes (*)(VTArg)>(func);
+        auto udfMatMat = reinterpret_cast<DenseMatrix<VTRes> *(*)(DenseMatrix<VTArg> *)>(func);
+        auto udfMatElem = reinterpret_cast<VTRes (*)(DenseMatrix<VTArg> *)>(func);
+
+        if (axis == 1) { // column-wise
+            size_t resNumRows = 1;
+            for (size_t c = 0; c < numCols; ++c) {
+                // Extract current column
+                DenseMatrix<VTArg> *currentCol =
+                    DataObjectFactory::create<DenseMatrix<VTArg>>(arg, 0, numRows, c, c + 1);
+                if (!udfReturnsScalar) { // Default: Matrix -> Matrix
+                    // Apply UDF and check shape of resulting column
+                    const DenseMatrix<VTRes> *resCol = udfMatMat(currentCol);
+                    if (resCol == nullptr || resCol->getNumCols() != 1)
+                        throw std::runtime_error("Matrix returned by UDF should be a single column!");
+                    // Set result matrix size in first iteration
+                    if (c == 0) {
+                        resNumRows = resCol->getNumRows();
+                        res = DataObjectFactory::create<DenseMatrix<VTRes>>(resNumRows, numCols, false);
+                    }
+                    // Set result row-wise
+                    valuesRes = res->getValues();
+                    const VTRes *valuesResCol = resCol->getValues();
+                    for (size_t r = 0; r < resNumRows; ++r) {
+                        valuesRes[c] = valuesResCol[0];
+                        valuesResCol += resCol->getRowSkip();
+                        valuesRes += res->getRowSkip();
+                    }
+                } else { // Matrix -> Scalar
+                    // Set result matrix size in first iteration
+                    if (c == 0) {
+                        res = DataObjectFactory::create<DenseMatrix<VTRes>>(resNumRows, numCols, false);
+                        valuesRes = res->getValues();
+                    }
+                    valuesRes[c] = udfMatElem(currentCol);
+                }
+                // TODO(#XXX): Call to destroy leads to res becoming a view -> segfault in print
+                // DataObjectFactory::destroy(currentCol);
+            }
+        } else { // row-wise or element-wise
+            const VTArg *valuesArg = arg->getValues();
+            size_t resNumCols = 1;
+            for (size_t r = 0; r < numRows; ++r) {
+                if (axis == 0) {
+                    // Extract current row
+                    DenseMatrix<VTArg> *currentRow = DataObjectFactory::create<DenseMatrix<VTArg>>(arg, r, r + 1);
+                    if (!udfReturnsScalar) { // Default: Matrix -> Matrix
+                        // Apply UDF and check shape of resulting row
+                        const DenseMatrix<VTRes> *resRow = udfMatMat(currentRow);
+                        if (resRow == nullptr || resRow->getNumRows() != 1)
+                            throw std::runtime_error("Matrix returned by UDF should be a single row!");
+                        // Set result matrix size in first iteration
+                        if (r == 0) {
+                            resNumCols = resRow->getNumCols();
+                            res = DataObjectFactory::create<DenseMatrix<VTRes>>(numRows, resNumCols, false);
+                            valuesRes = res->getValues();
+                        }
+                        // Set result by copying values
+                        const VTRes *valuesResRow = resRow->getValues();
+                        memcpy(valuesRes, valuesResRow, resNumCols * sizeof(VTRes));
+                    } else { // Matrix -> Scalar
+                        // Set result matrix size in first iteration
+                        if (r == 0) {
+                            res = DataObjectFactory::create<DenseMatrix<VTRes>>(numRows, resNumCols, false);
+                            valuesRes = res->getValues();
+                        }
+                        valuesRes[0] = udfMatElem(currentRow);
+                    }
+                    // TODO(#XXX): Call to destroy leads to res becoming a view -> segfault in print
+                    // DataObjectFactory::destroy(currentRow);
+                } else {
+                    for (size_t c = 0; c < numCols; ++c)
+                        valuesRes[c] = udfElem(valuesArg[c]);
+                    valuesArg += arg->getRowSkip();
+                }
+                valuesRes += res->getRowSkip();
+            }
         }
     }
 };
@@ -78,19 +149,91 @@ template <typename VTRes, typename VTArg> struct Map<DenseMatrix<VTRes>, DenseMa
 // ----------------------------------------------------------------------------
 
 template <typename VTRes, typename VTArg> struct Map<Matrix<VTRes>, Matrix<VTArg>> {
-    static void apply(Matrix<VTRes> *&res, const Matrix<VTArg> *arg, void *func, DCTX(ctx)) {
+    static void apply(Matrix<VTRes> *&res, const Matrix<VTArg> *arg, void *func, const int64_t axis,
+                      const bool udfReturnsScalar, DCTX(ctx)) {
         const size_t numRows = arg->getNumRows();
         const size_t numCols = arg->getNumCols();
 
-        if (res == nullptr)
+        if (axis != 0 && axis != 1) { // element-wise
             res = DataObjectFactory::create<DenseMatrix<VTRes>>(numRows, numCols, false);
+            res->prepareAppend();
+        }
 
-        auto udf = reinterpret_cast<VTRes (*)(VTArg)>(func);
-
-        res->prepareAppend();
-        for (size_t r = 0; r < numRows; ++r)
-            for (size_t c = 0; c < numCols; ++c)
-                res->append(r, c, udf(arg->get(r, c)));
+        auto udfElem = reinterpret_cast<VTRes (*)(VTArg)>(func);
+        auto udfMatMat = reinterpret_cast<Matrix<VTRes> *(*)(Matrix<VTArg> *)>(func);
+        auto udfMatElem = reinterpret_cast<VTRes (*)(Matrix<VTArg> *)>(func);
+
+        if (axis == 1) { // column-wise
+            size_t resNumRows = 1;
+            for (size_t c = 0; c < numCols; ++c) {
+                // Extract current column
+                Matrix<VTArg> *currentCol = DataObjectFactory::create<DenseMatrix<VTArg>>(
+                    dynamic_cast<const DenseMatrix<VTArg> *>(arg), 0, numRows, c, c + 1);
+                if (!udfReturnsScalar) { // Default: Matrix -> Matrix
+                    // TODO(#XXX): Append works only on strictly
+                    // increasing coordinates, therefore not supported
+                    throw std::runtime_error("Currently not supported!");
+                    // Apply UDF and check shape of resulting column
+                    const Matrix<VTRes> *resCol = udfMatMat(currentCol);
+                    if (resCol == nullptr || resCol->getNumCols() != 1)
+                        throw std::runtime_error("Matrix returned by UDF should be a single column!");
+                    // Set result matrix size in first iteration
+                    if (c == 0) {
+                        resNumRows = resCol->getNumRows();
+                        res = DataObjectFactory::create<DenseMatrix<VTRes>>(resNumRows, numCols, false);
+                        res->prepareAppend();
+                    }
+                    // Set result row-wise
+                    for (size_t r = 0; r < resNumRows; ++r)
+                        res->append(r, c, resCol->get(r, 0));
+                } else { // Matrix -> Scalar
+                    // Set result matrix size in first iteration
+                    if (c == 0) {
+                        res = DataObjectFactory::create<DenseMatrix<VTRes>>(resNumRows, numCols, false);
+                        res->prepareAppend();
+                    }
+                    res->append(0, c, udfMatElem(currentCol));
+                }
+                // TODO(#XXX): Call to destroy leads to res becoming a view -> segfault in print
+                // DataObjectFactory::destroy(currentCol);
+            }
+        } else { // row-wise or element-wise
+            size_t resNumCols = 1;
+            for (size_t r = 0; r < numRows; ++r) {
+                if (axis == 0) {
+                    // Extract current row
+                    Matrix<VTArg> *currentRow = DataObjectFactory::create<DenseMatrix<VTArg>>(
+                        dynamic_cast<const DenseMatrix<VTArg> *>(arg), r, r + 1);
+                    if (!udfReturnsScalar) { // Default: Matrix -> Matrix
+                        // Apply UDF and check shape of resulting row
+                        const Matrix<VTRes> *resRow = udfMatMat(currentRow);
+                        if (resRow == nullptr || resRow->getNumRows() != 1)
+                            throw std::runtime_error("Matrix returned by UDF should be a single row!");
+                        // Set result matrix size in first iteration
+                        if (r == 0) {
+                            resNumCols = resRow->getNumCols();
+                            res = DataObjectFactory::create<DenseMatrix<VTRes>>(numRows, resNumCols, false);
+                            res->prepareAppend();
+                        }
+                        // Set result by column-wise
+                        for (size_t c = 0; c < resNumCols; ++c)
+                            res->append(r, c, resRow->get(0, c));
+                    } else { // Matrix -> Scalar
+                        // Set result matrix size in first iteration
+                        if (r == 0) {
+                            res = DataObjectFactory::create<DenseMatrix<VTRes>>(numRows, resNumCols, false);
+                            res->prepareAppend();
+                        }
+                        res->append(r, 0, udfMatElem(currentRow));
+                    }
+                    // TODO(#XXX): Call to destroy leads to res becoming a view -> segfault in print
+                    // DataObjectFactory::destroy(currentRow);
+                } else {
+                    for (size_t c = 0; c < numCols; ++c)
+                        res->append(r, c, udfElem(arg->get(r, c)));
+                }
+            }
+        }
         res->finishAppend();
     }
 };
diff --git a/src/runtime/local/kernels/kernels.json b/src/runtime/local/kernels/kernels.json
index c663cc499..b97c54933 100644
--- a/src/runtime/local/kernels/kernels.json
+++ b/src/runtime/local/kernels/kernels.json
@@ -3183,6 +3183,14 @@
                 {
                     "type": "void *",
                     "name": "func"
+                },
+                {
+                    "type": "const int64_t",
+                    "name": "axis"
+                },
+                {
+                    "type": "const bool",
+                    "name": "udfReturnsScalar"
                 }
             ]
         },
@@ -3219,6 +3227,10 @@
                 ["DenseMatrix", "double"],
                 ["DenseMatrix", "uint8_t"]
             ],
+            [
+                ["DenseMatrix", "double"],
+                ["DenseMatrix", "bool"]
+            ],
 
             [
                 ["DenseMatrix", "float"],
@@ -3252,6 +3264,10 @@
                 ["DenseMatrix", "float"],
                 ["DenseMatrix", "uint8_t"]
             ],
+            [
+                ["DenseMatrix", "float"],
+                ["DenseMatrix", "bool"]
+            ],
 
             [
                 ["DenseMatrix", "int64_t"],
@@ -3285,6 +3301,10 @@
                 ["DenseMatrix", "int64_t"],
                 ["DenseMatrix", "uint8_t"]
             ],
+            [
+                ["DenseMatrix", "int64_t"],
+                ["DenseMatrix", "bool"]
+            ],
 
             [
                 ["DenseMatrix", "int32_t"],
@@ -3318,6 +3338,10 @@
                 ["DenseMatrix", "int32_t"],
                 ["DenseMatrix", "uint8_t"]
             ],
+            [
+                ["DenseMatrix", "int32_t"],
+                ["DenseMatrix", "bool"]
+            ],
 
             [
                 ["DenseMatrix", "int8_t"],
@@ -3351,6 +3375,10 @@
                 ["DenseMatrix", "int8_t"],
                 ["DenseMatrix", "uint8_t"]
             ],
+            [
+                ["DenseMatrix", "int8_t"],
+                ["DenseMatrix", "bool"]
+            ],
 
             [
                 ["DenseMatrix", "uint64_t"],
@@ -3384,6 +3412,10 @@
                 ["DenseMatrix", "uint64_t"],
                 ["DenseMatrix", "uint8_t"]
             ],
+            [
+                ["DenseMatrix", "uint64_t"],
+                ["DenseMatrix", "bool"]
+            ],
 
             [
                 ["DenseMatrix", "uint32_t"],
@@ -3417,6 +3449,10 @@
                 ["DenseMatrix", "uint32_t"],
                 ["DenseMatrix", "uint8_t"]
             ],
+            [
+                ["DenseMatrix", "uint32_t"],
+                ["DenseMatrix", "bool"]
+            ],
 
             [
                 ["DenseMatrix", "uint8_t"],
@@ -3449,6 +3485,47 @@
             [
                 ["DenseMatrix", "uint8_t"],
                 ["DenseMatrix", "uint8_t"]
+            ],
+            [
+                ["DenseMatrix", "uint8_t"],
+                ["DenseMatrix", "bool"]
+            ],
+
+            [
+                ["DenseMatrix", "bool"],
+                ["DenseMatrix", "double"]
+            ],
+            [
+                ["DenseMatrix", "bool"],
+                ["DenseMatrix", "float"]
+            ],
+            [
+                ["DenseMatrix", "bool"],
+                ["DenseMatrix", "int64_t"]
+            ],
+            [
+                ["DenseMatrix", "bool"],
+                ["DenseMatrix", "int32_t"]
+            ],
+            [
+                ["DenseMatrix", "bool"],
+                ["DenseMatrix", "int8_t"]
+            ],
+            [
+                ["DenseMatrix", "bool"],
+                ["DenseMatrix", "uint64_t"]
+            ],
+            [
+                ["DenseMatrix", "bool"],
+                ["DenseMatrix", "uint32_t"]
+            ],
+            [
+                ["DenseMatrix", "bool"],
+                ["DenseMatrix", "uint8_t"]
+            ],
+            [
+                ["DenseMatrix", "bool"],
+                ["DenseMatrix", "bool"]
             ]
         ]
     },
diff --git a/test/api/cli/secondorder/SecondOrderTest.cpp b/test/api/cli/secondorder/SecondOrderTest.cpp
index 4bd482f9a..ba73641cf 100644
--- a/test/api/cli/secondorder/SecondOrderTest.cpp
+++ b/test/api/cli/secondorder/SecondOrderTest.cpp
@@ -41,5 +41,5 @@ const std::string dirPath = "test/api/cli/secondorder/";
         }                                                                                                              \
     }
 
-MAKE_TEST_CASE("map", 3)
-MAKE_FAILURE_TEST_CASE("map", 5)
\ No newline at end of file
+MAKE_TEST_CASE("map", 7)
+MAKE_FAILURE_TEST_CASE("map", 8)
\ No newline at end of file
diff --git a/test/api/cli/secondorder/map_1.daphne b/test/api/cli/secondorder/map_1.daphne
index ff2f90ba2..316d495d0 100644
--- a/test/api/cli/secondorder/map_1.daphne
+++ b/test/api/cli/secondorder/map_1.daphne
@@ -18,4 +18,12 @@ print(map(X, timesThree));
 
 // matrix
 X = reshape([1, 0, 5, -2], 2, 2);
-print(map(X, timesThree));
\ No newline at end of file
+print(map(X, timesThree));
+
+// row-wise
+X = reshape([1, 0, 5, -2], 2, 2);
+print(map(X, timesThree, 0));
+
+// column-wise
+X = reshape([1, 0, 5, -2], 2, 2);
+print(map(X, timesThree, 1));
\ No newline at end of file
diff --git a/test/api/cli/secondorder/map_1.txt b/test/api/cli/secondorder/map_1.txt
index 1110f0dfd..0fb488e19 100644
--- a/test/api/cli/secondorder/map_1.txt
+++ b/test/api/cli/secondorder/map_1.txt
@@ -10,3 +10,9 @@ DenseMatrix(1x4, int64_t)
 DenseMatrix(2x2, int64_t)
 3 0
 15 -6
+DenseMatrix(2x2, int64_t)
+3 0
+15 -6
+DenseMatrix(2x2, int64_t)
+3 0
+15 -6
diff --git a/test/api/cli/secondorder/map_2.daphne b/test/api/cli/secondorder/map_2.daphne
index ead65c580..ad1cb9715 100644
--- a/test/api/cli/secondorder/map_2.daphne
+++ b/test/api/cli/secondorder/map_2.daphne
@@ -18,4 +18,12 @@ print(map(X, timesThree));
 
 // matrix
 X = reshape([1.1, 0.0, 5.5, -2.2], 2, 2);
-print(map(X, timesThree));
\ No newline at end of file
+print(map(X, timesThree));
+
+// row-wise
+X = reshape([1.1, 0.0, 5.5, -2.2], 2, 2);
+print(map(X, timesThree, 0));
+
+// column-wise
+X = reshape([1.1, 0.0, 5.5, -2.2], 2, 2);
+print(map(X, timesThree, 1));
\ No newline at end of file
diff --git a/test/api/cli/secondorder/map_2.txt b/test/api/cli/secondorder/map_2.txt
index 8ba6ede95..2d275e7ec 100644
--- a/test/api/cli/secondorder/map_2.txt
+++ b/test/api/cli/secondorder/map_2.txt
@@ -10,3 +10,9 @@ DenseMatrix(1x4, double)
 DenseMatrix(2x2, double)
 3.3 0
 16.5 -6.6
+DenseMatrix(2x2, double)
+3.3 0
+16.5 -6.6
+DenseMatrix(2x2, double)
+3.3 0
+16.5 -6.6
diff --git a/test/api/cli/secondorder/map_3.daphne b/test/api/cli/secondorder/map_3.daphne
index 58bb551e6..493b7c64f 100644
--- a/test/api/cli/secondorder/map_3.daphne
+++ b/test/api/cli/secondorder/map_3.daphne
@@ -4,5 +4,14 @@ def timesThree(x) {
     return x * 3;
 }
 
+// element-wise
 X = fill(2, 1000, 1000);
-print(sum(map(X, timesThree)));
\ No newline at end of file
+print(sum(map(X, timesThree)));
+
+// row-wise
+X = fill(2, 1000, 1000);
+print(sum(map(X, timesThree, 0)));
+
+// column-wise
+X = fill(2, 1000, 1000);
+print(sum(map(X, timesThree, 0)));
\ No newline at end of file
diff --git a/test/api/cli/secondorder/map_3.txt b/test/api/cli/secondorder/map_3.txt
index f348d457b..20a8c6e32 100644
--- a/test/api/cli/secondorder/map_3.txt
+++ b/test/api/cli/secondorder/map_3.txt
@@ -1 +1,3 @@
 6000000
+6000000
+6000000
diff --git a/test/api/cli/secondorder/map_4.daphne b/test/api/cli/secondorder/map_4.daphne
new file mode 100644
index 000000000..8d4795b60
--- /dev/null
+++ b/test/api/cli/secondorder/map_4.daphne
@@ -0,0 +1,21 @@
+// Row-wise (Matrix -> Scalar)
+
+def sumOfFirstAndLast(row) {
+  return as.scalar<f64>(row[0, 0]) + as.scalar<f64>(row[0, ncol(row) - 1]);
+}
+
+// double
+X = reshape([1.1, 0.0, 3.0, 5.5, -2.2, -1.2], 2, 3);
+print(map(X, sumOfFirstAndLast, 0, true));
+
+// int64_t
+X = reshape([1, 2, 3, 4, 5, 6, 7, 8], 2, 4);
+print(map(X, sumOfFirstAndLast, 0, true));
+
+// one row
+X = reshape([1, 2, 3, 4, 5], 1, 5);
+print(map(X, sumOfFirstAndLast, 0, true));
+
+// one column
+X = reshape([1, 2, 3, 4, 5], 5, 1);
+print(map(X, sumOfFirstAndLast, 0, true));
\ No newline at end of file
diff --git a/test/api/cli/secondorder/map_4.txt b/test/api/cli/secondorder/map_4.txt
new file mode 100644
index 000000000..7ea0b55e4
--- /dev/null
+++ b/test/api/cli/secondorder/map_4.txt
@@ -0,0 +1,14 @@
+DenseMatrix(2x1, double)
+4.1
+4.3
+DenseMatrix(2x1, double)
+5
+13
+DenseMatrix(1x1, double)
+6
+DenseMatrix(5x1, double)
+2
+4
+6
+8
+10
diff --git a/test/api/cli/secondorder/map_5.daphne b/test/api/cli/secondorder/map_5.daphne
new file mode 100644
index 000000000..44f7d3012
--- /dev/null
+++ b/test/api/cli/secondorder/map_5.daphne
@@ -0,0 +1,21 @@
+// Row-wise (Matrix -> Matrix)
+
+def appendFirst(row) {
+  return cbind(row, row[0, 0]);
+}
+
+// double
+X = reshape([1.1, 0.0, 3.0, 5.5, -2.2, -1.2], 2, 3);
+print(map(X, appendFirst, 0));
+
+// int64_t
+X = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], 3, 4);
+print(map(X, appendFirst, 0));
+
+// one row
+X = reshape([1, 2, 3, 4, 5], 1, 5);
+print(map(X, appendFirst, 0));
+
+// one column
+X = reshape([1, 2, 3, 4, 5], 5, 1);
+print(map(X, appendFirst, 0));
\ No newline at end of file
diff --git a/test/api/cli/secondorder/map_5.txt b/test/api/cli/secondorder/map_5.txt
new file mode 100644
index 000000000..e617e09a6
--- /dev/null
+++ b/test/api/cli/secondorder/map_5.txt
@@ -0,0 +1,15 @@
+DenseMatrix(2x4, double)
+1.1 0 3 1.1
+5.5 -2.2 -1.2 5.5
+DenseMatrix(3x5, int64_t)
+1 2 3 4 1
+5 6 7 8 5
+9 10 11 12 9
+DenseMatrix(1x6, int64_t)
+1 2 3 4 5 1
+DenseMatrix(5x2, int64_t)
+1 1
+2 2
+3 3
+4 4
+5 5
diff --git a/test/api/cli/secondorder/map_6.daphne b/test/api/cli/secondorder/map_6.daphne
new file mode 100644
index 000000000..b1b2de5fc
--- /dev/null
+++ b/test/api/cli/secondorder/map_6.daphne
@@ -0,0 +1,17 @@
+// Column-wise (Matrix -> Matrix)
+
+def topThree(col) {
+  return order(col[:3, 0], 0, true, false);
+}
+
+// double
+X = reshape([1.1, 1.1, 1.1, 1.1, 0.3, 0.3, 0.3, 0.3, -2.2, -2.2, -2.2, -2.2, 5.5, 5.5, 5.5, 5.5], 4, 4);
+print(map(X, topThree, 1));
+
+// int64_t
+X = reshape([4, 3, 2, 1, 0, 4, 3, 2, 1, 0, 4, 3, 2, 1, 0, 4, 3, 2, 1, 0], 5, 4);
+print(map(X, topThree, 1));
+
+// one column
+X = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 10, 1);
+print(map(X, topThree, 1));
\ No newline at end of file
diff --git a/test/api/cli/secondorder/map_6.txt b/test/api/cli/secondorder/map_6.txt
new file mode 100644
index 000000000..40463e9b4
--- /dev/null
+++ b/test/api/cli/secondorder/map_6.txt
@@ -0,0 +1,12 @@
+DenseMatrix(3x4, double)
+-2.2 -2.2 -2.2 -2.2
+0.3 0.3 0.3 0.3
+1.1 1.1 1.1 1.1
+DenseMatrix(3x4, int64_t)
+0 0 2 1
+1 3 3 2
+4 4 4 3
+DenseMatrix(3x1, int64_t)
+1
+2
+3
diff --git a/test/api/cli/secondorder/map_7.daphne b/test/api/cli/secondorder/map_7.daphne
new file mode 100644
index 000000000..1d2705430
--- /dev/null
+++ b/test/api/cli/secondorder/map_7.daphne
@@ -0,0 +1,21 @@
+// Column-wise (Matrix -> Scalar)
+
+def isSumEven(col) {
+  return as.bool(as.scalar<si64>(sum(col) % 2 == 0));
+}
+
+// double
+X = reshape([1.1, 1.1, 1.1, 1.1, 0.3, 0.3, 0.3, 0.3, -2.2, -2.2, -2.2, -2.2, 5.5, 5.5, 5.5, 5.5], 4, 4);
+print(map(X, isSumEven, 1, true));
+
+// int64_t
+X = reshape([0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12], 3, 4);
+print(map(X, isSumEven, 1, true));
+
+// one column
+X = reshape([0, 1, 2], 3, 1);
+print(map(X, isSumEven, 1, true));
+
+// one row
+X = reshape([0, 1, 2], 1, 3);
+print(map(X, isSumEven, 1, true));
\ No newline at end of file
diff --git a/test/api/cli/secondorder/map_7.txt b/test/api/cli/secondorder/map_7.txt
new file mode 100644
index 000000000..625cc0eef
--- /dev/null
+++ b/test/api/cli/secondorder/map_7.txt
@@ -0,0 +1,8 @@
+DenseMatrix(1x4, bool)
+0 0 0 0
+DenseMatrix(1x4, bool)
+1 1 0 1
+DenseMatrix(1x1, bool)
+0
+DenseMatrix(1x3, bool)
+1 0 1
diff --git a/test/api/cli/secondorder/map_failure_6.daphne b/test/api/cli/secondorder/map_failure_6.daphne
new file mode 100644
index 000000000..2a91fa965
--- /dev/null
+++ b/test/api/cli/secondorder/map_failure_6.daphne
@@ -0,0 +1,8 @@
+// udfReturnsScalar doesn't match actual UDF result
+
+def sumOfFirstAndLast(row) {
+  return as.scalar<f64>(row[0, 0]) + as.scalar<f64>(row[0, ncol(row) - 1]);
+}
+
+X = reshape([1, 2, 3], 1, 3);
+print(map(X, sumOfFirstAndLast, 0, false));
\ No newline at end of file
diff --git a/test/api/cli/secondorder/map_failure_7.daphne b/test/api/cli/secondorder/map_failure_7.daphne
new file mode 100644
index 000000000..72846b26c
--- /dev/null
+++ b/test/api/cli/secondorder/map_failure_7.daphne
@@ -0,0 +1,8 @@
+// udfReturnsScalar doesn't match actual UDF result
+
+def appendFirst(row) {
+  return cbind(row, [row[0, 0]]);
+}
+
+X = reshape([1, 2, 3], 1, 3);
+print(map(X, sumOfFirstAndLast, 0, true));
\ No newline at end of file
diff --git a/test/api/cli/secondorder/map_failure_8.daphne b/test/api/cli/secondorder/map_failure_8.daphne
new file mode 100644
index 000000000..02658954c
--- /dev/null
+++ b/test/api/cli/secondorder/map_failure_8.daphne
@@ -0,0 +1,13 @@
+// UDF returns row matrix instead of column matrix and vice versa
+
+def transposeFunc(X) {
+  return transpose(X);
+}
+
+// row-wise
+X = reshape([1, 2, 3, 1, 2, 3, 1, 2, 3], 3, 3);
+print(map(X, transposeFunc, 0));
+
+// column-wise
+X = reshape([1, 2, 3, 1, 2, 3, 1, 2, 3], 3, 3);
+print(map(X, transposeFunc, 1));
\ No newline at end of file
diff --git a/test/runtime/local/kernels/MapTest.cpp b/test/runtime/local/kernels/MapTest.cpp
index a59613c7d..7cc3982a7 100644
--- a/test/runtime/local/kernels/MapTest.cpp
+++ b/test/runtime/local/kernels/MapTest.cpp
@@ -27,20 +27,21 @@
 
 #define TYPES double, float, int64_t, int32_t, int8_t, uint64_t, uint8_t
 
-template <class DTRes, class DTArg> void checkMap(const DTArg *arg, const DTRes *exp, void *func) {
+template <class DTRes, class DTArg>
+void checkMap(const DTArg *arg, const DTRes *exp, void *func, const int64_t axis, const bool udfReturnsScalar) {
     DTRes *res = nullptr;
-    map(res, arg, func, nullptr);
+    map(res, arg, func, axis, udfReturnsScalar, nullptr);
     CHECK(*res == *exp);
     DataObjectFactory::destroy(res);
 }
 
-template <typename VTarg, typename VTres> VTres mult3func(VTarg arg) { return static_cast<VTres>(arg) * 3; }
+template <typename VTArg, typename VTRes> VTRes mult3func(VTArg arg) { return static_cast<VTRes>(arg) * 3; }
 
-template <template <typename VT> class DT, class VTarg, class VTres> void checkMult3Map() {
-    using DTArg = DT<VTarg>;
-    using DTRes = DT<VTres>;
+template <template <typename VT> class DT, class VTArg, class VTRes> void checkMult3Map() {
+    using DTArg = DT<VTArg>;
+    using DTRes = DT<VTRes>;
 
-    void *mult3funcPtr = reinterpret_cast<void *>(&mult3func<VTarg, VTres>);
+    void *mult3funcPtr = reinterpret_cast<void *>(&mult3func<VTArg, VTRes>);
 
     auto m1 = genGivenVals<DTArg>(3, {
                                          0,
@@ -92,21 +93,389 @@ template <template <typename VT> class DT, class VTarg, class VTres> void checkM
 
     auto mult3_res3 = genGivenVals<DTRes>(3, {3, 15, 9});
 
-    checkMap(m1, mult3_res1, mult3funcPtr);
-    checkMap(m2, mult3_res2, mult3funcPtr);
-    checkMap(m3, mult3_res3, mult3funcPtr);
+    checkMap(m1, mult3_res1, mult3funcPtr, -1, true);
+    checkMap(m2, mult3_res2, mult3funcPtr, -1, true);
+    checkMap(m3, mult3_res3, mult3funcPtr, -1, true);
 
     DataObjectFactory::destroy(m1, m2, m3, mult3_res1, mult3_res2, mult3_res3);
 }
 
-template <template <typename VT> class DT, typename VTarg, typename VTres1, typename... VTresN>
-std::enable_if_t<(sizeof...(VTresN) > 0)> checkMult3Map() {
-    checkMult3Map<DT, VTarg, VTres1>();
-    checkMult3Map<DT, VTarg, VTresN...>();
+template <template <typename VT> class DT, typename VTArg, typename VTRes1, typename... VTResN>
+std::enable_if_t<(sizeof...(VTResN) > 0)> checkMult3Map() {
+    checkMult3Map<DT, VTArg, VTRes1>();
+    checkMult3Map<DT, VTArg, VTResN...>();
 }
 
 TEMPLATE_TEST_CASE("Map element-wise dense matrix", TAG_KERNELS, TYPES) {
     // Test all combination of types in TYPES
     checkMult3Map<DenseMatrix, TestType, TYPES>();
     checkMult3Map<Matrix, TestType, TYPES>();
+}
+
+// DenseMatrix
+template <typename VTArg, typename VTRes> VTRes sumOfFirstAndLastFunc(const DenseMatrix<VTArg> *row) {
+    return (static_cast<VTRes>(row->getValues()[0]) + static_cast<VTRes>(row->getValues()[row->getNumCols() - 1]));
+}
+
+// Matrix
+template <typename VTArg, typename VTRes> VTRes sumOfFirstAndLastFunc(const Matrix<VTArg> *row) {
+    return (static_cast<VTRes>(row->get(0, 0)) + static_cast<VTRes>(row->get(0, row->getNumCols() - 1)));
+}
+
+template <template <typename VT> class DT, class VTArg, class VTRes> void checkSumOfFirstAndLastMap() {
+    using DTArg = DT<VTArg>;
+    using DTRes = DT<VTRes>;
+
+    void *sumRowFuncPtr = reinterpret_cast<void *>((VTRes(*)(const DT<VTArg> *)) & sumOfFirstAndLastFunc<VTArg, VTRes>);
+
+    auto m1 = genGivenVals<DTArg>(3, {
+                                         0,
+                                         1,
+                                         2,
+                                         1,
+                                         2,
+                                         3,
+                                         3,
+                                         4,
+                                         5,
+                                     });
+
+    auto sumRow_res1 = genGivenVals<DTRes>(3, {
+                                                  2,
+                                                  4,
+                                                  8,
+                                              });
+
+    auto m2 = genGivenVals<DTArg>(2, {
+                                         1,
+                                         0,
+                                         2,
+                                         0,
+                                         3,
+                                         1,
+                                         2,
+                                         0,
+                                     });
+
+    auto sumRow_res2 = genGivenVals<DTRes>(2, {1, 3});
+
+    auto m3 = genGivenVals<DTArg>(1, {1, 5, 3});
+
+    auto sumRow_res3 = genGivenVals<DTRes>(1, {4});
+
+    auto m4 = genGivenVals<DTArg>(3, {1, 5, 3});
+
+    auto sumRow_res4 = genGivenVals<DTRes>(3, {2, 10, 6});
+
+    checkMap(m1, sumRow_res1, sumRowFuncPtr, 0, true);
+    checkMap(m2, sumRow_res2, sumRowFuncPtr, 0, true);
+    checkMap(m3, sumRow_res3, sumRowFuncPtr, 0, true);
+    checkMap(m4, sumRow_res4, sumRowFuncPtr, 0, true);
+
+    DataObjectFactory::destroy(m1, m2, m3, m4, sumRow_res1, sumRow_res2, sumRow_res3, sumRow_res4);
+}
+
+template <template <typename VT> class DT, typename VTArg, typename VTRes1, typename... VTResN>
+std::enable_if_t<(sizeof...(VTResN) > 0)> checkSumOfFirstAndLastMap() {
+    checkSumOfFirstAndLastMap<DT, VTArg, VTRes1>();
+    checkSumOfFirstAndLastMap<DT, VTArg, VTResN...>();
+}
+
+TEMPLATE_TEST_CASE("Map row-wise dense matrix (Matrix -> Scalar)", TAG_KERNELS, TYPES) {
+    // Test all combination of types in TYPES
+    checkSumOfFirstAndLastMap<DenseMatrix, TestType, TYPES>();
+    checkSumOfFirstAndLastMap<Matrix, TestType, TYPES>();
+}
+
+// DenseMatrix
+template <typename VTArg, typename VTRes> DenseMatrix<VTRes> *appendFirst(const DenseMatrix<VTArg> *row) {
+    auto numCols = row->getNumCols();
+    auto res = DataObjectFactory::create<DenseMatrix<VTRes>>(1, numCols + 1, true);
+    const VTArg *valuesRow = row->getValues();
+    VTRes *valuesRes = res->getValues();
+    for (size_t c = 0; c < numCols; c++)
+        valuesRes[c] = valuesRow[c];
+    valuesRes[numCols] = valuesRow[0];
+    return res;
+}
+
+// Matrix
+template <typename VTArg, typename VTRes> Matrix<VTRes> *appendFirst(const Matrix<VTArg> *row) {
+    auto numCols = row->getNumCols();
+    auto res = DataObjectFactory::create<DenseMatrix<VTRes>>(1, numCols + 1, true);
+    VTRes *valuesRes = res->getValues();
+    for (size_t c = 0; c < numCols; c++)
+        valuesRes[c] = row->get(0, c);
+    valuesRes[numCols] = row->get(0, 0);
+    return dynamic_cast<Matrix<VTRes> *>(res);
+}
+
+template <template <typename VT> class DT, class VTArg, class VTRes> void checkAppendFirstMap() {
+    using DTArg = DT<VTArg>;
+    using DTRes = DT<VTRes>;
+
+    void *appendFirstFuncPtr =
+        reinterpret_cast<void *>((DT<VTRes> * (*)(const DT<VTArg> *)) & appendFirst<VTArg, VTRes>);
+
+    auto m1 = genGivenVals<DTArg>(3, {
+                                         0,
+                                         1,
+                                         2,
+                                         1,
+                                         2,
+                                         3,
+                                         3,
+                                         4,
+                                         5,
+                                     });
+
+    auto appendFirst_res1 = genGivenVals<DTRes>(3, {
+                                                       0,
+                                                       1,
+                                                       2,
+                                                       0,
+                                                       1,
+                                                       2,
+                                                       3,
+                                                       1,
+                                                       3,
+                                                       4,
+                                                       5,
+                                                       3,
+                                                   });
+
+    auto m2 = genGivenVals<DTArg>(2, {
+                                         1,
+                                         0,
+                                         2,
+                                         0,
+                                         3,
+                                         1,
+                                         2,
+                                         0,
+                                     });
+
+    auto appendFirst_res2 = genGivenVals<DTRes>(2, {
+                                                       1,
+                                                       0,
+                                                       2,
+                                                       0,
+                                                       1,
+                                                       3,
+                                                       1,
+                                                       2,
+                                                       0,
+                                                       3,
+                                                   });
+
+    auto m3 = genGivenVals<DTArg>(1, {1, 5, 3});
+
+    auto appendFirst_res3 = genGivenVals<DTRes>(1, {1, 5, 3, 1});
+
+    auto m4 = genGivenVals<DTArg>(3, {1, 5, 3});
+
+    auto appendFirst_res4 = genGivenVals<DTRes>(3, {1, 1, 5, 5, 3, 3});
+
+    checkMap(m1, appendFirst_res1, appendFirstFuncPtr, 0, false);
+    checkMap(m2, appendFirst_res2, appendFirstFuncPtr, 0, false);
+    checkMap(m3, appendFirst_res3, appendFirstFuncPtr, 0, false);
+    checkMap(m4, appendFirst_res4, appendFirstFuncPtr, 0, false);
+
+    DataObjectFactory::destroy(m1, m2, m3, m4, appendFirst_res1, appendFirst_res2, appendFirst_res3, appendFirst_res4);
+}
+
+template <template <typename VT> class DT, typename VTArg, typename VTRes1, typename... VTResN>
+std::enable_if_t<(sizeof...(VTResN) > 0)> checkAppendFirstMap() {
+    checkAppendFirstMap<DT, VTArg, VTRes1>();
+    checkAppendFirstMap<DT, VTArg, VTResN...>();
+}
+
+TEMPLATE_TEST_CASE("Map row-wise dense matrix (Matrix -> Matrix)", TAG_KERNELS, TYPES) {
+    // Test all combination of types in TYPES
+    checkAppendFirstMap<DenseMatrix, TestType, TYPES>();
+    checkAppendFirstMap<Matrix, TestType, TYPES>();
+}
+
+// DenseMatrix
+template <typename VTArg, typename VTRes> DenseMatrix<VTRes> *topThreeFunc(const DenseMatrix<VTArg> *col) {
+    auto res = DataObjectFactory::create<DenseMatrix<VTRes>>(3, 1, true);
+    const VTArg *valuesCol = col->getValues();
+    VTRes *valuesRes = res->getValues();
+    for (size_t r = 0; r < 3; r++) {
+        valuesRes[0] = static_cast<VTRes>(valuesCol[0]);
+        valuesCol += col->getRowSkip();
+        valuesRes += res->getRowSkip();
+    }
+    return res;
+}
+
+// Matrix
+template <typename VTArg, typename VTRes> Matrix<VTRes> *topThreeFunc(const Matrix<VTArg> *col) {
+    auto res = DataObjectFactory::create<DenseMatrix<VTRes>>(3, 1, true);
+    res->prepareAppend();
+    for (size_t r = 0; r < 3; r++)
+        res->append(r, 0, col->get(r, 0));
+    res->finishAppend();
+    return dynamic_cast<Matrix<VTRes> *>(res);
+}
+
+template <template <typename VT> class DT, class VTArg, class VTRes> void checkTopThreeMap() {
+    using DTArg = DT<VTArg>;
+    using DTRes = DT<VTRes>;
+
+    void *topThreeFuncPtr = reinterpret_cast<void *>((DT<VTRes> * (*)(const DT<VTArg> *)) & topThreeFunc<VTArg, VTRes>);
+
+    auto m1 = genGivenVals<DTArg>(4, {
+                                         0,
+                                         1,
+                                         2,
+                                         3,
+                                     });
+
+    auto topThree_res1 = genGivenVals<DTRes>(3, {
+                                                    0,
+                                                    1,
+                                                    2,
+                                                });
+
+    auto m2 = genGivenVals<DTArg>(4, {
+                                         0,
+                                         1,
+                                         2,
+                                         1,
+                                         2,
+                                         3,
+                                         3,
+                                         4,
+                                         5,
+                                         5,
+                                         6,
+                                         7,
+                                     });
+
+    auto topThree_res2 = genGivenVals<DTRes>(3, {
+                                                    0,
+                                                    1,
+                                                    2,
+                                                    1,
+                                                    2,
+                                                    3,
+                                                    3,
+                                                    4,
+                                                    5,
+                                                });
+
+    auto m3 = genGivenVals<DTArg>(10, {1, 5, 3, 2, 4, 7, 2, 0, 6, 4});
+
+    auto topThree_res3 = genGivenVals<DTRes>(3, {1, 5, 3});
+
+    checkMap(m1, topThree_res1, topThreeFuncPtr, 1, false);
+    checkMap(m2, topThree_res2, topThreeFuncPtr, 1, false);
+    checkMap(m3, topThree_res3, topThreeFuncPtr, 1, false);
+
+    DataObjectFactory::destroy(m1, m2, m3, topThree_res1, topThree_res2, topThree_res3);
+}
+
+template <template <typename VT> class DT, typename VTArg, typename VTRes1, typename... VTResN>
+std::enable_if_t<(sizeof...(VTResN) > 0)> checkTopThreeMap() {
+    checkTopThreeMap<DT, VTArg, VTRes1>();
+    checkTopThreeMap<DT, VTArg, VTResN...>();
+}
+
+TEMPLATE_TEST_CASE("Map column-wise dense matrix (Matrix -> Matrix)", TAG_KERNELS, TYPES) {
+    // Test all combination of types in TYPES
+    checkTopThreeMap<DenseMatrix, TestType, TYPES>();
+    // checkTopThreeMap<Matrix, TestType, TYPES>();
+}
+
+// DenseMatrix
+template <typename VTArg, typename VTRes> VTRes isSumEven(const DenseMatrix<VTArg> *col) {
+    int64_t sum = 0;
+    const VTArg *valuesCol = col->getValues();
+    for (size_t r = 0; r < col->getNumRows(); r++) {
+        sum += static_cast<int64_t>(valuesCol[0]);
+        valuesCol += col->getRowSkip();
+    }
+    return static_cast<VTRes>(sum % 2 == 0);
+}
+
+// Matrix
+template <typename VTArg, typename VTRes> VTRes isSumEven(const Matrix<VTArg> *col) {
+    int64_t sum = 0;
+    for (size_t r = 0; r < col->getNumRows(); r++)
+        sum += static_cast<int64_t>(col->get(r, 0));
+    return static_cast<VTRes>(sum % 2 == 0);
+}
+
+template <template <typename VT> class DT, class VTArg, class VTRes> void checkIsSumEvenMap() {
+    using DTArg = DT<VTArg>;
+    using DTRes = DT<VTRes>;
+
+    void *isSumEvenFuncPtr = reinterpret_cast<void *>((VTRes(*)(const DT<VTArg> *)) & isSumEven<VTArg, VTRes>);
+
+    auto m1 = genGivenVals<DTArg>(3, {
+                                         0,
+                                         2,
+                                         2,
+                                         3,
+                                         4,
+                                         5,
+                                         7,
+                                         8,
+                                         9,
+                                     });
+
+    auto isSumEven_res1 = genGivenVals<DTRes>(1, {
+                                                     1,
+                                                     1,
+                                                     1,
+                                                 });
+
+    auto m2 = genGivenVals<DTArg>(4, {
+                                         0,
+                                         1,
+                                         2,
+                                         1,
+                                         2,
+                                         3,
+                                         3,
+                                         4,
+                                         5,
+                                         5,
+                                         6,
+                                         7,
+                                     });
+
+    auto isSumEven_res2 = genGivenVals<DTRes>(1, {
+                                                     0,
+                                                     0,
+                                                     0,
+                                                 });
+
+    auto m3 = genGivenVals<DTArg>(3, {1, 5, 3});
+
+    auto isSumEven_res3 = genGivenVals<DTRes>(1, {0});
+
+    auto m4 = genGivenVals<DTArg>(1, {1, 6, 3});
+
+    auto isSumEven_res4 = genGivenVals<DTRes>(1, {0, 1, 0});
+
+    checkMap(m1, isSumEven_res1, isSumEvenFuncPtr, 1, true);
+    checkMap(m2, isSumEven_res2, isSumEvenFuncPtr, 1, true);
+    checkMap(m3, isSumEven_res3, isSumEvenFuncPtr, 1, true);
+    checkMap(m4, isSumEven_res4, isSumEvenFuncPtr, 1, true);
+
+    DataObjectFactory::destroy(m1, m2, m3, m4, isSumEven_res1, isSumEven_res2, isSumEven_res3, isSumEven_res4);
+}
+
+template <template <typename VT> class DT, typename VTArg, typename VTRes1, typename... VTResN>
+std::enable_if_t<(sizeof...(VTResN) > 0)> checkIsSumEvenMap() {
+    checkIsSumEvenMap<DT, VTArg, VTRes1>();
+    checkIsSumEvenMap<DT, VTArg, VTResN...>();
+}
+
+TEMPLATE_TEST_CASE("Map column-wise dense matrix (Matrix -> Scalar)", TAG_KERNELS, TYPES) {
+    // Test all combination of types in TYPES
+    checkIsSumEvenMap<DenseMatrix, TestType, TYPES>();
+    checkIsSumEvenMap<Matrix, TestType, TYPES>();
 }
\ No newline at end of file