;;; Something like this
(define-test destructive-matrix-op
  (let* ((matrix (generate-matrix ...)))
    (destructive-matrix-op matrix ...)
    (assert-numerical-equal
     exact-matrix matrix)))

;;; As a design note, most destructive operations still return the
;;; modified object. So, it would be good if DESTRUCTIVE-MATRIX-OP
;;; return the matrix.

MONTE-CARLO-TEST> (run-tests)
MAP-MATRIX!-TEST: 1 assertions passed, 0 failed.
MAP-MATRIX-TEST: 1 assertions passed, 0 failed.
MATRIX-TEST: 8 assertions passed, 0 failed.
TOTAL: 10 assertions passed, 0 failed, 0 execution errors.
; No value
MONTE-CARLO-TEST> (run-tests)
MAP-MATRIX!-TEST: MATRIX failed: 
Expected #2A((0.84147096 0.9092974) (0.14112 -0.7568025)) 
but saw #2A((0.7456241 0.78907233) (0.14065208 -0.68660027))
MAP-MATRIX!-TEST: 0 assertions passed, 1 failed.
MAP-MATRIX-TEST: 1 assertions passed, 0 failed.
MATRIX-TEST: 8 assertions passed, 0 failed.
TOTAL: 9 assertions passed, 1 failed, 0 execution errors.

(define-test map-matrix!-test
  (let ((matrix #2A((1 2) (3 4))))
    (map-matrix! #'sin matrix)
    (assert-numerical-equal
     #2A((0.841470985 0.909297427)
         (0.141120001 -0.756802495))
     matrix)))

(define-test map-matrix!-test
  (let ((matrix (make-array '(2 2) :initial-contents
                            '((1 2)
                              (3 4)))))
    (map-matrix! #'sin matrix)
    (assert-numerical-equal
     #2A((0.841470985 0.909297427)
         (0.141120001 -0.756802495))
     matrix)))

(defun map-matrix (function first-matrix &rest more-matrixes)
  "Apply FUNCTION on the elements of matrix."
  (let ((result (make-array (array-dimensions first-matrix))))
    (handler-case
        (loop for i below (array-total-size first-matrix)
           do (setf (row-major-aref result i)
                    (apply function
                           (row-major-aref first-matrix i)
                           (mapcar (lambda (m)
                                     (row-major-aref m i))
                                   more-matrixes)))
           finally (return result))
      (error nil
        (format t "Matrix size does not match or too few arguments")
        (return-from map-matrix nil)))))

(defun map-matrix! (function result-matrix &rest more-matrixes)
  "Apply FUNCTION on the elements of matrix MAT."
  (handler-case
      (loop for i below (array-total-size result-matrix)
         do (setf (row-major-aref result-matrix i)
                  (apply function
                         (row-major-aref result-matrix i)
                         (mapcar (lambda (m)
                                   (row-major-aref m i))
                                 more-matrixes)))
         finally (return result-matrix))
    (error nil 
      (format t "Matrix size does not match or too few arguments")
      (return-from map-matrix! nil))))

;;; Lightly tested; needs to be factored; generally models MAP and MAP-INTO
(defun %validate-dimensions (matrix matrices)
  "Return the dimensions validating that all are equal."
  (if (= 2 (array-rank matrix))
      (let ((dimensions (array-dimensions matrix)))
        (dolist (mat matrices dimensions)
          (unless (equal dimensions (array-dimensions mat))
            (error "Unequal matrix dimensions."))))
      (error "Matrix rank, ~S, does not equal 2."
             (array-rank matrix))))

;;; FIXME : The result element type should be set to the most general
;;; type of all of the matrices according to UPGRADED-ARRAY-ELEMENT-TYPE
(defun map-matrix (function matrix &rest matrices)
  "Apply FUNCTION on the elements of matrix."
  (declare (optimize debug))
  (if matrices
      ;; The function accepts more than one matrix.
      (destructuring-bind (rows columns)
          (%validate-dimensions matrix matrices)
        (let ((arglist (make-list (length matrices)))
              (result  (make-array (list rows columns)
                                   :element-type
                                   (array-element-type matrix))))
          (dotimes (i0 rows result)
            (dotimes (i1 columns)
              (setf
               (aref result i0 i1)
               (apply function
                      (aref matrix i0 i1)
                      (map-into arglist
                                (lambda (mat) (aref mat i0 i1))
                                matrices)))))))
      ;; The function accepts one matrix.
      (destructuring-bind (rows columns)
          (array-dimensions matrix)
        (let ((result (make-array (list rows columns)
                                  :element-type
                                  (array-element-type matrix))))
          (dotimes (i0 rows result)
            (dotimes (i1 columns)
              (setf
               (aref result i0 i1)
               (funcall function (aref matrix i0 i1)))))))))

(defun map-into-matrix (result-matrix function &rest matrices)
  "Apply FUNCTION on the elements of matrix MAT."
  (if matrices
      ;; The result is built from other matrices.
      (destructuring-bind (rows columns)
          (%validate-dimensions result-matrix matrices)
        (let ((arglist (make-list (length matrices))))
          (dotimes (i0 rows result-matrix)
            (dotimes (i1 columns)
              (setf
               (aref result-matrix i0 i1)
               (apply function
                      (map-into arglist
                                (lambda (mat) (aref mat i0 i1))
                                matrices)))))))
      ;; The result is built from a function.
      (destructuring-bind (rows columns)
          (array-dimensions result-matrix)
        (dotimes (i0 rows result-matrix)
          (dotimes (i1 columns)
            (setf
             (aref result-matrix i0 i1)
             (funcall function)))))))

;;; FIXME : There are some opportunities for abstraction,
;;;         but you get what you pay for.