;;; 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.