;;; Lisplab, level1-dgt.lisp
;;; Tridiagonal double-float matrices

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

;;; Note: not optimzied, but I see no good reason to optimzized them either.

;;; TODO: 
;;; I should reimplement this with just on store and just move 
;;; the sap when calling the FFIs. Then the matrices would work 
;;; automatically with the vector operations.

(in-package :lisplab)

;;; Tridiagonal matrices

(defclass matrix-base-dgt 
    (structure-tridiagonal
     element-double-float
     implementation-base) 
  ((size
    :initarg :size
    :reader size
    :type type-blas-idx)
   (diagonal-store
    :initarg :diagonal-store
    :initform nil
    :type type-blas-store)
   (subdiagonal-store
    :initarg :subdiagonal-store
    :initform nil
    :type type-blas-store)
   (superdiagonal-store
    :initarg :superdiagonal-store
    :initform nil
    :type type-blas-store)))

(defmethod initialize-instance :after ((m matrix-base-dgt) &key dim (value 0))
  (with-slots (rowcols size diagonal-store 
		       subdiagonal-store superdiagonal-store) 
      m
    (setf rowcols (if (numberp dim) 
		      dim
		      (if (= (car dim) (cadr dim))
			  (car dim)
			  (error "You try to create a non-square diagonal matrix.")))
	  size (- (* 3 rowcols) 2))
    (unless diagonal-store
      (setf diagonal-store (allocate-real-store rowcols value)))
    (unless subdiagonal-store
      (setf subdiagonal-store (allocate-real-store (1- rowcols) value)))
    (unless superdiagonal-store
      (setf superdiagonal-store (allocate-real-store (1- rowcols) value)))))

(defclass matrix-lisp-dgt (implementation-lisp matrix-base-dgt) ())

(defclass matrix-foreign-dgt (implementation-foreign matrix-lisp-dgt) ())

(defclass matrix-dgt (matrix-foreign-dgt) ())

(defmethod make-matrix-class ((a (eql :d)) (b (eql :gt)) (c (eql :base)))
  (find-class 'matrix-base-dgt))

(defmethod make-matrix-class ((a (eql :d)) (b (eql :gt)) (c (eql :lisp)))
  (find-class 'matrix-lisp-dgt))

(defmethod make-matrix-class ((a (eql :d)) (b (eql :gt)) (c (eql :ffi)))
  (find-class 'matrix-foreign-dgt))

(defmethod make-matrix-class ((a (eql :d)) (b (eql :gt)) (c (eql :foreign)))
  (find-class 'matrix-foreign-dgt))

(defmethod make-matrix-class ((a (eql :d)) (b (eql :gt)) (c (eql :blas)))
  (find-class 'matrix-foreign-dgt))

(defmethod make-matrix-class ((a (eql :d)) (b (eql :gt)) (c (eql :any)))
  (find-class 'matrix-dgt))

;;; Methods spezilied for the tridiagnoal matrices

(defmethod mref ((matrix matrix-base-dgt) row col)
  (cond ((= row col)
	 (aref (slot-value matrix 'diagonal-store) row))
	((= (1- row) col) 
	 (aref (slot-value matrix 'subdiagonal-store) col))
	((= (1+ row) col) 8
	 (aref (slot-value matrix 'superdiagonal-store) row))
	(t 0d0)))

(defmethod (setf mref) (value (matrix  matrix-base-dgt) row col)
  (let ((val2 (coerce value 'double-float)))
    (cond ((= row col)
	   (setf (aref (slot-value matrix 'diagonal-store) row)
		 val2))
	  ((= (1- row) col)
	   (setf (aref (slot-value matrix 'subdiagonal-store) col)
		 val2))
	  ((= (1+ row) col)
	   (setf (aref (slot-value matrix 'superdiagonal-store) row)
		 val2))
	  #+nil (t (warn "Array out of bonds for tridiagonal matrix. Ignored.")))))

(defmethod vref ((matrix  matrix-base-dgt) idx)
  (let ((len (slot-value matrix 'rowcols)))
    (cond ((< idx len)	   
	   (aref (slot-value matrix 'diagonal-store) idx))
	  ((< idx (- (* 2 len) 1))
	   (aref (slot-value matrix 'superdiagonal-store) (- idx len)))
	  ((< idx (slot-value matrix 'size))
	   (aref (slot-value matrix 'subdiagonal-store) (- idx (- (* 2 len) 1))))
	  #+nil (t 
	   (warn "Array out of bonds for tridiagonal matrix. Ignored.")))))
	   
(defmethod (setf vref) (value (matrix matrix-base-dgt) idx)
  (let ((val2 (coerce value 'double-float))
	(len (slot-value matrix 'rowcols)))
    (cond ((< idx len)	   
	   (setf (aref (slot-value matrix 'diagonal-store) idx)
		 val2))
	  ((< idx (- (* 2 len) 1))
	   (setf (aref (slot-value matrix 'superdiagonal-store) (- idx len))
		 val2))
	  ((< idx (- (* 3 len) 2))
	   (setf (aref (slot-value matrix 'subdiagonal-store) (- idx (- (* 2 len) 1)))
		 val2))
	  #+nil (t 
	   (warn "Array out of bonds for tridiagonal matrix. Ignored.")))
    val2))