(require-extension srfi-1 srfi-4 blas test)

(define (zeros m n)
  (let ((size (* m n)))
    (list->f64vector (list-tabulate size (lambda (i) 0.0)))))

(define (ones m n)
  (let ((size (* m n)))
    (list->f64vector (list-tabulate size (lambda (i) 1.0)))))

(define order blas:RowMajor)

(define asym  (f64vector 1 2 3 2 1 3 3 3 1))
(define vsym  (f64vector 9 10 11))

(define j (f64vector 440.0 1112.0 1784.0 2456.0))
(define h (f64vector 110.0 278.0 446.0 614.0)) 
  

(test-group "BLAS Level 1 Test"
  (test "daxpy" 
	(f64vector 440.0 1112.0 1784.0 2456.0)
	(blas:daxpy 4 4.0 h (zeros 4 1)))

  (test "ddot"
	2661184.0
	(blas:ddot 4 j h))
  )


(define order blas:RowMajor)

(define asym  (f64vector 1 2 3 2 1 3 3 3 1))
(define vsym  (f64vector 9 10 11))
   
(define j (f64vector 440.0 1112.0 1784.0 2456.0))
(define g (f64vector 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16))
(define h (f64vector 110.0 278.0 446.0 614.0))
(define v (f64vector 9 10 11 12)) 

(test-group "BLAS Level 2 Test"

  (test "dgemv" 
	(f64vector 110.0 278.0 446.0 614.0)
	(blas:dgemv order blas:NoTrans 4 4 1 g v 0 (zeros 4 1))
	)

  (test "dsymv" 
	(f64vector 62.0 61.0 68.0)
	(blas:dsymv order blas:Upper 3 1.0 asym vsym 0.0 (zeros 3 1))
	)

  (test "dsyr" 
	(f64vector 81 90 99 108 0 100 110 120 0 0 121 132 0 0 0 144)
	(blas:dsyr order blas:Upper 4 1.0 v (zeros 4 4)) 
	)

  (test "dger"
	(f64vector 48401.0 122321.0 196241.0 270161.0 122321.0
			 309137.0 495953.0 682769.0 196241.0 495953.0 795665.0 1095377.0 
			 270161.0 682769.0 1095377.0 1507985.0)
	(blas:dger order 4 4 1.0 j h (ones 4 4))
	) 
  )

(define order blas:RowMajor)

(define asym  (f64vector 1 2 3 2 1 3 3 3 1))

(define a (f64vector 1 2 3 4 5 6 7 8))
(define b (f64vector 1 2 3 4 5 6 7 8 9 10 11 12))
(define c (f64vector 9 10 11 12 13 14 15 16))
(define d (f64vector 70.0 80.0 90.0 158.0 184.0 210.0)) 

(test-group "BLAS Level 3 Test"

   (test "dgemm" 
	 (f64vector 70.0 80.0 90.0 158.0 184.0 210.0)
	 (blas:dgemm order blas:NoTrans blas:NoTrans 2 3 4 1.0 a b 0.0 (zeros 2 3))
	 )
   
   (test "dgemm transpose" 
	 (f64vector 860.0 1088.0 1316.0 1544.0 1000.0 1264.0 1528.0 1792.0 1140.0 1440.0 1740.0 2040.0)
	 (blas:dgemm order blas:Trans blas:NoTrans 3 4 2 1.0 d a 0.0 (zeros 3 4))
	 )

   (test "dsymm" 
	 (f64vector 38.0 44.0 50.0 56.0 34.0 40.0 46.0 52.0 27.0 34.0 41.0 48.0)
	 (blas:dsymm order blas:Left blas:Upper 3 4 1.0 asym b 0.0 (zeros 3 4))
	 )

   (test "dsyrk" 
	 (f64vector 30.0 70.0 0.0 174.0)
	 (blas:dsyrk order blas:Upper blas:NoTrans 2 4 1.0 a 0.0 (zeros 2 2))
	 )

   (test "dsyrk transpose"
	 (f64vector 26.0 32.0 38.0 44.0 0.0 40.0 48.0 56.0 0.0 0.0 58.0 68.0 0.0 0.0 0.0 80.0)
	 (blas:dsyrk order blas:Upper blas:Trans 4 2 1.0 a 0.0 (zeros 4 4))
	 )

   (test "dsyr2k" 
	 (f64vector 220.0 428.0 0.0 764.0)
	 (blas:dsyr2k order blas:Upper blas:NoTrans 2 4 1.0 a c 0.0 (zeros 2 2))
	 )
)