From: Akshay S. <ak...@us...> - 2013-08-14 19:53:53
|
This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "matlisp". The branch, classy has been updated via f9bf6a61b1860942b520069596563b2db546f927 (commit) via be4148122456b5f7a6d4032fcba44e4652f4eb0b (commit) from 228188fe426f884dd6a1743578e879350b7050ec (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit f9bf6a61b1860942b520069596563b2db546f927 Author: Akshay Srinivasan <aks...@gm...> Date: Wed Aug 14 12:53:03 2013 -0700 Saving state. diff --git a/matlisp.asd b/matlisp.asd index 71aa515..8450621 100644 --- a/matlisp.asd +++ b/matlisp.asd @@ -147,11 +147,8 @@ :depends-on ("copy" "maker")) (:file "realimag" :depends-on ("copy")) - #+nil - ( (:file "trans" - :depends-on ("scal" "copy"))))) - #+nil + :depends-on ("scal" "copy")))) (:module "matlisp-level-2" :pathname "level-2" :depends-on ("matlisp-base" "matlisp-classes" "foreign-core" "matlisp-level-1") diff --git a/src/base/blas-helpers.lisp b/src/base/blas-helpers.lisp index 11ef0c0..de623a8 100644 --- a/src/base/blas-helpers.lisp +++ b/src/base/blas-helpers.lisp @@ -46,8 +46,8 @@ (cs (aref stds 1) :type index-type)) ;;Note that it is not required that (rs = nc * cs) or (cs = nr * rs) (cond - ((= cs 1) (values :row-major rs (fortran-nop op))) - ((= rs 1) (values :col-major cs op))))) + ((= cs 1) (values rs (fortran-nop op) :row-major)) + ((= rs 1) (values cs op :col-major))))) ;;Stride makers. (definline make-stride-rmj (dims) diff --git a/src/base/template.lisp b/src/base/template.lisp index bdf0d06..f027c83 100644 --- a/src/base/template.lisp +++ b/src/base/template.lisp @@ -145,8 +145,6 @@ :for eleb :in b :do (when (funcall func elea eleb) (return t)))) -) - ;;This one is hard to get one's brain around. (deft/generic (t/strict-coerce @@ -195,3 +193,5 @@ ;; (t/strict-coerce (fixnum real) x) -> (COERCE X 'REAL) ;; (t/strict-coerce (double-float t) x) -> X ;; (t/strict-coerce (fixnum (complex integer)) x) -> (COERCE X '(COMPLEX INTEGER)) + +) diff --git a/src/classes/numeric.lisp b/src/classes/numeric.lisp index aaeb41e..8605abb 100644 --- a/src/classes/numeric.lisp +++ b/src/classes/numeric.lisp @@ -99,3 +99,4 @@ (defleaf scomplex-tensor (complex-numeric-tensor) ()) (deft/method t/store-element-type (sym scomplex-tensor) () 'single-float)) + diff --git a/src/level-1/axpy.lisp b/src/level-1/axpy.lisp index 52aae00..d6dd42a 100644 --- a/src/level-1/axpy.lisp +++ b/src/level-1/axpy.lisp @@ -122,8 +122,6 @@ `(t/axpy! ,clx alpha x y)) y))) (axpy! alpha x y)) - ((coerceable? clx cly) - (axpy! alpha (coerce-tensor x cly) y)) (t (error "Don't know how to apply axpy! to classes ~a, ~a." clx cly))))) diff --git a/src/level-1/dot.lisp b/src/level-1/dot.lisp index 1866473..354c271 100644 --- a/src/level-1/dot.lisp +++ b/src/level-1/dot.lisp @@ -131,11 +131,5 @@ (t/dot ,clx x y t) (t/dot ,clx x y nil))))) (dot x y conjugate-p)) - ;;You pay the piper if you like mixing types. - ;;This is (or should be) a rare enough to not matter. - ((coerceable? clx cly) - (dot (coerce-tensor x cly) y conjugate-p)) - ((coerceable? cly clx) - (dot x (coerce-tensor y clx) conjugate-p)) (t (error "Don't know how to compute the dot product of ~a , ~a." clx cly))))) diff --git a/src/level-1/scal.lisp b/src/level-1/scal.lisp index 6a06b36..cedea4d 100644 --- a/src/level-1/scal.lisp +++ b/src/level-1/scal.lisp @@ -114,8 +114,6 @@ `(t/scdi! ,clx x y :scal? t :numx? nil)) y)) (scal! x y)) - ((coerceable? clx cly) - (scal! (coerce-tensor x cly) y)) (t (error "Don't know how to apply scal! to classes ~a, ~a." clx cly))))) @@ -170,8 +168,6 @@ `(t/scdi! ,clx x y :scal? nil :numx? nil)) y)) (div! x y)) - ((coerceable? clx cly) - (div! (coerce-tensor x cly) y)) (t (error "Don't know how to apply div! to classes ~a, ~a." clx cly))))) diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp index 28a82d0..177d7c8 100644 --- a/src/level-2/gemv.lisp +++ b/src/level-2/gemv.lisp @@ -17,8 +17,6 @@ (declare (type ,sym ,A ,x ,y) (type ,(field-type sym) ,alpha ,beta) (type index-type ,st-x ,st-y ,lda ,m ,n)) - (when (cl:char= (char-upcase ,transp) #\T) - (rotatef ,m ,n)) (,(macroexpand-1 `(t/blas-gemv-func ,sym)) ,transp ,m ,n ,alpha @@ -27,7 +25,7 @@ ,beta (the ,(store-type sym) (store ,y)) ,st-y (the index-type (head ,A)) (the index-type (head ,x)) (the index-type (head ,y))) - y)))) + ,y)))) ;; (deft/generic (t/gemv! #'subtypep) sym (alpha A x beta y transp)) @@ -37,20 +35,19 @@ (using-gensyms (decl (alpha A x beta y transp)) `(let (,@decl) (declare (type ,sym ,A ,x ,y) - (type ,(field-type sym) ,alpha ,beta)) + (type ,(field-type sym) ,alpha ,beta) + (type character ,transp)) (unless (t/f= ,(field-type sym) ,beta (t/fid* ,(field-type sym))) (t/scdi! ,sym ,beta ,y :scal? t :numx? t)) ;;These loops are optimized for column major matrices - (if ,transp - (einstein-sum ,sym (i j) (ref ,y i) (* ,alpha (ref ,A j i) (ref ,x j)) nil) - (einstein-sum ,sym (j i) (ref ,y i) (* ,alpha (ref ,A i j) (ref ,x j)) nil)) + (ecase (char-upcase ,transp) + (#\N (einstein-sum ,sym (j i) (ref ,y i) (* ,alpha (ref ,A i j) (ref ,x j)) nil)) + (#\T (einstein-sum ,sym (i j) (ref ,y i) (* ,alpha (ref ,A j i) (ref ,x j)) nil)) + (#\C (einstein-sum ,sym (j i) (ref ,y i) (* ,alpha (t/fc ,(field-type sym) (ref ,A i j)) (ref ,x j)) nil)) + (#\H (einstein-sum ,sym (i j) (ref ,y i) (* ,alpha (t/fc ,(field-type sym) (ref ,A j i)) (ref ,x j)) nil))) ,y))) - -;;Symbolic -#+maxima -(generate-typed-gemv! symbolic-base-typed-gemv! - (symbolic-tensor nil 0)) ;;---------------------------------------------------------------;; + (defgeneric gemv! (alpha A x beta y &optional job) (:documentation " @@ -79,74 +76,75 @@ :C alpha * conjugate(A) * x + beta * y :H alpha * transpose o conjugate(A) + beta * y ") - (:method :before ((alpha number) (A standard-matrix) (x standard-vector) - (beta number) (y standard-vector) + (:method :before (alpha (A standard-tensor) (x standard-tensor) + beta (y standard-tensor) &optional (job :n)) (assert (member job '(:n :t :c :h)) nil 'invalid-value :given job :expected `(member job '(:n :t :c :h)) - :message "Inside gemv!") + :message "GEMV!: Given an unknown job.") (assert (not (eq x y)) nil 'invalid-arguments :message "GEMV!: x and y cannot be the same vector") (assert (and - (= (aref (dimensions x) 0) - (aref (dimensions A) (if (eq job :t) 0 1))) - (= (aref (dimensions y) 0) - (aref (dimensions A) (if (eq job :t) 1 0)))) + (tensor-vectorp x) (tensor-vectorp y) (tensor-matrixp A) + (= (aref (the index-store-vector (dimensions x)) 0) + (aref (the index-store-vector (dimensions A)) (if (member job '(:t :h)) 0 1))) + (= (aref (the index-store-vector (dimensions y)) 0) + (aref (the index-store-vector (dimensions A)) (if (member job '(:t :h)) 1 0)))) nil 'tensor-dimension-mismatch))) -(defmethod gemv! ((alpha number) (A real-matrix) (x real-vector) - (beta number) (y real-vector) &optional (job :n)) - (real-typed-gemv! (coerce-real alpha) A x - (coerce-real beta) y job)) - -(defmethod gemv! ((alpha number) (A complex-matrix) (x complex-vector) - (beta number) (y complex-vector) &optional (job :n)) - (complex-typed-gemv! (coerce-complex alpha) A x - (coerce-complex beta) y job)) - -(defmethod gemv! ((alpha number) (A real-matrix) (x real-vector) - (beta number) (y complex-vector) &optional (job :n)) - (unless (= beta 1) - (complex-typed-scal! (coerce-complex beta) y)) - (unless (= alpha 0) - (if (not (zerop (imagpart alpha))) - (let ((A.x (make-real-tensor (aref (dimensions y) 0))) - (vw-y (tensor-realpart~ y))) - (real-typed-gemv! (coerce-real 1) A x (coerce-real 0) A.x job) - ;; - (real-typed-axpy! (coerce-real (realpart alpha)) A.x vw-y) - ;;Move view to the imaginary part - (incf (head vw-y)) - (real-typed-axpy! (coerce-real (imagpart alpha)) A.x vw-y)) - (real-typed-gemv! (coerce-real alpha) A x - (coerce-real 1) (tensor-realpart~ y) job))) - y) - -(defmethod gemv! ((alpha number) (A real-matrix) (x complex-vector) - (beta number) (y complex-matrix) &optional (job :n)) - (unless (= beta 1) - (complex-typed-scal! (coerce-complex beta) y)) - (unless (= alpha 0) - (let ((A.x (make-complex-tensor (aref (dimensions y) 0)))) - (let ((vw-x (tensor-realpart~ x)) - (vw-A.x (tensor-realpart~ x))) - ;;Re - (real-typed-gemv! (coerce-real 1) A vw-x (coerce-real 0) vw-A.x job) - ;;Im - (incf (head vw-x)) - (incf (head vw-A.x)) - (real-typed-gemv! (coerce-real 1) A vw-x (coerce-real 0) vw-A.x job)) - (complex-typed-axpy! (coerce-complex alpha) A.x y))) - y) - -(defmethod gemv! ((alpha number) (A complex-matrix) (x real-vector) - (beta number) (y complex-vector) &optional (job :n)) - (let ((cplx-x (make-complex-tensor (aref (dimensions x) 0)))) - (real-typed-copy! x (tensor-realpart~ cplx-x)) - (complex-typed-gemv! (coerce-complex alpha) A cplx-x - (coerce-complex beta) y job)) - y) - +(defun ieql (&rest args) + (loop :for ele :in (cdr args) + :do (unless (eql (car args) ele) + (return nil)) + :finally (return t))) + +(defun >class (cls) + (let ((clist (reverse (sort cls #'coerceable?)))) + (loop :for ele :in (cdr clist) + :do (unless (coerceable? ele (car clist)) + (return nil)) + :finally (return (car clist))))) + +(defun tensor-coerce (ten cls &optional (duplicate? t)) + (let ((clname (if (typep cls 'standard-class) (class-name cls) cls))) + (if (and (not duplicate?) (typep ten clname)) ten + (copy! ten (zeros (dimensions ten) clname))))) + +(defmethod gemv! (alpha (A standard-tensor) (x standard-tensor) beta (y standard-tensor) &optional (job :n)) + (let ((clx (class-name (class-of x))) + (cly (class-name (class-of y))) + (cla (class-name (class-of y)))) + (assert (and (member cla *tensor-type-leaves*) + (member clx *tensor-type-leaves*) + (member cly *tensor-type-leaves*)) + nil 'tensor-abstract-class :tensor-class (list cla clx cly)) + (cond + ((ieql clx cly cla) + (compile-and-eval + `(defmethod gemv! (alpha (A ,cla) (x ,clx) beta (y ,cly) &optional (job :n)) + (let ((alpha (t/coerce ,(field-type clx) alpha)) + (beta (t/coerce ,(field-type clx) beta)) + (cjob (aref (symbol-name job) 0))) + (declare (type ,(field-type clx) alpha beta) + (type character trans-op)) + ,(recursive-append + (when (subtypep clx 'blas-numeric-tensor) + `(if (call-fortran? A (t/l2-lb ,cla)) + (let ((A-copy (if (blas-matrix-compatiblep A cjob) A + (let ((*default-stride-ordering* :col-major)) + (t/copy! (,cla ,cla) A (t/zeros ,clx (dimensions A))))))) + (multiple-value-bind (lda op maj) (blas-matrix-compatiblep A-copy cjob) + (declare (ignore maj)) + (t/blas-gemv! ,cla alpha A-copy lda + x (aref (the index-store-vector (strides x)) 0) + beta + y (aref (the index-store-vector (strides y)) 0) + op))))) + `(t/gemv! ,cla alpha A x beta y cjob))) + y)) + (gemv! alpha A x beta y job)) + (t + (error "Don't know how to apply axpy! to classes ~a, ~a." clx cly))))) ;;---------------------------------------------------------------;; (defgeneric gemv (alpha A x beta y &optional job) (:documentation @@ -174,26 +172,6 @@ :H alpha * transpose o conjugate(A) + beta * y ")) -(defmethod gemv ((alpha number) (A standard-matrix) (x standard-vector) - (beta number) (y complex-vector) &optional (job :n)) - (let ((result (copy y))) - (gemv! alpha A x 1d0 result job))) - -(defmethod gemv ((alpha number) (A standard-matrix) (x standard-vector) - (beta number) (y real-vector) &optional (job :n)) - (let ((result (if (or (complexp alpha) (complexp beta) - (typep A 'complex-matrix) (typep x 'complex-vector)) - (make-complex-tensor (aref (dimensions y) 0)) - (make-real-tensor (aref (dimensions y) 0))))) - (scal! y result) - (gemv! alpha A x beta result job))) - -(defmethod gemv ((alpha number) (A standard-matrix) (x standard-vector) - (beta (eql nil)) (y (eql nil)) &optional (job :n)) - (let ((result (apply - (if (or (complexp alpha) (complexp beta) - (typep A 'complex-matrix) (typep x 'complex-vector)) - #'make-complex-tensor - #'make-real-tensor) - (list (ecase job ((:n :c) (nrows A)) ((:t :h) (ncols A))))))) - (gemv! alpha A x 0 result job))) +(defmethod gemv (alpha (A standard-tensor) (x standard-tensor) + beta (y standard-tensor) &optional (job :n)) + (gemv! alpha A x beta (copy y) job)) diff --git a/src/level-3/gemm.lisp b/src/level-3/gemm.lisp index b24fc76..776406f 100644 --- a/src/level-3/gemm.lisp +++ b/src/level-3/gemm.lisp @@ -1,260 +1,56 @@ -;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*- -;;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;;; -;;; Copyright (c) 2000 The Regents of the University of California. -;;; All rights reserved. -;;; -;;; Permission is hereby granted, without written agreement and without -;;; license or royalty fees, to use, copy, modify, and distribute this -;;; software and its documentation for any purpose, provided that the -;;; above copyright notice and the following two paragraphs appear in all -;;; copies of this software. -;;; -;;; IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY -;;; FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES -;;; ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF -;;; THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF -;;; SUCH DAMAGE. -;;; -;;; THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, -;;; INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF -;;; MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE -;;; PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF -;;; CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, -;;; ENHANCEMENTS, OR MODIFICATIONS. -;;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; - (in-package #:matlisp) -(defmacro generate-typed-gemm! (func (tensor-class blas-gemm-func fortran-lb-parameter)) - (let* ((opt (if-ret (get-tensor-class-optimization-hashtable tensor-class) - (error 'tensor-cannot-find-optimization :tensor-class tensor-class))) - (matrix-class (getf opt :matrix)) - (blas? blas-gemm-func)) - `(progn - (eval-when (:compile-toplevel :load-toplevel :execute) - (let ((opt (get-tensor-class-optimization-hashtable ',tensor-class))) - (assert opt nil 'tensor-cannot-find-optimization :tensor-class ',tensor-class) - (setf (getf opt :gemm) ',func - (get-tensor-class-optimization ',tensor-class) opt))) - (defun ,func (alpha A B beta C job) - (declare (type ,(getf opt :element-type) alpha beta) - (type ,matrix-class A B C) - (type symbol job)) - ,(let - ((lisp-routine - `(let-typed ((nr-C (nrows C) :type index-type) - (nc-C (ncols C) :type index-type) - (dotl (ecase job-A (:n (ncols A)) (:t (nrows A))) :type index-type) - ; - (rstp-A (row-stride A) :type index-type) - (cstp-A (col-stride A) :type index-type) - (hd-A (head A) :type index-type) - ; - (rstp-B (row-stride B) :type index-type) - (cstp-B (col-stride B) :type index-type) - (hd-B (head B) :type index-type) - ; - (rstp-C (row-stride C) :type index-type) - (cstp-C (col-stride C) :type index-type) - (hd-C (head C) :type index-type)) - ;; - (when (eq job-A :t) - (rotatef rstp-A cstp-A)) - (when (eq job-B :t) - (rotatef rstp-B cstp-B)) - ;; - (unless (,(getf opt :f=) beta (,(getf opt :fid*))) - (,(getf opt :num-scal) beta C)) - ;;Most of these loop orderings are borrowed from the Fortran reference - ;;implementation of BLAS. - (cond - ((and (= cstp-C 1) (= cstp-B 1)) - (let-typed ((of-A hd-A :type index-type) - (of-B hd-B :type index-type) - (of-C hd-C :type index-type) - (d.rstp-B (- rstp-B nc-C) :type index-type) - (d.rstp-A (- rstp-A (* cstp-A dotl)) :type index-type) - (sto-A (store A) :type ,(linear-array-type (getf opt :store-type))) - (sto-B (store B) :type ,(linear-array-type (getf opt :store-type))) - (sto-C (store C) :type ,(linear-array-type (getf opt :store-type)))) - (very-quickly - (loop :repeat nr-C - :do (progn - (loop :repeat dotl - :do (let-typed - ((ele-A (,(getf opt :f*) alpha (,(getf opt :reader) sto-A of-A)) :type ,(getf opt :element-type))) - (loop :repeat nc-C - :do (progn - (,(getf opt :value-incfer) (,(getf opt :f*) ele-A (,(getf opt :reader) sto-B of-B)) - sto-C of-C) - (incf of-C) - (incf of-B))) - (decf of-C nc-C) - (incf of-A cstp-A) - (incf of-B d.rstp-B))) - (incf of-C rstp-C) - (incf of-A d.rstp-A) - (setf of-B hd-B)))))) - ((and (= cstp-A 1) (= rstp-B 1)) - (let-typed ((of-A hd-A :type index-type) - (of-B hd-B :type index-type) - (of-C hd-C :type index-type) - (d.cstp-B (- cstp-B dotl) :type index-type) - (d.rstp-C (- rstp-C (* nc-C cstp-C)) :type index-type) - (sto-A (store A) :type ,(linear-array-type (getf opt :store-type))) - (sto-B (store B) :type ,(linear-array-type (getf opt :store-type))) - (sto-C (store C) :type ,(linear-array-type (getf opt :store-type))) - (dot (,(getf opt :fid+)) :type ,(getf opt :element-type))) - (very-quickly - (loop :repeat nr-C - :do (progn - (loop :repeat nc-C - :do (progn - (setf dot (,(getf opt :fid+))) - (loop :repeat dotl - :do (progn - (setf dot (,(getf opt :f+) dot (,(getf opt :f*) (,(getf opt :reader) sto-A of-A) (,(getf opt :reader) sto-B of-B)))) - (incf of-A) - (incf of-B))) - (,(getf opt :value-incfer) dot sto-C of-C) - (incf of-C cstp-C) - (decf of-A dotl) - (incf of-B d.cstp-B))) - (incf of-C d.rstp-C) - (incf of-A rstp-A) - (setf of-B hd-B)))))) - ((and (= cstp-A 1) (= rstp-B 1)) - (let-typed ((of-A hd-A :type index-type) - (of-B hd-B :type index-type) - (of-C hd-C :type index-type) - (d.cstp-B (- cstp-B dotl) :type index-type) - (d.rstp-C (- rstp-C (* nc-C cstp-C)) :type index-type) - (sto-A (store A) :type ,(linear-array-type (getf opt :store-type))) - (sto-B (store B) :type ,(linear-array-type (getf opt :store-type))) - (sto-C (store C) :type ,(linear-array-type (getf opt :store-type))) - (dot (,(getf opt :fid+)) :type ,(getf opt :element-type))) - (very-quickly - (loop :repeat nr-C - :do (progn - (loop :repeat nc-C - :do (progn - (setf dot (,(getf opt :fid+))) - (loop :repeat dotl - :do (progn - (setf dot (,(getf opt :f+) dot (,(getf opt :f*) (,(getf opt :reader) sto-A of-A) (,(getf opt :reader) sto-B of-B)))) - (incf of-A) - (incf of-B))) - (,(getf opt :value-incfer) dot sto-C of-C) - (incf of-C cstp-C) - (decf of-A dotl) - (incf of-B d.cstp-B))) - (incf of-C d.rstp-C) - (incf of-A rstp-A) - (setf of-B hd-B)))))) - ((and (= rstp-A 1) (= rstp-C 1)) - (let-typed ((of-A hd-A :type index-type) - (of-B hd-B :type index-type) - (of-C hd-C :type index-type) - (d.cstp-B (- cstp-B (* rstp-B dotl)) :type index-type) - (d.cstp-A (- cstp-A nr-C) :type index-type) - (sto-A (store A) :type ,(linear-array-type (getf opt :store-type))) - (sto-B (store B) :type ,(linear-array-type (getf opt :store-type))) - (sto-C (store C) :type ,(linear-array-type (getf opt :store-type)))) - (very-quickly - (loop :repeat nc-C - :do (progn - (loop :repeat dotl - :do (let-typed - ((ele-B (,(getf opt :f*) alpha (,(getf opt :reader) sto-B of-B)) :type ,(getf opt :element-type))) - (loop :repeat nr-C - :do (progn - (,(getf opt :value-incfer) (,(getf opt :f*) ele-B (,(getf opt :reader) sto-A of-A)) - sto-C of-C) - (incf of-C) - (incf of-A))) - (decf of-C nr-C) - (incf of-A d.cstp-A) - (incf of-B rstp-B))) - (incf of-C cstp-C) - (setf of-A hd-A) - (incf of-B d.cstp-B)))))) - (t - (let-typed ((of-A hd-A :type index-type) - (of-B hd-B :type index-type) - (of-C hd-C :type index-type) - (r.cstp-C (* cstp-C nc-C) :type index-type) - (d.rstp-B (- rstp-B (* cstp-B nc-C)) :type index-type) - (d.rstp-A (- rstp-A (* cstp-A dotl)) :type index-type) - (sto-A (store A) :type ,(linear-array-type (getf opt :store-type))) - (sto-B (store B) :type ,(linear-array-type (getf opt :store-type))) - (sto-C (store C) :type ,(linear-array-type (getf opt :store-type)))) - (very-quickly - (loop :repeat nr-C - :do (progn - (loop :repeat dotl - :do (let-typed - ((ele-A (,(getf opt :f*) alpha (,(getf opt :reader) sto-A of-A)) :type ,(getf opt :element-type))) - (loop :repeat nc-C - :do (progn - (,(getf opt :value-writer) - (,(getf opt :f+) - (,(getf opt :reader) sto-C of-C) - (,(getf opt :f*) ele-A (,(getf opt :reader) sto-B of-B))) - sto-C of-C) - (incf of-C cstp-C) - (incf of-B cstp-B))) - (decf of-C r.cstp-C) - (incf of-A cstp-A) - (incf of-B d.rstp-B))) - (incf of-C rstp-C) - (incf of-A d.rstp-A) - (setf of-B hd-B)))))))))) - ;;Tie together Fortran and lisp-routines. - `(mlet* (((job-A job-B) (ecase job - (:nn (values :n :n)) - (:nt (values :n :t)) - (:tn (values :t :n)) - (:tt (values :t :t))) - :type (symbol symbol)) - ,@(when blas? - `((call-fortran? (> (max (nrows C) (ncols C) (if (eq job-A :n) (ncols A) (nrows A))) - ,fortran-lb-parameter)) - ((maj-A ld-A fop-A) (blas-matrix-compatible-p A job-A) :type (symbol index-type (string 1))) - ((maj-B ld-B fop-B) (blas-matrix-compatible-p B job-B) :type (symbol index-type (string 1))) - ((maj-C ld-C fop-C) (blas-matrix-compatible-p C :n) :type (symbol index-type nil))))) - ,(if blas? - `(cond - (call-fortran? - (if (and maj-A maj-B maj-C) - (let-typed ((nr-C (nrows C) :type index-type) - (nc-C (ncols C) :type index-type) - (dotl (ecase job-A (:n (ncols A)) (:t (nrows A))) :type index-type)) - (when (eq maj-C :row-major) - (rotatef A B) - (rotatef ld-A ld-B) - (rotatef maj-A maj-B) - (rotatef nr-C nc-C) - (setf (values fop-A fop-B) - (values (fortran-snop fop-B) (fortran-snop fop-A)))) - (,blas-gemm-func fop-A fop-B nr-C nc-C dotl - alpha (store A) ld-A (store B) ld-B - beta (store C) ld-C - (head A) (head B) (head C))) - (let ((ret (,func alpha - (if maj-A A (,(getf opt :copy) A (,(getf opt :zero-maker) (dimensions A)))) - (if maj-B B (,(getf opt :copy) B (,(getf opt :zero-maker) (dimensions B)))) - beta - (if maj-C C (,(getf opt :copy) C (,(getf opt :zero-maker) (dimensions C)))) - job))) - (unless maj-C - (,(getf opt :copy) ret C))))) - (t - ,lisp-routine)) - lisp-routine))) - C)))) +(deft/generic (t/blas-gemm-func #'subtypep) sym ()) +(deft/method t/blas-gemm-func (sym real-tensor) () + 'dgemm) +(deft/method t/blas-gemm-func (sym complex-tensor) () + 'zgemm) +;; +(deft/generic (t/blas-gemm! #'subtypep) sym (alpha A lda B ldb beta C ldc transa transb)) + +(deft/method t/blas-gemm! (sym blas-numeric-tensor) (alpha A lda B ldb beta C ldc transa transb) + (using-gensyms (decl (alpha A lda B ldb beta C ldc transa transb)) + (with-gensyms (m n k) + `(let* (,@decl + (,m (aref (the index-store-vector (dimensions ,C)) 0)) + (,n (aref (the index-store-vector (dimensions ,C)) 1)) + (,k (aref (the index-store-vector (dimensions ,A)) (ecase (char-upcase ,transa) ((#\N #\C) 1) ((#\T #\H) 0))))) + (declare (type ,sym ,A ,B ,C) + (type ,(field-type sym) ,alpha ,beta) + (type index-type ,lda ,ldb ,ldc ,m ,n ,k) + (type character ,transa ,transb)) + (,(macroexpand-1 `(t/blas-gemm-func ,sym)) + ,transa ,transb + ,m ,n ,k + ,alpha + (the ,(store-type sym) (store ,A)) ,lda + (the ,(store-type sym) (store ,B)) ,ldb + ,beta + (the ,(store-type sym) (store ,C)) ,ldc + (the index-type (head ,A)) (the index-type (head ,B)) (the index-type (head ,C))) + ,C)))) + +;; +(deft/generic (t/gemm! #'subtypep) sym (alpha A B beta C transa transb)) + +;;Witness the power of macros, muggles! :) +(deft/method t/gemm! (sym standard-tensor) (alpha A B beta C transa transb) + (using-gensyms (decl (alpha A x beta y transp)) + `(let (,@decl) + (declare (type ,sym ,A ,x ,y) + (type ,(field-type sym) ,alpha ,beta) + (type character ,transp)) + (unless (t/f= ,(field-type sym) ,beta (t/fid* ,(field-type sym))) + (t/scdi! ,sym ,beta ,y :scal? t :numx? t)) + ;;These loops are optimized for column major matrices + (ecase (char-upcase ,transp) + (#\N (einstein-sum ,sym (j i) (ref ,y i) (* ,alpha (ref ,A i j) (ref ,x j)) nil)) + (#\T (einstein-sum ,sym (i j) (ref ,y i) (* ,alpha (ref ,A j i) (ref ,x j)) nil)) + (#\C (einstein-sum ,sym (j i) (ref ,y i) (* ,alpha (t/fc ,(field-type sym) (ref ,A i j)) (ref ,x j)) nil)) + (#\H (einstein-sum ,sym (i j) (ref ,y i) (* ,alpha (t/fc ,(field-type sym) (ref ,A j i)) (ref ,x j)) nil))) + ,y))) + + ;;Real (generate-typed-gemm! real-base-typed-gemm! (real-tensor dgemm *real-l3-fcall-lb*)) diff --git a/tests/tcomp.lisp b/tests/tcomp.lisp index 8f7e990..d7bd73d 100644 --- a/tests/tcomp.lisp +++ b/tests/tcomp.lisp @@ -47,9 +47,28 @@ t) (defun mv-test (A b c) - (t/gemv! real-tensor 1d0 A b 0d0 c nil)) - + (t/gemv! real-tensor 1d0 A b 2d0 c #\t)) + +(let ((A (zeros '(1000 1000))) + (x (zeros 1000)) + (y (zeros 1000))) + (let-typed ((sto-x (store x) :type (simple-array double-float)) + (sto-y (store y) :type (simple-array double-float)) + (sto-a (store A) :type (simple-array double-float))) + (loop :for i :from 0 :below (array-dimension sto-x 0) + :do (setf (aref sto-x i) (random 1d0) + (aref sto-y i) (random 1d0))) + (loop :for i :from 0 :below (array-dimension sto-a 0) + :do (setf (aref sto-a i) (random 1d0)))) + (time (let ((*real-l2-fcall-lb* (* 1000 2000))) (gemv! 1 A x 1 y))) + t) + (let ((A (copy! #2a((1 2) (3 4)) (zeros '(2 2)))) (b (copy! 1 (zeros 2))) (c (copy! #(1 2) (zeros 2)))) (time (dotimes (i 1000) (mv-test A b c)))) + +(let ((A (copy! #2a((1 2) (3 4)) (zeros '(2 2)))) + (b (copy! 1 (zeros 2))) + (c (copy! #(1 2) (zeros 2)))) + (time (dotimes (i 1000) (gemv! 1 A b 0 c :n)))) commit be4148122456b5f7a6d4032fcba44e4652f4eb0b Author: Akshay Srinivasan <aks...@gm...> Date: Sat Aug 10 09:51:28 2013 -0700 Added gemv template using einstein. diff --git a/src/base/einstein.lisp b/src/base/einstein.lisp index ac0ada7..64998a5 100644 --- a/src/base/einstein.lisp +++ b/src/base/einstein.lisp @@ -69,6 +69,7 @@ tmp))) (values refs tlist indices))) +;;Add options (allow function to compile the clause ?) for more compiler options. (defun loop-generator-base (type index-order place clause &key (testp t) (tight-iloop nil)) (multiple-value-bind (refs tlist indices) (parse-loopx type place clause) (let* ((tens (mapcar #'(lambda (x) (second (getf x :tensor))) tlist)) @@ -219,8 +220,8 @@ (defmacro einstein-sum (type idx-order place clause &optional (testp t)) (loop-generator type idx-order place clause :testp testp)) -;;Yes this is slow, but if you're *really* worried about computation then roll your custom loops -;;with einstein-sum-base. This is a super-adaptive on-the-fly loop generation function generation +;;Yes this has an overhead, but if you're *really* worried about efficiency then roll your custom loops +;;with einstein-sum-base. This is a super-adaptive on-the-fly-loop-generating-function generating ;;macro. You have the power now, without any of the tedium :) (defmacro define-einstein-sum (name args (type place clause &optional (testp t))) (multiple-value-bind (refs tlist indices) (parse-loopx type place clause) @@ -236,7 +237,7 @@ (let* ((code (loop-generator ',type idx-ord ',place ',clause :testp ,testp)) (funcnew (compile-and-eval (list 'lambda '(,@args) code)))) - (format t "Compiling code for index-order : ~a~%" idx-ord) + (format t "~a: Compiling code for index-order : ~a~%" ,(symbol-name name) idx-ord) (setf (gethash idx-ord ,functable) funcnew) funcnew)))) (apply func (list ,@args)))))))) diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp index 29c6cbb..28a82d0 100644 --- a/src/level-2/gemv.lisp +++ b/src/level-2/gemv.lisp @@ -9,194 +9,48 @@ (deft/generic (t/blas-gemv! #'subtypep) sym (alpha A lda x st-x beta y st-y transp)) (deft/method t/blas-gemv! (sym blas-numeric-tensor) (alpha A lda x st-x beta y st-y transp) - (using-gensyms (decl (alpha A lda x st-x beta y st-y transp)) - (with-gensyms (m n) - `(let* (,@decl - (,m (aref (the index-store-vector (dimensions ,A)) 0)) - (,n (aref (the index-store-vector (dimensions ,A)) 1))) - (declare (type ,sym ,A ,x ,y) - (type ,(field-type sym) ,alpha ,beta) - (type index-type ,st-x ,st-y ,lda ,m ,n)) - (when (cl:char= (char-upcase ,transp) #\T) - (rotatef ,m ,n)) - (,(macroexpand-1 `(t/blas-gemv-func ,sym)) - ,transp ,m ,n - ,alpha - (the ,(store-type sym) (store ,A)) ,lda - (the ,(store-type sym) (store ,x)) ,st-x - ,beta - (the ,(store-type sym) (store ,y)) ,st-y - (the index-type (head ,A)) (the index-type (head ,x)) (the index-type (head ,y))) - y)))) + (using-gensyms (decl (alpha A lda x st-x beta y st-y transp)) + (with-gensyms (m n) + `(let* (,@decl + (,m (aref (the index-store-vector (dimensions ,A)) 0)) + (,n (aref (the index-store-vector (dimensions ,A)) 1))) + (declare (type ,sym ,A ,x ,y) + (type ,(field-type sym) ,alpha ,beta) + (type index-type ,st-x ,st-y ,lda ,m ,n)) + (when (cl:char= (char-upcase ,transp) #\T) + (rotatef ,m ,n)) + (,(macroexpand-1 `(t/blas-gemv-func ,sym)) + ,transp ,m ,n + ,alpha + (the ,(store-type sym) (store ,A)) ,lda + (the ,(store-type sym) (store ,x)) ,st-x + ,beta + (the ,(store-type sym) (store ,y)) ,st-y + (the index-type (head ,A)) (the index-type (head ,x)) (the index-type (head ,y))) + y)))) ;; (deft/generic (t/gemv! #'subtypep) sym (alpha A x beta y transp)) -(deft/method t/gemv! (sym standard-tensor) (alpha A x beta y transp)) - -(let ((A (copy! #2a((2 1) (3 4)) (zeros '(2 2 )))) - (x (copy! #(1 2) (zeros 2))) - (y (copy! #(0 1) (zeros 2)))) - (t/blas-gemv! real-tensor 1d0 A 2 x 1 0d0 y 1 t)) - - -(deft/method t/blas-axpy! (sym blas-numeric-tensor) (a x st-x y st-y) - (let ((apy? (null x))) - (using-gensyms (decl (a x y)) - (with-gensyms (sto-x stp-x) - `(let (,@decl) - (declare (type ,sym ,@(unless apy? `(,x)) ,y) - ,@(when apy? `((ignore ,x)))) - (let ((,sto-x ,(if apy? `(t/store-allocator ,sym 1) `(store ,x))) - (,stp-x ,(if apy? 0 st-x))) - (declare (type ,(store-type sym) ,sto-x) - (type index-type ,stp-x)) - ,@(when apy? - `((t/store-set ,sym (t/fid* ,(field-type sym)) ,sto-x 0))) - (,(macroexpand-1 `(t/blas-axpy-func ,sym)) - (the index-type (size ,y)) - (the ,(field-type sym) ,a) - ,sto-x ,stp-x - (the ,(store-type sym) (store ,y)) (the index-type ,st-y) - ,(if apy? 0 `(head ,x)) (head ,y)) - ,y)))))) - -(deft/generic (t/axpy! #'subtypep) sym (a x y)) -(deft/method t/axpy! (sym standard-tensor) (a x y) - (let ((apy? (null x))) - (using-gensyms (decl (a x y)) - (with-gensyms (idx sto-x sto-y of-x of-y) - `(let (,@decl) - (declare (type ,sym ,@(unless apy? `(,x)) ,y) - (type ,(field-type sym) ,a) - ,@(when apy? `((ignore ,x)))) - (let (,@(unless apy? `((,sto-x (store ,x)))) - (,sto-y (store ,y))) - (declare (type ,(store-type sym) ,@(unless apy? `(,sto-x)) ,sto-y)) - (very-quickly - (mod-dotimes (,idx (dimensions ,y)) - :with (linear-sums - ,@(unless apy? `((,of-x (strides ,x) (head ,x)))) - (,of-y (strides ,y) (head ,y))) - :do (t/store-set ,sym (t/f+ ,(field-type sym) - ,@(if apy? - `(,a) - `((t/f* ,(field-type sym) - ,a (t/store-ref ,sym ,sto-x ,of-x)))) - (t/store-ref ,sym ,sto-y ,of-y)) - ,sto-y ,of-y))) - ,y)))))) - - -(defmacro generate-typed-gemv! (func - (tensor-class blas-gemv-func - fortran-call-lb)) - (let* ((opt (if-ret (get-tensor-class-optimization-hashtable tensor-class) - (error 'tensor-cannot-find-optimization :tensor-class tensor-class))) - (matrix-class (getf opt :matrix)) - (vector-class (getf opt :vector))) - `(progn - (eval-when (:compile-toplevel :load-toplevel :execute) - (let ((opt (get-tensor-class-optimization-hashtable ',tensor-class))) - (assert opt nil 'tensor-cannot-find-optimization :tensor-class ',tensor-class) - (setf (getf opt :gemv) ',func - (get-tensor-class-optimization ',tensor-class) opt))) - (defun ,func (alpha A x beta y job) - (declare (type ,(getf opt :element-type) alpha beta) - (type ,matrix-class A) - (type ,vector-class x y) - (type list job)) - ,(let - ((lisp-routine - `(let-typed ((nr-A (nrows A) :type index-type) - (nc-A (ncols A) :type index-type) - (rs-A (row-stride A) :type index-type) - (cs-A (col-stride A) :type index-type) - (sto-A (store A) :type ,(linear-array-type (getf opt :store-type))) - ; - (stp-x (aref (strides x) 0) :type index-type) - (sto-x (store x) :type ,(linear-array-type (getf opt :store-type))) - (hd-x (head x) :type index-type) - ; - (stp-y (aref (strides y) 0) :type index-type) - (sto-y (store y) :type ,(linear-array-type (getf opt :store-type))) - ; - (job (car job) :type character)) - (when (char= job #\T) - (rotatef nr-A nc-A) - (rotatef rs-A cs-A)) - (very-quickly - (loop :repeat nr-A - :for of-y :of-type index-type := (head y) :then (+ of-y stp-y) - :for rof-A :of-type index-type := (head A) :then (+ rof-A rs-A) - :do (let-typed ((val (,(getf opt :f*) beta (,(getf opt :reader) sto-y of-y)) :type ,(getf opt :element-type))) - (loop :repeat nc-A - :for of-x :of-type index-type := hd-x :then (+ of-x stp-x) - :for of-A :of-type index-type := rof-A :then (+ of-A cs-A) - :with dot :of-type ,(getf opt :element-type) = (,(getf opt :fid+)) - :do (let-typed ((xval (,(getf opt :reader) sto-x of-x) :type ,(getf opt :element-type)) - (Aval (,(getf opt :reader) sto-A of-A) :type ,(getf opt :element-type))) - (setf dot (,(getf opt :f+) dot (,(getf opt :f*) xval Aval)))) - :finally (,(getf opt :value-writer) (,(getf opt :f+) (,(getf opt :f*) alpha dot) val) sto-y of-y)))))))) - (if blas-gemv-func - `(mlet* - ((call-fortran? (> (max (nrows A) (ncols A)) ,fortran-call-lb)) - ((maj-A ld-A fop-A) (blas-matrix-compatible-p A job) :type (symbol index-type character))) - (cond - (call-fortran? - (if maj-A - (let-typed ((nr-A (nrows A) :type index-type) - (nc-A (ncols A) :type index-type)) - (when (eq maj-A :row-major) - (rotatef nr-A nc-A)) - (,blas-gemv-func fop-a nr-A nc-A - alpha (store A) ld-A - (store x) (aref (strides x) 0) - beta - (store y) (aref (strides y) 0) - (head A) (head x) (head y))) - (,func alpha (,(getf opt :copy) A (,(getf opt :zero-maker) (dimensions A))) x beta y job))) - (t - ,lisp-routine))) - lisp-routine)) - y)))) - -;;Real -(generate-typed-gemv! real-base-typed-gemv! - (real-tensor dgemv *real-l2-fcall-lb*)) - -(definline real-typed-gemv! (alpha A x beta y job) - (real-base-typed-gemv! alpha A x beta y (ecase job ((:n :t) job) (:h :t) (:c :n)))) - -;;Complex -(generate-typed-gemv! complex-base-typed-gemv! - (complex-tensor zgemv *complex-l2-fcall-lb*)) - -(definline complex-typed-gemv! (alpha A x beta y job) - (declare (type complex-matrix A) - (type complex-vector x y) - (type complex-type alpha beta) - (type symbol job)) - (if (member job '(:n :t)) - (complex-base-typed-gemv! alpha A x beta y job) - ;; - (let-typed ((cx (let ((ret (complex-typed-copy! x (complex-typed-zeros (dimensions x))))) - (real-typed-num-scal! -1d0 (tensor-imagpart~ ret)) - ret) :type complex-vector) - (y.view (tensor-imagpart~ y))) - (real-typed-num-scal! -1d0 y.view) - (complex-base-typed-gemv! (cl:conjugate alpha) A cx - (cl:conjugate beta) y (ecase job (:h :t) (:c :n))) - (real-typed-num-scal! -1d0 y.view))) - y) +;;Witness the power of macros, muggles! :) +(deft/method t/gemv! (sym standard-tensor) (alpha A x beta y transp) + (using-gensyms (decl (alpha A x beta y transp)) + `(let (,@decl) + (declare (type ,sym ,A ,x ,y) + (type ,(field-type sym) ,alpha ,beta)) + (unless (t/f= ,(field-type sym) ,beta (t/fid* ,(field-type sym))) + (t/scdi! ,sym ,beta ,y :scal? t :numx? t)) + ;;These loops are optimized for column major matrices + (if ,transp + (einstein-sum ,sym (i j) (ref ,y i) (* ,alpha (ref ,A j i) (ref ,x j)) nil) + (einstein-sum ,sym (j i) (ref ,y i) (* ,alpha (ref ,A i j) (ref ,x j)) nil)) + ,y))) ;;Symbolic #+maxima (generate-typed-gemv! symbolic-base-typed-gemv! (symbolic-tensor nil 0)) - ;;---------------------------------------------------------------;; - (defgeneric gemv! (alpha A x beta y &optional job) (:documentation " diff --git a/src/utilities/macros.lisp b/src/utilities/macros.lisp index 738ab74..f19868f 100644 --- a/src/utilities/macros.lisp +++ b/src/utilities/macros.lisp @@ -462,6 +462,14 @@ Macro which encloses @arg{forms} inside (declare (optimize (speed 3) (safety 0) (space 0))) " + #+matlisp-debug + `(with-optimization + #+lispworks + (:safety 3) + #-lispworks + (:safety 3) + ,@forms) + #-matlisp-debug `(with-optimization #+lispworks (:safety 0 :space 0 :speed 3 :float 0 :fixnum-safety 0) diff --git a/tests/tcomp.lisp b/tests/tcomp.lisp index 36c0683..8f7e990 100644 --- a/tests/tcomp.lisp +++ b/tests/tcomp.lisp @@ -45,3 +45,11 @@ (aref sto-y i) (random 1d0)))) (time (mm-test x y z)) t) + +(defun mv-test (A b c) + (t/gemv! real-tensor 1d0 A b 0d0 c nil)) + +(let ((A (copy! #2a((1 2) (3 4)) (zeros '(2 2)))) + (b (copy! 1 (zeros 2))) + (c (copy! #(1 2) (zeros 2)))) + (time (dotimes (i 1000) (mv-test A b c)))) ----------------------------------------------------------------------- Summary of changes: matlisp.asd | 5 +- src/base/blas-helpers.lisp | 4 +- src/base/einstein.lisp | 7 +- src/base/template.lisp | 4 +- src/classes/numeric.lisp | 1 + src/level-1/axpy.lisp | 2 - src/level-1/dot.lisp | 6 - src/level-1/scal.lisp | 4 - src/level-2/gemv.lisp | 362 ++++++++++++-------------------------------- src/level-3/gemm.lisp | 306 ++++++------------------------------- src/utilities/macros.lisp | 8 + tests/tcomp.lisp | 27 ++++ 12 files changed, 193 insertions(+), 543 deletions(-) hooks/post-receive -- matlisp |