(require-extension srfi-4 blas )

(define matrix-size f64vector-length)


(define read-matrix
  (lambda (port)
    (let ((m (read port))
	  (n (read port)))
        (let ((size (* m n))
	      (mat (make-f64vector (* m n))))
          (do ((i 0 (+ 1 i)))
              ((= i size) mat)
	    (f64vector-set! mat i (read port)))))))

;; match-error is called to complain when mul receives a pair of
;; incompatible arguments.
(define match-error
  (lambda (what1 what2)
    (error 'mul
           "~s and ~s are incompatible operands"
           what1
           what2)))

(define (zeros m n)
  (let ((size (* m n)))
    (make-f64vector size 0.0)))

(define bench-multiply
  (lambda (a b times)
    (let ((c #f))
      (do ((i 0 (+ i 1)))
          ((> i times))
          (set! c (blas:dgemm blas:RowMajor blas:NoTrans blas:NoTrans 100 100 4 1.0 a b 0.0 (zeros 100 100))))
      c)))

(define main
  (lambda (times)
    (let ((port-a (open-input-file "x.mat"))
          (port-b (open-input-file "y.mat")))
      (let ((a (read-matrix port-a))
            (b (read-matrix port-b))
            (c #f))
        (display (format "------------"))
        (newline)
        (time (set! c (bench-multiply a b times)))
        (display (format "------------"))
        (newline)
        (if (< (matrix-size c) 100)
            (display c))
        (newline)
        0))))

(newline)
(main (read))