;;; Lisplab, level2-matrix-zge.lisp
;;; Optimizations for complex matrices.

;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU General Public License as published by
;;; the Free Software Foundation; either version 2 of the License, or
;;; (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License along
;;; with this program; if not, write to the Free Software Foundation, Inc.,
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

(in-package :lisplab)

(defmethod mfill ((a vector-z) value)
  (let ((rx (coerce (realpart value) 'double-float))
	(cx (coerce (imagpart value) 'double-float))
	(store (vector-store a)))
    (declare (type type-blas-store store))
    (loop for i from 0 below (length store) by 2 do 
      (setf (aref store i) rx
	    (aref store (1+ i)) cx))))

(defmethod copy  ((a vector-z))
  (make-instance (class-name (class-of a))
		 :store (copy-seq (the type-blas-store (vector-store a)))
		 :dim (dim a)))

(defmethod copy-contents ((a vector-z) (b vector-z) 
			  &optional (converter nil))
  (let ((store-a (vector-store a))
	(store-b (vector-store b)))
    (if converter 
	(map-into store-b converter store-a)
	(copy-matrix-stores store-a store-b)))    
    b)

(defmethod copy-contents ((from vector-d) (to vector-z) 
			  &optional (converter nil))
  (if converter 
      (call-next-method) ;; Could have some testes here to improve performance
      (let* ((store-a (vector-store from))
	     (store-b (vector-store to))
	     (len (length store-a)))
	(declare (type type-blas-store store-a store-b)
		 (type type-blas-idx len))
	(dotimes (i len)
	  (setf (aref store-b (truly-the type-blas-idx (* 2 i))) (aref store-a i)))
	to))) 

(defmethod msum ((m vector-z))
  (let ((sum-r 0d0)
        (sum-i 0d0)
        (m0 (vector-store m)))
    (declare (type double-float sum-r sum-i)
             (type type-blas-store m0))
    (loop for i from 0 below (length m0) by 2 do
          (incf sum-r (aref m0 i))
          (incf sum-i (aref m0 (1+ i))))
    (complex sum-r sum-i)))


;;; Cration of complex vector

(defun dge-to-zge (a)
  (let ((b (mcreate* a :element-type :z)))
    (copy_dfa-to-cdfa (vector-store a) (vector-store b))
    b))

(defmethod .complex ((a vector-d) (b real))
  (let ((c (mcreate* a :element-type :z)))
    (complex_dfa-df (vector-store a) (coerce b 'double-float) (vector-store c))
    c))

(defmethod .complex ((a real) (b vector-d))
  (let ((c (mcreate* b :element-type :z)))
    (complex_dfa-df (coerce a 'double-float) (vector-store b) (vector-store c))
    c))

(defmethod .complex ((a vector-d) (b vector-d))
  (let ((c (mcreate* a :element-type :z)))
    (complex_dfa-dfa (vector-store a) (vector-store b) (vector-store c))
    c))

 ;;; Vector and complex

(define-constant +generic-function-cdfa-cdf-map+ 
    '((.add . +_cdfa-cdf)
      (.sub . -_cdfa-cdf)
      (.mul . *_cdfa-cdf)
      (.div . /_cdfa-cdf)
      (.expt . ^_cdfa-cdf)))

(defmacro defmethod-cdfa-cdf (name underlying-function)  
  (let ((a (gensym "a"))
	(b (gensym "b"))
	(c (gensym "c")))
  `(defmethod ,name ((,a vector-z) (,b number))
     (let ((,c (mcreate ,a)))
       (,underlying-function (vector-store ,a) 
			     (coerce ,b '(complex double-float)) 
			     (vector-store ,c))
       ,c))))

(defmacro expand-generic-function-cdfa-cdf-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defmethod-cdfa-cdf ,(car x) ,(cdr x)))
	      +generic-function-cdfa-cdf-map+)))

(expand-generic-function-cdfa-cdf-map)

;;; Complex and vector 

(define-constant +generic-function-cdf-cdfa-map+ 
    '((.add . +_cdf-cdfa)
      (.sub . -_cdf-cdfa)
      (.mul . *_cdf-cdfa)
      (.div . /_cdf-cdfa)
      (.expt . ^_cdf-cdfa)))

(defmacro defmethod-cdf-cdfa (name underlying-function)  
  (let ((a (gensym "a"))
	(b (gensym "b"))
	(c (gensym "c")))
  `(defmethod ,name ((,a number) (,b vector-z))
     (let ((,c (mcreate ,b)))
       (,underlying-function (coerce ,a '(complex double-float)) 
			     (vector-store ,b) 
			     (vector-store ,c))
       ,c))))

(defmacro expand-generic-function-cdf-cdfa-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defmethod-cdf-cdfa ,(car x) ,(cdr x)))
	      +generic-function-cdf-cdfa-map+)))

(expand-generic-function-cdf-cdfa-map)

;;;; Vector and vector

(define-constant +generic-function-cdfa-cdfa-map+ 
    '((.add . +_cdfa-cdfa)
      (.sub . -_cdfa-cdfa)
      (.mul . *_cdfa-cdfa)
      (.div . /_cdfa-cdfa)
      (.expt . ^_cdfa-cdfa) ))

(defmacro defmethod-cdfa-cdfa (name underlying-function)  
  (let ((a (gensym "a"))
	(b (gensym "b"))
	(c (gensym "c")))
  `(defmethod ,name ((,a vector-z) (,b vector-z))
     (let ((,c (mcreate ,a)))
       (,underlying-function (vector-store ,a) (vector-store ,b) (vector-store ,c))
       ,c))))

(defmacro expand-generic-function-cdfa-cdfa-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defmethod-cdfa-cdfa ,(car x) ,(cdr x)))
	      +generic-function-cdfa-cdfa-map+)))

(expand-generic-function-cdfa-cdfa-map)

;;;; The cross terms, where one input is a real vector, the other complex

;;; Add

(defmethod .add ((a vector-d) (b complex))
  (.add (dge-to-zge a) (coerce b '(complex double-float))))

(defmethod .add ((a complex) (b vector-d))
  (.add (coerce a '(complex double-float)) (dge-to-zge b)))

(defmethod .add ((a vector-z) (b real))
  (let ((c (mcreate a)))
    (+_cdfa-df (vector-store a) (coerce b 'double-float) (vector-store c))
    c))

(defmethod .add ((b real) (a vector-z) )
  (.add a b))

(defmethod .add ((a vector-z) (b vector-d))
  (let ((c (mcreate a)))
    (+_cdfa-dfa (vector-store a) (vector-store b) (vector-store c))
    c))

(defmethod .add ((b vector-d) (a vector-z))
  (.add a b))

;;; Sub

(defmethod .sub ((a vector-d) (b complex))
  (.sub (dge-to-zge a) (coerce b '(complex double-float))))

(defmethod .add ((a complex) (b vector-d))
  (.sub (coerce a '(complex double-float)) (dge-to-zge b)))

(defmethod .sub ((a vector-z) (b real))
  (let ((c (mcreate a)))
    (-_cdfa-df (vector-store a) (coerce b 'double-float) (vector-store c))
    c))

(defmethod .sub ((a real) (b vector-z))
  (let ((c (mcreate b)))
    (-_df-cdfa (coerce a 'double-float) (vector-store b) (vector-store c))
    c))

(defmethod .sub ((a vector-z) (b vector-d))
  (let ((c (mcreate a)))
    (-_cdfa-dfa (vector-store a) (vector-store b) (vector-store c))
    c))

(defmethod .sub ((a vector-d) (b vector-z))
  (let ((c (mcreate b)))
    (-_dfa-cdfa (vector-store a) (vector-store b) (vector-store c))
    c))

;;; Mul

(defmethod .mul ((a vector-d) (b complex))
  (.mul (dge-to-zge a) (coerce b '(complex double-float))))

(defmethod .mul ((a complex) (b vector-d))
  (.mul (coerce a '(complex double-float)) (dge-to-zge b)))

(defmethod .mul ((a vector-z) (b real))
  (let ((c (mcreate a)))
    (*_cdfa-df (vector-store a) (coerce b 'double-float) (vector-store c))
    c))

(defmethod .mul ((b real) (a vector-z) )
  (.mul a b))

(defmethod .mul ((a vector-z) (b vector-d))
  (let ((c (mcreate a)))
    (*_cdfa-dfa (vector-store a) (vector-store b) (vector-store c))
    c))

(defmethod .mul ((b vector-d) (a vector-z))
  (.mul a b))

;;; Div

(defmethod .div ((a vector-d) (b complex))
  (.div (dge-to-zge a) (coerce b '(complex double-float))))

(defmethod .div ((a complex) (b vector-d))
  (.div (coerce a '(complex double-float)) (dge-to-zge b)))

(defmethod .div ((a vector-z) (b real))
  (let ((c (mcreate a)))
    (/_cdfa-df (vector-store a) (coerce b 'double-float) (vector-store c))
    c))

(defmethod .div ((a vector-z) (b vector-d))
  (let ((c (mcreate a)))
    (/_cdfa-dfa (vector-store a) (vector-store b) (vector-store c))
    c))

(defmethod .div ((a vector-d) (b vector-z))
  (.div (dge-to-zge a) b))

#+nil (def-cross-complex-real-method .div vector-d vector-z)

;;; Expt

(defmethod .expt ((a vector-d) (b complex))
  (.expt (dge-to-zge a) b))

(defmethod .expt ((a complex) (b vector-d))
  (.expt a (dge-to-zge b)))

(defmethod .expt ((a vector-d) (b vector-z))
  (.expt (dge-to-zge a) b))

(defmethod .expt ((a vector-z) (b vector-d))
  (.expt a (dge-to-zge b)))

#+nil (def-cross-complex-real-method .expt vector-z vector-d)
#+nil (def-cross-complex-real-method .expt vector-d vector-z)


;;;; Ordinary functions

(define-constant +generic-function-cdfa-to-cdfa-map+ ;really bad name    
    '((.sin . sin_cdfa-to-cdfa) 
      (.cos . cos_cdfa-to-cdfa) 
      (.tan . tan_cdfa-to-cdfa) 
      (.asin . asin_cdfa-to-cdfa)   
      (.acos . acos_cdfa-to-cdfa)  
      (.atan . atan_cdfa-to-cdfa) 
      (.sinh . sinh_cdfa-to-cdfa) 
      (.cosh . cosh_cdfa-to-cdfa) 
      (.tanh . tanh_cdfa-to-cdfa) 
      (.asinh . asinh_cdfa-to-cdfa) 
      (.acosh . acosh_cdfa-to-cdfa) 
      (.atanh . atanh_cdfa-to-cdfa)
      (.exp . exp_cdfa-to-cdfa) 
      (.ln . log_cdfa-to-cdfa) 
      (.sqrt . sqrt_cdfa-to-cdfa) 
      (.conj . conjugate_cdfa-to-cdfa)))

(defmacro defmethod-cdfa-to-cdfa (name underlying-function)  
  (let ((a (gensym "a"))
	(b (gensym "b")))
  `(defmethod ,name ((,a vector-z))
     (let ((,b (mcreate ,a)))
       (,underlying-function (vector-store ,a) (vector-store ,b) )
       ,b))))

(defmacro expand-generic-function-cdfa-to-cdfa-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defmethod-cdfa-to-cdfa ,(car x) ,(cdr x)))
	      +generic-function-cdfa-to-cdfa-map+)))

(expand-generic-function-cdfa-to-cdfa-map)

;;;; Exceptions, in that output is real for complex input

(defmethod .im ((a vector-z))
  (let ((out (mcreate* a :element-type :d)))
    (imagpart_cdfa (vector-store a) (vector-store out))
    out))

(defmethod .re ((a vector-z))
  (let ((out (mcreate* a :element-type :d)))
    (realpart_cdfa (vector-store a) (vector-store out))
    out))

(defmethod .abs ((a vector-z))
  (let ((out (mcreate* a :element-type :d)))
    (abs_cdfa (vector-store a) (vector-store out))
    out))