;;; Lisplab, level1-util.lisp
;;; Level1, utility functions for matrix definitions.
;;; 

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

;;; TODO: change name of this to something about blas store
;;;
;;; 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)

(defun allocate-real-store (size &optional (initial-element 0d0))
  ;; All matrix double and complex double constructors 
  ;; should call this one
  (let ((x (coerce initial-element 'double-float)))
    (declare (type double-float x)
	     (type type-blas-idx size))
    ;; Stupid efficiency hack for SBCL. Allocations of arrays with zeros 
    ;; is significantly faster than others!
    (if (= x 0d0) 	
	(make-array size
		    :element-type 'double-float
		    :initial-element 0d0)	
	(make-array size
		    :element-type 'double-float
		    :initial-element x))))

;;; Hopfully helpful iterators

(defmacro with-indices-df-1 (a m idx &body body)
  "Iterats over all the indices of one array"
  `(let ((,m ,a))
     (declare (type type-blas-store ,m))
     (dotimes (,idx (length ,m))
       (declare (type type-blas-idx ,idx)) 
       ,@body)))

(defmacro with-elements-df-1 (a elm &body body)
  "Iterats over all the elements of one array"
  (let ((m (gensym))
	(idx (gensym)))
    `(let ((,m ,a))
       (declare (type type-blas-store ,m))
       (dotimes (,idx (length ,m))
	 (declare (type type-blas-idx ,idx))
	 (let ((,elm (aref ,m ,idx)))
	   (declare (type double-float ,elm))       
	   ,@body)))))

;;;; The complex store

(defun allocate-complex-store (size &optional (value 0d0))
  (let* ((2size (* 2 size))
	 (rv (coerce (realpart value) 'double-float))
	 (iv (coerce (imagpart value) 'double-float))
	 (store (allocate-real-store 2size iv)))
    (declare (type type-blas-idx 2size)
	     (type double-float rv iv))
    (when (/= rv iv)    
      (loop for i from 0 below 2size by 2 do
	   (setf (aref store i) rv)))
    store))

;;; The unsigend-byte 1 store

(defun allocate-ub1-store (size &optional (initial-element 0))
  (let ((x (coerce initial-element '(unsigned-byte 1))))
    (declare (type (unsigned-byte 1) x)
	     (type type-blas-idx size))
    ;; Stupid efficiency hack for SBCL. Allocations of arrays with zeros 
    ;; is significantly faster than others!
    (if (= x 0) 	
	(make-array size
		    :element-type '(unsigned-byte 1)
		    :initial-element 0)	
	(make-array size
		    :element-type '(unsigned-byte 1)
		    :initial-element x))))

;;; The unsigend-byte 8 store

(defun allocate-idx-store (size &optional (initial-element 0))
  (let ((x (coerce initial-element 'type-blas-idx)))
    (declare (type type-blas-idx x)
	     (type type-blas-idx size))
    ;; Stupid efficiency hack for SBCL. Allocations of arrays with zeros 
    ;; is significantly faster than others!
    (if (= x 0) 	
	(make-array size
		    :element-type 'type-blas-idx
		    :initial-element 0)	
	(make-array size
		    :element-type 'type-blas-idx
		    :initial-element x))))

;;; The unsigend-byte 8 store

(defun allocate-ub8-store (size &optional (initial-element 0))
  (let ((x (coerce initial-element '(unsigned-byte 8))))
    (declare (type (unsigned-byte 8) x)
	     (type type-blas-idx size))
    ;; Stupid efficiency hack for SBCL. Allocations of arrays with zeros 
    ;; is significantly faster than others!
    (if (= x 0) 	
	(make-array size
		    :element-type '(unsigned-byte 8)
		    :initial-element 0)	
	(make-array size
		    :element-type '(unsigned-byte 8)
		    :initial-element x))))

;;; The unsigend-byte 16 store

(defun allocate-ub16-store (size &optional (initial-element 0))
  (let ((x (coerce initial-element '(unsigned-byte 16))))
    (declare (type (unsigned-byte 16) x)
	     (type type-blas-idx size))
    ;; Stupid efficiency hack for SBCL. Allocations of arrays with zeros 
    ;; is significantly faster than others!
    (if (= x 0) 	
	(make-array size
		    :element-type '(unsigned-byte 16)
		    :initial-element 0)	
	(make-array size
		    :element-type '(unsigned-byte 16)
		    :initial-element x))))

;;; The unsigend-byte 32 store

(defun allocate-ub32-store (size &optional (initial-element 0))
  (let ((x (coerce initial-element '(unsigned-byte 32))))
    (declare (type (unsigned-byte 32) x)
	     (type type-blas-idx size))
    ;; Stupid efficiency hack for SBCL. Allocations of arrays with zeros 
    ;; is significantly faster than others!
    (if (= x 0) 	
	(make-array size
		    :element-type '(unsigned-byte 32)
		    :initial-element 0)	
	(make-array size
		    :element-type '(unsigned-byte 32)
		    :initial-element x))))