;;; Lisplab, ext-store-operators.lisp
;;; Double float and complex double float operators (such as +,-,*, etc) on
;;; simple arrays. 
;;; Extended parameters list with offset and step.
;;; Used by the matrix implementations. 

;;; Copyright (C) 2009 Joern Inge Vestgaarden
;;;
;;; 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.

;;; This file contains manipulations of simple double-float arrays 
;;; and should be called by the spesialized matrix methods. 
;;; The purpose of this layer is that it can be used by 
;;; many classes such as matrix-base-dge and matrix-base-ddi, etc. 
;;; 
;;; The content of this file must be highly optimized 
;;; and should not depend anything exept Common Lisp itself.

(in-package :lisplab)

;;; TODO: there must be some easier way to generate the code in this file,
;;;       but I have not the energy to do it. I do, however, think that
;;;       the basic idea of having a layer of ordinary functions is a good one. 

;;; The reason for generating ordinary functions and not using methods,
;;; is that the real and complex stores have the same type. The fortran-compatible
;;; complex arrays are just subsequent real and complex double-floats.

;;; The reason for having both real and complex in the same file is that 
;;; not all operators function on both real and complex arguments. Care must 
;;; be taken. This is also the reason why it's hard to generate more code 
;;; automatically. 

;;; The below code generates ordinary lisp functions 
;;; for elementwise operations on simple double-float arrays. 
;;; They use a naming conventions, which should be pretty easy to 
;;; guess, such as df = double float and cdfa = complex double float array. 
;;;
;;; (The convention for complex should for consistnt naming be changed to zdfa, 
;;; but its not really important)
;;; 
;;; I use map-into when its performance is equal or better to the iterations. 
;;; The iterative version for all operations are still in the file, since other lisps  
;;; than sbcl might have a slower map-into, so that they can be needed later.
;;; For real numbers, map-into can be used for all operations, while for complex 
;;; number only + and - (*, / and expt mix the real and complex parts)



(declaim (inline double-float-simple-array-ref-ext))
(declaim (ftype (function 
		 (type-blas-store
		  type-blas-idx
		  type-blas-idx
		  type-blas-idx)
		 double-float)
		double-float-simple-array-ref-ext))
(defun double-float-simple-array-ref-ext (a i off step)
  (declare (type type-blas-idx i off step)
	   (type type-blas-store a))
  (aref a (column-major-idx off i step)))

(declaim (inline (setf double-float-simple-array-ref-ext)))
(declaim (ftype (function 
		 (double-float
		  type-blas-store
		  type-blas-idx
		  type-blas-idx
		  type-blas-idx)
		 double-float)
		(setf double-float-simple-array-ref-ext)))
(defun (setf double-float-simple-array-ref-ext) (value a i off step)
  (declare (type type-blas-idx i off step)
	   (type double-float value)
	   (type type-blas-store a))
  (setf (aref a (column-major-idx off i step))
	value)
  value)

;;;  Array and number
	   
(defmacro defun-ext-dfa-df (name op)
  `(defun ,name (len a aoff astep b c coff cstep) 
     (declare (type type-blas-store a c)
	      (type double-float b))
     (dotimes (i len)
       (setf (double-float-simple-array-ref-ext c i coff cstep)
	     (,op (double-float-simple-array-ref-ext a i aoff astep)
		  b)))
     c))

(defun-ext-dfa-df +_ext-dfa-df +)
(defun-ext-dfa-df -_ext-dfa-df -)
(defun-ext-dfa-df *_ext-dfa-df *)
(defun-ext-dfa-df /_ext-dfa-df /)
(defun-ext-dfa-df ^_ext-dfa-df expt)
(defun-ext-dfa-df max_ext-dfa-df max)
(defun-ext-dfa-df min_ext-dfa-df min)
(defun-ext-dfa-df log_ext-dfa-df log)

;;; Number and array

(defmacro defun-ext-df-dfa (name op)
  `(defun ,name (len a b boff bstep c coff cstep) 
     (declare (type type-blas-store b c)
	      (type double-float a))
     (dotimes (i len)
       (setf (double-float-simple-array-ref-ext c i coff cstep)
	     (,op a 
		  (double-float-simple-array-ref-ext b i boff bstep))))
     c))

(defun-ext-df-dfa +_ext-df-dfa +)
(defun-ext-df-dfa -_ext-df-dfa -)
(defun-ext-df-dfa *_ext-df-dfa *)
(defun-ext-df-dfa /_ext-df-dfa /)
(defun-ext-df-dfa ^_ext-df-dfa expt)
(defun-ext-df-dfa max_ext-df-dfa max)
(defun-ext-df-dfa min_ext-df-dfa min)
(defun-ext-df-dfa log_ext-df-dfa log)

;;; Array and array

(defmacro defun-ext-dfa-dfa (name op)
  `(defun ,name (len a aoff astep b boff bstep c coff cstep) 
     (declare (type type-blas-store b a c))
     (dotimes (i len)
       (setf (double-float-simple-array-ref-ext c i coff cstep)
	     (,op (double-float-simple-array-ref-ext a i aoff astep)
		  (double-float-simple-array-ref-ext b i boff bstep))))
     c))

(defun-ext-dfa-dfa +_ext-dfa-dfa +)
(defun-ext-dfa-dfa -_ext-dfa-dfa -)
(defun-ext-dfa-dfa *_ext-dfa-dfa *)
(defun-ext-dfa-dfa /_ext-dfa-dfa /)
(defun-ext-dfa-dfa ^_ext-dfa-dfa expt)
(defun-ext-dfa-dfa max_ext-dfa-dfa max)
(defun-ext-dfa-dfa min_ext-dfa-dfa min)
(defun-ext-dfa-dfa log_ext-dfa-dfa log)