You can subscribe to this list here.
2012 |
Jan
|
Feb
|
Mar
(34) |
Apr
(4) |
May
(2) |
Jun
(11) |
Jul
(22) |
Aug
(9) |
Sep
|
Oct
|
Nov
|
Dec
(4) |
---|---|---|---|---|---|---|---|---|---|---|---|---|
2013 |
Jan
(15) |
Feb
(17) |
Mar
(3) |
Apr
|
May
|
Jun
(3) |
Jul
(1) |
Aug
(5) |
Sep
(5) |
Oct
(2) |
Nov
(1) |
Dec
(1) |
2014 |
Jan
|
Feb
(1) |
Mar
(1) |
Apr
|
May
(1) |
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
2016 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
(1) |
Sep
|
Oct
|
Nov
|
Dec
|
From: Akshay S. <ak...@us...> - 2013-02-26 06:33:59
|
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, tensor has been updated via 8d1294ce927280b8f2633b1e8636bd8aacfbf45d (commit) from 36eb4f85eaf7042a8f2bcf7a1af1f48dddd7bb45 (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 8d1294ce927280b8f2633b1e8636bd8aacfbf45d Author: Akshay Srinivasan <aks...@gm...> Date: Mon Feb 25 22:33:15 2013 -0800 Minor changes to gnuplot.lisp. diff --git a/lib-src/gnuplot/gnuplot.lisp b/lib-src/gnuplot/gnuplot.lisp index e895bb4..4894a43 100644 --- a/lib-src/gnuplot/gnuplot.lisp +++ b/lib-src/gnuplot/gnuplot.lisp @@ -8,24 +8,30 @@ ccl:run-program gnuplot-binary nil :input :stream :wait nil :output t)) -(defun plot2d (data &key (color (list "#FF0000")) (stream (#+:sbcl - sb-ext:process-input - #+:ccl - ccl:external-process-input-stream - *current-gnuplot-process*))) - (with-open-file (s "/tmp/matlisp-gnuplot.out" :direction :output :if-exists :supersede :if-does-not-exist :create) - (loop :for i :from 0 :below (loop :for x :in data :minimizing (number-of-elements x)) - :do (loop :for x :in data :do (format s "~a " (tensor-ref x i)) :finally (format s "~%")))) - (format stream "plot '/tmp/matlisp-gnuplot.out' with lines linecolor rgb ~s~%" color) - (finish-output stream)) +(defun plot2d (data &key (color (list "#FF0000"))) + (unless *current-gnuplot-process* + (setf *current-gnuplot-process* (open-gnuplot-stream))) + (let ((stream (#+:sbcl + sb-ext:process-input + #+:ccl + ccl:external-process-input-stream + *current-gnuplot-process*))) + (with-open-file (s "/tmp/matlisp-gnuplot.out" :direction :output :if-exists :supersede :if-does-not-exist :create) + (loop :for i :from 0 :below (loop :for x :in data :minimizing (number-of-elements x)) + :do (loop :for x :in data :do (format s "~a " (tensor-ref x i)) :finally (format s "~%")))) + (format stream "plot '/tmp/matlisp-gnuplot.out' with lines linecolor rgb ~s~%" color) + (finish-output stream))) -(defun gnuplot-send (str &key (stream (#+:sbcl - sb-ext:process-input - #+:ccl - ccl:external-process-input-stream - *current-gnuplot-process*))) - (format stream "~a~%" str) - (finish-output stream)) +(defun gnuplot-send (str) + (unless *current-gnuplot-process* + (setf *current-gnuplot-process* (open-gnuplot-stream))) + (let (stream (#+:sbcl + sb-ext:process-input + #+:ccl + ccl:external-process-input-stream + *current-gnuplot-process*)) + (format stream "~a~%" str) + (finish-output stream))) ;; (defclass gnuplot-plot-info () ----------------------------------------------------------------------- Summary of changes: lib-src/gnuplot/gnuplot.lisp | 40 +++++++++++++++++++++++----------------- 1 files changed, 23 insertions(+), 17 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-25 18:00:07
|
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, tensor has been updated via 36eb4f85eaf7042a8f2bcf7a1af1f48dddd7bb45 (commit) from 6d64374a148b18ceea55608f39b0bc0d5f8cca75 (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 36eb4f85eaf7042a8f2bcf7a1af1f48dddd7bb45 Author: Akshay Srinivasan <aks...@gm...> Date: Mon Feb 25 09:54:59 2013 -0800 Added seq.lisp to the asdf load file. diff --git a/matlisp.asd b/matlisp.asd index e83b2ee..09ec278 100644 --- a/matlisp.asd +++ b/matlisp.asd @@ -157,11 +157,13 @@ :pathname "lapack" :depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3") :components ((:file "getrf"))) - #+nil (:module "matlisp-sugar" :pathname "sugar" :depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3") - :components ((:file "mplusminus") + :components ((:file "seq") + #+nil + (:file "mplusminus") + #+nil (:file "mtimesdivide"))) #+nil (:module "matlisp-reader" diff --git a/src/sugar/seq.lisp b/src/sugar/seq.lisp index 9da94ca..9c96efa 100644 --- a/src/sugar/seq.lisp +++ b/src/sugar/seq.lisp @@ -3,14 +3,14 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; ;;; Copyright (c) 2000 The Regents of the University of California. -;;; All rights reserved. -;;; +;;; 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 @@ -38,60 +38,19 @@ :do (setf (aref sto-r i) ori)) ret)))) -(defun linspace (start end &optional (num-points (1+ (abs (- start end))))) +(defun alinspace (start end &optional (num-points (1+ (abs (- start end))))) (let ((h (/ (- end start) (1- num-points)))) (arange start (+ h end) (abs h)))) -#+nil (export '(seq)) - -(if (not (fboundp '%push-on-end%)) -(defmacro %push-on-end% (value location) - `(setf ,location (nconc ,location (list ,value))))) - - -(defun seq (start step &optional end) - " - Syntax - ====== - (SEQ start [step] end) - - Purpose - ======= - Creates a list containing the sequence - - START, START+STEP, ..., START+(N-1)*STEP - - where | END-START | - N = | --------- | + 1 - | STEP | - -- -- - - The representations of START,STEP,END are - assumed to be exact (i.e. the arguments - are rationalized. - - The type of the elements in the - sequence are of the same type as STEP. - - The optional argument STEP defaults to 1. -" - (when (not end) - (setq end step) - (setq step 1)) - - (let ((start (rationalize start)) - (type (if (typep step 'integer) - 'integer - 'rational)) - (step (rationalize step)) - (end (rationalize end)) - (seq nil)) +(defun range (start end &optional (h 1)) + (declare (type real start end h)) + (let ((quo (ceiling (if (> start end) (- start end) (- end start)) h))) + (if (= quo 0) nil + (let ((h (if (> start end) (- h) h))) + (loop :for i :from 0 :below quo + :for ori := start :then (+ ori h) + :collect ori))))) - (when (zerop step) - (error "STEP equal to 0")) - (do ((x start (+ x step))) - ((if (> step 0) - (> x end) - (< x end)) - (nreverse seq)) - (push (coerce x type) seq)))) +(defun linspace (start end &optional (num-points (1+ (abs (- start end))))) + (let ((h (/ (- end start) (1- num-points)))) + (range start (+ h end) (abs h)))) ----------------------------------------------------------------------- Summary of changes: matlisp.asd | 6 +++- src/sugar/seq.lisp | 71 +++++++++++----------------------------------------- 2 files changed, 19 insertions(+), 58 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-25 09:41:25
|
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, tensor has been updated via 6d64374a148b18ceea55608f39b0bc0d5f8cca75 (commit) from 7338eb40c61c8fcab6a6180f68c0db7e013ad9ee (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 6d64374a148b18ceea55608f39b0bc0d5f8cca75 Author: Akshay Srinivasan <aks...@gm...> Date: Mon Feb 25 01:40:50 2013 -0800 Added CCL support to gnuplot.lisp diff --git a/lib-src/gnuplot/gnuplot.lisp b/lib-src/gnuplot/gnuplot.lisp index 906a981..e895bb4 100644 --- a/lib-src/gnuplot/gnuplot.lisp +++ b/lib-src/gnuplot/gnuplot.lisp @@ -1,15 +1,18 @@ (in-package :matlisp) -(defvar *current-gnuplot-stream* nil) -(defvar *gnuplot-binary* "/usr/bin/gnuplot") +(defvar *current-gnuplot-process* nil) -(defun open-gnuplot-stream () +(defun open-gnuplot-stream (&optional (gnuplot-binary (pathname "/usr/bin/gnuplot"))) (#+:sbcl sb-ext:run-program - *gnuplot-binary* nil :input :stream :wait nil :output t)) + #+:ccl + ccl:run-program + gnuplot-binary nil :input :stream :wait nil :output t)) (defun plot2d (data &key (color (list "#FF0000")) (stream (#+:sbcl - sb-ext:process-input - *current-gnuplot-stream*))) + sb-ext:process-input + #+:ccl + ccl:external-process-input-stream + *current-gnuplot-process*))) (with-open-file (s "/tmp/matlisp-gnuplot.out" :direction :output :if-exists :supersede :if-does-not-exist :create) (loop :for i :from 0 :below (loop :for x :in data :minimizing (number-of-elements x)) :do (loop :for x :in data :do (format s "~a " (tensor-ref x i)) :finally (format s "~%")))) @@ -18,7 +21,9 @@ (defun gnuplot-send (str &key (stream (#+:sbcl sb-ext:process-input - *current-gnuplot-stream*))) + #+:ccl + ccl:external-process-input-stream + *current-gnuplot-process*))) (format stream "~a~%" str) (finish-output stream)) ----------------------------------------------------------------------- Summary of changes: lib-src/gnuplot/gnuplot.lisp | 19 ++++++++++++------- 1 files changed, 12 insertions(+), 7 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-19 04:13:00
|
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, tensor has been updated via 7338eb40c61c8fcab6a6180f68c0db7e013ad9ee (commit) from 039adc9fd4f65be5beea4d7f077a54892c3310ca (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 7338eb40c61c8fcab6a6180f68c0db7e013ad9ee Author: Akshay Srinivasan <aks...@gm...> Date: Mon Feb 18 20:12:29 2013 -0800 First stab at getting gnuplot to work. diff --git a/lib-src/gnuplot/gnuplot.lisp b/lib-src/gnuplot/gnuplot.lisp index a329ccb..906a981 100644 --- a/lib-src/gnuplot/gnuplot.lisp +++ b/lib-src/gnuplot/gnuplot.lisp @@ -1,60 +1,71 @@ +(in-package :matlisp) (defvar *current-gnuplot-stream* nil) - -(defvar *gnuplot-binary* "gnuplot") - -(defclass gnuplot-plot-info () - ((title - :initform "GNU PLOT" - :accessor gnuplot-title) - (x-label - :initform "X" - :accessor gnuplot-x-label) - (y-label - :initform "Y" - :accessor gnuplot-y-label) - (x-data - :accessor gnuplot-x-data) - (y-data - :accessor gnuplot-y-data) - (z-data - :accessor gnuplot-z-data))) +(defvar *gnuplot-binary* "/usr/bin/gnuplot") (defun open-gnuplot-stream () - (#-:sbcl - ext:run-program - #+:sbcl + (#+:sbcl sb-ext:run-program *gnuplot-binary* nil :input :stream :wait nil :output t)) -(defun gnuplot-plot (info &key (stream (#-:sbcl - ext:process-input - #+:sbcl - sb-ext:process-input - *current-gnuplot-stream*))) - (with-accessors ((title gnuplot-title) - (x-label gnuplot-x-label) - (y-label gnuplot-y-label) - (x-data gnuplot-x-data) - (y-data gnuplot-y-data) - (z-data gnuplot-z-data)) - info - (format stream "~&set title '~S'~%" title) - (format stream "~&set xlabel '~S'~%" x-label) - (format stream "~&set ylabel '~S'~%" y-label) - (finish-output stream) - (map nil #'(lambda (x y z) - (with-open-file (s "/tmp/gnuplot.out" :direction :output - :if-exists :overwrite) - (map nil #'(lambda (xs ys zs) - (if zs - (format s "~A ~A ~A~%" xs ys zs) - (format s "~A ~A~%" xs ys))) - x y z) - (format stream "~A '/tmp/gnuplot.out'~%" - (if z "splot" "plot")) - (finish-output stream) - (sleep 5))) - x-data y-data z-data) - )) - - +(defun plot2d (data &key (color (list "#FF0000")) (stream (#+:sbcl + sb-ext:process-input + *current-gnuplot-stream*))) + (with-open-file (s "/tmp/matlisp-gnuplot.out" :direction :output :if-exists :supersede :if-does-not-exist :create) + (loop :for i :from 0 :below (loop :for x :in data :minimizing (number-of-elements x)) + :do (loop :for x :in data :do (format s "~a " (tensor-ref x i)) :finally (format s "~%")))) + (format stream "plot '/tmp/matlisp-gnuplot.out' with lines linecolor rgb ~s~%" color) + (finish-output stream)) + +(defun gnuplot-send (str &key (stream (#+:sbcl + sb-ext:process-input + *current-gnuplot-stream*))) + (format stream "~a~%" str) + (finish-output stream)) + + +;; (defclass gnuplot-plot-info () +;; ((title +;; :initform "GNU PLOT" +;; :accessor gnuplot-title) +;; (x-label +;; :initform "X" +;; :accessor gnuplot-x-label) +;; (y-label +;; :initform "Y" +;; :accessor gnuplot-y-label) +;; (x-data +;; :accessor gnuplot-x-data) +;; (y-data +;; :accessor gnuplot-y-data) +;; (z-data +;; :accessor gnuplot-z-data))) + + +;; (defun gnuplot-plot (info &key (stream (#+:sbcl +;; sb-ext:process-input +;; *current-gnuplot-stream*))) +;; (with-accessors ((title gnuplot-title) +;; (x-label gnuplot-x-label) +;; (y-label gnuplot-y-label) +;; (x-data gnuplot-x-data) +;; (y-data gnuplot-y-data) +;; (z-data gnuplot-z-data)) +;; info +;; (format stream "~&set title '~S'~%" title) +;; (format stream "~&set xlabel '~S'~%" x-label) +;; (format stream "~&set ylabel '~S'~%" y-label) +;; (finish-output stream) +;; (map nil #'(lambda (x y z) +;; (with-open-file (s "/tmp/gnuplot.out" :direction :output +;; :if-exists :overwrite) +;; (map nil #'(lambda (xs ys zs) +;; (if zs +;; (format s "~A ~A ~A~%" xs ys zs) +;; (format s "~A ~A~%" xs ys))) +;; x y z) +;; (format stream "~A '/tmp/gnuplot.out'~%" +;; (if z "splot" "plot")) +;; (finish-output stream) +;; (sleep 5))) +;; x-data y-data z-data) +;; )) diff --git a/src/base/print.lisp b/src/base/print.lisp index 7200306..f8ad402 100644 --- a/src/base/print.lisp +++ b/src/base/print.lisp @@ -90,7 +90,7 @@ of a matrix (default 0) (1 (format stream (format-to-string "~~~AT" *print-indent*)) (dotimes (i (aref dims 0)) - (if (< i *print-max-len*) + (if (or (eq *print-max-len* t) (< i *print-max-len*)) (progn (print-element tensor (tensor-ref tensor `(,i)) stream) (format stream "~,4T")) diff --git a/src/old/seq.lisp b/src/sugar/seq.lisp similarity index 67% rename from src/old/seq.lisp rename to src/sugar/seq.lisp index 96eb780..9da94ca 100644 --- a/src/old/seq.lisp +++ b/src/sugar/seq.lisp @@ -25,41 +25,22 @@ ;;; ENHANCEMENTS, OR MODIFICATIONS. ;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;;; -;;; Originally written by Tunc Simsek, Univ. of California, Berkeley, -;;; 2000, si...@ee... -;;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;;; -;;; $Id: seq.lisp,v 1.5 2004/12/03 18:01:38 rtoy Exp $ -;;; -;;; $Log: seq.lisp,v $ -;;; Revision 1.5 2004/12/03 18:01:38 rtoy -;;; We were doing (type-of step) and doing (coerce x type). But (type-of -;;; 1) can be (integer 1 1), and that doesn't work so well with coerce. -;;; So make type be either integer or rational, as appropriate. -;;; -;;; But the real question is why we need to do this at all. -;;; -;;; Revision 1.4 2002/06/24 18:05:30 rtoy -;;; Modified SEQ to run much faster by pushing the values onto the -;;; beginning of the list and reversing the result instead of -;;; destructively appending to the end, for an O(n^2) process instead of -;;; O(n). -;;; -;;; Revision 1.3 2000/07/11 18:02:03 simsek -;;; o Added credits -;;; -;;; Revision 1.2 2000/07/11 02:11:56 simsek -;;; o Added support for Allegro CL -;;; -;;; Revision 1.1 2000/04/14 00:12:48 simsek -;;; Initial revision. -;;; -;;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +(in-package #:matlisp) + +(defun arange (start end &optional (h 1d0)) + (let ((quo (ceiling (if (> start end) (- start end) (- end start)) h))) + (if (= quo 0) nil + (let*-typed ((ret (real-typed-zeros (idxv quo)) :type real-tensor) + (sto-r (store ret) :type real-store-vector) + (h (coerce-real-unforgiving (if (> start end) (- h) h)) :type real-type)) + (loop :for i :from 0 :below quo + :for ori := (coerce-real-unforgiving start) :then (+ ori h) + :do (setf (aref sto-r i) ori)) + ret)))) -(in-package "MATLISP") +(defun linspace (start end &optional (num-points (1+ (abs (- start end))))) + (let ((h (/ (- end start) (1- num-points)))) + (arange start (+ h end) (abs h)))) #+nil (export '(seq)) ----------------------------------------------------------------------- Summary of changes: lib-src/gnuplot/gnuplot.lisp | 119 +++++++++++++++++++++++------------------- src/base/print.lisp | 2 +- src/{old => sugar}/seq.lisp | 49 +++++------------ 3 files changed, 81 insertions(+), 89 deletions(-) rename src/{old => sugar}/seq.lisp (67%) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-15 02:54:12
|
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, master has been updated via c213febdfa60e0b1a9a11c796911eb5b93fef90e (commit) from 9980ae3686cf6361c2e8d8dec95d85f355b3a5d8 (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 c213febdfa60e0b1a9a11c796911eb5b93fef90e Author: Akshay Srinivasan <aks...@gm...> Date: Thu Feb 14 18:49:07 2013 -0800 Fixed path bug in configure.ac (line 380). diff --git a/configure.ac b/configure.ac index 9bbb440..035ac06 100644 --- a/configure.ac +++ b/configure.ac @@ -377,7 +377,7 @@ int main() EOF $CC $CFLAGS -c conftest.c $F77 $FFLAGS -o a.out conftest.o -L${ATLAS_DIR} -latlas -lcblas -lf77blas -llapack - if a.out; then + if ./a.out; then AC_MSG_RESULT([yes]) F2C=-ff2c else ----------------------------------------------------------------------- Summary of changes: configure.ac | 2 +- 1 files changed, 1 insertions(+), 1 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-11 07:44:16
|
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, tensor has been updated via 039adc9fd4f65be5beea4d7f077a54892c3310ca (commit) from f01e8cb0e33110ab72e6e828f052e7c534c84ebc (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 039adc9fd4f65be5beea4d7f077a54892c3310ca Author: Akshay Srinivasan <aks...@gm...> Date: Sun Feb 10 23:39:21 2013 -0800 Added a usability test. Code still very much resembles Fortran/C (with quirks), must get the Inline parser to work. diff --git a/tests/usability.lisp b/tests/usability.lisp new file mode 100644 index 0000000..8c0df45 --- /dev/null +++ b/tests/usability.lisp @@ -0,0 +1,39 @@ +(defpackage #:mpc + (:use (:cl :matlisp))) + +(defun pcart-field (time y &optional (ydot (real-typed-zeros (dimensions y)))) + (declare (ignore time)) + (assert (= 4 (length (store y)) (length (store ydot))) nil "nooo!") + (let*-typed ((sto-y (store y) :type real-store-vector) + (sto-yd (store ydot) :type real-store-vector) + (theta (aref sto-y 1) :type double-float) + (xdot (aref sto-y 2) :type double-float) + (thetadot (aref sto-y 3) :type double-float)) + (very-quickly + (setf (aref sto-yd 0) xdot + (aref sto-yd 1) thetadot + (aref sto-yd 2) (/ (+ (* (cos theta) (sin theta)) (* (sin theta) (expt thetadot 2))) (- 2 (expt (cos theta) 2))) + (aref sto-yd 3) (/ (+ (* 2 (sin theta)) (* (cos theta) (sin theta) (expt thetadot 2))) (- (expt (cos theta) 2) 2))))) + ydot) + +(defun rk4-stepper (field dim) + (compile-and-eval + `(let ((k1 (make-real-tensor ,dim)) + (k2 (make-real-tensor ,dim)) + (k3 (make-real-tensor ,dim)) + (k4 (make-real-tensor ,dim)) + (xtmp (make-real-tensor ,dim))) + (lambda (dtm tm x0) + (declare (type double-float dtm tm) + (type real-tensor x0)) + (scal! dtm (funcall ,field tm x0 k1)) + (scal! dtm (funcall ,field (+ tm (/ dtm 2d0)) (axpy! 0.5d0 k1 (copy! x0 xtmp)) k2)) + (scal! dtm (funcall ,field (+ tm (/ dtm 2d0)) (axpy! 0.5d0 k2 (copy! x0 xtmp)) k3)) + (scal! dtm (funcall ,field (+ tm dtm) (axpy! 1d0 k3 (copy! x0 xtmp)) k4)) + (axpy! (/ 1d0 6d0) k1 x0) + (axpy! (/ 1d0 3d0) k2 x0) + (axpy! (/ 1d0 3d0) k3 x0) + (axpy! (/ 1d0 6d0) k4 x0))))) + +;;This is far too verbose. Probably okay for performance intensive stuff, but annoying +;;for quick prototyping ----------------------------------------------------------------------- Summary of changes: tests/usability.lisp | 39 +++++++++++++++++++++++++++++++++++++++ 1 files changed, 39 insertions(+), 0 deletions(-) create mode 100644 tests/usability.lisp hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-11 06:59:38
|
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, tensor has been updated via f01e8cb0e33110ab72e6e828f052e7c534c84ebc (commit) via 18563e60200ee394dc974dd81ec63a48e2b96ea5 (commit) from 2bf1a9674740ca3829a86c45f0dc14a2ee3157c4 (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 f01e8cb0e33110ab72e6e828f052e7c534c84ebc Merge: 18563e6 2bf1a96 Author: Akshay Srinivasan <aks...@gm...> Date: Sun Feb 10 22:58:47 2013 -0800 Merge branch 'tensor' of ssh://matlisp.git.sourceforge.net/gitroot/matlisp/matlisp into tensor commit 18563e60200ee394dc974dd81ec63a48e2b96ea5 Author: Akshay Srinivasan <aks...@gm...> Date: Sun Feb 10 22:58:21 2013 -0800 Brought the &rest argument back to tensor-ref. Changed function names for static methods for tensor classes. diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp index 4e8e80a..8ef5435 100644 --- a/src/base/standard-tensor.lisp +++ b/src/base/standard-tensor.lisp @@ -294,7 +294,7 @@ (setf (number-of-elements tensor) sz))))))))) ;; -(defgeneric tensor-ref (tensor subscripts) +(defgeneric tensor-ref (tensor &rest subscripts) (:documentation " Syntax ====== @@ -312,7 +312,7 @@ of the store.")) -(defgeneric (setf tensor-ref) (value tensor subscripts)) +(defgeneric (setf tensor-ref) (value tensor &rest subscripts)) ;; (defgeneric tensor-store-ref (tensor store-idx) @@ -362,12 +362,12 @@ (defclass ,vector (standard-vector ,tensor-class) ()) ;;Store refs - (defmethod tensor-ref ((tensor ,tensor-class) subs) - (let-typed ((lidx (store-indexing subs tensor) :type index-type) + (defmethod tensor-ref ((tensor ,tensor-class) &rest subs) + (let-typed ((lidx (store-indexing (if (typep (car subs) '(or cons vector)) (car subs) subs) tensor) :type index-type) (sto-x (store tensor) :type ,(linear-array-type store-element-type))) (,reader sto-x lidx))) - (defmethod (setf tensor-ref) (value (tensor ,tensor-class) subs) - (let-typed ((lidx (store-indexing subs tensor) :type index-type) + (defmethod (setf tensor-ref) (value (tensor ,tensor-class) &rest subs) + (let-typed ((lidx (store-indexing (if (typep (car subs) '(or cons vector)) (car subs) subs) tensor) :type index-type) (sto-x (store tensor) :type ,(linear-array-type store-element-type))) (,value-writer (,coercer-unforgiving value) sto-x lidx))) (defmethod tensor-store-ref ((tensor ,tensor-class) lidx) diff --git a/src/classes/complex-tensor.lisp b/src/classes/complex-tensor.lisp index 382641a..df2c0a2 100644 --- a/src/classes/complex-tensor.lisp +++ b/src/classes/complex-tensor.lisp @@ -14,41 +14,41 @@ '(cl:complex complex-base-type)) ;;Field operations -(definline complex-type.f+ (a b) +(definline complex-typed.f+ (a b) (declare (type complex-type a b)) (+ a b)) -(definline complex-type.f- (a b) +(definline complex-typed.f- (a b) (declare (type complex-type a b)) (- a b)) -(definline complex-type.finv+ (a) +(definline complex-typed.finv+ (a) (declare (type complex-type a)) (- a)) -(definline complex-type.fid+ () +(definline complex-typed.fid+ () #c(0.0d0 0.0d0)) -(definline complex-type.f* (a b) +(definline complex-typed.f* (a b) (declare (type complex-type a b)) (* a b)) -(definline complex-type.f/ (a b) +(definline complex-typed.f/ (a b) (declare (type complex-type a b)) (/ a b)) -(definline complex-type.finv* (a) +(definline complex-typed.finv* (a) (declare (type complex-type a)) (/ a)) -(definline complex-type.fid* () +(definline complex-typed.fid* () #c(1.0d0 0.0d0)) -(definline complex-type.fconj (a) +(definline complex-typed.fconj (a) (declare (type complex-type a)) (conjugate a)) -(definline complex-type.f= (a b) +(definline complex-typed.f= (a b) (declare (type complex-type a b)) (= a b)) @@ -76,33 +76,33 @@ (restart-case (coerce-complex-base-unforgiving x) (use-value (value) (coerce-complex-base value)))) ;; -(definline complex-type.reader (tstore idx) +(definline complex-typed.reader (tstore idx) (declare (type complex-store-vector tstore) (type index-type idx)) (complex (aref tstore (* 2 idx)) (aref tstore (1+ (* 2 idx))))) -(definline complex-type.value-writer (value store idx) +(definline complex-typed.value-writer (value store idx) (declare (type complex-store-vector store) (type index-type idx) (type complex-type value)) (setf (aref store (* 2 idx)) (realpart value) (aref store (1+ (* 2 idx))) (imagpart value))) -(definline complex-type.value-incfer (value store idx) +(definline complex-typed.value-incfer (value store idx) (declare (type complex-store-vector store) (type index-type idx) (type complex-type value)) (incf (aref store (* 2 idx)) (realpart value)) (incf (aref store (1+ (* 2 idx))) (imagpart value))) -(definline complex-type.reader-writer (fstore fidx tstore tidx) +(definline complex-typed.reader-writer (fstore fidx tstore tidx) (declare (type complex-store-vector fstore tstore) (type index-type fidx tidx)) (setf (aref tstore (* 2 tidx)) (aref fstore (* 2 fidx)) (aref tstore (1+ (* 2 tidx))) (aref fstore (1+ (* 2 fidx))))) -(definline complex-type.swapper (fstore fidx tstore tidx) +(definline complex-typed.swapper (fstore fidx tstore tidx) (declare (type complex-store-vector fstore tstore) (type index-type fidx tidx)) (rotatef (aref tstore (* 2 tidx)) (aref fstore (* 2 fidx))) @@ -113,16 +113,16 @@ (:documentation "Tensor class with complex elements.")) :matrix complex-matrix :vector complex-vector ;; - :f+ complex-type.f+ - :f- complex-type.f- - :finv+ complex-type.finv+ - :fid+ complex-type.fid+ - :f* complex-type.f* - :f/ complex-type.f/ - :finv* complex-type.finv* - :fid* complex-type.fid* - :f= complex-type.f= - :fconj complex-type.fconj + :f+ complex-typed.f+ + :f- complex-typed.f- + :finv+ complex-typed.finv+ + :fid+ complex-typed.fid+ + :f* complex-typed.f* + :f/ complex-typed.f/ + :finv* complex-typed.finv* + :fid* complex-typed.fid* + :f= complex-typed.f= + :fconj complex-typed.fconj ;; :store-allocator allocate-complex-store :coercer coerce-complex @@ -130,11 +130,11 @@ ;; :matrix complex-matrix :vector complex-vector ;; - :reader complex-type.reader - :value-writer complex-type.value-writer - :value-incfer complex-type.value-incfer - :reader-writer complex-type.reader-writer - :swapper complex-type.swapper) + :reader complex-typed.reader + :value-writer complex-typed.value-writer + :value-incfer complex-typed.value-incfer + :reader-writer complex-typed.reader-writer + :swapper complex-typed.swapper) ;; #+nil diff --git a/src/classes/real-tensor.lisp b/src/classes/real-tensor.lisp index ac5d144..9224e61 100644 --- a/src/classes/real-tensor.lisp +++ b/src/classes/real-tensor.lisp @@ -9,64 +9,64 @@ `(simple-array real-type (,size))) ;;Field definitions -(definline real-type.f+ (a b) +(definline real-typed.f+ (a b) (declare (type real-type a b)) (+ a b)) -(definline real-type.f- (a b) +(definline real-typed.f- (a b) (declare (type real-type a b)) (- a b)) -(definline real-type.finv+ (a) +(definline real-typed.finv+ (a) (declare (type real-type a)) (- a)) -(definline real-type.fid+ () +(definline real-typed.fid+ () 0.0d0) -(definline real-type.f* (a b) +(definline real-typed.f* (a b) (declare (type real-type a b)) (* a b)) -(definline real-type.f/ (a b) +(definline real-typed.f/ (a b) (declare (type real-type a b)) (/ a b)) -(definline real-type.finv* (a) +(definline real-typed.finv* (a) (declare (type real-type a)) (/ a)) -(definline real-type.fid* () +(definline real-typed.fid* () 1.0d0) -(definline real-type.fid= (a b) +(definline real-typed.fid= (a b) (declare (type real-type a b)) (= a b)) ;;Store definitions -(definline real-type.reader (tstore idx) +(definline real-typed.reader (tstore idx) (declare (type index-type idx) (type real-store-vector tstore)) (aref tstore idx)) -(definline real-type.value-writer (value store idx) +(definline real-typed.value-writer (value store idx) (declare (type index-type idx) (type real-store-vector store) (type real-type value)) (setf (aref store idx) value)) -(definline real-type.value-incfer (value store idx) +(definline real-typed.value-incfer (value store idx) (declare (type index-type idx) (type real-store-vector store) (type real-type value)) (incf (aref store idx) value)) -(definline real-type.reader-writer (fstore fidx tstore tidx) +(definline real-typed.reader-writer (fstore fidx tstore tidx) (declare (type index-type fidx tidx) (type real-store-vector fstore tstore)) (setf (aref tstore tidx) (aref fstore fidx))) -(definline real-type.swapper (fstore fidx tstore tidx) +(definline real-typed.swapper (fstore fidx tstore tidx) (declare (type index-type fidx tidx) (type real-store-vector fstore tstore)) (rotatef (aref tstore tidx) (aref fstore fidx))) @@ -88,26 +88,26 @@ Allocates real storage. Default initial-element = 0d0.") (:documentation "Tensor class with real double elements.")) :matrix real-matrix :vector real-vector ;; - :f+ real-type.f+ - :f- real-type.f- - :finv+ real-type.finv+ - :fid+ real-type.fid+ - :f* real-type.f* - :f/ real-type.f/ - :finv* real-type.finv* - :fid* real-type.fid* - :f= real-type.fid= + :f+ real-typed.f+ + :f- real-typed.f- + :finv+ real-typed.finv+ + :fid+ real-typed.fid+ + :f* real-typed.f* + :f/ real-typed.f/ + :finv* real-typed.finv* + :fid* real-typed.fid* + :f= real-typed.fid= :fconj nil ;; :store-allocator allocate-real-store :coercer coerce-real :coercer-unforgiving coerce-real-unforgiving ;; - :reader real-type.reader - :value-writer real-type.value-writer - :value-incfer real-type.value-incfer - :reader-writer real-type.reader-writer - :swapper real-type.swapper) + :reader real-typed.reader + :value-writer real-typed.value-writer + :value-incfer real-typed.value-incfer + :reader-writer real-typed.reader-writer + :swapper real-typed.swapper) (defmethod print-element ((tensor real-tensor) element stream) diff --git a/src/classes/symbolic-tensor.lisp b/src/classes/symbolic-tensor.lisp index 0365929..4e2c7e2 100644 --- a/src/classes/symbolic-tensor.lisp +++ b/src/classes/symbolic-tensor.lisp @@ -9,44 +9,44 @@ `(simple-array symbolic-type (,size))) ;;Field definitions -(definline symbolic-type.f+ (a b) +(definline symbolic-typed.f+ (a b) (declare (type symbolic-type a b)) (maxima::add a b)) -(definline symbolic-type.f- (a b) +(definline symbolic-typed.f- (a b) (declare (type symbolic-type a b)) (maxima::sub a b)) -(definline symbolic-type.finv+ (a) +(definline symbolic-typed.finv+ (a) (declare (type symbolic-type a)) (maxima::mul -1 a)) -(definline symbolic-type.fid+ () +(definline symbolic-typed.fid+ () 0) -(definline symbolic-type.f* (a b) +(definline symbolic-typed.f* (a b) (declare (type symbolic-type a b)) (maxima::mul a b)) -(definline symbolic-type.f/ (a b) +(definline symbolic-typed.f/ (a b) (declare (type symbolic-type a b)) (maxima::div a b)) -(definline symbolic-type.finv* (a) +(definline symbolic-typed.finv* (a) (declare (type symbolic-type a)) (maxima::div 1 a)) -(definline symbolic-type.fid* () +(definline symbolic-typed.fid* () 1) -(definline symbolic-type.f= (a b) +(definline symbolic-typed.f= (a b) (declare (type symbolic-type a b)) (maxima::equal a b)) -(definline symbolic-type.fconj (a) +(definline symbolic-typed.fconj (a) (maxima::meval `((maxima::$conjugate maxima::simp) ,a))) -(definline symbolic-type.diff (a x) +(definline symbolic-typed.diff (a x) (etypecase a (symbolic-type (maxima::$diff a x)) @@ -56,29 +56,29 @@ :store (map 'symbolic-store-vector #'(lambda (f) (maxima::$diff f x)) (store a)))))) ;; ;;Store definitions -(definline symbolic-type.reader (tstore idx) +(definline symbolic-typed.reader (tstore idx) (declare (type index-type idx) (type symbolic-store-vector tstore)) (aref tstore idx)) -(definline symbolic-type.value-writer (value store idx) +(definline symbolic-typed.value-writer (value store idx) (declare (type index-type idx) (type symbolic-store-vector store) (type symbolic-type value)) (setf (aref store idx) value)) -(definline symbolic-type.value-incfer (value store idx) +(definline symbolic-typed.value-incfer (value store idx) (declare (type index-type idx) (type symbolic-store-vector store) (type symbolic-type value)) - (setf (aref store idx) (symbolic-type.f+ (aref store idx) value))) + (setf (aref store idx) (symbolic-typed.f+ (aref store idx) value))) -(definline symbolic-type.reader-writer (fstore fidx tstore tidx) +(definline symbolic-typed.reader-writer (fstore fidx tstore tidx) (declare (type index-type fidx tidx) (type symbolic-store-vector fstore tstore)) (setf (aref tstore tidx) (aref fstore fidx))) -(definline symbolic-type.swapper (fstore fidx tstore tidx) +(definline symbolic-typed.swapper (fstore fidx tstore tidx) (declare (type index-type fidx tidx) (type symbolic-store-vector fstore tstore)) (rotatef (aref tstore tidx) (aref fstore fidx))) @@ -99,26 +99,26 @@ Allocates symbolic storage. Default initial-element = 0.") (:documentation "Tensor class with symbolic double elements.")) :matrix symbolic-matrix :vector symbolic-vector ;; - :f+ symbolic-type.f+ - :f- symbolic-type.f- - :finv+ symbolic-type.finv+ - :fid+ symbolic-type.fid+ - :f* symbolic-type.f* - :f/ symbolic-type.f/ - :finv* symbolic-type.finv* - :fid* symbolic-type.fid* - :f= symbolic-type.f= - :fconj symbolic-type.fconj + :f+ symbolic-typed.f+ + :f- symbolic-typed.f- + :finv+ symbolic-typed.finv+ + :fid+ symbolic-typed.fid+ + :f* symbolic-typed.f* + :f/ symbolic-typed.f/ + :finv* symbolic-typed.finv* + :fid* symbolic-typed.fid* + :f= symbolic-typed.f= + :fconj symbolic-typed.fconj ;; :store-allocator allocate-symbolic-store :coercer coerce-symbolic :coercer-unforgiving coerce-symbolic-unforgiving ;; - :reader symbolic-type.reader - :value-writer symbolic-type.value-writer - :value-incfer symbolic-type.value-incfer - :reader-writer symbolic-type.reader-writer - :swapper symbolic-type.swapper) + :reader symbolic-typed.reader + :value-writer symbolic-typed.value-writer + :value-incfer symbolic-typed.value-incfer + :reader-writer symbolic-typed.reader-writer + :swapper symbolic-typed.swapper) (defmethod initialize-instance ((tensor symbolic-tensor) &rest initargs) (if (getf initargs :store) diff --git a/src/ffi/foreign-vector.lisp b/src/ffi/foreign-vector.lisp index f70c459..1521986 100644 --- a/src/ffi/foreign-vector.lisp +++ b/src/ffi/foreign-vector.lisp @@ -18,24 +18,14 @@ (defun fv-ref (x n) (declare (type foreign-vector x) (type fixnum n)) - (let ((sap (fv-pointer x)) - (ss (fv-size x)) - (sty (fv-type x))) - (unless (< -1 n ss) - (error 'out-of-bounds-error :requested n :bound ss - :message "From inside fv-ref.")) - (cffi:mem-aref sap sty n))) + (assert (< -1 n (fv-size x)) nil 'out-of-bounds-error :requested n :bound (fv-size x) :message "From inside fv-ref.") + (cffi:mem-aref (fv-pointer x) (fv-type x) n)) (defun (setf fv-ref) (value x n) (declare (type foreign-vector x) (type fixnum n)) - (let ((sap (fv-pointer x)) - (ss (fv-size x)) - (sty (fv-type x))) - (unless (< -1 n ss) - (error 'out-of-bounds-error :requested n :bound ss - :message "From inside (setf fv-ref).")) - (setf (cffi:mem-aref sap sty n) value))) + (assert (< -1 n (fv-size x)) nil 'out-of-bounds-error :requested n :bound (fv-size x) :message "From inside fv-ref.") + (setf (cffi:mem-aref (fv-pointer x) (fv-type x) n) value)) ;;; Rudimentary support for making it a bit easier to deal with Fortran ;;; arrays in callbacks. ----------------------------------------------------------------------- Summary of changes: src/base/standard-tensor.lisp | 12 +++--- src/classes/complex-tensor.lisp | 60 ++++++++++++++++++------------------ src/classes/real-tensor.lisp | 56 ++++++++++++++++---------------- src/classes/symbolic-tensor.lisp | 64 +++++++++++++++++++------------------- src/ffi/foreign-vector.lisp | 18 ++-------- 5 files changed, 100 insertions(+), 110 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-05 06:41:57
|
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, tensor has been updated via 2bf1a9674740ca3829a86c45f0dc14a2ee3157c4 (commit) via 9253264f16cb969a88e0411ba4e8441a18406395 (commit) from b8b3f1fdeb9425c9e36e096428fc04e044f720f6 (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 2bf1a9674740ca3829a86c45f0dc14a2ee3157c4 Merge: 9253264 b8b3f1f Author: Akshay Srinivasan <aks...@gm...> Date: Mon Feb 4 22:35:47 2013 -0800 Merge branch 'tensor' of ssh://matlisp.git.sourceforge.net/gitroot/matlisp/matlisp into tensor commit 9253264f16cb969a88e0411ba4e8441a18406395 Author: Akshay Srinivasan <aks...@gm...> Date: Mon Feb 4 22:35:04 2013 -0800 Used unwind-protect instead of multiple-value-prog1 inside the with-foreign-objects-heaped macro. diff --git a/src/ffi/ffi-cffi.lisp b/src/ffi/ffi-cffi.lisp index 5698605..8ac6db3 100644 --- a/src/ffi/ffi-cffi.lisp +++ b/src/ffi/ffi-cffi.lisp @@ -59,7 +59,7 @@ declarations)))) ;; Store result and then free foreign-objects (when declarations - `(multiple-value-prog1)) + `(unwind-protect)) `((progn ,@body) ;;Free C objects ----------------------------------------------------------------------- Summary of changes: src/ffi/ffi-cffi.lisp | 2 +- 1 files changed, 1 insertions(+), 1 deletions(-) hooks/post-receive -- matlisp |
From: Raymond T. <rt...@us...> - 2013-02-05 03:50:23
|
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, tensor has been updated via b8b3f1fdeb9425c9e36e096428fc04e044f720f6 (commit) from 8c2e30aafd1d5481e28d62cbf5382687fe0b4fe7 (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 b8b3f1fdeb9425c9e36e096428fc04e044f720f6 Author: Raymond Toy <toy...@gm...> Date: Mon Feb 4 19:50:01 2013 -0800 Get rid of bad style warnings for ecase by converting them to case. diff --git a/src/ffi/f77-ffi.lisp b/src/ffi/f77-ffi.lisp index 31031c2..3fad009 100644 --- a/src/ffi/f77-ffi.lisp +++ b/src/ffi/f77-ffi.lisp @@ -48,7 +48,7 @@ `(:pointer ,(%f77.cffi-type :single-float))) ((eq type :complex-double-float) `(:pointer ,(%f77.cffi-type :double-float))) - (t (ecase type + (t (case type (:void :void) (:integer :int32) (:character :char) diff --git a/src/level-1/realimag.lisp b/src/level-1/realimag.lisp index db9566a..c09d67a 100644 --- a/src/level-1/realimag.lisp +++ b/src/level-1/realimag.lisp @@ -43,7 +43,10 @@ " (etypecase tensor (real-tensor tensor) - (complex-tensor (make-instance (ecase (rank tensor) (2 'real-matrix) (1 'real-vector) (t 'real-tensor)) + (complex-tensor (make-instance (case (rank tensor) + (2 'real-matrix) + (1 'real-vector) + (t 'real-tensor)) :parent-tensor tensor :store (store tensor) :store-size (length (store tensor)) :dimensions (dimensions tensor) :strides (map 'index-store-vector #'(lambda (x) (* 2 x)) (strides tensor)) @@ -65,7 +68,10 @@ " (etypecase tensor (real-tensor tensor) - (complex-tensor (make-instance (ecase (rank tensor) (2 'real-matrix) (1 'real-vector) (t 'real-tensor)) + (complex-tensor (make-instance (case (rank tensor) + (2 'real-matrix) + (1 'real-vector) + (t 'real-tensor)) :parent-tensor tensor :store (store tensor) :store-size (length (store tensor)) :dimensions (dimensions tensor) :strides (map 'index-store-vector #'(lambda (x) (* 2 x)) (strides tensor)) ----------------------------------------------------------------------- Summary of changes: src/ffi/f77-ffi.lisp | 2 +- src/level-1/realimag.lisp | 10 ++++++++-- 2 files changed, 9 insertions(+), 3 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-04 06:38:51
|
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, tensor has been updated via 8c2e30aafd1d5481e28d62cbf5382687fe0b4fe7 (commit) from bff689247f6e345a4ff7505a2973d70a3e30acea (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 8c2e30aafd1d5481e28d62cbf5382687fe0b4fe7 Author: Akshay Srinivasan <aks...@gm...> Date: Sun Feb 3 22:31:44 2013 -0800 Removed type declarations in store-indexing. CMUCL does not seem to like them :) diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp index 815164d..81c33c6 100644 --- a/src/base/standard-tensor.lisp +++ b/src/base/standard-tensor.lisp @@ -256,9 +256,7 @@ /_ i i i = 0 " - (declare (type standard-tensor tensor) - (type (or index-store-vector cons) idx)) - (typecase idx + (etypecase idx (cons (store-indexing-lst idx (head tensor) (strides tensor) (dimensions tensor))) (vector (store-indexing-vec idx (head tensor) (strides tensor) (dimensions tensor))))) ;; ----------------------------------------------------------------------- Summary of changes: src/base/standard-tensor.lisp | 4 +--- 1 files changed, 1 insertions(+), 3 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-04 06:01:19
|
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, tensor has been updated via bff689247f6e345a4ff7505a2973d70a3e30acea (commit) from d5850b327492c870f2fba052f997dd9cea2a488f (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 bff689247f6e345a4ff7505a2973d70a3e30acea Author: Akshay Srinivasan <aks...@gm...> Date: Sun Feb 3 21:56:26 2013 -0800 Moved method definitions out of the eval-when macro inside define-tensor. diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp index 4e8e80a..815164d 100644 --- a/src/base/standard-tensor.lisp +++ b/src/base/standard-tensor.lisp @@ -352,7 +352,7 @@ ;;Error checking (assert (and f+ f- finv+ fid+ f* f/ finv* fid* f= store-allocator coercer coercer-unforgiving matrix vector reader value-writer value-incfer reader-writer swapper)) ;; - `(eval-when (:compile-toplevel :load-toplevel :execute) + `(progn ;;Class definitions (defclass ,tensor-class (standard-tensor) ((store :type ,store-type)) @@ -379,34 +379,35 @@ (let-typed ((sto-x (store tensor) :type ,(linear-array-type store-element-type))) (,value-writer (,coercer-unforgiving value) sto-x lidx))) ;; - (let ((hst (list - :tensor ',tensor-class - :matrix ',matrix - :vector ',vector - :element-type ',element-type - :f+ ',f+ - :f- ',f- - :finv+ ',finv+ - :fid+ ',fid+ - :f* ',f* - :f/ ',f/ - :finv* ',finv* - :fid* ',fid* - :f= ',f= - :fconj ',fconj - :reader ',reader - :value-writer ',value-writer - :value-incfer ',value-incfer - :reader-writer ',reader-writer - :swapper ',swapper - :store-allocator ',store-allocator - :coercer ',coercer - :coercer-unforgiving ',coercer-unforgiving - :store-type ',store-element-type))) - (setf (get-tensor-class-optimization ',tensor-class) hst - (get-tensor-class-optimization ',matrix) ',tensor-class - (get-tensor-class-optimization ',vector) ',tensor-class) - (setf (symbol-plist ',tensor-class) hst)))) + (eval-when (:compile-toplevel :load-toplevel :execute) + (let ((hst (list + :tensor ',tensor-class + :matrix ',matrix + :vector ',vector + :element-type ',element-type + :f+ ',f+ + :f- ',f- + :finv+ ',finv+ + :fid+ ',fid+ + :f* ',f* + :f/ ',f/ + :finv* ',finv* + :fid* ',fid* + :f= ',f= + :fconj ',fconj + :reader ',reader + :value-writer ',value-writer + :value-incfer ',value-incfer + :reader-writer ',reader-writer + :swapper ',swapper + :store-allocator ',store-allocator + :coercer ',coercer + :coercer-unforgiving ',coercer-unforgiving + :store-type ',store-element-type))) + (setf (get-tensor-class-optimization ',tensor-class) hst + (get-tensor-class-optimization ',matrix) ',tensor-class + (get-tensor-class-optimization ',vector) ',tensor-class) + (setf (symbol-plist ',tensor-class) hst))))) ;; (defun tensor-typep (tensor subscripts) ----------------------------------------------------------------------- Summary of changes: src/base/standard-tensor.lisp | 59 +++++++++++++++++++++-------------------- 1 files changed, 30 insertions(+), 29 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-04 03:55:56
|
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, tensor has been updated via d5850b327492c870f2fba052f997dd9cea2a488f (commit) from 926c7dc0991dbe277f1802ab9074b89599cc5008 (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 d5850b327492c870f2fba052f997dd9cea2a488f Author: Akshay Srinivasan <aks...@gm...> Date: Sun Feb 3 19:55:34 2013 -0800 Saving changes to getrs.lisp. diff --git a/src/foreign-core/lapack.lisp b/src/foreign-core/lapack.lisp index ad75f1b..d368258 100644 --- a/src/foreign-core/lapack.lisp +++ b/src/foreign-core/lapack.lisp @@ -307,7 +307,7 @@ The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0. - A (input/output) DOUBLE PRECISION array, dimension (LDA,N) + A (input) DOUBLE PRECISION array, dimension (LDA,N) The factors L and U from the factorization A = P*L*U as computed by DGETRF. diff --git a/src/lapack/getrf.lisp b/src/lapack/getrf.lisp index 7e32538..df4a244 100644 --- a/src/lapack/getrf.lisp +++ b/src/lapack/getrf.lisp @@ -43,17 +43,23 @@ (type permutation-pivot-flip ipiv)) (mlet* (((maj-A ld-A fop-A) (blas-matrix-compatible-p A :n) :type (symbol index-type nil))) (if (eq maj-A :col-major) - (multiple-value-bind (n-A n-ipiv info) (,lapack-func - (nrows A) (ncols A) (store A) - ld-A (store ipiv) 0) - (declare (ignore n-A n-ipiv)) - (assert (= info 0) nil 'invalid-arguments :argnum (1- (- info)) :message (format-to-string "GETRF returned INFO: ~a." info)) - ;;Convert back to 0-based indexing. + (progn + ;;Convert to 1-based indexing. (let-typed ((pidv (store ipiv) :type pindex-store-vector)) (very-quickly (loop :for i :of-type index-type :from 0 :below (length pidv) - :do (decf (aref pidv i))))) - (values A ipiv info)) + :do (incf (aref pidv i))))) + (multiple-value-bind (n-A n-ipiv info) (,lapack-func + (nrows A) (ncols A) (store A) + ld-A (store ipiv) 0) + (declare (ignore n-A n-ipiv)) + (assert (= info 0) nil 'invalid-arguments :argnum (1- (- info)) :message (format-to-string "GETRF returned INFO: ~a." info)) + ;;Convert back to 0-based indexing. + (let-typed ((pidv (store ipiv) :type pindex-store-vector)) + (very-quickly + (loop :for i :of-type index-type :from 0 :below (length pidv) + :do (decf (aref pidv i))))) + (values A ipiv info))) (let ((tmp (,(getf opt :copy) A (with-order :col-major (,(getf opt :zero-maker) (dimensions A)))))) (multiple-value-bind (n-tmp n-ipiv info) (,func-name tmp ipiv) (declare (ignore n-tmp n-ipiv)) diff --git a/src/lapack/getrs.lisp b/src/lapack/getrs.lisp index d25f6f1..45a2bbe 100644 --- a/src/lapack/getrs.lisp +++ b/src/lapack/getrs.lisp @@ -28,7 +28,73 @@ (in-package #:matlisp) -(defgeneric getrs! (A B &optional job-a) +(defmacro generate-typed-getrs! (func-name (tensor-class lapack-func)) + (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 :getrs) ',func-name + (get-tensor-class-optimization ',tensor-class) opt))) + (defun ,func-name (A B ipiv job-A) + (declare (type ,matrix-class A) + (type (or ,matrix-class ,vector-class) B) + (type permutation-pivot-flip ipiv)) + (mlet* (((maj-A ld-A fop-A) (blas-matrix-compatible-p A job-A) :type (symbol index-type string)) + ((maj-B ld-B fop-B) (etypecase B + (,matrix-class (blas-matrix-compatible-p B :n)) + (,vector-class (if (= (aref (strides B) 0) 1) + (values :col-major (aref (dimensions B) 0) nil) + (values nil 0 nil)))) + :type (symbol index-type nil))) + (if (and (eq maj-B :col-major) (eq maj-A :col-major)) + (progn + ;;Convert to 1-based indexing. + (let-typed ((pidv (store ipiv) :type pindex-store-vector)) + (very-quickly + (loop :for i :of-type index-type :from 0 :below (length pidv) + :do (incf (aref pidv i))))) + (multiple-value-bind (n-B info) (,lapack-func + fop-A (nrows A) (etypecase B (,matrix-class (ncols B)) (,vector-class 1)) + (store A) ld-A + (store ipiv) + (store B) ld-b + 0) + (declare (ignore n-B)) + (assert (= info 0) nil 'invalid-arguments :argnum (1- (- info)) :message (format-to-string "GETRF returned INFO: ~a." info)) + ;;Convert back to 0-based indexing. + (let-typed ((pidv (store ipiv) :type pindex-store-vector)) + (very-quickly + (loop :for i :of-type index-type :from 0 :below (length pidv) + :do (decf (aref pidv i))))) + (values B info))) + (let ((n-A (if (eq maj-A :col-major) A (,(getf opt :copy) A (with-order :col-major (,(getf opt :zero-maker) (dimensions A)))))) + (n-B (if (eq maj-B :col-major) B (,(getf opt :copy) B (with-order :col-major (,(getf opt :zero-maker) (dimensions B))))))) + (multiple-value-bind (B.ret info) (,func-name n-A n-B ipiv job-A) + (declare (ignore B.ret)) + (unless (eq maj-B :col-major) + (,(getf opt :copy) n-B B)) + (values B info))))))))) + +(generate-typed-getrs! real-base-typed-getrs! (real-tensor dgetrs)) +(definline real-typed-getrs! (A B ipiv job-A) + (real-base-typed-getrs! A B ipiv (ecase job-A (:c :n) (:h :t) ((:n :t) job-A)))) + +(generate-typed-getrs! complex-base-typed-getrs! (complex-tensor zgetrs)) +(definline complex-typed-getrs! (A B ipiv job-A) + (let ((A (ecase job-A + ((:n :t) A) + ((:h :c) (let ((ret (complex-typed-copy! A (complex-typed-zeros (dimensions A))))) + (real-typed-num-scal! -1d0 (tensor-imagpart~ ret)) + ret))))) + (complex-base-typed-getrs! A B ipiv job-A))) + + +;; +(defgeneric getrs! (ipiv A B &optional job-a) (:documentation " Syntax @@ -55,87 +121,24 @@ but the factor U is exactly singular. Solution could not be computed. ") - (:method :before ((A standard-matrix) (B standard-matrix) &optional job-a) + (:method :before ((ipiv permutation) (A standard-matrix) (B standard-matrix) &optional job-a) (declare (ignore job-a)) - (assert (= (nrows A) (ncols A) - (nrows B)) nil 'tensor-dimension-mismatch) - ;;Well, we can't re-implement every algorithm in LAPACK to work - ;;with two strided matrices. Threw in the towel after BLAS. - (assert (and (consecutive-store-p A) - (consecutive-store-p B)) - nil 'tensor-store-not-consecutive))) - -(defmacro generate-typed-getrs! (func-name (matrix-class lapack-func)) - (let* ((opt (get-tensor-class-optimization matrix-class))) - (assert opt nil 'tensor-cannot-find-optimization :tensor-class matrix-class) - `(defun ,func-name (A B job-A) - (declare (type ,matrix-class A B) - (type symbol job-A)) - (mlet* (((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 :n) :type (symbol index-type (string 1)))) - (if (eq (maj-B - - (cond - ((eq maj-A :col-major) - -(defmethod getrs! ((A real-matrix) (B real-matrix) &optional (job-A :n)) - (mlet* (((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 :n) :type (symbol index-type (string 1)))) - - (let* ((n (nrows a)) - (m (ncols b))) - - (declare (type fixnum n m) - (type (simple-array (unsigned-byte 32) (*)) ipiv)) - - (multiple-value-bind (x info) - (dgetrs (case trans - (:C "C") - (:T "T") - (t "N")) - n - m - (store a) - n - ipiv - (store b) - n - 0) - - (values - (make-instance 'real-matrix :nrows n :ncols m :store x) - (if (zerop info) - t - info))))) - -(defmethod getrs! ((a complex-matrix) ipiv (b complex-matrix) &key trans) - - (let* ((n (nrows a)) - (m (ncols b))) - - (declare (type fixnum n m) - (type (simple-array (unsigned-byte 32) (*)) ipiv)) - - (multiple-value-bind (x info) - (zgetrs (case trans - (:C "C") - (:T "T") - (t "N")) - n - m - (store a) - n - ipiv - (store b) - n - 0) - - (values - (make-instance 'complex-matrix :nrows n :ncols m :store x) - (if (zerop info) - t - info))))) + (assert (and (= (nrows A) (ncols A) (nrows B)) (>= (permutation-size ipiv) (nrows A))) + nil 'tensor-dimension-mismatch))) +(defmethod getrs! ((ipiv permutation-pivot-flip) (A real-matrix) (B real-matrix) &optional (job-A :n)) + (real-typed-getrs! A B ipiv job-A)) + +(defmethod getrs! ((ipiv permutation-pivot-flip) (A real-matrix) (B real-vector) &optional (job-A :n)) + (real-typed-getrs! A B ipiv job-A)) + +(defmethod getrs! ((ipiv permutation-pivot-flip) (A complex-matrix) (B complex-matrix) &optional (job-A :n)) + (complex-typed-getrs! A B ipiv job-A)) + +(defmethod getrs! ((ipiv permutation-pivot-flip) (A complex-matrix) (B complex-vector) &optional (job-A :n)) + (complex-typed-getrs! A B ipiv job-A)) + +;; (defmethod getrs! ((a standard-matrix) ipiv (b standard-matrix) &key trans) (let ((a (typecase a (real-matrix (copy! a (make-complex-matrix-dim (nrows a) (ncols a)))) ----------------------------------------------------------------------- Summary of changes: src/foreign-core/lapack.lisp | 2 +- src/lapack/getrf.lisp | 22 ++++-- src/lapack/getrs.lisp | 155 +++++++++++++++++++++-------------------- 3 files changed, 94 insertions(+), 85 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-03 22:29:04
|
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, tensor has been updated via 926c7dc0991dbe277f1802ab9074b89599cc5008 (commit) from 2f7115f30d58288fdc9b160ce621ab19d1277772 (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 926c7dc0991dbe277f1802ab9074b89599cc5008 Author: Akshay Srinivasan <aks...@gm...> Date: Sun Feb 3 14:28:44 2013 -0800 Moved function definitions out of (eval-when ..). diff --git a/src/ffi/f77-ffi.lisp b/src/ffi/f77-ffi.lisp index b03b2c5..31031c2 100644 --- a/src/ffi/f77-ffi.lisp +++ b/src/ffi/f77-ffi.lisp @@ -546,7 +546,7 @@ ,@pars)) (setq hack-return-type :void))) - `(eval-when (:compile-toplevel :load-toplevel :execute) + `(progn (cffi:defcfun (,fortran-name ,lisp-name) ,(%f77.get-return-type hack-return-type) ,@(%f77.parse-fortran-parameters hack-body)) ,@(%f77.def-fortran-interface name hack-return-type hack-body hidden-var-name))))) diff --git a/src/lapack/getrf.lisp b/src/lapack/getrf.lisp index d5a17ae..7e32538 100644 --- a/src/lapack/getrf.lisp +++ b/src/lapack/getrf.lisp @@ -32,11 +32,12 @@ (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))) - `(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 :getrf) ',func-name - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :getrf) ',func-name + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func-name (A ipiv) (declare (type ,matrix-class A) (type permutation-pivot-flip ipiv)) @@ -136,13 +137,12 @@ (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))) - `(eval-when (:compile-toplevel :load-toplevel :execute) - (defmethod lu ((A ,matrix-class) &optional (split-lu? t)) - (multiple-value-bind (lu ipiv info) - (getrf! (with-order :col-major - (,(getf opt :copy) A (,(getf opt :zero-maker) (dimensions A))))) - (declare (ignore info)) - (let* ((n (nrows a)) + `(defmethod lu ((A ,matrix-class) &optional (split-lu? t)) + (multiple-value-bind (lu ipiv info) + (getrf! (with-order :col-major + (,(getf opt :copy) A (,(getf opt :zero-maker) (dimensions A))))) + (declare (ignore info)) + (let* ((n (nrows a)) (m (ncols a)) (p (min n m))) (declare (type fixnum n m p)) @@ -197,7 +197,7 @@ (incf lu.of (- lu.cstd (the index-type (* (- n j 2) lu.rstd)))) (incf l.of (- l.cstd (the index-type (* (- n j 2) l.rstd))))))) (values lmat lu ipiv))) - (values lu ipiv)))))))) + (values lu ipiv))))))) (make-lu real-tensor) (make-lu complex-tensor) diff --git a/src/level-1/axpy.lisp b/src/level-1/axpy.lisp index d161d97..25966de 100644 --- a/src/level-1/axpy.lisp +++ b/src/level-1/axpy.lisp @@ -34,11 +34,12 @@ ;;Use only after checking the arguments for compatibility. (let* ((opt (get-tensor-class-optimization-hashtable tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) - `(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 :axpy) ',func - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :axpy) ',func + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func (alpha from to) (declare (type ,tensor-class from to) (type ,(getf opt :element-type) alpha)) @@ -80,11 +81,12 @@ ;;Use only after checking the arguments for compatibility. (let* ((opt (get-tensor-class-optimization tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) - `(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 :num-axpy) ',func - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :num-axpy) ',func + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func (num-from to) (declare (type ,tensor-class to) (type ,(getf opt :element-type) num-from)) diff --git a/src/level-1/copy.lisp b/src/level-1/copy.lisp index 6680e8a..880ccf5 100644 --- a/src/level-1/copy.lisp +++ b/src/level-1/copy.lisp @@ -34,11 +34,12 @@ ;;Use only after checking the arguments for compatibility. (let* ((opt (get-tensor-class-optimization-hashtable tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) - `(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 :copy) ',func - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :copy) ',func + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func (from to) (declare (type ,tensor-class from to)) ,(let @@ -76,11 +77,12 @@ ;;Use only after checking the arguments for compatibility. (let* ((opt (get-tensor-class-optimization-hashtable tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) - `(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 :num-copy) ',func - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :num-copy) ',func + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func (num-from to) (declare (type ,tensor-class to) (type ,(getf opt :element-type) num-from)) diff --git a/src/level-1/dot.lisp b/src/level-1/dot.lisp index 435cff1..8e1fe89 100644 --- a/src/level-1/dot.lisp +++ b/src/level-1/dot.lisp @@ -34,11 +34,12 @@ (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) (setf (getf opt :dot) func (get-tensor-class-optimization tensor-class) opt) - `(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 :axpy) ',func - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :axpy) ',func + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func (x y conjugate-p) (declare (type ,tensor-class x y) ,(if conj? diff --git a/src/level-1/scal.lisp b/src/level-1/scal.lisp index ad2749d..72a7962 100644 --- a/src/level-1/scal.lisp +++ b/src/level-1/scal.lisp @@ -31,11 +31,12 @@ (defmacro generate-typed-scal! (func (tensor-class fortran-func fortran-lb)) (let* ((opt (get-tensor-class-optimization-hashtable tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) - `(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 :scal) ',func - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :scal) ',func + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func (from to) (declare (type ,tensor-class from to)) ,(let @@ -69,11 +70,12 @@ (defmacro generate-typed-num-scal! (func (tensor-class blas-func fortran-lb)) (let ((opt (get-tensor-class-optimization-hashtable tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) - `(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 :num-scal) ',func - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :num-scal) ',func + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func (alpha to) (declare (type ,tensor-class to) (type ,(getf opt :element-type) alpha)) @@ -101,11 +103,12 @@ (defmacro generate-typed-div! (func (tensor-class fortran-func fortran-lb)) (let* ((opt (get-tensor-class-optimization-hashtable tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) - `(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 :div) ',func - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :div) ',func + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func (from to) (declare (type ,tensor-class from to)) ,(let @@ -139,11 +142,12 @@ (defmacro generate-typed-num-div! (func (tensor-class fortran-func fortran-lb)) (let ((opt (get-tensor-class-optimization tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) - `(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 :num-div) ',func - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :num-div) ',func + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func (alpha to) (declare (type ,tensor-class to) (type ,(getf opt :element-type) alpha)) diff --git a/src/level-1/swap.lisp b/src/level-1/swap.lisp index 9464bc8..d9febc0 100644 --- a/src/level-1/swap.lisp +++ b/src/level-1/swap.lisp @@ -34,11 +34,12 @@ ;;Use only after checking the arguments for compatibility. (let* ((opt (get-tensor-class-optimization-hashtable tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) - `(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 :swap) ',func - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :swap) ',func + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func (x y) (declare (type ,tensor-class x y)) ,(let diff --git a/src/level-1/tensor-maker.lisp b/src/level-1/tensor-maker.lisp index f4e66bc..1620666 100644 --- a/src/level-1/tensor-maker.lisp +++ b/src/level-1/tensor-maker.lisp @@ -3,11 +3,12 @@ (defmacro make-tensor-maker (func-name (tensor-class)) (let ((opt (get-tensor-class-optimization-hashtable tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) - `(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 :maker) ',func-name - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :maker) ',func-name + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func-name (&rest args) (labels ((make-dims (dims) (declare (type cons dims)) @@ -67,11 +68,12 @@ (defmacro make-zeros-dims (func-name (tensor-class)) (let ((opt (get-tensor-class-optimization-hashtable tensor-class))) (assert opt nil 'tensor-cannot-find-optimization :tensor-class tensor-class) - `(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 :zero-maker) ',func-name - (get-tensor-class-optimization ',tensor-class) opt)) + `(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 :zero-maker) ',func-name + (get-tensor-class-optimization ',tensor-class) opt))) (defun ,func-name (dims) (declare (type (or cons index-store-vector) dims)) (let*-typed ((dims (if (consp dims) (make-index-store dims) (copy-seq dims)) :type index-store-vector) diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp index 033d36e..9307caa 100644 --- a/src/level-2/gemv.lisp +++ b/src/level-2/gemv.lisp @@ -7,11 +7,12 @@ (error 'tensor-cannot-find-optimization :tensor-class tensor-class))) (matrix-class (getf opt :matrix)) (vector-class (getf opt :vector))) - `(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)) + `(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) diff --git a/src/level-3/gemm.lisp b/src/level-3/gemm.lisp index bd59f29..b24fc76 100644 --- a/src/level-3/gemm.lisp +++ b/src/level-3/gemm.lisp @@ -33,11 +33,12 @@ (error 'tensor-cannot-find-optimization :tensor-class tensor-class))) (matrix-class (getf opt :matrix)) (blas? blas-gemm-func)) - `(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)) + `(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) ----------------------------------------------------------------------- Summary of changes: src/ffi/f77-ffi.lisp | 2 +- src/lapack/getrf.lisp | 26 ++++++++++++------------ src/level-1/axpy.lisp | 22 +++++++++++--------- src/level-1/copy.lisp | 22 +++++++++++--------- src/level-1/dot.lisp | 11 +++++---- src/level-1/scal.lisp | 44 ++++++++++++++++++++++------------------ src/level-1/swap.lisp | 11 +++++---- src/level-1/tensor-maker.lisp | 22 +++++++++++--------- src/level-2/gemv.lisp | 11 +++++---- src/level-3/gemm.lisp | 11 +++++---- 10 files changed, 98 insertions(+), 84 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-03 09:29:14
|
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, tensor has been updated via 2f7115f30d58288fdc9b160ce621ab19d1277772 (commit) from 468ad7e97a5c82ae3ed1c39d47085811b719e8b0 (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 2f7115f30d58288fdc9b160ce621ab19d1277772 Author: Akshay Srinivasan <aks...@gm...> Date: Sun Feb 3 01:24:22 2013 -0800 Moved lapack code from old/ to lapack/. diff --git a/src/old/geev.lisp b/src/lapack/geev.lisp similarity index 100% rename from src/old/geev.lisp rename to src/lapack/geev.lisp diff --git a/src/old/gels.lisp b/src/lapack/gels.lisp similarity index 100% rename from src/old/gels.lisp rename to src/lapack/gels.lisp diff --git a/src/old/geqr.lisp b/src/lapack/geqr.lisp similarity index 100% rename from src/old/geqr.lisp rename to src/lapack/geqr.lisp diff --git a/src/old/gesv.lisp b/src/lapack/gesv.lisp similarity index 100% rename from src/old/gesv.lisp rename to src/lapack/gesv.lisp diff --git a/src/old/potrf.lisp b/src/lapack/potrf.lisp similarity index 100% rename from src/old/potrf.lisp rename to src/lapack/potrf.lisp diff --git a/src/old/potrs.lisp b/src/lapack/potrs.lisp similarity index 100% rename from src/old/potrs.lisp rename to src/lapack/potrs.lisp diff --git a/src/old/svd.lisp b/src/lapack/svd.lisp similarity index 100% rename from src/old/svd.lisp rename to src/lapack/svd.lisp ----------------------------------------------------------------------- Summary of changes: src/{old => lapack}/geev.lisp | 0 src/{old => lapack}/gels.lisp | 0 src/{old => lapack}/geqr.lisp | 0 src/{old => lapack}/gesv.lisp | 0 src/{old => lapack}/potrf.lisp | 0 src/{old => lapack}/potrs.lisp | 0 src/{old => lapack}/svd.lisp | 0 7 files changed, 0 insertions(+), 0 deletions(-) rename src/{old => lapack}/geev.lisp (100%) rename src/{old => lapack}/gels.lisp (100%) rename src/{old => lapack}/geqr.lisp (100%) rename src/{old => lapack}/gesv.lisp (100%) rename src/{old => lapack}/potrf.lisp (100%) rename src/{old => lapack}/potrs.lisp (100%) rename src/{old => lapack}/svd.lisp (100%) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-03 09:15:30
|
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, tensor has been updated via 468ad7e97a5c82ae3ed1c39d47085811b719e8b0 (commit) from d43798f53f1a103b043af6fa6742ea347c777a2f (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 468ad7e97a5c82ae3ed1c39d47085811b719e8b0 Author: Akshay Srinivasan <aks...@gm...> Date: Sun Feb 3 01:10:31 2013 -0800 Fixed some bugs. Finally ported getrf(LU) of LAPACK. diff --git a/matlisp.asd b/matlisp.asd index 520d66c..e83b2ee 100644 --- a/matlisp.asd +++ b/matlisp.asd @@ -153,10 +153,10 @@ :pathname "level-3" :depends-on ("matlisp-base" "matlisp-classes" "foreign-core" "matlisp-level-1" "matlisp-level-2") :components ((:file "gemm"))) - #+nil(:module "matlisp-lapack" + (:module "matlisp-lapack" :pathname "lapack" :depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3") - :components ((:file "gesv"))) + :components ((:file "getrf"))) #+nil (:module "matlisp-sugar" :pathname "sugar" diff --git a/packages.lisp b/packages.lisp index 5f02a07..64e09af 100644 --- a/packages.lisp +++ b/packages.lisp @@ -62,6 +62,7 @@ #:tensor-cannot-find-optimization #:tensor-class #:tensor-dimension-mismatch #:tensor-store-not-consecutive + #:tensor-method-does-not-exist )) (defpackage "MATLISP-UTILITIES" diff --git a/src/base/permutation.lisp b/src/base/permutation.lisp index a0d1eb0..778ceba 100644 --- a/src/base/permutation.lisp +++ b/src/base/permutation.lisp @@ -106,12 +106,10 @@ (declare (ignore arg)) (let ((len (length seq))) (assert (>= len (permutation-size perm)) nil - 'permutation-permute-error :seq-len len :per-size (permutation-size perm))))) - -;; (:method :before ((ten standard-tensor) (perm permutation) &optional (arg 0)) -;; (let ((len (aref (dimensions ten) arg))) -;; (assert (>= len (permutation-size perm)) nil -;; 'permutation-permute-error :seq-len len :permutation-size (permutation-size perm))))) + 'permutation-permute-error :seq-len len :per-size (permutation-size perm)))) + (:method :before ((ten standard-tensor) (perm permutation) &optional (arg 0)) + (assert (>= (aref (dimensions ten) arg) (permutation-size perm)) nil + 'permutation-permute-error :seq-len (aref (dimensions ten) arg) :permutation-size (permutation-size perm)))) (definline permute (thing perm &optional (arg 0)) (permute! (copy thing) perm arg)) @@ -136,9 +134,8 @@ :do (setf (aref seq i) (aref cseq (aref act i))) :finally (return seq)))) -;; (defmethod permute! ((ten standard-tensor) (perm permutation-action) &optional (arg 0)) -;; (let ((cyc (action->cycle perm))) -;; (permute! ten cyc arg))) +(defmethod permute! ((ten standard-tensor) (perm permutation-action) &optional (arg 0)) + (permute! ten (action->pivot-flip perm) arg)) ;;Cycle (definline apply-cycle! (seq pcyc) @@ -168,29 +165,8 @@ :do (apply-cycle! seq cyc))) seq) -;; (defmethod permute! ((A standard-tensor) (perm permutation-cycle) &optional (arg 0)) -;; (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) -;; (rplaca (nthcdr arg slst) 0) -;; (values (sub-tensor~ A slst) (sub-tensor~ A slst))) -;; (let-typed ((cyclst (store perm) :type cons) -;; (cp-ten (make-instance (class-of tone) -;; :dimensions (copy-seq (dimensions tone)))) -;; (std-arg (aref (strides A) arg) :type index-type) -;; (hd-sl (head ttwo) :type index-type)) -;; (dolist (cyc cyclst) -;; (declare (type pindex-store-vector cyc)) -;; (setf (head tone) (+ hd-sl (* std-arg (aref cyc (1- (length cyc)))))) -;; (copy! tone cp-ten) -;; (loop :for i :of-type index-type :downfrom (1- (length cyc)) :to 0 -;; :do (progn -;; (setf (head tone) (+ hd-sl (* std-arg (aref cyc i)))) -;; (copy! -;; (if (= i 0) cp-ten -;; (progn -;; (setf (head ttwo) (+ hd-sl (* std-arg (aref cyc (1- i))))) -;; ttwo)) -;; tone)))))) -;; A) +(defmethod permute! ((A standard-tensor) (perm permutation-cycle) &optional (arg 0)) + (permute! A (action->pivot-flip (cycle->action perm)) arg)) ;Pivot idx (definline apply-flips! (seq pflip) @@ -213,20 +189,21 @@ (copy-n cseq seq size)) seq) -;; (defmethod permute! ((A standard-tensor) (perm permutation-pivot-flip) &optional (arg 0)) -;; (let ((idiv (store perm))) -;; (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) -;; (rplaca (nthcdr arg slst) 0) -;; (values (sub-tensor~ A slst nil) (sub-tensor~ A slst nil))) -;; (let ((argstd (aref (strides A) arg)) -;; (hd-sl (head ttwo))) -;; (loop :for i :from 0 :below (length idiv) -;; :do (progn -;; (unless (= i (aref idiv i)) -;; (setf (head ttwo) (+ hd-sl (* (aref idiv i) argstd))) -;; (swap! tone ttwo)) -;; (incf (head tone) argstd)))))) -;; A) +(defmethod permute! ((A standard-tensor) (perm permutation-pivot-flip) &optional (arg 0)) + (multiple-value-bind (t1 t2) (let ((slst (make-list (rank A) :initial-element '(* * *)))) + (rplaca (nthcdr arg slst) (list 0 '* 1)) + (values (sub-tensor~ A slst nil) (sub-tensor~ A slst nil))) + (let-typed ((argstd (aref (strides A) arg) :type index-type) + (hd-sl (head t2) :type index-type) + (idiv (store perm) :type pindex-store-vector)) + (very-quickly + (loop :for i :from 0 :below (length idiv) + :do (progn + (unless (= i (aref idiv i)) + (setf (head t2) (the index-type (+ hd-sl (the index-type (* (aref idiv i) argstd))))) + (swap! t1 t2)) + (setf (head t1) (the index-type (+ argstd (the index-type (head t1)))))))))) + A) ;;Conversions----------------------------------------------------;; (defun action->cycle (act) diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp index 12eb8d2..4e8e80a 100644 --- a/src/base/standard-tensor.lisp +++ b/src/base/standard-tensor.lisp @@ -374,7 +374,7 @@ (declare (type index-type lidx)) (let-typed ((sto-x (store tensor) :type ,(linear-array-type store-element-type))) (,reader sto-x lidx))) - (defmethod (setf tensor-ref) (value (tensor ,tensor-class) lidx) + (defmethod (setf tensor-store-ref) (value (tensor ,tensor-class) lidx) (declare (type index-type lidx)) (let-typed ((sto-x (store tensor) :type ,(linear-array-type store-element-type))) (,value-writer (,coercer-unforgiving value) sto-x lidx))) diff --git a/src/conditions.lisp b/src/conditions.lisp index 862dab0..43bdacc 100644 --- a/src/conditions.lisp +++ b/src/conditions.lisp @@ -224,3 +224,11 @@ (:report (lambda (c stream) (declare (ignore c)) (format stream "The strides of the store, of the given tensor are not conscutive.")))) + +(define-condition tensor-method-does-not-exist (tensor-error) + ((method-name :reader method-name :initarg :method-name) + (tensor-class :reader tensor-class :initarg :tensor-class)) + (:documentation "The method is not defined for the given tensor class.") + (:report (lambda (c stream) + (when (slots-boundp c 'method 'tensor-class) + (format stream "The tensor class: ~a, does not have the method: ~a defined as a specialization." (method-name c) (tensor-class c)))))) diff --git a/src/lapack/getrf.lisp b/src/lapack/getrf.lisp index da4bea0..d5a17ae 100644 --- a/src/lapack/getrf.lisp +++ b/src/lapack/getrf.lisp @@ -3,14 +3,14 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; ;;; Copyright (c) 2000 The Regents of the University of California. -;;; All rights reserved. -;;; +;;; 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 @@ -44,11 +44,11 @@ (if (eq maj-A :col-major) (multiple-value-bind (n-A n-ipiv info) (,lapack-func (nrows A) (ncols A) (store A) - ld-A (repr ipiv) 0) + ld-A (store ipiv) 0) (declare (ignore n-A n-ipiv)) (assert (= info 0) nil 'invalid-arguments :argnum (1- (- info)) :message (format-to-string "GETRF returned INFO: ~a." info)) ;;Convert back to 0-based indexing. - (let-typed ((pidv (repr ipiv) :type perrepr-vector)) + (let-typed ((pidv (store ipiv) :type pindex-store-vector)) (very-quickly (loop :for i :of-type index-type :from 0 :below (length pidv) :do (decf (aref pidv i))))) @@ -74,54 +74,58 @@ Given an NxM matrix A, compute its LU factorization using partial pivoting, row or column interchanges: - A = P * L * U (if A is row-major ordered) - A = L * U * P' (if A is col-major ordered) + A = P * L * U (if A is row-major ordered) + A = L * U * P' (if A is col-major ordered) where: - P: permutation matrix - L: lower triangular with unit diagonal elements - (lower trapezoidal when N>M) - U: upper triangular - (upper trapezoidal when N<M) + P: permutation matrix + L: lower triangular with unit diagonal elements + (lower trapezoidal when N>M) + U: upper triangular + (upper trapezoidal when N<M) Return Values ============= - [1] The factors L and U from the factorization A = P*L*U where the + [1] The factors L and U from the factorization A = P*L*U where the unit diagonal elements of L are not stored. (overwriting A) [2] IPIV [3] INFO = T: successful - i: U(i,i) is exactly zero. + i: U(i,i) is exactly zero. ")) (defmethod getrf! ((A real-matrix)) - (let ((ipiv (make-instance 'permutation-pivot-flip - :repr (perrepr-id-action (lvec-min (dimensions A)))))) + (let ((ipiv (let* ((*check-after-initializing?* nil) + (ret (make-instance 'permutation-pivot-flip :store (pindex-id (lvec-min (dimensions A)))))) + (setf (permutation-size ret) (length (store ret))) + ret))) (real-typed-getrf! A ipiv))) (defmethod getrf! ((A complex-matrix)) - (let ((ipiv (make-instance 'permutation-pivot-flip - :repr (perrepr-id-action (lvec-min (dimensions A)))))) + (let ((ipiv (let* ((*check-after-initializing?* nil) + (ret (make-instance 'permutation-pivot-flip :store (pindex-id (lvec-min (dimensions A)))))) + (setf (permutation-size ret) (length (store ret))) + ret))) (complex-typed-getrf! A ipiv))) ;; -(defgeneric lu (a &key with-l with-u with-p) +(defgeneric lu (a &optional split-lu?) (:documentation " Syntax ====== (LU a [:WITH-P with-p] [:WITH-L with-l] [:WITH-U with-u]) - + Purpose ======= - Computes the LU decomposition of A. + Computes the LU decomposition of A. This functions is an interface to GETRF! Return Values ============= [1] the factors L,U from the factorization in a single matrix, - where the unit diagonal elements of L are not stored + where the unit diagonal elements of L are not stored [2]-[4] If WITH-X then X, in the order L,U,P By default WITH-L,WITH-U,WITH-P. @@ -130,81 +134,70 @@ ;;Sure I can do this with a defmethod (slowly), but where's the fun in that ? :) (defmacro make-lu (tensor-class) (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))) + (error 'tensor-cannot-find-optimization :tensor-class tensor-class))) + (matrix-class (getf opt :matrix))) `(eval-when (:compile-toplevel :load-toplevel :execute) - (defmethod lu ((A ,matrix-class) &key with-l with-u with-p) + (defmethod lu ((A ,matrix-class) &optional (split-lu? t)) (multiple-value-bind (lu ipiv info) (getrf! (with-order :col-major (,(getf opt :copy) A (,(getf opt :zero-maker) (dimensions A))))) (declare (ignore info)) - (let* ((ret (list lu)) - (n (nrows a)) + (let* ((n (nrows a)) (m (ncols a)) (p (min n m))) (declare (type fixnum n m p)) ;; Extract the lower triangular part, if requested - (when with-l - (let*-typed ((lmat (,(getf opt :zero-maker) (list n p)) :type ,matrix-class) - ;; - (l.rstd (row-stride lmat) :type index-type) - (l.cstd (col-stride lmat) :type index-type) - (l.of (head lmat) :type index-type) - (l.sto (store lmat) :type ,(linear-array-type (getf opt :element-type))) - ;; - (lu.rstd (row-stride lu) :type index-type) - (lu.cstd (col-stride lu) :type index-type) - (lu.of (head lu) :type index-type) - (lu.sto (store lu) :type ,(linear-array-type (getf opt :element-type)))) - (very-quickly - (loop :for j :of-type index-type :from 0 :below p - :do (progn - (,(getf opt :value-writer) (,(getf opt :fid*)) l.sto l.of) - (loop :repeat (- n j 1) - :do (progn - (incf lu.of lu.rstd) - (incf l.of l.rstd) - (,(getf opt :reader-writer) lu.sto lu.of l.sto l.of))) - (incf lu.of (- lu.cstd (the index-type (* (- n j 2) lu.rstd)))) - (incf l.of (- l.cstd (the index-type (* (- n j 2) l.rstd)))))) - (push lmat ret)))) - ;; Extract the upper triangular part, if requested - (when with-u - (let*-typed ((umat (,(getf opt :zero-maker) (list p m)) :type ,matrix-class) - ;; - (u.rstd (row-stride umat) :type index-type) - (u.cstd (col-stride umat) :type index-type) - (u.of (head umat) :type index-type) - (u.sto (store umat) :type ,(linear-array-type (getf opt :element-type))) - ;; - (lu.rstd (row-stride lu) :type index-type) - (lu.cstd (col-stride lu) :type index-type) - (lu.of (head lu) :type index-type) - (lu.sto (store lu) :type ,(linear-array-type (getf opt :element-type)))) - (progn - (loop :for i :of-type index-type :from 0 :below p - :do (progn - (loop :repeat (- m i) - :do (progn - (,(getf opt :reader-writer) lu.sto lu.of u.sto u.of) - (incf lu.of lu.cstd) - (incf u.of u.cstd))) - (incf lu.of (- lu.rstd (the index-type (* (- m i 1) lu.cstd)))) - (incf u.of (- u.rstd (the index-type (* (- m i 1) u.cstd)))))) - (push umat ret)))) - ;; Extract the permutation matrix, if requested - (if with-p - (let* ((npiv (length ipiv)) - (pmat (make-real-matrix-dim n n)) - (pidx (make-array n :element-type '(unsigned-byte 32)))) - ;; Compute the P matrix from the pivot vector - (dotimes (k n) - (setf (aref pidx k) k)) - (dotimes (k npiv) - (rotatef (aref pidx k) (aref pidx (1- (aref ipiv k))))) - (dotimes (k n) - (setf (matrix-ref pmat (aref pidx k) k) 1)) - (push pmat result))) - ;; Return the final result - (values-list (nreverse ret)))))))) + (if split-lu? + (if (= p m) + (let*-typed ((umat (,(getf opt :zero-maker) (list p m)) :type ,matrix-class) + ;; + (u.rstd (row-stride umat) :type index-type) + (u.cstd (col-stride umat) :type index-type) + (u.of (head umat) :type index-type) + (u.sto (store umat) :type ,(linear-array-type (getf opt :store-type))) + ;; + (lu.rstd (row-stride lu) :type index-type) + (lu.cstd (col-stride lu) :type index-type) + (lu.of (head lu) :type index-type) + (lu.sto (store lu) :type ,(linear-array-type (getf opt :store-type)))) + (very-quickly + (loop :for i :of-type index-type :from 0 :below p + :do (let-typed ((lu.of-ii lu.of :type index-type)) + (loop :repeat (- m i) + :do (progn + (,(getf opt :reader-writer) lu.sto lu.of u.sto u.of) + (,(getf opt :value-writer) (,(getf opt :fid+)) lu.sto lu.of) + (incf lu.of lu.cstd) + (incf u.of u.cstd))) + (,(getf opt :value-writer) (,(getf opt :fid*)) lu.sto lu.of-ii) + (incf lu.of (- lu.rstd (the index-type (* (- m i 1) lu.cstd)))) + (incf u.of (- u.rstd (the index-type (* (- m i 1) u.cstd))))))) + (values lu umat ipiv)) + (let*-typed ((lmat (,(getf opt :zero-maker) (list n p)) :type ,matrix-class) + ;; + (l.rstd (row-stride lmat) :type index-type) + (l.cstd (col-stride lmat) :type index-type) + (l.of (head lmat) :type index-type) + (l.sto (store lmat) :type ,(linear-array-type (getf opt :store-type))) + ;; + (lu.rstd (row-stride lu) :type index-type) + (lu.cstd (col-stride lu) :type index-type) + (lu.of (head lu) :type index-type) + (lu.sto (store lu) :type ,(linear-array-type (getf opt :store-type)))) + (very-quickly + (loop :for j :of-type index-type :from 0 :below p + :do (progn + (,(getf opt :value-writer) (,(getf opt :fid*)) l.sto l.of) + (loop :repeat (- n j 1) + :do (progn + (incf lu.of lu.rstd) + (incf l.of l.rstd) + (,(getf opt :reader-writer) lu.sto lu.of l.sto l.of) + (,(getf opt :value-writer) (,(getf opt :fid+)) lu.sto lu.of))) + (incf lu.of (- lu.cstd (the index-type (* (- n j 2) lu.rstd)))) + (incf l.of (- l.cstd (the index-type (* (- n j 2) l.rstd))))))) + (values lmat lu ipiv))) + (values lu ipiv)))))))) (make-lu real-tensor) +(make-lu complex-tensor) ----------------------------------------------------------------------- Summary of changes: matlisp.asd | 4 +- packages.lisp | 1 + src/base/permutation.lisp | 69 ++++++----------- src/base/standard-tensor.lisp | 2 +- src/conditions.lisp | 8 ++ src/lapack/getrf.lisp | 171 ++++++++++++++++++++--------------------- 6 files changed, 117 insertions(+), 138 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-03 02:23:05
|
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, tensor has been updated via d43798f53f1a103b043af6fa6742ea347c777a2f (commit) from 7a8ab5c0db938424bfc8d2ebca6022e2673b6a9a (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 d43798f53f1a103b043af6fa6742ea347c777a2f Author: Akshay Srinivasan <aks...@gm...> Date: Sat Feb 2 18:22:45 2013 -0800 Reworked a lot of ugly code in permutation.lisp (I really was a newbie 6 months ago\!). diff --git a/packages.lisp b/packages.lisp index 87b3afc..5f02a07 100644 --- a/packages.lisp +++ b/packages.lisp @@ -66,7 +66,8 @@ (defpackage "MATLISP-UTILITIES" (:use #:common-lisp #:matlisp-conditions) - (:export #:ensure-list #:id ;;#:make-ring + (:export #:ensure-list #:id + #:vectorify #:copy-n #:zip #:zip-eq #:cut-cons-chain! #:slot-values diff --git a/src/base/permutation.lisp b/src/base/permutation.lisp index 9f5bbda..a0d1eb0 100644 --- a/src/base/permutation.lisp +++ b/src/base/permutation.lisp @@ -70,16 +70,16 @@ (when *check-after-initializing?* (if (null (store per)) (setf (permutation-size per) 1) - (very-quickly - (loop - :for cyc :of-type pindex-store-vector :in (store per) - :with ss :of-type pindex-type := 0 - :do (loop + (loop + :for cyc :of-type pindex-store-vector :in (store per) + :with ss :of-type pindex-type := 0 + :do (very-quickly + (loop :for i :of-type index-type :from 1 :below (length cyc) :with scyc :of-type pindex-store-vector := (sort (copy-seq cyc) #'<) :do (assert (/= (aref scyc (1- i)) (aref scyc i)) nil 'permutation-invalid-error) - :finally (setf ss (max ss (aref scyc (1- (length scyc)))))) - :finally (setf (permutation-size per) (1+ ss))))))) + :finally (setf ss (max ss (aref scyc (1- (length scyc))))))) + :finally (setf (permutation-size per) (1+ ss)))))) ;; (defclass permutation-pivot-flip (permutation) @@ -93,21 +93,7 @@ (very-quickly (loop :for i :of-type index-type :from 0 :below len :do (assert (< -1 (aref repr i) len) nil 'permutation-invalid-error))) - (setf (store-size per) len)))) - -;; -(defclass permutation-matrix (permutation) - ((store :type pindex-store-vector))) - -(defmethod initialize-instance :after ((perm permutation-matrix) &rest initargs) - (declare (ignore initargs)) - (when *check-after-initializing?* - (let-typed ((repr (store perm) :type pindex-store-vector)) - (very-quickly - (loop :for i :of-type index-type :from 0 :below (length repr) - :with srepr :of-type pindex-store-vector := (sort (copy-seq repr) #'<) - :do (assert (= (aref srepr i) i) nil 'permutation-invalid-error))) - (setf (permutation-size perm) (length repr))))) + (setf (permutation-size per) len)))) ;;Generic permute! method. (defgeneric permute! (thing permutation &optional argument) @@ -119,80 +105,74 @@ (:method :before ((seq sequence) (perm permutation) &optional (arg 0)) (declare (ignore arg)) (let ((len (length seq))) - (assert (>= len (group-rank perm)) nil - 'permutation-permute-error :seq-len len :group-rank (group-rank perm)))) - (:method :before ((ten standard-tensor) (perm permutation) &optional (arg 0)) - (let ((len (aref (dimensions ten) arg))) - (assert (>= len (group-rank perm)) nil - 'permutation-permute-error :seq-len len :group-rank (group-rank perm))))) + (assert (>= len (permutation-size perm)) nil + 'permutation-permute-error :seq-len len :per-size (permutation-size perm))))) + +;; (:method :before ((ten standard-tensor) (perm permutation) &optional (arg 0)) +;; (let ((len (aref (dimensions ten) arg))) +;; (assert (>= len (permutation-size perm)) nil +;; 'permutation-permute-error :seq-len len :permutation-size (permutation-size perm))))) (definline permute (thing perm &optional (arg 0)) (permute! (copy thing) perm arg)) -;; ;;Action -;; (defmethod permute! ((seq cons) (perm permutation-action) &optional arg) -;; (declare (ignore arg)) -;; (let ((cseq (make-array (length seq) :initial-contents seq)) -;; (act (repr perm)) -;; (glen (group-rank perm))) -;; (mapl -;; (let ((i 0)) -;; (declare (type fixnum i)) -;; (lambda (x) -;; (when (< i glen) -;; (rplaca x (aref cseq (aref act i))) -;; (incf i)))) seq))) - -;; (defmethod permute! ((seq vector) (perm permutation-action) &optional arg) -;; (declare (ignore arg)) -;; (let ((cseq (make-array (length seq) :initial-contents seq)) -;; (act (repr perm))) -;; (loop -;; :for i :from 0 :below (group-rank perm) -;; :do (unless (= i (aref act i)) -;; (setf (aref seq i) (aref cseq (aref act i)))) -;; :finally (return seq)))) +;;Action +(defmethod permute! ((seq cons) (perm permutation-action) &optional arg) + (declare (ignore arg)) + (let* ((size (permutation-size perm)) + (cseq (vectorify seq size)) + (act (store perm))) + (loop :for i :from 0 :below size + :for lst := seq :then (cdr lst) + :do (setf (car lst) (aref cseq (aref act i))) + :finally (return seq)))) + +(defmethod permute! ((seq vector) (perm permutation-action) &optional arg) + (declare (ignore arg)) + (let* ((size (permutation-size perm)) + (cseq (vectorify seq size)) + (act (store perm))) + (loop :for i :from 0 :below size + :do (setf (aref seq i) (aref cseq (aref act i))) + :finally (return seq)))) ;; (defmethod permute! ((ten standard-tensor) (perm permutation-action) &optional (arg 0)) ;; (let ((cyc (action->cycle perm))) ;; (permute! ten cyc arg))) -;; ;;Cycle -;; ;;Might be useful ? -;; (defun apply-cycle! (seq pcyc) -;; (declare (type pindex-store-vector pcyc) -;; (type vector seq)) -;; (let ((xl (aref seq (aref pcyc (1- (length pcyc)))))) -;; (loop :for i :of-type index-type :downfrom (1- (length pcyc)) :to 0 -;; :do (setf (aref seq (aref pcyc i)) -;; (if (= i 0) xl -;; (aref seq (aref pcyc (1- i)))))))) - -;; (defmethod permute! ((seq cons) (perm permutation-cycle) &optional arg) -;; (declare (ignore arg)) -;; (let ((cseq (make-array (length seq) :initial-contents seq)) -;; (glen (group-rank perm))) -;; (dolist (cyc (repr perm)) -;; (declare (type pindex-store-vector cyc)) -;; (apply-cycle! cseq cyc)) -;; (mapl -;; (let ((i 0)) -;; (lambda (x) -;; (when (< i glen) -;; (rplaca x (aref cseq i)) -;; (incf i)))) seq))) - -;; (defmethod permute! ((seq vector) (perm permutation-cycle) &optional arg) -;; (declare (ignore arg)) -;; (dolist (cyc (repr perm) seq) -;; (declare (type pindex-store-vector cyc)) -;; (apply-cycle! seq cyc))) +;;Cycle +(definline apply-cycle! (seq pcyc) + (declare (type pindex-store-vector pcyc) + (type vector seq)) + (loop :for i :of-type index-type :downfrom (1- (length pcyc)) :to 1 + :with xl := (aref seq (aref pcyc (1- (length pcyc)))) + :do (setf (aref seq (aref pcyc i)) (aref seq (aref pcyc (1- i)))) + :finally (progn + (setf (aref seq (aref pcyc 0)) xl) + (return seq)))) + +(defmethod permute! ((seq cons) (perm permutation-cycle) &optional arg) + (declare (ignore arg)) + (unless (= (permutation-size perm) 1) + (let* ((size (permutation-size perm)) + (cseq (vectorify seq size))) + (loop :for cyc :of-type pindex-store-vector :in (store perm) + :do (apply-cycle! cseq cyc)) + (copy-n cseq seq size))) + seq) + +(defmethod permute! ((seq vector) (perm permutation-cycle) &optional arg) + (declare (ignore arg)) + (unless (= (permutation-size perm) 1) + (loop :for cyc :of-type pindex-store-vector :in (store perm) + :do (apply-cycle! seq cyc))) + seq) ;; (defmethod permute! ((A standard-tensor) (perm permutation-cycle) &optional (arg 0)) ;; (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) ;; (rplaca (nthcdr arg slst) 0) ;; (values (sub-tensor~ A slst) (sub-tensor~ A slst))) -;; (let-typed ((cyclst (repr perm) :type cons) +;; (let-typed ((cyclst (store perm) :type cons) ;; (cp-ten (make-instance (class-of tone) ;; :dimensions (copy-seq (dimensions tone)))) ;; (std-arg (aref (strides A) arg) :type index-type) @@ -212,29 +192,29 @@ ;; tone)))))) ;; A) -;; ;;Pivot idx -;; (defmethod permute! ((seq vector) (perm permutation-pivot-flip) &optional arg) -;; (declare (ignore arg)) -;; (let-typed ((pidx (repr perm) :type pindex-store-vector)) -;; (loop :for i :of-type index-type :from 0 :below (group-rank perm) -;; :unless (= i (aref pidx i)) -;; :do (rotatef (aref seq i) (aref seq (aref pidx i))) -;; :finally (return seq)))) - -;; (defmethod permute! ((seq cons) (perm permutation-pivot-flip) &optional arg) -;; (declare (ignore arg)) -;; (let ((cseq (make-array (length seq) :initial-contents seq)) -;; (glen (group-rank perm))) -;; (permute! cseq perm) -;; (mapl -;; (let ((i 0)) -;; (lambda (x) -;; (when (< i glen) -;; (rplaca x (aref cseq i)) -;; (incf i)))) seq))) +;Pivot idx +(definline apply-flips! (seq pflip) + (declare (type pindex-store-vector pflip) + (type vector seq)) + (loop :for i :of-type index-type :from 0 :below (length pflip) + :unless (= i (aref pflip i)) + :do (rotatef (aref seq i) (aref seq (aref pflip i))) + :finally (return seq))) + +(defmethod permute! ((seq vector) (perm permutation-pivot-flip) &optional arg) + (declare (ignore arg)) + (apply-flips! seq (store perm))) + +(defmethod permute! ((seq cons) (perm permutation-pivot-flip) &optional arg) + (declare (ignore arg)) + (let* ((size (permutation-size perm)) + (cseq (vectorify seq size))) + (apply-flips! cseq (store perm)) + (copy-n cseq seq size)) + seq) ;; (defmethod permute! ((A standard-tensor) (perm permutation-pivot-flip) &optional (arg 0)) -;; (let ((idiv (repr perm))) +;; (let ((idiv (store perm))) ;; (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) ;; (rplaca (nthcdr arg slst) 0) ;; (values (sub-tensor~ A slst nil) (sub-tensor~ A slst nil))) @@ -248,186 +228,151 @@ ;; (incf (head tone) argstd)))))) ;; A) -;; ;;Conversions----------------------------------------------------;; -;; (defun action->cycle (act) -;; " -;; (action->cycle act) - -;; This function obtains the canonical cycle representation -;; of a permutation. The first argument \"act\" is the action of the -;; permutation on the array #(0 1 2 3 ..): an object of the class -;; permutation-action. - -;; \"Canonical\" may be a bit of an overstatement; this is the way -;; S_n was presented in Van der Waerden's book. -;; " -;; (declare (type permutation-action act)) -;; (mlet* -;; ((arr (repr act) :type pindex-store-vector)) -;; (labels ((find-cycle (x0) -;; ;; This function obtains the cycle starting from x_0. -;; (declare (type pindex-type x0)) -;; (if (= (aref arr x0) x0) (values 0 nil) -;; (very-quickly -;; (loop -;; :for x :of-type pindex-type := (aref arr x0) :then (aref arr x) -;; :and ret :of-type cons := (list x0) :then (cons x ret) -;; :counting t :into i :of-type index-type -;; :when (= x x0) -;; :do (return (values i ret)))))) -;; (cycle-walk (cyc ignore) -;; ;; Finds all cycles -;; (let ((x0 (find-if-not #'(lambda (x) (member x ignore)) arr))) -;; (if (null x0) -;; cyc -;; (multiple-value-bind (clen clst) (find-cycle x0) -;; (declare (type index-type clen) -;; (type list clst)) -;; (cycle-walk -;; (if (= clen 0) cyc -;; (cons (make-array clen :element-type 'pindex-type :initial-contents clst) cyc)) -;; (nconc ignore (if (= clen 0) (list x0) clst)))))))) -;; (let ((cyc-lst (cycle-walk nil nil))) -;; (make-instance 'permutation-cycle -;; :repr cyc-lst))))) - -;; (defun cycle->action (cyc) -;; " -;; (cycle->action cyc) - -;; This function obtains the action representation of a permutation -;; from the cyclic one. The first argument \"cyc\" is the cyclic -;; representation of the permutation: an object of the class -;; permutation-cycle. -;; " -;; (declare (type permutation-cycle cyc)) -;; (let ((act-repr (pindex-id-action (group-rank cyc))) -;; (cycs-repr (repr cyc))) -;; (declare (type pindex-store-vector act-repr)) -;; (dolist (cyc cycs-repr) -;; (declare (type pindex-store-vector cyc)) -;; (let ((xl (aref act-repr (aref cyc (1- (length cyc)))))) -;; (very-quickly -;; (loop -;; :for i :of-type index-type :downfrom (1- (length cyc)) :to 0 -;; :do (setf (aref act-repr (aref cyc i)) -;; (if (= i 0) xl -;; (aref act-repr (aref cyc (1- i))))))))) -;; (make-instance 'permutation-action :repr act-repr))) - -;; (defun pivot-flip->action (pflip) -;; (declare (type permutation-pivot-flip pflip)) -;; (let* ((idiv (repr pflip)) -;; (len (length idiv))) -;; (declare (type pindex-store-vector idiv) -;; (type index-type len)) -;; (let ((act (pindex-id-action len))) -;; (declare (type pindex-store-vector act)) -;; (very-quickly -;; (loop :for i :from 0 :below len -;; :do (let ((val (aref idiv i))) -;; (unless (= val i) -;; (rotatef (aref act i) (aref act val)))))) -;; (make-instance 'permutation-action :repr act)))) - -;; (defun mod-max (seq lidx uidx) -;; (declare (type pindex-store-vector seq)) -;; (let ((len (length seq))) -;; (very-quickly -;; (loop :for idx :of-type index-type :downfrom uidx :above lidx -;; :with max :of-type pindex-type := (aref seq uidx) -;; :with max-idx :of-type index-type := uidx -;; :do (let ((ele (aref seq (mod idx len)))) -;; (when (> ele max) -;; (setf max ele -;; max-idx idx))) -;; :finally (return (values max max-idx)))))) - -;; #| - -;; (defun cycle->pivot-flip (cyc) -;; (let ((cp-cyc (copy-seq cyc))) -;; (let -;; (labels ((mod-max (seq lidx uidx) -;; (declare (type pindex-store-vector seq)) -;; (let ((len (length cyc))) -;; (very-quickly -;; (loop :for idx :of-type index-type :downfrom uidx :above lidx -;; :with max :of-type pindex-type := (aref seq uidx) -;; :with max-idx :of-type index-type := uidx -;; :do (let ((ele (aref seq (mod idx len)))) -;; (when (> ele max) -;; (setf max ele -;; max-idx idx))) -;; :finally (return (values max max-idx)))))) -;; (get-flip (lidx uidx) -;; (multiple-value-bind (max max-idx) (mod-max cyc lidx uidx) -;; (let ((ele-0 (aref cyc (mod max-idx len))) -;; (ele-1 (aref cyc (mod (1- max-idx) len)))) -;; (setf (aref pidx (max ele-0 ele-1)) -;; (min ele-0 ele-1)) -;; |# - -;; #+nil -;; (defun permute-argument (func-symbol perm) -;; (declare (type symbol func-symbol) -;; (type permutation perm)) -;; (let* ((glen (group-rank perm)) -;; (args (loop :for i :from 0 :below glen -;; :collect (gensym)))) -;; (eval `(lambda (,@args &rest rest) -;; (apply ',func-symbol (append (list ,@(permute! args perm)) rest)))))) - -;; (defun argument-permute (func perm) -;; (declare (type function func) -;; (type permutation perm)) -;; (lambda (&rest args) -;; (apply func (permute! args perm)))) - -;; (defun curry (func perm &rest curried-args) -;; (declare (type function func) -;; (type permutation perm)) -;; (lambda (&rest args) -;; (apply func (permute! (append curried-args args) perm)))) - -;; (defun compose (func-a func-b perm) -;; (declare (type function func-a func-b) -;; (type permutation perm)) -;; (lambda (&rest args) -;; (apply func-a (permute! (multiple-value-list (funcall func-b args)) perm)))) -;; ;; - - -;; (defstruct (ring (:constructor nil)) -;; (circular-list nil :type list) -;; (end-list nil :type list)) - -;; (defun make-ring (n) -;; (let ((ret (cons nil nil))) -;; (loop :for i :from 0 :below (1- n) -;; :with tail :of-type cons := ret -;; :do (setf (cdr tail) (cons nil nil) -;; tail (cdr tail)) -;; :finally (setf (cdr tail) ret)) -;; ret)) - - +;;Conversions----------------------------------------------------;; +(defun action->cycle (act) + " + (action->cycle act) + + This function obtains the canonical cycle representation + of a permutation. The first argument \"act\" is the action of the + permutation on the array #(0 1 2 3 ..): an object of the class + permutation-action. + + \"Canonical\" may be a bit of an overstatement; this is the way + S_n was presented in Van der Waerden's book. +" + (declare (type permutation-action act)) + (let-typed ((arr (store act) :type pindex-store-vector)) + (labels ((find-cycle (x0) + ;; This function obtains the cycle starting from x_0. + (declare (type pindex-type x0)) + (if (= (aref arr x0) x0) (values 0 nil) + (very-quickly + (loop + :for x :of-type pindex-type := (aref arr x0) :then (aref arr x) + :and ret :of-type cons := (list x0) :then (cons x ret) + :counting t :into i :of-type index-type + :when (= x x0) + :do (return (values i ret)))))) + (cycle-walk (cyc ignore) + ;; Finds all cycles + (let ((x0 (find-if-not #'(lambda (x) (member x ignore)) arr))) + (if (null x0) + cyc + (multiple-value-bind (clen clst) (find-cycle x0) + (declare (type index-type clen) + (type list clst)) + (cycle-walk + (if (= clen 0) cyc + (cons (make-array clen :element-type 'pindex-type :initial-contents clst) cyc)) + (nconc ignore (if (= clen 0) (list x0) clst)))))))) + (make-instance 'permutation-cycle :store (cycle-walk nil nil))))) + +(defun action->pivot-flip (act) + (declare (type permutation-action act)) + (let*-typed ((size (permutation-size act) :type index-type) + (actr (store act) :type pindex-store-vector) + (ret (pindex-id size) :type pindex-store-vector) + (inv (pindex-id size) :type pindex-store-vector) + (for (pindex-id size) :type pindex-store-vector)) + (very-quickly + (loop :for i :of-type index-type :from 0 :below size + :do (let ((flip (aref inv (aref actr i)))) + (setf (aref ret i) flip + (aref inv (aref for i)) flip + (aref for flip) (aref for i))))) + (make-instance 'permutation-pivot-flip :store ret))) + +(defun cycle->action (cyc) + " + (cycle->action cyc) + + This function obtains the action representation of a permutation + from the cyclic one. The first argument \"cyc\" is the cyclic + representation of the permutation: an object of the class + permutation-cycle. +" + (declare (type permutation-cycle cyc)) + (let-typed ((act-repr (pindex-id (permutation-size cyc)) :type pindex-store-vector) + (cycs (store cyc))) + (very-quickly + (loop :for cyc :of-type pindex-store-vector :in cycs + :do (apply-cycle! act-repr cyc))) + (make-instance 'permutation-action :store act-repr))) + +(defun pivot-flip->action (pflip) + (declare (type permutation-pivot-flip pflip)) + (let*-typed ((idiv (store pflip) :type pindex-store-vector) + (len (permutation-size pflip) :type index-type) + (ret (pindex-id len) :type pindex-store-vector)) + (make-instance 'permutation-action :store (very-quickly (apply-flips! ret idiv))))) + +;;Uber-functional stuff +;;None of these are ever useful (I've found), neat things for showing off though :] +(defun permute-arguments-and-compile (func perm) + (declare (type function func) + (type permutation perm)) + (let ((args (loop :for i :from 0 :below (permutation-size perm) + :collect (gensym)))) + (compile-and-eval `(lambda (,@args &rest rest) + (apply ,func (append (list ,@(permute! args perm)) rest)))))) + +(defun permute-arguments (func perm) + (declare (type function func) + (type permutation perm)) + (lambda (&rest args) + (apply func (permute! args perm)))) + +(defun curry (func perm &rest curried-args) + (declare (type function func) + (type permutation perm)) + (lambda (&rest args) + (apply func (permute! (append curried-args args) perm)))) + +(defun curry-and-compile (func perm &rest curried-args) + (declare (type function func) + (type permutation perm)) + (let ((args (loop :for i :from 0 :below (permutation-size perm) + :collect (gensym)))) + (compile-and-eval + `(let (,@(mapcar #'(lambda (a b) `(,a ,b)) args curried-args)) + (lambda (,@(nthcdr (length curried-args) args) &rest rest) + (apply ,func (append (list ,@(permute! args perm)) rest))))))) + +(defun compose (func-a func-b perm) + (declare (type function func-a func-b) + (type permutation perm)) + (lambda (&rest args) + (apply func-a (permute! (multiple-value-list (funcall func-b args)) perm)))) + +(defun compose-and-compile (func-a func-b perm) + (declare (type function func-a func-b) + (type permutation perm)) + (let ((syms (loop :for i :from 0 :below (permutation-size perm) + :collect (gensym)))) + (compile-and-eval + `(lambda (&rest args) + (destructuring-bind (,@syms &rest rest) (multiple-value-list (apply ,func-b args)) + (apply ,func-a (append (list ,@(permute! syms perm)) rest))))))) + +;;Back to practical matters. ;;This function is ugly of-course, but is also very very quick! (definline sort-permute (seq predicate &key (key #'matlisp-utilities:id)) " Sorts a lisp-vector in-place, by using the function @arg{predicate} as the - order. Also computes the permutation which would sort the original sequence - @arg{seq}. + order. Also computes the permutation action which would sort the original + sequence @arg{seq} when applied. " (declare (type vector seq)) (let*-typed ((len (length seq) :type fixnum) (perm (pindex-id len) :type pindex-store-vector) - (jobs (list `(0 ,len)))) + (jobs (make-array len :adjustable t :fill-pointer 0))) (loop - :for bounds := (car jobs) :then (pop jobs) + :for bounds := (cons 0 len) :then (unless (zerop (length jobs)) + (vector-pop jobs)) :until (null bounds) - :do (let*-typed ((below-idx (first bounds) :type fixnum) - (above-idx (second bounds) :type fixnum) + :do (let*-typed ((below-idx (car bounds) :type fixnum) + (above-idx (cdr bounds) :type fixnum) (piv (+ below-idx (floor (- above-idx below-idx) 2)) :type fixnum)) (loop :with ele := (funcall key (aref seq piv)) @@ -438,14 +383,14 @@ :until (or (= i piv) (funcall predicate ele (funcall key (aref seq i)))) :finally (setq lbound i)) (loop :for i :of-type fixnum :downfrom ubound :to piv - :until (or (= i piv) (funcall predicate ele (funcall key (aref seq i)))) + :until (or (= i piv) (funcall predicate (funcall key (aref seq i)) ele)) :finally (setq ubound i)) (cond ((= ubound lbound piv) (when (> (- piv below-idx) 1) - (push `(,below-idx ,piv) jobs)) + (vector-push-extend (cons below-idx piv) jobs)) (when (> (- above-idx (1+ piv)) 1) - (push `(,(1+ piv) ,above-idx) jobs)) + (vector-push-extend (cons (1+ piv) above-idx) jobs)) t) ((< lbound piv ubound) (rotatef (aref seq lbound) (aref seq ubound)) diff --git a/src/conditions.lisp b/src/conditions.lisp index 87a526a..862dab0 100644 --- a/src/conditions.lisp +++ b/src/conditions.lisp @@ -122,12 +122,12 @@ (define-condition permutation-permute-error (permutation-error) ((sequence-length :reader seq-len :initarg :seq-len) - (group-rank :reader group-rank :initarg :group-rank)) + (permutation-size :reader per-size :initarg :per-size)) (:documentation "Cannot permute sequence.") - (:report (lambda (c stream) + (:report (lambda (c stream) (format stream "Cannot permute sequence.~%") - (when (slots-boundp c 'sequence-length 'group-rank) - (format stream "~%sequence-length : ~a group-rank: ~a" (seq-len c) (group-rank c)))))) + (when (slots-boundp c 'sequence-length 'permutation-size) + (format stream "~%sequence-length : ~a, permutation size: ~a" (seq-len c) (per-size c)))))) ;;Tensor conditions----------------------------------------------;; (define-condition tensor-error (error) diff --git a/src/utilities/functions.lisp b/src/utilities/functions.lisp index 3c5c40a..aefe273 100644 --- a/src/utilities/functions.lisp +++ b/src/utilities/functions.lisp @@ -6,15 +6,32 @@ (declaim (inline id)) (defun id (x) x) - #+nil - (defun make-ring (n) - (let ((ret (cons nil nil))) - (loop :for i :from 0 :below (1- n) - :with tail :of-type cons := ret - :do (setf (cdr tail) (cons nil nil) - tail (cdr tail)) - :finally (setf (cdr tail) ret)) - ret)) + (declaim (inline vectorify)) + (defun vectorify (seq n &optional (element-type t)) + (declare (type (or vector list) seq)) + (etypecase seq + (cons + (let ((ret (make-array n :element-type element-type))) + (loop :for i :of-type fixnum :from 0 :below n + :for lst := seq :then (cdr lst) + :do (setf (aref ret i) (car lst)) + :finally (return ret)))) + (vector + (let ((ret (make-array n :element-type element-type))) + (loop :for i :of-type fixnum :from 0 :below n + :for ele :across seq + :do (setf (aref ret i) ele) + :finally (return ret)))))) + + (declaim (inline copy-n)) + (defun copy-n (vec lst n) + (declare (type vector vec) + (type list lst) + (type fixnum n)) + (loop :for i :of-type fixnum :from 0 :below n + :for vlst := lst :then (cdr vlst) + :do (setf (car vlst) (aref vec i))) + lst) (declaim (inline slot-values)) (defun slot-values (obj slots) ----------------------------------------------------------------------- Summary of changes: packages.lisp | 3 +- src/base/permutation.lisp | 497 +++++++++++++++++++----------------------- src/conditions.lisp | 8 +- src/utilities/functions.lisp | 35 +++- 4 files changed, 253 insertions(+), 290 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-02-02 18:30:51
|
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, tensor has been updated via 7a8ab5c0db938424bfc8d2ebca6022e2673b6a9a (commit) via 28177ada352e525b622e38a5c7baf97986854288 (commit) from 702e41559e7f9cc9d84d108e467a1ec132e6f50f (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 7a8ab5c0db938424bfc8d2ebca6022e2673b6a9a Author: Akshay Srinivasan <aks...@gm...> Date: Sat Feb 2 10:25:32 2013 -0800 Saving changes in permutation.lisp diff --git a/packages.lisp b/packages.lisp index 6e70c2a..87b3afc 100644 --- a/packages.lisp +++ b/packages.lisp @@ -66,7 +66,7 @@ (defpackage "MATLISP-UTILITIES" (:use #:common-lisp #:matlisp-conditions) - (:export #:ensure-list + (:export #:ensure-list #:id ;;#:make-ring #:zip #:zip-eq #:cut-cons-chain! #:slot-values diff --git a/src/base/permutation.lisp b/src/base/permutation.lisp index 9439689..9f5bbda 100644 --- a/src/base/permutation.lisp +++ b/src/base/permutation.lisp @@ -1,35 +1,33 @@ (in-package #:matlisp) ;;This must match the type used in LAPACK -(deftype perrepr-type () +(deftype pindex-type () '(unsigned-byte 32)) -(deftype perrepr-vector (&optional (size '*)) - `(simple-array perrepr-type (,size))) +(deftype pindex-store-vector (&optional (size '*)) + `(simple-array pindex-type (,size))) -(make-array-allocator allocate-perrepr-store 'perrepr-type 0 - " +(make-array-allocator allocate-pindex-store 'pindex-type 0 + " Syntax ====== - (ALLOCATE-PERREPR-STORE SIZE [INITIAL-ELEMENT 0]) + (ALLOCATE-PINDEX-STORE SIZE [INITIAL-ELEMENT 0]) Purpose ======= Allocates integer4 (32-bits) storage.") ;; -(defun perrepr-id-action (n) +(definline pindex-id (n) (declare (type fixnum n)) - (let ((ret (allocate-perrepr-store n))) - (declare (type perrepr-vector ret)) - (very-quickly - (loop - :for i :of-type perrepr-type :from 0 :below n - :do (setf (aref ret i) i))) - ret)) + (let-typed ((ret (allocate-pindex-store n) :type pindex-store-vector)) + (very-quickly + (loop :for i :of-type pindex-type :from 0 :below n + :do (setf (aref ret i) i))) + ret)) -(definline perv (&rest contents) - (make-array (length contents) :element-type 'perrepr-type :initial-contents contents)) +(definline pidxv (&rest contents) + (make-array (length contents) :element-type 'pindex-type :initial-contents contents)) ;;Write a uniform randomiser (defun seqrnd (seq) @@ -39,141 +37,77 @@ ;;Class definitions----------------------------------------------;; (defclass permutation () - ((representation :accessor repr - :initarg :repr) - (group-rank :accessor group-rank - :type index-type))) + ((store :accessor store :initarg :store) + (permutation-size :accessor permutation-size :type index-type))) (defmethod print-object ((per permutation) stream) (print-unreadable-object (per stream :type t) - (format stream "S_~a~%" (group-rank per)) - (if (zerop (group-rank per)) + (format stream "S_~a~%" (permutation-size per)) + (if (= (permutation-size per) 1) (format stream "ID~%") - (format stream "~a~%" (repr per))))) -;; + (format stream "~a~%" (store per))))) +;; (defclass permutation-action (permutation) - ((representation :type perrepr-vector))) + ((store :type pindex-store-vector))) - -(defmethod initialize-instance :after ((per permutation-action) &rest initargs) +(defmethod initialize-instance :after ((perm permutation-action) &rest initargs) (declare (ignore initargs)) - (labels ((action-repr-p (act) - " - Checks if ARR is a possible permutation vector. A permutation pi - is characterized by a vector containing the indices from 0,..., - @function{length}(@arg{perm})-1 in some order. -" - (if (not (typep act 'perrepr-vector)) nil - (locally - (declare (type perrepr-vector act)) - (let* ((len (length act)) - (sort (very-quickly (sort (copy-seq act) #'<)))) - (declare (type perrepr-vector sort) - (type index-type len)) - (very-quickly - (loop :for i :of-type index-type :from 0 :below len - :unless (= (aref sort i) i) - :do (return nil) - :finally (return t)))))))) - (assert (action-repr-p (repr per)) nil 'permutation-invalid-error) - (setf (group-rank per) (length (repr per))))) -;; + (when *check-after-initializing?* + (let-typed ((repr (store perm) :type pindex-store-vector)) + (very-quickly + (loop :for i :of-type index-type :from 0 :below (length repr) + :with srepr :of-type pindex-store-vector := (sort (copy-seq repr) #'<) + :do (assert (= (aref srepr i) i) nil 'permutation-invalid-error))) + (setf (permutation-size perm) (length repr))))) +;; (defclass permutation-cycle (permutation) - ((representation :type list))) + ((store :type list))) (defmethod initialize-instance :after ((per permutation-cycle) &rest initargs) (declare (ignore initargs)) - (labels ((cycle-repr-p (perm) - " - Does a sorting operation to check for duplicate elements in - the cycle representation of a permutation. -" - (if (not (typep perm 'perrepr-vector)) nil - (locally - (declare (type perrepr-vector perm)) - (let ((len (length perm))) - (declare (type index-type len)) - (if (<= len 1) nil - (let ((sort (very-quickly (sort (copy-seq perm) #'<)))) - (declare (type perrepr-vector sort)) - (very-quickly - (loop :for i :of-type index-type :from 1 :below len - :when (= (aref sort i) (aref sort (1- i))) - :do (return nil) - :finally (return t)))))))))) - (if (null (repr per)) (setf (group-rank per) 0) + (when *check-after-initializing?* + (if (null (store per)) + (setf (permutation-size per) 1) (very-quickly (loop - :for cyc :of-type perrepr-vector :in (repr per) - :unless (cycle-repr-p cyc) - :do (error 'permutation-invalid-error) - :maximizing (lvec-max cyc) into g-rnk of-type index-type - :finally (setf (group-rank per) (the index-type (1+ g-rnk)))))))) + :for cyc :of-type pindex-store-vector :in (store per) + :with ss :of-type pindex-type := 0 + :do (loop + :for i :of-type index-type :from 1 :below (length cyc) + :with scyc :of-type pindex-store-vector := (sort (copy-seq cyc) #'<) + :do (assert (/= (aref scyc (1- i)) (aref scyc i)) nil 'permutation-invalid-error) + :finally (setf ss (max ss (aref scyc (1- (length scyc)))))) + :finally (setf (permutation-size per) (1+ ss))))))) ;; (defclass permutation-pivot-flip (permutation) - ((representation :type perrepr-vector))) + ((store :type pindex-store-vector))) (defmethod initialize-instance :after ((per permutation-pivot-flip) &rest initargs) (declare (ignore initargs)) - (labels ((pivot-flip-p (idiv) - " - Checks if ARR is a possible permutation pivot-flip. A pivot - flip is more algorithmic in its representation. If a sequence - is given, apply a pivot-flip on it is equivalent to running from - the left to the right of the permutation by flipping (pi(i), i) - sequentially. This is the representation used in LAPACK. -" - (if (not (typep idiv 'perrepr-vector)) nil - (let ((len (length idiv))) - (declare (type perrepr-vector idiv) - (type index-type len)) - (very-quickly - (loop :for i :of-type index-type :from 0 :below len - :do (unless (< -1 (aref idiv i) len) - (return nil)) - :finally (return t))))))) - (assert (pivot-flip-p (repr per)) nil 'permutation-invalid-error) - (setf (group-rank per) (length (repr per))))) - -;; -;; (defclass permutation-matrix (permutation) -;; ((representation :type perrepr-vector))) - -;; (defmethod initialize-instance :after ((per permutation-pivot-flip) &rest initargs) -;; (declare (ignore initargs)) -;; (labels ((pivot-flip-p (idiv) -;; " -;; Checks if ARR is a possible permutation pivot-flip. A pivot -;; flip is more algorithmic in its representation. If a sequence -;; is given, apply a pivot-flip on it is equivalent to running from -;; the left to the right of the permutation by flipping (pi(i), i) -;; sequentially. This is the representation used in LAPACK. -;; " -;; (if (not (typep idiv 'perrepr-vector)) nil -;; (let ((len (length idiv))) -;; (declare (type perrepr-vector idiv) -;; (type index-type len)) -;; (very-quickly -;; (loop :for i :of-type index-type :from 0 :below len -;; :do (unless (< -1 (aref idiv i) len) -;; (return nil)) -;; :finally (return t))))))) -;; (assert (pivot-flip-p (repr per)) nil 'permutation-invalid-error) -;; (setf (group-rank per) (length (repr per))))) - + (when *check-after-initializing?* + (let*-typed ((repr (store per) :type pindex-store-vector) + (len (length repr) :type index-type)) + (very-quickly + (loop :for i :of-type index-type :from 0 :below len + :do (assert (< -1 (aref repr i) len) nil 'permutation-invalid-error))) + (setf (store-size per) len)))) ;; -(definline make-pcycle (&rest args) - (make-instance 'permutation-cycle :repr args)) +(defclass permutation-matrix (permutation) + ((store :type pindex-store-vector))) -(definline make-paction (pact) - (make-instance 'permutation-action :repr pact)) - -(definline make-pidx (pact) - (make-instance 'permutation-pivot-flip :repr pact)) +(defmethod initialize-instance :after ((perm permutation-matrix) &rest initargs) + (declare (ignore initargs)) + (when *check-after-initializing?* + (let-typed ((repr (store perm) :type pindex-store-vector)) + (very-quickly + (loop :for i :of-type index-type :from 0 :below (length repr) + :with srepr :of-type pindex-store-vector := (sort (copy-seq repr) #'<) + :do (assert (= (aref srepr i) i) nil 'permutation-invalid-error))) + (setf (permutation-size perm) (length repr))))) ;;Generic permute! method. (defgeneric permute! (thing permutation &optional argument) @@ -195,311 +129,323 @@ (definline permute (thing perm &optional (arg 0)) (permute! (copy thing) perm arg)) -;;Action -(defmethod permute! ((seq cons) (perm permutation-action) &optional arg) - (declare (ignore arg)) - (let ((cseq (make-array (length seq) :initial-contents seq)) - (act (repr perm)) - (glen (group-rank perm))) - (mapl - (let ((i 0)) - (declare (type fixnum i)) - (lambda (x) - (when (< i glen) - (rplaca x (aref cseq (aref act i))) - (incf i)))) seq))) - -(defmethod permute! ((seq vector) (perm permutation-action) &optional arg) - (declare (ignore arg)) - (let ((cseq (make-array (length seq) :initial-contents seq)) - (act (repr perm))) - (loop - :for i :from 0 :below (group-rank perm) - :do (unless (= i (aref act i)) - (setf (aref seq i) (aref cseq (aref act i)))) - :finally (return seq)))) - -(defmethod permute! ((ten standard-tensor) (perm permutation-action) &optional (arg 0)) - (let ((cyc (action->cycle perm))) - (permute! ten cyc arg))) - -;;Cycle -;;Might be useful ? -(defun apply-cycle! (seq pcyc) - (declare (type perrepr-vector pcyc) - (type vector seq)) - (let ((xl (aref seq (aref pcyc (1- (length pcyc)))))) - (loop :for i :of-type index-type :downfrom (1- (length pcyc)) :to 0 - :do (setf (aref seq (aref pcyc i)) - (if (= i 0) xl - (aref seq (aref pcyc (1- i)))))))) - -(defmethod permute! ((seq cons) (perm permutation-cycle) &optional arg) - (declare (ignore arg)) - (let ((cseq (make-array (length seq) :initial-contents seq)) - (glen (group-rank perm))) - (dolist (cyc (repr perm)) - (declare (type perrepr-vector cyc)) - (apply-cycle! cseq cyc)) - (mapl - (let ((i 0)) - (lambda (x) - (when (< i glen) - (rplaca x (aref cseq i)) - (incf i)))) seq))) - -(defmethod permute! ((seq vector) (perm permutation-cycle) &optional arg) - (declare (ignore arg)) - (dolist (cyc (repr perm) seq) - (declare (type perrepr-vector cyc)) - (apply-cycle! seq cyc))) - -(defmethod permute! ((A standard-tensor) (perm permutation-cycle) &optional (arg 0)) - (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) - (rplaca (nthcdr arg slst) 0) - (values (sub-tensor~ A slst) (sub-tensor~ A slst))) - (let-typed ((cyclst (repr perm) :type cons) - (cp-ten (make-instance (class-of tone) - :dimensions (copy-seq (dimensions tone)))) - (std-arg (aref (strides A) arg) :type index-type) - (hd-sl (head ttwo) :type index-type)) - (dolist (cyc cyclst) - (declare (type perrepr-vector cyc)) - (setf (head tone) (+ hd-sl (* std-arg (aref cyc (1- (length cyc)))))) - (copy! tone cp-ten) - (loop :for i :of-type index-type :downfrom (1- (length cyc)) :to 0 - :do (progn - (setf (head tone) (+ hd-sl (* std-arg (aref cyc i)))) - (copy! - (if (= i 0) cp-ten - (progn - (setf (head ttwo) (+ hd-sl (* std-arg (aref cyc (1- i))))) - ttwo)) - tone)))))) - A) - -;;Pivot idx -(defmethod permute! ((seq vector) (perm permutation-pivot-flip) &optional arg) - (declare (ignore arg)) - (let-typed ((pidx (repr perm) :type perrepr-vector)) - (loop :for i :of-type index-type :from 0 :below (group-rank perm) - :unless (= i (aref pidx i)) - :do (rotatef (aref seq i) (aref seq (aref pidx i))) - :finally (return seq)))) - -(defmethod permute! ((seq cons) (perm permutation-pivot-flip) &optional arg) - (declare (ignore arg)) - (let ((cseq (make-array (length seq) :initial-contents seq)) - (glen (group-rank perm))) - (permute! cseq perm) - (mapl - (let ((i 0)) - (lambda (x) - (when (< i glen) - (rplaca x (aref cseq i)) - (incf i)))) seq))) - -(defmethod permute! ((A standard-tensor) (perm permutation-pivot-flip) &optional (arg 0)) - (let ((idiv (repr perm))) - (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) - (rplaca (nthcdr arg slst) 0) - (values (sub-tensor~ A slst nil) (sub-tensor~ A slst nil))) - (let ((argstd (aref (strides A) arg)) - (hd-sl (head ttwo))) - (loop :for i :from 0 :below (length idiv) - :do (progn - (unless (= i (aref idiv i)) - (setf (head ttwo) (+ hd-sl (* (aref idiv i) argstd))) - (swap! tone ttwo)) - (incf (head tone) argstd)))))) - A) - -;;Conversions----------------------------------------------------;; -(defun action->cycle (act) - " - (action->cycle act) - - This function obtains the canonical cycle representation - of a permutation. The first argument \"act\" is the action of the - permutation on the array #(0 1 2 3 ..): an object of the class - permutation-action. - - \"Canonical\" may be a bit of an overstatement; this is the way - S_n was presented in Van der Waerden's book. -" - (declare (type permutation-action act)) - (mlet* - ((arr (repr act) :type perrepr-vector)) - (labels ((find-cycle (x0) - ;; This function obtains the cycle starting from x_0. - (declare (type perrepr-type x0)) - (if (= (aref arr x0) x0) (values 0 nil) - (very-quickly - (loop - :for x :of-type perrepr-type := (aref arr x0) :then (aref arr x) - :and ret :of-type cons := (list x0) :then (cons x ret) - :counting t :into i :of-type index-type - :when (= x x0) - :do (return (values i ret)))))) - (cycle-walk (cyc ignore) - ;; Finds all cycles - (let ((x0 (find-if-not #'(lambda (x) (member x ignore)) arr))) - (if (null x0) - cyc - (multiple-value-bind (clen clst) (find-cycle x0) - (declare (type index-type clen) - (type list clst)) - (cycle-walk - (if (= clen 0) cyc - (cons (make-array clen :element-type 'perrepr-type :initial-contents clst) cyc)) - (nconc ignore (if (= clen 0) (list x0) clst)))))))) - (let ((cyc-lst (cycle-walk nil nil))) - (make-instance 'permutation-cycle - :repr cyc-lst))))) - -(defun cycle->action (cyc) - " - (cycle->action cyc) - - This function obtains the action representation of a permutation - from the cyclic one. The first argument \"cyc\" is the cyclic - representation of the permutation: an object of the class - permutation-cycle. -" - (declare (type permutation-cycle cyc)) - (let ((act-repr (perrepr-id-action (group-rank cyc))) - (cycs-repr (repr cyc))) - (declare (type perrepr-vector act-repr)) - (dolist (cyc cycs-repr) - (declare (type perrepr-vector cyc)) - (let ((xl (aref act-repr (aref cyc (1- (length cyc)))))) - (very-quickly - (loop - :for i :of-type index-type :downfrom (1- (length cyc)) :to 0 - :do (setf (aref act-repr (aref cyc i)) - (if (= i 0) xl - (aref act-repr (aref cyc (1- i))))))))) - (make-instance 'permutation-action :repr act-repr))) - -(defun pivot-flip->action (pflip) - (declare (type permutation-pivot-flip pflip)) - (let* ((idiv (repr pflip)) - (len (length idiv))) - (declare (type perrepr-vector idiv) - (type index-type len)) - (let ((act (perrepr-id-action len))) - (declare (type perrepr-vector act)) - (very-quickly - (loop :for i :from 0 :below len - :do (let ((val (aref idiv i))) - (unless (= val i) - (rotatef (aref act i) (aref act val)))))) - (make-instance 'permutation-action :repr act)))) - -(defun mod-max (seq lidx uidx) - (declare (type perrepr-vector seq)) - (let ((len (length seq))) - (very-quickly - (loop :for idx :of-type index-type :downfrom uidx :above lidx - :with max :of-type perrepr-type := (aref seq uidx) - :with max-idx :of-type index-type := uidx - :do (let ((ele (aref seq (mod idx len)))) - (when (> ele max) - (setf max ele - max-idx idx))) - :finally (return (values max max-idx)))))) - -#| - -(defun cycle->pivot-flip (cyc) - (let ((cp-cyc (copy-seq cyc))) - (let - (labels ((mod-max (seq lidx uidx) - (declare (type perrepr-vector seq)) - (let ((len (length cyc))) - (very-quickly - (loop :for idx :of-type index-type :downfrom uidx :above lidx - :with max :of-type perrepr-type := (aref seq uidx) - :with max-idx :of-type index-type := uidx - :do (let ((ele (aref seq (mod idx len)))) - (when (> ele max) - (setf max ele - max-idx idx))) - :finally (return (values max max-idx)))))) - (get-flip (lidx uidx) - (multiple-value-bind (max max-idx) (mod-max cyc lidx uidx) - (let ((ele-0 (aref cyc (mod max-idx len))) - (ele-1 (aref cyc (mod (1- max-idx) len)))) - (setf (aref pidx (max ele-0 ele-1)) - (min ele-0 ele-1)) - |# - -#+nil -(defun permute-argument (func-symbol perm) - (declare (type symbol func-symbol) - (type permutation perm)) - (let* ((glen (group-rank perm)) - (args (loop :for i :from 0 :below glen - :collect (gensym)))) - (eval `(lambda (,@args &rest rest) - (apply ',func-symbol (append (list ,@(permute! args perm)) rest)))))) - -(defun argument-permute (func perm) - (declare (type function func) - (type permutation perm)) - (lambda (&rest args) - (apply func (permute! args perm)))) - -(defun curry (func perm &rest curried-args) - (declare (type function func) - (type permutation perm)) - (lambda (&rest args) - (apply func (permute! (append curried-args args) perm)))) - -(defun compose (func-a func-b perm) - (declare (type function func-a func-b) - (type permutation perm)) - (lambda (&rest args) - (apply func-a (permute! (multiple-value-list (funcall func-b args)) perm)))) -;; - -(definline sort-permute (seq predicate) +;; ;;Action +;; (defmethod permute! ((seq cons) (perm permutation-action) &optional arg) +;; (declare (ignore arg)) +;; (let ((cseq (make-array (length seq) :initial-contents seq)) +;; (act (repr perm)) +;; (glen (group-rank perm))) +;; (mapl +;; (let ((i 0)) +;; (declare (type fixnum i)) +;; (lambda (x) +;; (when (< i glen) +;; (rplaca x (aref cseq (aref act i))) +;; (incf i)))) seq))) + +;; (defmethod permute! ((seq vector) (perm permutation-action) &optional arg) +;; (declare (ignore arg)) +;; (let ((cseq (make-array (length seq) :initial-contents seq)) +;; (act (repr perm))) +;; (loop +;; :for i :from 0 :below (group-rank perm) +;; :do (unless (= i (aref act i)) +;; (setf (aref seq i) (aref cseq (aref act i)))) +;; :finally (return seq)))) + +;; (defmethod permute! ((ten standard-tensor) (perm permutation-action) &optional (arg 0)) +;; (let ((cyc (action->cycle perm))) +;; (permute! ten cyc arg))) + +;; ;;Cycle +;; ;;Might be useful ? +;; (defun apply-cycle! (seq pcyc) +;; (declare (type pindex-store-vector pcyc) +;; (type vector seq)) +;; (let ((xl (aref seq (aref pcyc (1- (length pcyc)))))) +;; (loop :for i :of-type index-type :downfrom (1- (length pcyc)) :to 0 +;; :do (setf (aref seq (aref pcyc i)) +;; (if (= i 0) xl +;; (aref seq (aref pcyc (1- i)))))))) + +;; (defmethod permute! ((seq cons) (perm permutation-cycle) &optional arg) +;; (declare (ignore arg)) +;; (let ((cseq (make-array (length seq) :initial-contents seq)) +;; (glen (group-rank perm))) +;; (dolist (cyc (repr perm)) +;; (declare (type pindex-store-vector cyc)) +;; (apply-cycle! cseq cyc)) +;; (mapl +;; (let ((i 0)) +;; (lambda (x) +;; (when (< i glen) +;; (rplaca x (aref cseq i)) +;; (incf i)))) seq))) + +;; (defmethod permute! ((seq vector) (perm permutation-cycle) &optional arg) +;; (declare (ignore arg)) +;; (dolist (cyc (repr perm) seq) +;; (declare (type pindex-store-vector cyc)) +;; (apply-cycle! seq cyc))) + +;; (defmethod permute! ((A standard-tensor) (perm permutation-cycle) &optional (arg 0)) +;; (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) +;; (rplaca (nthcdr arg slst) 0) +;; (values (sub-tensor~ A slst) (sub-tensor~ A slst))) +;; (let-typed ((cyclst (repr perm) :type cons) +;; (cp-ten (make-instance (class-of tone) +;; :dimensions (copy-seq (dimensions tone)))) +;; (std-arg (aref (strides A) arg) :type index-type) +;; (hd-sl (head ttwo) :type index-type)) +;; (dolist (cyc cyclst) +;; (declare (type pindex-store-vector cyc)) +;; (setf (head tone) (+ hd-sl (* std-arg (aref cyc (1- (length cyc)))))) +;; (copy! tone cp-ten) +;; (loop :for i :of-type index-type :downfrom (1- (length cyc)) :to 0 +;; :do (progn +;; (setf (head tone) (+ hd-sl (* std-arg (aref cyc i)))) +;; (copy! +;; (if (= i 0) cp-ten +;; (progn +;; (setf (head ttwo) (+ hd-sl (* std-arg (aref cyc (1- i))))) +;; ttwo)) +;; tone)))))) +;; A) + +;; ;;Pivot idx +;; (defmethod permute! ((seq vector) (perm permutation-pivot-flip) &optional arg) +;; (declare (ignore arg)) +;; (let-typed ((pidx (repr perm) :type pindex-store-vector)) +;; (loop :for i :of-type index-type :from 0 :below (group-rank perm) +;; :unless (= i (aref pidx i)) +;; :do (rotatef (aref seq i) (aref seq (aref pidx i))) +;; :finally (return seq)))) + +;; (defmethod permute! ((seq cons) (perm permutation-pivot-flip) &optional arg) +;; (declare (ignore arg)) +;; (let ((cseq (make-array (length seq) :initial-contents seq)) +;; (glen (group-rank perm))) +;; (permute! cseq perm) +;; (mapl +;; (let ((i 0)) +;; (lambda (x) +;; (when (< i glen) +;; (rplaca x (aref cseq i)) +;; (incf i)))) seq))) + +;; (defmethod permute! ((A standard-tensor) (perm permutation-pivot-flip) &optional (arg 0)) +;; (let ((idiv (repr perm))) +;; (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) +;; (rplaca (nthcdr arg slst) 0) +;; (values (sub-tensor~ A slst nil) (sub-tensor~ A slst nil))) +;; (let ((argstd (aref (strides A) arg)) +;; (hd-sl (head ttwo))) +;; (loop :for i :from 0 :below (length idiv) +;; :do (progn +;; (unless (= i (aref idiv i)) +;; (setf (head ttwo) (+ hd-sl (* (aref idiv i) argstd))) +;; (swap! tone ttwo)) +;; (incf (head tone) argstd)))))) +;; A) + +;; ;;Conversions----------------------------------------------------;; +;; (defun action->cycle (act) +;; " +;; (action->cycle act) + +;; This function obtains the canonical cycle representation +;; of a permutation. The first argument \"act\" is the action of the +;; permutation on the array #(0 1 2 3 ..): an object of the class +;; permutation-action. + +;; \"Canonical\" may be a bit of an overstatement; this is the way +;; S_n was presented in Van der Waerden's book. +;; " +;; (declare (type permutation-action act)) +;; (mlet* +;; ((arr (repr act) :type pindex-store-vector)) +;; (labels ((find-cycle (x0) +;; ;; This function obtains the cycle starting from x_0. +;; (declare (type pindex-type x0)) +;; (if (= (aref arr x0) x0) (values 0 nil) +;; (very-quickly +;; (loop +;; :for x :of-type pindex-type := (aref arr x0) :then (aref arr x) +;; :and ret :of-type cons := (list x0) :then (cons x ret) +;; :counting t :into i :of-type index-type +;; :when (= x x0) +;; :do (return (values i ret)))))) +;; (cycle-walk (cyc ignore) +;; ;; Finds all cycles +;; (let ((x0 (find-if-not #'(lambda (x) (member x ignore)) arr))) +;; (if (null x0) +;; cyc +;; (multiple-value-bind (clen clst) (find-cycle x0) +;; (declare (type index-type clen) +;; (type list clst)) +;; (cycle-walk +;; (if (= clen 0) cyc +;; (cons (make-array clen :element-type 'pindex-type :initial-contents clst) cyc)) +;; (nconc ignore (if (= clen 0) (list x0) clst)))))))) +;; (let ((cyc-lst (cycle-walk nil nil))) +;; (make-instance 'permutation-cycle +;; :repr cyc-lst))))) + +;; (defun cycle->action (cyc) +;; " +;; (cycle->action cyc) + +;; This function obtains the action representation of a permutation +;; from the cyclic one. The first argument \"cyc\" is the cyclic +;; representation of the permutation: an object of the class +;; permutation-cycle. +;; " +;; (declare (type permutation-cycle cyc)) +;; (let ((act-repr (pindex-id-action (group-rank cyc))) +;; (cycs-repr (repr cyc))) +;; (declare (type pindex-store-vector act-repr)) +;; (dolist (cyc cycs-repr) +;; (declare (type pindex-store-vector cyc)) +;; (let ((xl (aref act-repr (aref cyc (1- (length cyc)))))) +;; (very-quickly +;; (loop +;; :for i :of-type index-type :downfrom (1- (length cyc)) :to 0 +;; :do (setf (aref act-repr (aref cyc i)) +;; (if (= i 0) xl +;; (aref act-repr (aref cyc (1- i))))))))) +;; (make-instance 'permutation-action :repr act-repr))) + +;; (defun pivot-flip->action (pflip) +;; (declare (type permutation-pivot-flip pflip)) +;; (let* ((idiv (repr pflip)) +;; (len (length idiv))) +;; (declare (type pindex-store-vector idiv) +;; (type index-type len)) +;; (let ((act (pindex-id-action len))) +;; (declare (type pindex-store-vector act)) +;; (very-quickly +;; (loop :for i :from 0 :below len +;; :do (let ((val (aref idiv i))) +;; (unless (= val i) +;; (rotatef (aref act i) (aref act val)))))) +;; (make-instance 'permutation-action :repr act)))) + +;; (defun mod-max (seq lidx uidx) +;; (declare (type pindex-store-vector seq)) +;; (let ((len (length seq))) +;; (very-quickly +;; (loop :for idx :of-type index-type :downfrom uidx :above lidx +;; :with max :of-type pindex-type := (aref seq uidx) +;; :with max-idx :of-type index-type := uidx +;; :do (let ((ele (aref seq (mod idx len)))) +;; (when (> ele max) +;; (setf max ele +;; max-idx idx))) +;; :finally (return (values max max-idx)))))) + +;; #| + +;; (defun cycle->pivot-flip (cyc) +;; (let ((cp-cyc (copy-seq cyc))) +;; (let +;; (labels ((mod-max (seq lidx uidx) +;; (declare (type pindex-store-vector seq)) +;; (let ((len (length cyc))) +;; (very-quickly +;; (loop :for idx :of-type index-type :downfrom uidx :above lidx +;; :with max :of-type pindex-type := (aref seq uidx) +;; :with max-idx :of-type index-type := uidx +;; :do (let ((ele (aref seq (mod idx len)))) +;; (when (> ele max) +;; (setf max ele +;; max-idx idx))) +;; :finally (return (values max max-idx)))))) +;; (get-flip (lidx uidx) +;; (multiple-value-bind (max max-idx) (mod-max cyc lidx uidx) +;; (let ((ele-0 (aref cyc (mod max-idx len))) +;; (ele-1 (aref cyc (mod (1- max-idx) len)))) +;; (setf (aref pidx (max ele-0 ele-1)) +;; (min ele-0 ele-1)) +;; |# + +;; #+nil +;; (defun permute-argument (func-symbol perm) +;; (declare (type symbol func-symbol) +;; (type permutation perm)) +;; (let* ((glen (group-rank perm)) +;; (args (loop :for i :from 0 :below glen +;; :collect (gensym)))) +;; (eval `(lambda (,@args &rest rest) +;; (apply ',func-symbol (append (list ,@(permute! args perm)) rest)))))) + +;; (defun argument-permute (func perm) +;; (declare (type function func) +;; (type permutation perm)) +;; (lambda (&rest args) +;; (apply func (permute! args perm)))) + +;; (defun curry (func perm &rest curried-args) +;; (declare (type function func) +;; (type permutation perm)) +;; (lambda (&rest args) +;; (apply func (permute! (append curried-args args) perm)))) + +;; (defun compose (func-a func-b perm) +;; (declare (type function func-a func-b) +;; (type permutation perm)) +;; (lambda (&rest args) +;; (apply func-a (permute! (multiple-value-list (funcall func-b args)) perm)))) +;; ;; + + +;; (defstruct (ring (:constructor nil)) +;; (circular-list nil :type list) +;; (end-list nil :type list)) + +;; (defun make-ring (n) +;; (let ((ret (cons nil nil))) +;; (loop :for i :from 0 :below (1- n) +;; :with tail :of-type cons := ret +;; :do (setf (cdr tail) (cons nil nil) +;; tail (cdr tail)) +;; :finally (setf (cdr tail) ret)) +;; ret)) + + +;;This function is ugly of-course, but is also very very quick! +(definline sort-permute (seq predicate &key (key #'matlisp-utilities:id)) " Sorts a lisp-vector in-place, by using the function @arg{predicate} as the order. Also computes the permutation which would sort the original sequence @arg{seq}. " (declare (type vector seq)) - ;;This function is ugly of-course, but is also very very quick! (let*-typed ((len (length seq) :type fixnum) - (perm (perrepr-id-action len) :type perrepr-vector) + (perm (pindex-id len) :type pindex-store-vector) (jobs (list `(0 ,len)))) - (loop ;;:repeat 10 - :for bounds := (pop jobs) :then (pop jobs) + (loop + :for bounds := (car jobs) :then (pop jobs) :until (null bounds) :do (let*-typed ((below-idx (first bounds) :type fixnum) (above-idx (second bounds) :type fixnum) (piv (+ below-idx (floor (- above-idx below-idx) 2)) :type fixnum)) - (loop ;;:repeat 10 - :with ele := (aref seq piv) + (loop + :with ele := (funcall key (aref seq piv)) :with lbound :of-type fixnum := below-idx :with ubound :of-type fixnum := (1- above-idx) :until (progn - ;;(format t "~%~a ~%" (list lbound piv ubound)) (loop :for i :of-type fixnum :from lbound :to piv - :until (or (= i piv) (funcall predicate ele (aref seq i))) + :until (or (= i piv) (funcall predicate ele (funcall key (aref seq i)))) :finally (setq lbound i)) (loop :for i :of-type fixnum :downfrom ubound :to piv - :until (or (= i piv) (funcall predicate (aref seq i) ele)) + :until (or (= i piv) (funcall predicate ele (funcall key (aref seq i)))) :finally (setq ubound i)) - ;;(format t "~a ~%" (list lbound piv ubound)) (cond ((= ubound lbound piv) (when (> (- piv below-idx) 1) (push `(,below-idx ,piv) jobs)) (when (> (- above-idx (1+ piv)) 1) (push `(,(1+ piv) ,above-idx) jobs)) - ;;(format t "~a~%" jobs) t) ((< lbound piv ubound) (rotatef (aref seq lbound) (aref seq ubound)) @@ -525,4 +471,4 @@ (decf piv) (decf ubound) nil))))) - :finally (return (values seq (make-instance 'permutation-action :repr perm)))))) + :finally (return (values seq (make-instance 'permutation-action :store perm)))))) diff --git a/src/utilities/functions.lisp b/src/utilities/functions.lisp index 079c6ba..3c5c40a 100644 --- a/src/utilities/functions.lisp +++ b/src/utilities/functions.lisp @@ -2,6 +2,19 @@ ;;These functions are used all over the place inside Matlisp's macros. (eval-when (:compile-toplevel :load-toplevel :execute) + + (declaim (inline id)) + (defun id (x) x) + + #+nil + (defun make-ring (n) + (let ((ret (cons nil nil))) + (loop :for i :from 0 :below (1- n) + :with tail :of-type cons := ret + :do (setf (cdr tail) (cons nil nil) + tail (cdr tail)) + :finally (setf (cdr tail) ret)) + ret)) (declaim (inline slot-values)) (defun slot-values (obj slots) diff --git a/src/utilities/macros.lisp b/src/utilities/macros.lisp index bad1160..259fb95 100644 --- a/src/utilities/macros.lisp +++ b/src/utilities/macros.lisp @@ -63,11 +63,12 @@ ,@body)))) (defmacro make-array-allocator (allocator-name type init &optional doc) - `(definline ,allocator-name (size &optional (initial-element ,init)) - ,@(unless (null doc) - `(,doc)) - (make-array size - :element-type ,type :initial-element initial-element))) + `(eval-when (:compile-toplevel :load-toplevel :execute) + (definline ,allocator-name (size &optional (initial-element ,init)) + ,@(unless (null doc) + `(,doc)) + (make-array size + :element-type ,type :initial-element initial-element)))) (defmacro let-typed (bindings &rest body) " commit 28177ada352e525b622e38a5c7baf97986854288 Author: Akshay Srinivasan <aks...@gm...> Date: Fri Feb 1 22:40:07 2013 -0800 Cosmetic changes to permutation.lisp, dlsode.lisp. diff --git a/src/base/permutation.lisp b/src/base/permutation.lisp index 0a69941..9439689 100644 --- a/src/base/permutation.lisp +++ b/src/base/permutation.lisp @@ -477,7 +477,6 @@ (loop ;;:repeat 10 :for bounds := (pop jobs) :then (pop jobs) :until (null bounds) - :finally (return (values seq (make-instance 'permutation-action :repr perm))) :do (let*-typed ((below-idx (first bounds) :type fixnum) (above-idx (second bounds) :type fixnum) (piv (+ below-idx (floor (- above-idx below-idx) 2)) :type fixnum)) @@ -525,4 +524,5 @@ (rotatef (aref perm lbound) (aref perm piv))) (decf piv) (decf ubound) - nil)))))))) + nil))))) + :finally (return (values seq (make-instance 'permutation-action :repr perm)))))) diff --git a/src/packages/odepack/dlsode.lisp b/src/packages/odepack/dlsode.lisp index a18b957..dab15bb 100644 --- a/src/packages/odepack/dlsode.lisp +++ b/src/packages/odepack/dlsode.lisp @@ -1,24 +1,11 @@ (in-package #:matlisp) (cffi:define-foreign-library libodepack - #+nil(:unix #.(translate-logical-pathname - (merge-pathnames "matlisp:lib;libodepack" - *shared-library-pathname-extension*))) (t (:default "libodepack"))) (cffi:use-foreign-library libodepack) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -#+nil -(def-fortran-routine testde :void - (field (:callback :void - (c-neq :integer :input) - (c-t :double-float :input) - (c-y (* :double-float :size c-neq) :input) - (c-ydot (* :double-float :size c-neq) :output))) - (neq :integer :input) - (y (* :double-float) :input-output)) - (def-fortran-routine dlsode :void "DLSODE in ODEPACK" (field (:callback dlsode-field-callback ----------------------------------------------------------------------- Summary of changes: packages.lisp | 2 +- src/base/permutation.lisp | 764 ++++++++++++++++++-------------------- src/packages/odepack/dlsode.lisp | 13 - src/utilities/functions.lisp | 13 + src/utilities/macros.lisp | 11 +- 5 files changed, 375 insertions(+), 428 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-01-31 01:53:52
|
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, tensor has been updated via 702e41559e7f9cc9d84d108e467a1ec132e6f50f (commit) from 13671606829449440d65e08dd7ec8b54560c8f23 (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 702e41559e7f9cc9d84d108e467a1ec132e6f50f Author: Akshay Srinivasan <aks...@gm...> Date: Wed Jan 30 17:44:04 2013 -0800 Removed comments from mm.c. diff --git a/tests/mm.c b/tests/mm.c index d7a8bb3..4fa56c7 100644 --- a/tests/mm.c +++ b/tests/mm.c @@ -18,9 +18,6 @@ int main(int argc, char *argv[]){ a = calloc(n * n, sizeof(double)); b = calloc(n * n, sizeof(double)); c = calloc(n * n, sizeof(double));; - /* double a[] = {1, 2, 3, 4}, */ - /* b[] = {4, 5, 6, 7}, */ - /* c[] = {0, 0, 0, 0}; */ int i, j, k, l; printf("N: %d\n", n); @@ -35,38 +32,7 @@ int main(int argc, char *argv[]){ int of_a = 0, of_b = 0, of_c = 0; clock_t co, ci; - ci = clock(); - - /* for(i = 0; i < n; i++){ */ - /* for(j = 0; j < n; j++){ */ - /* tmp = 0.; */ - /* for(k = 0; k < n; k++){ */ - /* /\* if((of_c > n * n) || (of_c < 0) || \ *\/ */ - /* /\* (of_b > n * n) || (of_b < 0) || \ *\/ */ - /* /\* (of_a > n * n) || (of_a < 0)){ *\/ */ - /* /\* fprintf(stderr, "%d %d %d\n", i, j, k); *\/ */ - /* /\* fprintf(stderr, "of_a: %d\nof_b: %d\nof_c: %d\n", of_a, of_b, of_c); *\/ */ - /* /\* exit(2); *\/ */ - /* /\* } *\/ */ - /* tmp += a[of_a] * b[of_b]; */ - /* of_a += 1; // k */ - /* of_b += n; // k */ - /* } */ - /* c[of_c] += tmp; */ - /* of_b -= n * n; */ - /* of_a -= 1 * n; */ - /* // */ - /* of_c += 1; // j */ - /* of_b += 1; // j */ - /* } */ - /* of_c -= 1 * n; */ - /* // */ - /* of_c += n;//i j */ - /* of_a += n;//i k */ - /* of_b = 0;// k j */ - /* } */ - of_a = 0; of_b = 0; @@ -79,24 +45,12 @@ int main(int argc, char *argv[]){ of_c += 1; of_b += 1; } - /* for(l = 0; l < n * n; l++) */ - /* fprintf(stdout, "%lf\t", c[l]); */ - /* fprintf(stdout, "\n", c[l]); */ of_c -= n; of_a += 1; } of_c += n; of_b = 0; - } - - /* double err = 0.0; */ - /* for(i = 0; i < n * n; i++) */ - /* err += fabs(c[i]); */ - /* fprintf(stdout, "err: %lf\n", err); */ - - /* for(i = 0; i < n * n; i++) */ - /* fprintf(stdout, "%lf\n", c[i]); */ - + } co = clock(); //GCC is lazy! @@ -105,23 +59,5 @@ int main(int argc, char *argv[]){ fclose(fdump); fprintf(stdout, "Time: %lf\n", (double) (co - ci) / CLOCKS_PER_SEC); - /* FILE *out; */ - - /* out = fopen("a", "w"); */ - /* for(i = 0; i < n; i++){ */ - /* for(j = 0; j < n; j++) */ - /* fprintf(out, "%lf\t", a[n * i + j]); */ - /* fprintf(out, "\n"); */ - /* } */ - /* fclose(out); */ - - /* out = fopen("c", "w"); */ - /* for(i = 0; i < n; i++){ */ - /* for(j = 0; j < n; j++) */ - /* fprintf(out, "%lf\t", c[n * i + j]); */ - /* fprintf(out, "\n"); */ - /* } */ - /* fclose(out); */ - free(a); free(b); free(c); } ----------------------------------------------------------------------- Summary of changes: tests/mm.c | 66 +----------------------------------------------------------- 1 files changed, 1 insertions(+), 65 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-01-29 17:57:28
|
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, tensor has been updated via 13671606829449440d65e08dd7ec8b54560c8f23 (commit) via 10d713510529c4cc86dee7c000eb400b6ae81b06 (commit) from b0e4abdd67877cc2e4a2880542028d0d15c8b3eb (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 13671606829449440d65e08dd7ec8b54560c8f23 Author: Akshay Srinivasan <aks...@gm...> Date: Tue Jan 29 00:54:36 2013 -0800 Added C version of mm (for benchmarking later). diff --git a/tests/mm.c b/tests/mm.c new file mode 100644 index 0000000..d7a8bb3 --- /dev/null +++ b/tests/mm.c @@ -0,0 +1,127 @@ +#include <stdio.h> +#include <stdlib.h> +#include <time.h> +#include <math.h> + +int main(int argc, char *argv[]){ + int n; + if(argc < 2){ + fprintf(stderr, "Expected argument!\n"); + exit(2); + } + n = atoi(argv[1]); + + //This doesn't really help. + /* double *restrict a, *restrict b, *restrict c; */ + double *a, *b, *c; + + a = calloc(n * n, sizeof(double)); + b = calloc(n * n, sizeof(double)); + c = calloc(n * n, sizeof(double));; + /* double a[] = {1, 2, 3, 4}, */ + /* b[] = {4, 5, 6, 7}, */ + /* c[] = {0, 0, 0, 0}; */ + + int i, j, k, l; + printf("N: %d\n", n); + + for(i = 0; i < n; i++) + for(j = 0; j < n; j++){ + a[n * i + j] = ((double) i + j) * 1e-3; + b[n * i + j] = ((double) i + j) * 1e-3; + } + + double tmp; + int of_a = 0, of_b = 0, of_c = 0; + + clock_t co, ci; + + ci = clock(); + + /* for(i = 0; i < n; i++){ */ + /* for(j = 0; j < n; j++){ */ + /* tmp = 0.; */ + /* for(k = 0; k < n; k++){ */ + /* /\* if((of_c > n * n) || (of_c < 0) || \ *\/ */ + /* /\* (of_b > n * n) || (of_b < 0) || \ *\/ */ + /* /\* (of_a > n * n) || (of_a < 0)){ *\/ */ + /* /\* fprintf(stderr, "%d %d %d\n", i, j, k); *\/ */ + /* /\* fprintf(stderr, "of_a: %d\nof_b: %d\nof_c: %d\n", of_a, of_b, of_c); *\/ */ + /* /\* exit(2); *\/ */ + /* /\* } *\/ */ + /* tmp += a[of_a] * b[of_b]; */ + /* of_a += 1; // k */ + /* of_b += n; // k */ + /* } */ + /* c[of_c] += tmp; */ + /* of_b -= n * n; */ + /* of_a -= 1 * n; */ + /* // */ + /* of_c += 1; // j */ + /* of_b += 1; // j */ + /* } */ + /* of_c -= 1 * n; */ + /* // */ + /* of_c += n;//i j */ + /* of_a += n;//i k */ + /* of_b = 0;// k j */ + /* } */ + + + of_a = 0; + of_b = 0; + of_c = 0; + for(i = 0; i < n; i++){ + for(k = 0; k < n; k++){ + tmp = a[of_a]; + for(j = 0; j < n; j++){ + c[of_c] += tmp * b[of_b]; + of_c += 1; + of_b += 1; + } + /* for(l = 0; l < n * n; l++) */ + /* fprintf(stdout, "%lf\t", c[l]); */ + /* fprintf(stdout, "\n", c[l]); */ + of_c -= n; + of_a += 1; + } + of_c += n; + of_b = 0; + } + + /* double err = 0.0; */ + /* for(i = 0; i < n * n; i++) */ + /* err += fabs(c[i]); */ + /* fprintf(stdout, "err: %lf\n", err); */ + + /* for(i = 0; i < n * n; i++) */ + /* fprintf(stdout, "%lf\n", c[i]); */ + + co = clock(); + + //GCC is lazy! + FILE *fdump = fopen("/dev/null", "w"); + fprintf(fdump, "%lf\n", c[n * (n - 1) + (n - 1)]); + fclose(fdump); + + fprintf(stdout, "Time: %lf\n", (double) (co - ci) / CLOCKS_PER_SEC); + /* FILE *out; */ + + /* out = fopen("a", "w"); */ + /* for(i = 0; i < n; i++){ */ + /* for(j = 0; j < n; j++) */ + /* fprintf(out, "%lf\t", a[n * i + j]); */ + /* fprintf(out, "\n"); */ + /* } */ + /* fclose(out); */ + + /* out = fopen("c", "w"); */ + /* for(i = 0; i < n; i++){ */ + /* for(j = 0; j < n; j++) */ + /* fprintf(out, "%lf\t", c[n * i + j]); */ + /* fprintf(out, "\n"); */ + /* } */ + /* fclose(out); */ + + free(a); free(b); free(c); +} commit 10d713510529c4cc86dee7c000eb400b6ae81b06 Author: Akshay Srinivasan <aks...@gm...> Date: Tue Jan 29 00:49:02 2013 -0800 o Cosmetic changes to loops, whitespace cleanup. o Work on getrf.lisp is in progress. diff --git a/matlisp.asd b/matlisp.asd index 287657f..520d66c 100644 --- a/matlisp.asd +++ b/matlisp.asd @@ -36,7 +36,7 @@ (asdf:defsystem matlisp-packages :depends-on (#:cffi) :pathname #.(translate-logical-pathname "matlisp:srcdir;") - :components + :components ((:file "packages"))) (asdf:defsystem matlisp-config @@ -84,7 +84,7 @@ "matlisp-utilities" "fortran-names") :components ((:module "foreign-interface" - :pathname "ffi" + :pathname "ffi" :components ((:file "ffi-cffi") (:file "ffi-cffi-implementation-specific") (:file "foreign-vector") @@ -103,7 +103,9 @@ (:module "matlisp-base" :depends-on ("foreign-core") :pathname "base" - :components ((:file "standard-tensor") + :components ((:file "tweakable") + (:file "standard-tensor" + :depends-on ("tweakable")) ;; (:file "loopy" :depends-on ("standard-tensor")) @@ -116,9 +118,7 @@ (:file "blas-helpers" :depends-on ("standard-tensor" "permutation")) (:file "print" - :depends-on ("standard-tensor")) - ;;Probably not the right place, but should do. - (:file "tweakable"))) + :depends-on ("standard-tensor")))) (:module "matlisp-classes" :pathname "classes" :depends-on ("matlisp-base") diff --git a/src/base/blas-helpers.lisp b/src/base/blas-helpers.lisp index 07dd214..f30164e 100644 --- a/src/base/blas-helpers.lisp +++ b/src/base/blas-helpers.lisp @@ -13,18 +13,18 @@ (perm-b-dims (permute (dimensions ten-b) std-a-perm) :type index-store-vector)) (very-quickly (loop - for i of-type index-type from 0 below (rank ten-a) - for sost-a across sort-std-a - for a-aoff of-type index-type = (aref sort-std-a 0) then (the index-type (* a-aoff (aref perm-a-dims (1- i)))) + :for i :of-type index-type :from 0 :below (rank ten-a) + :for sost-a :across sort-std-a + :for a-aoff :of-type index-type := (aref sort-std-a 0) :then (the index-type (* a-aoff (aref perm-a-dims (1- i)))) ;; - for sost-b across sort-std-b - for b-aoff of-type index-type = (aref sort-std-b 0) then (the index-type (* b-aoff (aref perm-b-dims (1- i)))) + :for sost-b :across sort-std-b + :for b-aoff :of-type index-type := (aref sort-std-b 0) :then (the index-type (* b-aoff (aref perm-b-dims (1- i)))) ;; - do (progn - (unless (and (= sost-a a-aoff) - (= sost-b b-aoff)) - (return nil))) - finally (return (list (aref sort-std-a 0) (aref sort-std-b 0))))))) + :do (progn + (unless (and (= sost-a a-aoff) + (= sost-b b-aoff)) + (return nil))) + :finally (return (list (aref sort-std-a 0) (aref sort-std-b 0))))))) (defun consecutive-store-p (tensor) (declare (type standard-tensor tensor)) @@ -34,11 +34,11 @@ (perm-dims (permute (dimensions tensor) std-perm) :type index-store-vector)) (very-quickly (loop - for so-st across sort-std - for so-di across perm-dims - and accumulated-off = (aref sort-std 0) then (the index-type (* accumulated-off so-di)) - unless (= so-st accumulated-off) do (return nil) - finally (return (aref sort-std 0)))))) + :for so-st :across sort-std + :for so-di :across perm-dims + :and accumulated-off := (aref sort-std 0) :then (the index-type (* accumulated-off so-di)) + :unless (= so-st accumulated-off) :do (return nil) + :finally (return (aref sort-std 0)))))) (defun blas-matrix-compatible-p (matrix op) (declare (type standard-matrix matrix)) @@ -100,3 +100,6 @@ (assert (> st 0) nil 'tensor-invalid-dimension-value :argument i :dimension (aref dims i)) (setf (aref stds i) st)) :finally (return (values stds st)))))) + +(defun make-stride (dims) + (ecase *default-stride-ordering* (:row-major (make-stride-rmj dims)) (:col-major (make-stride-cmj dims)))) diff --git a/src/base/permutation.lisp b/src/base/permutation.lisp index 81a5281..0a69941 100644 --- a/src/base/permutation.lisp +++ b/src/base/permutation.lisp @@ -8,7 +8,7 @@ `(simple-array perrepr-type (,size))) (make-array-allocator allocate-perrepr-store 'perrepr-type 0 - " + " Syntax ====== (ALLOCATE-PERREPR-STORE SIZE [INITIAL-ELEMENT 0]) @@ -24,8 +24,8 @@ (declare (type perrepr-vector ret)) (very-quickly (loop - for i of-type perrepr-type from 0 below n - do (setf (aref ret i) i))) + :for i :of-type perrepr-type :from 0 :below n + :do (setf (aref ret i) i))) ret)) (definline perv (&rest contents) @@ -72,10 +72,10 @@ (declare (type perrepr-vector sort) (type index-type len)) (very-quickly - (loop for i of-type index-type from 0 below len - unless (= (aref sort i) i) - do (return nil) - finally (return t)))))))) + (loop :for i :of-type index-type :from 0 :below len + :unless (= (aref sort i) i) + :do (return nil) + :finally (return t)))))))) (assert (action-repr-p (repr per)) nil 'permutation-invalid-error) (setf (group-rank per) (length (repr per))))) ;; @@ -99,21 +99,20 @@ (let ((sort (very-quickly (sort (copy-seq perm) #'<)))) (declare (type perrepr-vector sort)) (very-quickly - (loop for i of-type index-type from 1 below len - when (= (aref sort i) (aref sort (1- i))) - do (return nil) - finally (return t)))))))))) + (loop :for i :of-type index-type :from 1 :below len + :when (= (aref sort i) (aref sort (1- i))) + :do (return nil) + :finally (return t)))))))))) (if (null (repr per)) (setf (group-rank per) 0) (very-quickly (loop - for cyc of-type perrepr-vector in (repr per) - unless (cycle-repr-p cyc) - do (error 'permutation-invalid-error) - maximizing (lvec-max cyc) into g-rnk of-type index-type - finally (setf (group-rank per) (the index-type (1+ g-rnk)))))))) + :for cyc :of-type perrepr-vector :in (repr per) + :unless (cycle-repr-p cyc) + :do (error 'permutation-invalid-error) + :maximizing (lvec-max cyc) into g-rnk of-type index-type + :finally (setf (group-rank per) (the index-type (1+ g-rnk)))))))) ;; - (defclass permutation-pivot-flip (permutation) ((representation :type perrepr-vector))) @@ -132,14 +131,41 @@ (declare (type perrepr-vector idiv) (type index-type len)) (very-quickly - (loop for i of-type index-type from 0 below len - do (unless (< -1 (aref idiv i) len) - (return nil)) - finally (return t))))))) + (loop :for i :of-type index-type :from 0 :below len + :do (unless (< -1 (aref idiv i) len) + (return nil)) + :finally (return t))))))) (assert (pivot-flip-p (repr per)) nil 'permutation-invalid-error) (setf (group-rank per) (length (repr per))))) ;; +;; (defclass permutation-matrix (permutation) +;; ((representation :type perrepr-vector))) + +;; (defmethod initialize-instance :after ((per permutation-pivot-flip) &rest initargs) +;; (declare (ignore initargs)) +;; (labels ((pivot-flip-p (idiv) +;; " +;; Checks if ARR is a possible permutation pivot-flip. A pivot +;; flip is more algorithmic in its representation. If a sequence +;; is given, apply a pivot-flip on it is equivalent to running from +;; the left to the right of the permutation by flipping (pi(i), i) +;; sequentially. This is the representation used in LAPACK. +;; " +;; (if (not (typep idiv 'perrepr-vector)) nil +;; (let ((len (length idiv))) +;; (declare (type perrepr-vector idiv) +;; (type index-type len)) +;; (very-quickly +;; (loop :for i :of-type index-type :from 0 :below len +;; :do (unless (< -1 (aref idiv i) len) +;; (return nil)) +;; :finally (return t))))))) +;; (assert (pivot-flip-p (repr per)) nil 'permutation-invalid-error) +;; (setf (group-rank per) (length (repr per))))) + + +;; (definline make-pcycle (&rest args) (make-instance 'permutation-cycle :repr args)) @@ -165,7 +191,7 @@ (let ((len (aref (dimensions ten) arg))) (assert (>= len (group-rank perm)) nil 'permutation-permute-error :seq-len len :group-rank (group-rank perm))))) - + (definline permute (thing perm &optional (arg 0)) (permute! (copy thing) perm arg)) @@ -182,16 +208,16 @@ (when (< i glen) (rplaca x (aref cseq (aref act i))) (incf i)))) seq))) - + (defmethod permute! ((seq vector) (perm permutation-action) &optional arg) (declare (ignore arg)) (let ((cseq (make-array (length seq) :initial-contents seq)) (act (repr perm))) - (loop - for i from 0 below (group-rank perm) - do (unless (= i (aref act i)) - (setf (aref seq i) (aref cseq (aref act i)))) - finally (return seq)))) + (loop + :for i :from 0 :below (group-rank perm) + :do (unless (= i (aref act i)) + (setf (aref seq i) (aref cseq (aref act i)))) + :finally (return seq)))) (defmethod permute! ((ten standard-tensor) (perm permutation-action) &optional (arg 0)) (let ((cyc (action->cycle perm))) @@ -203,10 +229,10 @@ (declare (type perrepr-vector pcyc) (type vector seq)) (let ((xl (aref seq (aref pcyc (1- (length pcyc)))))) - (loop for i of-type index-type downfrom (1- (length pcyc)) to 0 - do (setf (aref seq (aref pcyc i)) - (if (= i 0) xl - (aref seq (aref pcyc (1- i)))))))) + (loop :for i :of-type index-type :downfrom (1- (length pcyc)) :to 0 + :do (setf (aref seq (aref pcyc i)) + (if (= i 0) xl + (aref seq (aref pcyc (1- i)))))))) (defmethod permute! ((seq cons) (perm permutation-cycle) &optional arg) (declare (ignore arg)) @@ -230,37 +256,37 @@ (defmethod permute! ((A standard-tensor) (perm permutation-cycle) &optional (arg 0)) (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) - (rplaca (nthcdr arg slst) 0) - (values (sub-tensor~ A slst) (sub-tensor~ A slst))) + (rplaca (nthcdr arg slst) 0) + (values (sub-tensor~ A slst) (sub-tensor~ A slst))) (let-typed ((cyclst (repr perm) :type cons) (cp-ten (make-instance (class-of tone) :dimensions (copy-seq (dimensions tone)))) (std-arg (aref (strides A) arg) :type index-type) (hd-sl (head ttwo) :type index-type)) - (dolist (cyc cyclst) - (declare (type perrepr-vector cyc)) - (setf (head tone) (+ hd-sl (* std-arg (aref cyc (1- (length cyc)))))) - (copy! tone cp-ten) - (loop for i of-type index-type downfrom (1- (length cyc)) to 0 - do (progn - (setf (head tone) (+ hd-sl (* std-arg (aref cyc i)))) - (copy! - (if (= i 0) cp-ten - (progn - (setf (head ttwo) (+ hd-sl (* std-arg (aref cyc (1- i))))) - ttwo)) - tone)))))) + (dolist (cyc cyclst) + (declare (type perrepr-vector cyc)) + (setf (head tone) (+ hd-sl (* std-arg (aref cyc (1- (length cyc)))))) + (copy! tone cp-ten) + (loop :for i :of-type index-type :downfrom (1- (length cyc)) :to 0 + :do (progn + (setf (head tone) (+ hd-sl (* std-arg (aref cyc i)))) + (copy! + (if (= i 0) cp-ten + (progn + (setf (head ttwo) (+ hd-sl (* std-arg (aref cyc (1- i))))) + ttwo)) + tone)))))) A) ;;Pivot idx (defmethod permute! ((seq vector) (perm permutation-pivot-flip) &optional arg) (declare (ignore arg)) (let-typed ((pidx (repr perm) :type perrepr-vector)) - (loop for i of-type index-type from 0 below (group-rank perm) - unless (= i (aref pidx i)) - do (rotatef (aref seq i) (aref seq (aref pidx i))) - finally (return seq)))) - + (loop :for i :of-type index-type :from 0 :below (group-rank perm) + :unless (= i (aref pidx i)) + :do (rotatef (aref seq i) (aref seq (aref pidx i))) + :finally (return seq)))) + (defmethod permute! ((seq cons) (perm permutation-pivot-flip) &optional arg) (declare (ignore arg)) (let ((cseq (make-array (length seq) :initial-contents seq)) @@ -280,12 +306,12 @@ (values (sub-tensor~ A slst nil) (sub-tensor~ A slst nil))) (let ((argstd (aref (strides A) arg)) (hd-sl (head ttwo))) - (loop for i from 0 below (length idiv) - do (progn - (unless (= i (aref idiv i)) - (setf (head ttwo) (+ hd-sl (* (aref idiv i) argstd))) - (swap! tone ttwo)) - (incf (head tone) argstd)))))) + (loop :for i :from 0 :below (length idiv) + :do (progn + (unless (= i (aref idiv i)) + (setf (head ttwo) (+ hd-sl (* (aref idiv i) argstd))) + (swap! tone ttwo)) + (incf (head tone) argstd)))))) A) ;;Conversions----------------------------------------------------;; @@ -297,7 +323,7 @@ of a permutation. The first argument \"act\" is the action of the permutation on the array #(0 1 2 3 ..): an object of the class permutation-action. - + \"Canonical\" may be a bit of an overstatement; this is the way S_n was presented in Van der Waerden's book. " @@ -310,11 +336,11 @@ (if (= (aref arr x0) x0) (values 0 nil) (very-quickly (loop - for x of-type perrepr-type = (aref arr x0) then (aref arr x) - and ret of-type cons = (list x0) then (cons x ret) - counting t into i of-type index-type - when (= x x0) - do (return (values i ret)))))) + :for x :of-type perrepr-type := (aref arr x0) :then (aref arr x) + :and ret :of-type cons := (list x0) :then (cons x ret) + :counting t :into i :of-type index-type + :when (= x x0) + :do (return (values i ret)))))) (cycle-walk (cyc ignore) ;; Finds all cycles (let ((x0 (find-if-not #'(lambda (x) (member x ignore)) arr))) @@ -349,10 +375,10 @@ (let ((xl (aref act-repr (aref cyc (1- (length cyc)))))) (very-quickly (loop - for i of-type index-type downfrom (1- (length cyc)) to 0 - do (setf (aref act-repr (aref cyc i)) - (if (= i 0) xl - (aref act-repr (aref cyc (1- i))))))))) + :for i :of-type index-type :downfrom (1- (length cyc)) :to 0 + :do (setf (aref act-repr (aref cyc i)) + (if (= i 0) xl + (aref act-repr (aref cyc (1- i))))))))) (make-instance 'permutation-action :repr act-repr))) (defun pivot-flip->action (pflip) @@ -364,10 +390,10 @@ (let ((act (perrepr-id-action len))) (declare (type perrepr-vector act)) (very-quickly - (loop for i from 0 below len - do (let ((val (aref idiv i))) - (unless (= val i) - (rotatef (aref act i) (aref act val)))))) + (loop :for i :from 0 :below len + :do (let ((val (aref idiv i))) + (unless (= val i) + (rotatef (aref act i) (aref act val)))))) (make-instance 'permutation-action :repr act)))) (defun mod-max (seq lidx uidx) @@ -375,46 +401,46 @@ (let ((len (length seq))) (very-quickly (loop :for idx :of-type index-type :downfrom uidx :above lidx - with max of-type perrepr-type = (aref seq uidx) - with max-idx of-type index-type = uidx - do (let ((ele (aref seq (mod idx len)))) - (when (> ele max) - (setf max ele - max-idx idx))) - finally (return (values max max-idx)))))) + :with max :of-type perrepr-type := (aref seq uidx) + :with max-idx :of-type index-type := uidx + :do (let ((ele (aref seq (mod idx len)))) + (when (> ele max) + (setf max ele + max-idx idx))) + :finally (return (values max max-idx)))))) #| (defun cycle->pivot-flip (cyc) (let ((cp-cyc (copy-seq cyc))) (let - (labels ((mod-max (seq lidx uidx) - (declare (type perrepr-vector seq)) - (let ((len (length cyc))) - (very-quickly - (loop :for idx :of-type index-type :downfrom uidx :above lidx - with max of-type perrepr-type = (aref seq uidx) - with max-idx of-type index-type = uidx - do (let ((ele (aref seq (mod idx len)))) - (when (> ele max) - (setf max ele - max-idx idx))) - finally (return (values max max-idx)))))) - (get-flip (lidx uidx) - (multiple-value-bind (max max-idx) (mod-max cyc lidx uidx) - (let ((ele-0 (aref cyc (mod max-idx len))) - (ele-1 (aref cyc (mod (1- max-idx) len)))) - (setf (aref pidx (max ele-0 ele-1)) - (min ele-0 ele-1)) -|# + (labels ((mod-max (seq lidx uidx) + (declare (type perrepr-vector seq)) + (let ((len (length cyc))) + (very-quickly + (loop :for idx :of-type index-type :downfrom uidx :above lidx + :with max :of-type perrepr-type := (aref seq uidx) + :with max-idx :of-type index-type := uidx + :do (let ((ele (aref seq (mod idx len)))) + (when (> ele max) + (setf max ele + max-idx idx))) + :finally (return (values max max-idx)))))) + (get-flip (lidx uidx) + (multiple-value-bind (max max-idx) (mod-max cyc lidx uidx) + (let ((ele-0 (aref cyc (mod max-idx len))) + (ele-1 (aref cyc (mod (1- max-idx) len)))) + (setf (aref pidx (max ele-0 ele-1)) + (min ele-0 ele-1)) + |# #+nil (defun permute-argument (func-symbol perm) (declare (type symbol func-symbol) (type permutation perm)) (let* ((glen (group-rank perm)) - (args (loop for i from 0 below glen - collect (gensym)))) + (args (loop :for i :from 0 :below glen + :collect (gensym)))) (eval `(lambda (,@args &rest rest) (apply ',func-symbol (append (list ,@(permute! args perm)) rest)))))) @@ -436,7 +462,7 @@ (lambda (&rest args) (apply func-a (permute! (multiple-value-list (funcall func-b args)) perm)))) ;; - + (definline sort-permute (seq predicate) " Sorts a lisp-vector in-place, by using the function @arg{predicate} as the @@ -448,55 +474,55 @@ (let*-typed ((len (length seq) :type fixnum) (perm (perrepr-id-action len) :type perrepr-vector) (jobs (list `(0 ,len)))) - (loop ;;:repeat 10 - :for bounds := (pop jobs) :then (pop jobs) - :until (null bounds) - :finally (return (values seq (make-instance 'permutation-action :repr perm))) - :do (let*-typed ((below-idx (first bounds) :type fixnum) - (above-idx (second bounds) :type fixnum) - (piv (+ below-idx (floor (- above-idx below-idx) 2)) :type fixnum)) - (loop ;;:repeat 10 - :with ele := (aref seq piv) - :with lbound :of-type fixnum := below-idx - :with ubound :of-type fixnum := (1- above-idx) - :until (progn - ;;(format t "~%~a ~%" (list lbound piv ubound)) - (loop :for i :of-type fixnum :from lbound :to piv - :until (or (= i piv) (funcall predicate ele (aref seq i))) - :finally (setq lbound i)) - (loop :for i :of-type fixnum :downfrom ubound :to piv - :until (or (= i piv) (funcall predicate (aref seq i) ele)) - :finally (setq ubound i)) - ;;(format t "~a ~%" (list lbound piv ubound)) - (cond - ((= ubound lbound piv) - (when (> (- piv below-idx) 1) - (push `(,below-idx ,piv) jobs)) - (when (> (- above-idx (1+ piv)) 1) - (push `(,(1+ piv) ,above-idx) jobs)) - ;;(format t "~a~%" jobs) - t) - ((< lbound piv ubound) - (rotatef (aref seq lbound) (aref seq ubound)) - (rotatef (aref perm lbound) (aref perm ubound)) - (incf lbound) - (decf ubound) - nil) - ((= lbound piv) - (rotatef (aref seq piv) (aref seq (1+ piv))) - (rotatef (aref perm piv) (aref perm (1+ piv))) - (unless (= ubound (1+ piv)) - (rotatef (aref seq piv) (aref seq ubound)) - (rotatef (aref perm piv) (aref perm ubound))) - (incf piv) - (incf lbound) - nil) - ((= ubound piv) - (rotatef (aref seq (1- piv)) (aref seq piv)) - (rotatef (aref perm (1- piv)) (aref perm piv)) - (unless (= lbound (1- piv)) - (rotatef (aref seq lbound) (aref seq piv)) - (rotatef (aref perm lbound) (aref perm piv))) - (decf piv) - (decf ubound) - nil)))))))) + (loop ;;:repeat 10 + :for bounds := (pop jobs) :then (pop jobs) + :until (null bounds) + :finally (return (values seq (make-instance 'permutation-action :repr perm))) + :do (let*-typed ((below-idx (first bounds) :type fixnum) + (above-idx (second bounds) :type fixnum) + (piv (+ below-idx (floor (- above-idx below-idx) 2)) :type fixnum)) + (loop ;;:repeat 10 + :with ele := (aref seq piv) + :with lbound :of-type fixnum := below-idx + :with ubound :of-type fixnum := (1- above-idx) + :until (progn + ;;(format t "~%~a ~%" (list lbound piv ubound)) + (loop :for i :of-type fixnum :from lbound :to piv + :until (or (= i piv) (funcall predicate ele (aref seq i))) + :finally (setq lbound i)) + (loop :for i :of-type fixnum :downfrom ubound :to piv + :until (or (= i piv) (funcall predicate (aref seq i) ele)) + :finally (setq ubound i)) + ;;(format t "~a ~%" (list lbound piv ubound)) + (cond + ((= ubound lbound piv) + (when (> (- piv below-idx) 1) + (push `(,below-idx ,piv) jobs)) + (when (> (- above-idx (1+ piv)) 1) + (push `(,(1+ piv) ,above-idx) jobs)) + ;;(format t "~a~%" jobs) + t) + ((< lbound piv ubound) + (rotatef (aref seq lbound) (aref seq ubound)) + (rotatef (aref perm lbound) (aref perm ubound)) + (incf lbound) + (decf ubound) + nil) + ((= lbound piv) + (rotatef (aref seq piv) (aref seq (1+ piv))) + (rotatef (aref perm piv) (aref perm (1+ piv))) + (unless (= ubound (1+ piv)) + (rotatef (aref seq piv) (aref seq ubound)) + (rotatef (aref perm piv) (aref perm ubound))) + (incf piv) + (incf lbound) + nil) + ((= ubound piv) + (rotatef (aref seq (1- piv)) (aref seq piv)) + (rotatef (aref perm (1- piv)) (aref perm piv)) + (unless (= lbound (1- piv)) + (rotatef (aref seq lbound) (aref seq piv)) + (rotatef (aref perm lbound) (aref perm piv))) + (decf piv) + (decf ubound) + nil)))))))) diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp index 5113980..12eb8d2 100644 --- a/src/base/standard-tensor.lisp +++ b/src/base/standard-tensor.lisp @@ -261,40 +261,37 @@ (typecase idx (cons (store-indexing-lst idx (head tensor) (strides tensor) (dimensions tensor))) (vector (store-indexing-vec idx (head tensor) (strides tensor) (dimensions tensor))))) +;; -;;You have to love lexical scoping :) -(defparameter *default-stride-ordering* :row-major - " - Determines whether strides are row or column major by default. - Doing: - > (let ((*default-stride-ordering* :col-major)) - (make-real-tensor 10 10)) - returns a 10x10 matrix with Column major order. -") +(defmacro with-order (order &rest code) + `(let ((*default-stride-ordering* ,order)) + ,@code)) +;; (defmethod initialize-instance :after ((tensor standard-tensor) &rest initargs) (declare (ignore initargs)) (let-typed ((dims (dimensions tensor) :type index-store-vector)) (setf (rank tensor) (length dims)) - (assert (>= (head tensor) 0) nil 'tensor-invalid-head-value :head (head tensor) :tensor tensor) - (if (not (slot-boundp tensor 'strides)) - (multiple-value-bind (stds size) (ecase *default-stride-ordering* (:row-major (make-stride-rmj dims)) (:col-major (make-stride-cmj dims))) - (declare (type index-store-vector stds) - (type index-type size)) - (setf (number-of-elements tensor) size - (strides tensor) stds) - (assert (<= (+ (head tensor) (1- (number-of-elements tensor))) (store-size tensor)) nil 'tensor-insufficient-store :store-size (store-size tensor) :max-idx (+ (head tensor) (1- (number-of-elements tensor))) :tensor tensor)) - (very-quickly - (let-typed ((stds (strides tensor) :type index-store-vector)) - (loop :for i :of-type index-type :from 0 :below (rank tensor) - :for sz :of-type index-type := (aref dims 0) :then (the index-type (* sz (aref dims i))) - :for lidx :of-type index-type := (the index-type (* (aref stds 0) (1- (aref dims 0)))) :then (the index-type (+ lidx (the index-type (* (aref stds i) (1- (aref dims i)))))) - :do (progn - (assert (> (aref stds i) 0) nil 'tensor-invalid-stride-value :argument i :stride (aref stds i) :tensor tensor) - (assert (> (aref dims i) 0) nil 'tensor-invalid-dimension-value :argument i :dimension (aref dims i) :tensor tensor)) - :finally (progn - (assert (>= (the index-type (store-size tensor)) (the index-type (+ (the index-type (head tensor)) lidx))) nil 'tensor-insufficient-store :store-size (store-size tensor) :max-idx lidx :tensor tensor) - (setf (number-of-elements tensor) sz)))))))) + (when *check-after-initializing?* + (assert (>= (head tensor) 0) nil 'tensor-invalid-head-value :head (head tensor) :tensor tensor) + (if (not (slot-boundp tensor 'strides)) + (multiple-value-bind (stds size) (make-stride dims) + (declare (type index-store-vector stds) + (type index-type size)) + (setf (number-of-elements tensor) size + (strides tensor) stds) + (assert (<= (+ (head tensor) (1- (number-of-elements tensor))) (store-size tensor)) nil 'tensor-insufficient-store :store-size (store-size tensor) :max-idx (+ (head tensor) (1- (number-of-elements tensor))) :tensor tensor)) + (very-quickly + (let-typed ((stds (strides tensor) :type index-store-vector)) + (loop :for i :of-type index-type :from 0 :below (rank tensor) + :for sz :of-type index-type := (aref dims 0) :then (the index-type (* sz (aref dims i))) + :for lidx :of-type index-type := (the index-type (* (aref stds 0) (1- (aref dims 0)))) :then (the index-type (+ lidx (the index-type (* (aref stds i) (1- (aref dims i)))))) + :do (progn + (assert (> (aref stds i) 0) nil 'tensor-invalid-stride-value :argument i :stride (aref stds i) :tensor tensor) + (assert (> (aref dims i) 0) nil 'tensor-invalid-dimension-value :argument i :dimension (aref dims i) :tensor tensor)) + :finally (progn + (assert (>= (the index-type (store-size tensor)) (the index-type (+ (the index-type (head tensor)) lidx))) nil 'tensor-insufficient-store :store-size (store-size tensor) :max-idx lidx :tensor tensor) + (setf (number-of-elements tensor) sz))))))))) ;; (defgeneric tensor-ref (tensor subscripts) @@ -521,7 +518,7 @@ (end (if (eq end '*) (aref dims i) (progn (assert (and (typep end 'index-type) (<= 0 end (aref dims i))) nil 'tensor-index-out-of-bounds :argument i :index start :dimension (aref dims i)) - end)))) + end)))) (declare (type index-type start step end)) ;; (let-typed ((dim (ceiling (the index-type (- end start)) step) :type index-type)) diff --git a/src/base/tweakable.lisp b/src/base/tweakable.lisp index 536af42..2b203ea 100644 --- a/src/base/tweakable.lisp +++ b/src/base/tweakable.lisp @@ -1,13 +1,43 @@ (in-package #:matlisp) +;;These describe some global defaults governing the +;;runtime behavior of Matlisp. Although you can change +;;these by doing a setf during runtime, it is recommended +;;that you use lexical scoping to affect local changes to +;;code (global variables are only bad if you overwrite them :) + +;;Default ordering of strides +(defparameter *default-stride-ordering* :row-major + " + Determines whether strides are row or column major by default. + Doing: + > (let ((*default-stride-ordering* :col-major)) + (make-real-tensor 10 10)) + returns a 10x10 matrix with Column major order. +") + +(defparameter *check-after-initializing?* t + " + If t, then check for invalid values in the field of + the class in the :after specialized method (if defined), + else do nothing. One ought to be very carful when doing, + much of Matlisp's code is written on the assumption that + the fields of a tensor don't take invalid values; failing + which case, may lead to memory error. Use at your own risk. +") + ;;Level 1--------------------------------------------------------;; (defparameter *real-l1-fcall-lb* 20000 "If the size of the array is less than this parameter, the - lisp version of axpy is called in order to avoid FFI overheads") + lisp version of axpy is called in order to avoid FFI overheads. + The Fortran function is not called if the tensor does not have + a consecutive store (see blas-helpers.lisp/consecutive-store-p).") (defparameter *complex-l1-fcall-lb* 10000 "If the size of the array is less than this parameter, the - lisp version of axpy is called in order to avoid FFI overheads") + lisp version of axpy is called in order to avoid FFI overheads. + The Fortran function is not called if the tensor does not have + a consecutive store (see blas-helpers.lisp/consecutive-store-p).") ;;Level 2--------------------------------------------------------;; (defparameter *real-l2-fcall-lb* 1000 @@ -15,7 +45,10 @@ If the maximum dimension in the MV is lower than this parameter, then the lisp code is used by default, instead of calling BLAS. Used to avoid the FFI overhead when calling - MM with small matrices. + MM with small matrices. Note that if the dimensions do exceed + this lower bound, then the Fortran function is called even if + the matrix has a BLAS incompatible stride (by doing a copy). + Default set with SBCL on x86-64 linux. A reasonable value is something between 800 and 2000.") @@ -24,7 +57,10 @@ If the maximum dimension in the MV is lower than this parameter, then the lisp code is used by default, instead of calling BLAS. Used to avoid the FFI overhead when calling - MM with small matrices. + MM with small matrices. Note that if the dimensions do exceed + this lower bound, then the Fortran function is called even when + the matrices have a BLAS incompatible stride (by using a copy). + Default set with SBCL on x86-64 linux. A reasonable value is something between 400 and 1000.") ;;Level 3--------------------------------------------------------;; diff --git a/src/lapack/getrf.lisp b/src/lapack/getrf.lisp index 7e5f0eb..da4bea0 100644 --- a/src/lapack/getrf.lisp +++ b/src/lapack/getrf.lisp @@ -27,44 +27,37 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (in-package #:matlisp) +;;Possibly add Lisp routine in the future. (defmacro generate-typed-getrf! (func-name (tensor-class lapack-func)) (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))) - `(defun ,func-name (A ipiv) - (declare (type ,matrix-class A) - (type permutation-pivot-flip ipiv)) - (mlet* (((maj-A ld-A fop-A) (blas-matrix-compatible-p A :n) :type (symbol index-type nil))) - (if (not (eq maj-A :col-major)) - (let*-typed ((dims (dimensions A) :type index-store-vector) - ;;Column major - (stds (let*-typed ((rank (rank A) :type fixnum) - (stds (allocate-index-store rank) :type index-store-vector)) - (very-quickly - (loop - :for i :from 0 :below rank - :and st = 1 :then (the index-type (* st (aref dims i))) - :do (setf (aref stds i) st))) - stds) - :type index-store-vector) - (tmp (,(get tensor-class :copy) A (make-instance (class-of A) :dimensions dims :strides stds :store (,(get tensor-class :store-allocator) (lvec-foldr #'* dims)))))) - (mlet* (((maj-tmp ld-tmp fop-tmp) (blas-matrix-compatible-p tmp :n) :type (nil index-type nil)) - ((new-tmp new-ipiv info) (,lapack-func - (nrows tmp) (ncols tmp) (store tmp) - ld-tmp (repr ipiv) 0) :type (nil nil integer))) - (assert (= info 0) nil 'invalid-arguments :argnum (1- (- info)) :message (format-to-string "GETRF returned INFO: ~a." info)) - (,(get tensor-class :copy) tmp A))) - (mlet* (((new-A new-ipiv info) (,lapack-func - (nrows A) (ncols A) (store A) - ld-A (repr ipiv) 0) :type (nil nil integer))) - (assert (= info 0) nil 'invalid-arguments :argnum (1- (- info)) :message (format-to-string "GETRF returned INFO: ~a." info)))) - ;;Convert from 1-based indexing to 0-based indexing, and fix - ;;other Fortran-ic quirks - (let-typed ((pidv (repr ipiv) :type perrepr-vector)) - (very-quickly - (loop for i from 0 below (length pidv) - do (decf (aref pidv i))))) - (values A ipiv))))) + `(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 :getrf) ',func-name + (get-tensor-class-optimization ',tensor-class) opt)) + (defun ,func-name (A ipiv) + (declare (type ,matrix-class A) + (type permutation-pivot-flip ipiv)) + (mlet* (((maj-A ld-A fop-A) (blas-matrix-compatible-p A :n) :type (symbol index-type nil))) + (if (eq maj-A :col-major) + (multiple-value-bind (n-A n-ipiv info) (,lapack-func + (nrows A) (ncols A) (store A) + ld-A (repr ipiv) 0) + (declare (ignore n-A n-ipiv)) + (assert (= info 0) nil 'invalid-arguments :argnum (1- (- info)) :message (format-to-string "GETRF returned INFO: ~a." info)) + ;;Convert back to 0-based indexing. + (let-typed ((pidv (repr ipiv) :type perrepr-vector)) + (very-quickly + (loop :for i :of-type index-type :from 0 :below (length pidv) + :do (decf (aref pidv i))))) + (values A ipiv info)) + (let ((tmp (,(getf opt :copy) A (with-order :col-major (,(getf opt :zero-maker) (dimensions A)))))) + (multiple-value-bind (n-tmp n-ipiv info) (,func-name tmp ipiv) + (declare (ignore n-tmp n-ipiv)) + (,(getf opt :copy) tmp A) + (values A ipiv info))))))))) (generate-typed-getrf! real-typed-getrf! (real-tensor dgetrf)) (generate-typed-getrf! complex-typed-getrf! (complex-tensor zgetrf)) @@ -92,16 +85,6 @@ U: upper triangular (upper trapezoidal when N<M) - If the optional argument IPIV is provided it must - be a (SIMPLE-ARRAY (UNSIGNED-BYTE 32) (*)) of dimension >= (MIN N M) - - IPIV is filled with the pivot indices that define the permutation - matrix P: - - row i of the matrix was interchanged with row IPIV(i). - - If IPIV is not provided, it is allocated by GESV. - Return Values ============= [1] The factors L and U from the factorization A = P*L*U where the @@ -109,37 +92,21 @@ [2] IPIV [3] INFO = T: successful i: U(i,i) is exactly zero. -") - (:method :before ((A standard-matrix) (ipiv permutation-pivot-flip)) - (assert (>= (group-rank ipiv) (idx-min (dimensions A))) nil 'invalid-value - :given (group-rank ipiv) :expected '(>= (group-rank ipiv) (idx-min (dimensions A)))))) - -(defmethod getrf! ((A real-matrix) (ipiv permutation-pivot-flip)) - (let* ((copy? (not (consecutive-store-p A))) - (cp-A (if copy? (copy A) A)) - (ret (multiple-value-list (real-typed-getrf! cp-A ipiv)))) - (when copy? - (copy! (first ret) A) - (rplaca ret A)) - (values-list ret))) - -(defmethod getrf! ((A complex-matrix) (ipiv permutation-pivot-flip)) - (let* ((copy? (not (consecutive-store-p A))) - (cp-A (if copy? (copy A) A)) - (ret (multiple-value-list (complex-typed-getrf! cp-A ipiv)))) - (when copy? - (copy! (first ret) A) - (rplaca ret A)) - (values-list ret))) - -(defun split-lu (A op-info) - (declare (type standard-matrix A)) - (destructuring-bind (&key decomposition-type row-permutation col-permutation) op-info - (assert (member decomposition-type '(:|U_ii=1| :|L_ii=1|)) nil 'invalid-arguments :message "Bad decomposition-type") - +")) + +(defmethod getrf! ((A real-matrix)) + (let ((ipiv (make-instance 'permutation-pivot-flip + :repr (perrepr-id-action (lvec-min (dimensions A)))))) + (real-typed-getrf! A ipiv))) + +(defmethod getrf! ((A complex-matrix)) + (let ((ipiv (make-instance 'permutation-pivot-flip + :repr (perrepr-id-action (lvec-min (dimensions A)))))) + (complex-typed-getrf! A ipiv))) + ;; (defgeneric lu (a &key with-l with-u with-p) - (:documentation + (:documentation " Syntax ====== @@ -160,57 +127,84 @@ By default WITH-L,WITH-U,WITH-P. ")) - -(defmethod lu ((a standard-matrix) &key (with-l t) (with-u t) (with-p t)) - (multiple-value-bind (lu ipiv info) - (getrf! (copy a)) - (declare (ignore info)) - - (let* ((result (list lu)) - (n (nrows a)) - (m (ncols a)) - (p (min n m))) - - (declare (type fixnum n m p)) - - ;; Extract the lower triangular part, if requested - (when with-l - (let ((lmat (typecase lu - (real-matrix (make-real-matrix-dim n p)) - (complex-matrix (make-complex-matrix-dim n p))))) - (dotimes (i p) - (setf (matrix-ref lmat i i) 1.0d0)) - (dotimes (i n) - (dotimes (j (min i p)) - (setf (matrix-ref lmat i j) (matrix-ref lu i j)))) - - (push lmat result))) - - - ;; Extract the upper triangular part, if requested - (when with-u - (let ((umat (typecase lu - (real-matrix (make-real-matrix-dim p m)) - (complex-matrix (make-complex-matrix-dim p m))))) - (dotimes (i p) - (loop for j from i to (1- m) - do (setf (matrix-ref umat i j) (matrix-ref lu i j)))) - - (push umat result))) - - ;; Extract the permutation matrix, if requested - (when with-p - (let* ((npiv (length ipiv)) - (pmat (make-real-matrix-dim n n)) - (pidx (make-array n :element-type '(unsigned-byte 32)))) - ;; Compute the P matrix from the pivot vector - (dotimes (k n) - (setf (aref pidx k) k)) - (dotimes (k npiv) - (rotatef (aref pidx k) (aref pidx (1- (aref ipiv k))))) - (dotimes (k n) - (setf (matrix-ref pmat (aref pidx k) k) 1)) - (push pmat result))) - - ;; Return the final result - (values-list (nreverse result))))) +;;Sure I can do this with a defmethod (slowly), but where's the fun in that ? :) +(defmacro make-lu (tensor-class) + (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))) + `(eval-when (:compile-toplevel :load-toplevel :execute) + (defmethod lu ((A ,matrix-class) &key with-l with-u with-p) + (multiple-value-bind (lu ipiv info) + (getrf! (with-order :col-major + (,(getf opt :copy) A (,(getf opt :zero-maker) (dimensions A))))) + (declare (ignore info)) + (let* ((ret (list lu)) + (n (nrows a)) + (m (ncols a)) + (p (min n m))) + (declare (type fixnum n m p)) + ;; Extract the lower triangular part, if requested + (when with-l + (let*-typed ((lmat (,(getf opt :zero-maker) (list n p)) :type ,matrix-class) + ;; + (l.rstd (row-stride lmat) :type index-type) + (l.cstd (col-stride lmat) :type index-type) + (l.of (head lmat) :type index-type) + (l.sto (store lmat) :type ,(linear-array-type (getf opt :element-type))) + ;; + (lu.rstd (row-stride lu) :type index-type) + (lu.cstd (col-stride lu) :type index-type) + (lu.of (head lu) :type index-type) + (lu.sto (store lu) :type ,(linear-array-type (getf opt :element-type)))) + (very-quickly + (loop :for j :of-type index-type :from 0 :below p + :do (progn + (,(getf opt :value-writer) (,(getf opt :fid*)) l.sto l.of) + (loop :repeat (- n j 1) + :do (progn + (incf lu.of lu.rstd) + (incf l.of l.rstd) + (,(getf opt :reader-writer) lu.sto lu.of l.sto l.of))) + (incf lu.of (- lu.cstd (the index-type (* (- n j 2) lu.rstd)))) + (incf l.of (- l.cstd (the index-type (* (- n j 2) l.rstd)))))) + (push lmat ret)))) + ;; Extract the upper triangular part, if requested + (when with-u + (let*-typed ((umat (,(getf opt :zero-maker) (list p m)) :type ,matrix-class) + ;; + (u.rstd (row-stride umat) :type index-type) + (u.cstd (col-stride umat) :type index-type) + (u.of (head umat) :type index-type) + (u.sto (store umat) :type ,(linear-array-type (getf opt :element-type))) + ;; + (lu.rstd (row-stride lu) :type index-type) + (lu.cstd (col-stride lu) :type index-type) + (lu.of (head lu) :type index-type) + (lu.sto (store lu) :type ,(linear-array-type (getf opt :element-type)))) + (progn + (loop :for i :of-type index-type :from 0 :below p + :do (progn + (loop :repeat (- m i) + :do (progn + (,(getf opt :reader-writer) lu.sto lu.of u.sto u.of) + (incf lu.of lu.cstd) + (incf u.of u.cstd))) + (incf lu.of (- lu.rstd (the index-type (* (- m i 1) lu.cstd)))) + (incf u.of (- u.rstd (the index-type (* (- m i 1) u.cstd)))))) + (push umat ret)))) + ;; Extract the permutation matrix, if requested + (if with-p + (let* ((npiv (length ipiv)) + (pmat (make-real-matrix-dim n n)) + (pidx (make-array n :element-type '(unsigned-byte 32)))) + ;; Compute the P matrix from the pivot vector + (dotimes (k n) + (setf (aref pidx k) k)) + (dotimes (k npiv) + (rotatef (aref pidx k) (aref pidx (1- (aref ipiv k))))) + (dotimes (k n) + (setf (matrix-ref pmat (aref pidx k) k) 1)) + (push pmat result))) + ;; Return the final result + (values-list (nreverse ret)))))))) + +(make-lu real-tensor) diff --git a/src/level-1/tensor-maker.lisp b/src/level-1/tensor-maker.lisp index 37c0b20..f4e66bc 100644 --- a/src/level-1/tensor-maker.lisp +++ b/src/level-1/tensor-maker.lisp @@ -14,9 +14,13 @@ (let*-typed ((vdim (make-index-store dims) :type index-store-vector) (ss (very-quickly (lvec-foldl #'(lambda (x y) (the index-type (* x y))) vdim))) (store (,(getf opt :store-allocator) ss)) - (rnk (length vdim))) - (make-instance (case rnk (2 ',(getf opt :matrix)) (1 ',(getf opt :vector)) (t ',tensor-class)) - :store store :store-size ss :dimensions vdim))) + (rnk (length vdim)) + (ret (let ((*check-after-initializing?* nil)) + (make-instance (case rnk (2 ',(getf opt :matrix)) (1 ',(getf opt :vector)) (t ',tensor-class)) + :strides (make-stride vdim) + :store store :store-size ss :dimensions vdim)))) + (setf (number-of-elements ret) ss) + ret)) (make-from-array (arr) (declare (type (array * *) arr)) (let* ((ret (make-dims (array-dimensions arr))) @@ -69,11 +73,12 @@ (setf (getf opt :zero-maker) ',func-name (get-tensor-class-optimization ',tensor-class) opt)) (defun ,func-name (dims) - (declare (type index-store-vector dims)) - (let-typed ((rnk (length dims) :type index-type) - (size (very-quickly (lvec-foldl #'(lambda (a b) (declare (type index-type a b)) (the index-type (* a b))) dims)))) - (make-instance (case rnk (2 ',(getf opt :matrix)) (1 ',(getf opt :vector)) (t ',tensor-class)) - :dimensions (copy-seq dims) :store (,(getf opt :store-allocator) size) :store-size size)))))) + (declare (type (or cons index-store-vector) dims)) + (let*-typed ((dims (if (consp dims) (make-index-store dims) (copy-seq dims)) :type index-store-vector) + (rnk (length dims) :type index-type) + (size (very-quickly (lvec-foldl #'(lambda (a b) (declare (type index-type a b)) (the index-type (* a b))) dims)))) + (make-instance (case rnk (2 ',(getf opt :matrix)) (1 ',(getf opt :vector)) (t ',tensor-class)) + :dimensions dims :store (,(getf opt :store-allocator) size) :store-size size)))))) (make-zeros-dims real-typed-zeros (real-tensor)) (make-zeros-dims complex-typed-zeros (complex-tensor)) ----------------------------------------------------------------------- Summary of changes: matlisp.asd | 12 +- src/base/blas-helpers.lisp | 33 +++-- src/base/permutation.lisp | 332 ++++++++++++++++++++++------------------- src/base/standard-tensor.lisp | 55 ++++---- src/base/tweakable.lisp | 44 +++++- src/lapack/getrf.lisp | 248 +++++++++++++++---------------- src/level-1/tensor-maker.lisp | 21 ++- tests/mm.c | 127 ++++++++++++++++ 8 files changed, 530 insertions(+), 342 deletions(-) create mode 100644 tests/mm.c hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-01-28 02:10:47
|
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, tensor has been updated via b0e4abdd67877cc2e4a2880542028d0d15c8b3eb (commit) from bee3634bf96916ccb63248bd8fa61c16d7c25da3 (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 b0e4abdd67877cc2e4a2880542028d0d15c8b3eb Author: Akshay Srinivasan <aks...@gm...> Date: Sun Jan 27 18:10:23 2013 -0800 Made gemm use the <type>-zeros function instead of make-<type>-tensor method. diff --git a/src/level-3/gemm.lisp b/src/level-3/gemm.lisp index d6597e0..bd59f29 100644 --- a/src/level-3/gemm.lisp +++ b/src/level-3/gemm.lisp @@ -448,10 +448,11 @@ (defmethod gemm ((alpha number) (a standard-matrix) (b standard-matrix) (beta number) (c real-matrix) &optional (job :nn)) - (let ((result (if (or (complexp alpha) (complexp beta) - (typep a 'complex-matrix) (typep b 'complex-matrix)) - (make-complex-tensor (nrows C) (ncols C)) - (make-real-tensor (nrows C) (ncols C))))) + (let ((result (funcall (if (or (complexp alpha) (complexp beta) + (typep a 'complex-matrix) (typep b 'complex-matrix)) + #'complex-typed-zeros + #'real-typed-zeros) + (dimensions C)))) (copy! C result) (gemm! alpha A B beta result job))) @@ -459,11 +460,10 @@ (beta (eql nil)) (c (eql nil)) &optional (job :nn)) (multiple-value-bind (job-A job-B) (split-job job) - (let ((result (apply - (if (or (complexp alpha) (complexp beta) - (typep a 'complex-matrix) (typep b 'complex-matrix)) - #'make-complex-tensor - #'make-real-tensor) - (list (if (member job-A '(:n :c)) (nrows A) (ncols A)) - (if (member job-B '(:n :c)) (ncols B) (nrows B)))))) + (let ((result (funcall (if (or (complexp alpha) (complexp beta) + (typep a 'complex-matrix) (typep b 'complex-matrix)) + #'complex-typed-zeros + #'real-typed-zeros) + (make-index-store (list (if (member job-A '(:n :c)) (nrows A) (ncols A)) + (if (member job-B '(:n :c)) (ncols B) (nrows B))))))) (gemm! alpha A B 0 result job)))) ----------------------------------------------------------------------- Summary of changes: src/level-3/gemm.lisp | 22 +++++++++++----------- 1 files changed, 11 insertions(+), 11 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-01-27 09:20:13
|
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, tensor has been updated via bee3634bf96916ccb63248bd8fa61c16d7c25da3 (commit) from 0c9fe59c75989013c9154326a382abd22fc056f0 (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 bee3634bf96916ccb63248bd8fa61c16d7c25da3 Author: Akshay Srinivasan <aks...@gm...> Date: Sun Jan 27 01:15:23 2013 -0800 Added macrology to ignore multiple outputs inside the callback. diff --git a/src/ffi/f77-ffi.lisp b/src/ffi/f77-ffi.lisp index 1870a0d..b03b2c5 100644 --- a/src/ffi/f77-ffi.lisp +++ b/src/ffi/f77-ffi.lisp @@ -252,9 +252,6 @@ ,@(mapcar #'second return-vars))))))))) ;;TODO: Outputs are messed up inside the callback - ;;TODO: Define callbacks outside the function call and lexically bind functions inside the - ;; call. Callbacks allocate memory in some non-GC'ed part of the heap. Runs out of memory - ;; quite quickly. (defun %f77.def-fortran-callback (func callback-name return-type parm) (let* ((hack-return-type `,return-type) (hack-parm `(,@parm)) @@ -333,8 +330,9 @@ `(let (,@array-vars))) ;; `(multiple-value-bind (,retvar ,@(mapcar #'car return-vars)) (funcall ,func ,@func-pars) - ,@(when (eq hack-return-type :void) - `((declare (ignore ,retvar)))) + (declare (ignore ,@(mapcar #'car return-vars) + ,@(when (eq hack-return-type :void) + `(,retvar)))) ,@(mapcar #'(lambda (decl) (destructuring-bind (func-var ffi-var type) decl (if (member type '(:complex-single-float :complex-double-float)) ----------------------------------------------------------------------- Summary of changes: src/ffi/f77-ffi.lisp | 8 +++----- 1 files changed, 3 insertions(+), 5 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-01-27 09:09:18
|
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, tensor has been updated via 0c9fe59c75989013c9154326a382abd22fc056f0 (commit) from 23ae23fc820468a4dc5149850f7f85bd524ebd69 (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 0c9fe59c75989013c9154326a382abd22fc056f0 Author: Akshay Srinivasan <aks...@gm...> Date: Sun Jan 27 01:03:01 2013 -0800 o Fixed callback issue. Callbacks are now declared once for each FFI. o Added example to dlsode.lisp. diff --git a/src/ffi/f77-ffi.lisp b/src/ffi/f77-ffi.lisp index 301520e..1870a0d 100644 --- a/src/ffi/f77-ffi.lisp +++ b/src/ffi/f77-ffi.lisp @@ -59,7 +59,7 @@ ;; Pass a pointer to the function. (:callback :void) (t (error 'unknown-token :token type - :message "Don't know the given Fortran type.")))))) + :message "Don't know the given Fortran type.")))))) (defun %f77.get-return-type (type) " @@ -68,7 +68,7 @@ (if (or (%f77.cast-as-array-p type) (%f77.callback-type-p type)) (error 'invalid-type :given type :expected '(not (or (%f77.cast-as-array-p type) (%f77.callback-type-p type))) - :message "A Fortran function cannot return the given type.") + :message "A Fortran function cannot return the given type.") (%f77.cffi-type type))) (definline %f77.output-p (style) @@ -85,7 +85,7 @@ " Get the input type to be passed to CFFI." (assert (member style +ffi-styles+) nil 'unknown-token :token style - :message "Don't know how to handle style.") + :message "Don't know how to handle style.") (cond ;; Can't do much else if type is an array/complex or input is passed-by-value. ((or (%f77.callback-type-p type) @@ -106,242 +106,248 @@ (let* ((aux-pars nil) (new-pars - (mapcar #'(lambda (decl) - (destructuring-bind (name type &optional (style :input-reference)) decl - (case type - (:string - ;; String lengths are appended to the function arguments, - ;; passed by value. - (nconsc aux-pars `((,(scat "LEN-" name) ,(%f77.cffi-type :integer)))) - `(,name ,(%f77.cffi-type :string))) - (t - `(,name ,(%f77.get-read-in-type type style)))))) - pars))) + (mapcar #'(lambda (decl) + (destructuring-bind (name type &optional (style :input-reference)) decl + (case type + (:string + ;; String lengths are appended to the function arguments, + ;; passed by value. + (nconsc aux-pars `((,(scat "LEN-" name) ,(%f77.cffi-type :integer)))) + `(,name ,(%f77.cffi-type :string))) + (t + `(,name ,(%f77.get-read-in-type type style)))))) + pars))) `( ;; don't want documentation for direct interface, not useful ;; ,@doc ,@new-pars ,@aux-pars)))) ;; Create a form specifying a simple Lisp function that calls the -;; underlying Fortran routine of the same name. -(defun %f77.def-fortran-interface (name return-type body hidden-var-name) - (multiple-value-bind (doc pars) - (parse-doc-&-parameters body) - (let ((ffi-fn (make-fortran-ffi-name name)) - (return-vars nil) - (array-vars nil) - (ref-vars nil) - (callback-code nil) - ;; - (defun-args nil) - (defun-keyword-args nil) - ;; - (aux-args nil) - ;; - (ffi-args nil) - (aux-ffi-args nil)) - (dolist (decl pars) - (destructuring-bind (var type &optional style) decl - (let ((ffi-var nil) - (aux-var nil)) - (cond - ;; Callbacks are tricky. - ((%f77.callback-type-p type) - (let* ((callback-name (gensym (symbol-name var))) - (c-callback-code (%f77.def-fortran-callback var callback-name (second type) (cddr type)))) - (nconsc callback-code c-callback-code) - (setq ffi-var `(cffi:callback ,callback-name)))) - ;; Can't really enforce "style" when given an array. - ;; Complex numbers do not latch onto this case, they - ;; are passed by value. - ((%f77.array-p type) - (setq ffi-var (scat "ADDR-" var)) - (nconsc array-vars `((,ffi-var ,var))) - ;; - (when-let (arg (getf type :inc)) - (nconsc defun-keyword-args - `((,arg 0))) - (nconc (car (last array-vars)) `(:inc-type ,(cadr type) :inc ,arg)))) - ;; Strings - ((%f77.string-p type) - (setq ffi-var var) - (setq aux-var (scat "LEN-" var)) - (nconsc aux-args `((,aux-var (length (the string ,var)))))) - ;; Pass-by-value variables - ((eq style :input-value) - (setq ffi-var var)) - ;; Pass-by-reference variables - (t - (cond - ;; Makes more sense to copy complex numbers into - ;; arrays, rather than twiddling around with lisp - ;; memory internals. - ((member type '(:complex-single-float :complex-double-float)) - (setq ffi-var (scat "ADDR-REAL-CAST-" var)) - (nconsc ref-vars - `((,ffi-var ,(second (%f77.cffi-type type)) :count 2 :initial-contents (list (realpart ,var) (imagpart ,var)))))) - (t - (setq ffi-var (scat "REF-" var)) - (nconsc ref-vars - `((,ffi-var ,(%f77.cffi-type type) :initial-element ,var))))))) - ;; Output variables - (when (and (%f77.output-p style) (not (eq type :string))) - (nconsc return-vars - `((,ffi-var ,var ,type)))) - ;; Arguments for the lisp wrapper - (unless (eq var hidden-var-name) - (nconsc defun-args - `(,var))) - ;; Arguments for the FFI function - (nconsc ffi-args - `(,ffi-var)) - ;; Auxillary arguments for FFI - (unless (null aux-var) - (nconsc aux-ffi-args - `(,aux-var)))))) - ;;Complex returns through hidden variable. - (unless (null hidden-var-name) - (nconsc aux-args `((,hidden-var-name ,(ecase (second (first pars)) - (:complex-single-float #c(0e0 0e0)) - (:complex-double-float #c(0d0 0d0))))))) - ;;Keyword argument list - (unless (null defun-keyword-args) - (setq defun-keyword-args (cons '&optional defun-keyword-args))) - ;;Return the function definition - (let ((retvar (gensym))) - `( - ,(recursive-append - `(defun ,name ,(append defun-args defun-keyword-args) - ,@doc) + ;; underlying Fortran routine of the same name. + (defun %f77.def-fortran-interface (name return-type body hidden-var-name) + (multiple-value-bind (doc pars) + (parse-doc-&-parameters body) + (let ((ffi-fn (make-fortran-ffi-name name)) + (return-vars nil) + (array-vars nil) + (ref-vars nil) + (callback-code nil) ;; - (unless (null aux-args) - `(let (,@aux-args))) - ;;Don't use with-foreign.. if ref-vars is nil - (unless (null ref-vars) - `(with-foreign-objects-stacked (,@ref-vars))) - ;;Don't use with-vector-dat.. if array-vars is nil - (unless (null array-vars) - `(with-vector-data-addresses (,@array-vars))) + (defun-args nil) + (defun-keyword-args nil) + ;; + (aux-args nil) + ;; + (ffi-args nil) + (aux-ffi-args nil) + (callback-args nil)) + (dolist (decl pars) + (destructuring-bind (var type &optional style) decl + (let ((ffi-var nil) + (aux-var nil)) + (cond + ;; Callbacks are tricky. + ((%f77.callback-type-p type) + (let* ((callback-name (second type)) + (field-gvar (intern (string+ "*" (symbol-name (gensym (symbol-name var))) "*"))) + (c-callback-code (%f77.def-fortran-callback field-gvar callback-name (third type) (cdddr type)))) + (nconsc callback-code `((defvar ,field-gvar nil) ,@c-callback-code)) + (nconsc callback-args `((,field-gvar ,var))) + (setq ffi-var `(cffi:callback ,callback-name)))) + ;; Can't really enforce "style" when given an array. + ;; Complex numbers do not latch onto this case, they + ;; are passed by value. + ((%f77.array-p type) + (setq ffi-var (scat "ADDR-" var)) + (nconsc array-vars `((,ffi-var ,var))) + ;; + (when-let (arg (getf type :inc)) + (nconsc defun-keyword-args + `((,arg 0))) + (nconc (car (last array-vars)) `(:inc-type ,(cadr type) :inc ,arg)))) + ;; Strings + ((%f77.string-p type) + (setq ffi-var var) + (setq aux-var (scat "LEN-" var)) + (nconsc aux-args `((,aux-var (length (the string ,var)))))) + ;; Pass-by-value variables + ((eq style :input-value) + (setq ffi-var var)) + ;; Pass-by-reference variables + (t + (cond + ;; Makes more sense to copy complex numbers into + ;; arrays, rather than twiddling around with lisp + ;; memory internals. + ((member type '(:complex-single-float :complex-double-float)) + (setq ffi-var (scat "ADDR-REAL-CAST-" var)) + (nconsc ref-vars + `((,ffi-var ,(second (%f77.cffi-type type)) :count 2 :initial-contents (list (realpart ,var) (imagpart ,var)))))) + (t + (setq ffi-var (scat "REF-" var)) + (nconsc ref-vars + `((,ffi-var ,(%f77.cffi-type type) :initial-element ,var))))))) + ;; Output variables + (when (and (%f77.output-p style) (not (eq type :string))) + (nconsc return-vars + `((,ffi-var ,var ,type)))) + ;; Arguments for the lisp wrapper + (unless (eq var hidden-var-name) + (nconsc defun-args + `(,var))) + ;; Arguments for the FFI function + (nconsc ffi-args + `(,ffi-var)) + ;; Auxillary arguments for FFI + (unless (null aux-var) + (nconsc aux-ffi-args + `(,aux-var)))))) + ;;Complex returns through hidden variable. + (unless (null hidden-var-name) + (nconsc aux-args `((,hidden-var-name ,(ecase (second (first pars)) + (:complex-single-float #c(0e0 0e0)) + (:complex-double-float #c(0d0 0d0))))))) + ;;Keyword argument list + (unless (null defun-keyword-args) + (setq defun-keyword-args (cons '&optional defun-keyword-args))) + ;;Return the function definition + (let ((retvar (gensym))) + `( ;;Declare callbacks - callback-code - ;;Call the foreign-function - `(let ((,retvar (,ffi-fn ,@ffi-args ,@aux-ffi-args))) - ;;Ignore return if type is :void - ,@(when (eq return-type :void) - `((declare (ignore ,retvar)))) - ;; Copy values in reference pointers back to local - ;; variables. Lisp has local scope; its safe to - ;; modify variables in parameter lists. - ,@(mapcar #'(lambda (decl) - (destructuring-bind (ffi-var var type) decl - (if (member type '(:complex-single-float :complex-double-float)) - `(setq ,var (complex (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 0) - (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 1))) - `(setq ,var (cffi:mem-aref ,ffi-var ,(%f77.cffi-type type)))))) - (remove-if-not #'(lambda (x) - (member (first x) ref-vars :key #'car)) - return-vars)) - (values - ,@(unless (eq return-type :void) - `(,retvar)) - ,@(mapcar #'second return-vars))))))))) - -;;TODO: Outputs are messed up inside the callback -;;TODO: Define callbacks outside the function call and lexically bind functions inside the -;; call. Callbacks allocate memory in some non-GC'ed part of the heap. Runs out of memory -;; quite quickly. -(defun %f77.def-fortran-callback (func callback-name return-type parm) - (let* ((hack-return-type `,return-type) - (hack-parm `(,@parm)) - (hidden-var-name nil)) - ;; - (when (member hack-return-type '(:complex-single-float :complex-double-float)) - (setq hidden-var-name (gensym "HIDDEN-COMPLEX-RETURN-")) - (setq hack-parm `((,hidden-var-name ,hack-return-type :output) - ,@parm)) - (setq hack-return-type :void)) - ;; - (let* ((new-pars nil) - (aux-pars nil) - (func-pars nil) - (array-vars nil) - (return-vars nil) - (ref-vars nil)) - (dolist (decl hack-parm) - (destructuring-bind (var type &optional (style :input)) decl - (let ((ffi-var nil) - (func-var nil)) - (cond - ;; Callbacks are tricky. - ((%f77.callback-type-p type) - (setq ffi-var var) - (setq func-var var)) - ;; - ((%f77.array-p type) - (setq ffi-var (scat "ADDR-" var)) - (setq func-var var) - (nconsc array-vars `((,func-var (make-foreign-vector :pointer ,ffi-var :type ,(second (%f77.cffi-type type)) - :size ,(if-let (size (getf type :size)) - size - 1)))))) + ,@callback-code + ,(recursive-append + `(defun ,name ,(append defun-args defun-keyword-args) + ,@doc) ;; - ((%f77.string-p type) - (setq ffi-var var) - (setq func-var var) - (nconsc aux-pars - `((,(scat "LEN-" var) ,(%f77.cffi-type :integer))))) + (unless (null aux-args) + `(let (,@aux-args))) + ;;Don't use with-foreign.. if ref-vars is nil + (unless (null ref-vars) + `(with-foreign-objects-stacked (,@ref-vars))) + ;;Don't use with-vector-dat.. if array-vars is nil + (unless (null array-vars) + `(with-vector-data-addresses (,@array-vars))) + ;;Point the the dummy global variables to the proper functions + (unless (null callback-args) + `(let (,@callback-args))) + ;;Call the foreign-function + `(let ((,retvar (,ffi-fn ,@ffi-args ,@aux-ffi-args))) + ;;Ignore return if type is :void + ,@(when (eq return-type :void) + `((declare (ignore ,retvar)))) + ;; Copy values in reference pointers back to local + ;; variables. Lisp has local scope; its safe to + ;; modify variables in parameter lists. + ,@(mapcar #'(lambda (decl) + (destructuring-bind (ffi-var var type) decl + (if (member type '(:complex-single-float :complex-double-float)) + `(setq ,var (complex (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 0) + (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 1))) + `(setq ,var (cffi:mem-aref ,ffi-var ,(%f77.cffi-type type)))))) + (remove-if-not #'(lambda (x) + (member (first x) ref-vars :key #'car)) + return-vars)) + (values + ,@(unless (eq return-type :void) + `(,retvar)) + ,@(mapcar #'second return-vars))))))))) + + ;;TODO: Outputs are messed up inside the callback + ;;TODO: Define callbacks outside the function call and lexically bind functions inside the + ;; call. Callbacks allocate memory in some non-GC'ed part of the heap. Runs out of memory + ;; quite quickly. + (defun %f77.def-fortran-callback (func callback-name return-type parm) + (let* ((hack-return-type `,return-type) + (hack-parm `(,@parm)) + (hidden-var-name nil)) + ;; + (when (member hack-return-type '(:complex-single-float :complex-double-float)) + (setq hidden-var-name (gensym "HIDDEN-COMPLEX-RETURN-")) + (setq hack-parm `((,hidden-var-name ,hack-return-type :output) + ,@parm)) + (setq hack-return-type :void)) + ;; + (let* ((new-pars nil) + (aux-pars nil) + (func-pars nil) + (array-vars nil) + (return-vars nil) + (ref-vars nil)) + (dolist (decl hack-parm) + (destructuring-bind (var type &optional (style :input)) decl + (let ((ffi-var nil) + (func-var nil)) + (cond + ;; Callbacks are tricky. + ((%f77.callback-type-p type) + (setq ffi-var var) + (setq func-var var)) + ;; + ((%f77.array-p type) + (setq ffi-var (scat "ADDR-" var)) + (setq func-var var) + (nconsc array-vars `((,func-var (make-foreign-vector :pointer ,ffi-var :type ,(second (%f77.cffi-type type)) + :size ,(if-let (size (getf type :size)) + size + 1)))))) + ;; + ((%f77.string-p type) + (setq ffi-var var) + (setq func-var var) + (nconsc aux-pars + `((,(scat "LEN-" var) ,(%f77.cffi-type :integer))))) + ;; + ((eq style :input-value) + (setq ffi-var var) + (setq func-var var)) + ;; Pass-by-reference variables + (t + (cond + ((member type '(:complex-single-float :complex-double-float)) + (setq ffi-var (scat "ADDR-REAL-CAST-" var)) + (setq func-var var) + (nconsc ref-vars + `((,func-var (complex (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 0) + (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 1)))))) + (t + (setq ffi-var (scat "REF-" var)) + (setq func-var var) + (nconsc ref-vars + `((,func-var (cffi:mem-aref ,ffi-var ,(%f77.cffi-type type))))))))) ;; - ((eq style :input-value) - (setq ffi-var var) - (setq func-var var)) - ;; Pass-by-reference variables - (t - (cond - ((member type '(:complex-single-float :complex-double-float)) - (setq ffi-var (scat "ADDR-REAL-CAST-" var)) - (setq func-var var) - (nconsc ref-vars - `((,func-var (complex (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 0) - (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 1)))))) - (t - (setq ffi-var (scat "REF-" var)) - (setq func-var var) - (nconsc ref-vars - `((,func-var (cffi:mem-aref ,ffi-var ,(%f77.cffi-type type))))))))) - ;; - (nconsc new-pars `((,ffi-var ,(%f77.get-read-in-type type style)))) - (nconsc func-pars `(,func-var)) - (when (and (%f77.output-p style) (not (eq type :string))) - (nconsc return-vars - `((,func-var ,ffi-var ,type))))))) - - (let ((retvar (gensym))) - `( - ,(recursive-append - `(cffi:defcallback ,callback-name ,(%f77.get-return-type hack-return-type) + (nconsc new-pars `((,ffi-var ,(%f77.get-read-in-type type style)))) + (nconsc func-pars `(,func-var)) + (when (and (%f77.output-p style) (not (eq type :string))) + (nconsc return-vars + `((,func-var ,ffi-var ,type))))))) + + (let ((retvar (gensym))) + `( + ,(recursive-append + `(cffi:defcallback ,callback-name ,(%f77.get-return-type hack-return-type) (,@new-pars ,@aux-pars)) - ;; - (when ref-vars - `(let (,@ref-vars))) - ;; - (when array-vars - `(let (,@array-vars))) - ;; - `(multiple-value-bind (,retvar ,@(mapcar #'car return-vars)) (funcall ,func ,@func-pars) - ,@(when (eq hack-return-type :void) - `((declare (ignore ,retvar)))) - ,@(mapcar #'(lambda (decl) - (destructuring-bind (func-var ffi-var type) decl - (if (member type '(:complex-single-float :complex-double-float)) - `(setf (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 0) (realpart ,func-var) - (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 1) (imagpart ,func-var)) - `(setf (cffi:mem-aref ,ffi-var ,(%f77.cffi-type type)) ,func-var)))) - (remove-if-not #'(lambda (x) - (member (first x) ref-vars :key #'car)) - return-vars)) - ,(if (eq hack-return-type :void) - nil - retvar)))))))) -) + ;; + (when ref-vars + `(let (,@ref-vars))) + ;; + (when array-vars + `(let (,@array-vars))) + ;; + `(multiple-value-bind (,retvar ,@(mapcar #'car return-vars)) (funcall ,func ,@func-pars) + ,@(when (eq hack-return-type :void) + `((declare (ignore ,retvar)))) + ,@(mapcar #'(lambda (decl) + (destructuring-bind (func-var ffi-var type) decl + (if (member type '(:complex-single-float :complex-double-float)) + `(setf (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 0) (realpart ,func-var) + (cffi:mem-aref ,ffi-var ,(second (%f77.cffi-type type)) 1) (imagpart ,func-var)) + `(setf (cffi:mem-aref ,ffi-var ,(%f77.cffi-type type)) ,func-var)))) + (remove-if-not #'(lambda (x) + (member (first x) ref-vars :key #'car)) + return-vars)) + ,(if (eq hack-return-type :void) + nil + retvar)))))))) + ) (defmacro def-fortran-routine (name-and-options return-type &rest body) " @@ -542,7 +548,7 @@ ,@pars)) (setq hack-return-type :void))) - `(progn + `(eval-when (:compile-toplevel :load-toplevel :execute) (cffi:defcfun (,fortran-name ,lisp-name) ,(%f77.get-return-type hack-return-type) ,@(%f77.parse-fortran-parameters hack-body)) ,@(%f77.def-fortran-interface name hack-return-type hack-body hidden-var-name))))) diff --git a/src/packages/odepack/dlsode.lisp b/src/packages/odepack/dlsode.lisp index 323c76e..a18b957 100644 --- a/src/packages/odepack/dlsode.lisp +++ b/src/packages/odepack/dlsode.lisp @@ -19,14 +19,14 @@ (neq :integer :input) (y (* :double-float) :input-output)) - (def-fortran-routine dlsode :void "DLSODE in ODEPACK" - (field (:callback :void - (c-neq :integer :input) - (c-t :double-float :input) - (c-y (* :double-float :size c-neq) :input) - (c-ydot (* :double-float :size c-neq) :output))) + (field (:callback dlsode-field-callback + :void + (c-neq :integer :input) + (c-t :double-float :input) + (c-y (* :double-float :size c-neq) :input) + (c-ydot (* :double-float :size c-neq) :output))) (neq :integer :input) (y (* :double-float) :input-output) (ts :double-float :input-output) @@ -41,14 +41,15 @@ (lrw :integer :input) (iwork (* :integer) :input-output) (liw :integer :input) - (jacobian (:callback :void - (c-neq :integer :input) - (c-t :double-float :input) - (c-y (* :double-float :size c-neq) :input) - (c-upper-bandwidth :integer :input) - (c-lower-bandwidth :integer :input) - (c-pd (* :double-float :size (* c-neq c-neq)) :output) - (c-nrowpd :integer :input))) + (jacobian (:callback dlsode-jacobian-callback + :void + (c-neq :integer :input) + (c-t :double-float :input) + (c-y (* :double-float :size c-neq) :input) + (c-upper-bandwidth :integer :input) + (c-lower-bandwidth :integer :input) + (c-pd (* :double-float :size (* c-neq c-neq)) :output) + (c-nrowpd :integer :input))) (mf :integer :input)) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; @@ -61,7 +62,7 @@ (ts (aref t-array 0)) (tout (aref t-array 0)) (itol 1) - (atol (make-array 1 :element-type 'double-float :initial-element 1d-8)) + (atol (make-array 1 :element-type 'double-float :initial-element 1d-12)) (rtol (make-array 1 :element-type 'double-float :initial-element 1d-12)) (itask 1) (istate 1) @@ -98,19 +99,29 @@ (fv-ref ydot 3) (/ (+ (* 2 (sin theta)) (* (cos theta) (sin theta) (expt thetadot 2))) (- (expt (cos theta) 2) 2)))))) (defun pcart-report (ts y) - (declare (ignore ts y)) - #+nil (format t "~A ~A ~A ~A ~A ~%" ts (aref y 0) (aref y 1) (aref y 2) (aref y 3))) #+nil (let ((y (make-array 2 :element-type 'double-float :initial-contents `(,(/ pi 2) 0d0)))) (lsode-evolve #'pend-field y #(0d0 1d0 2d0) #'pend-report)) +;; Should return +;; 1.0d0 1.074911802207049d0 -0.975509986605856d0 +;; 2.0d0 -0.20563950412081608d0 -1.3992359518735706d0 + +#+nil (let ((y (make-array 4 :element-type 'double-float :initial-contents `(0d0 ,(/ pi 3) 0d0 0d0))) - (ts (make-array 200 :element-type 'double-float :initial-contents (loop :for i :from 0 :below 200 + (ts (make-array 11 :element-type 'double-float :initial-contents (loop :for i :from 0 :to 10 :collect (coerce i 'double-float))))) (time (lsode-evolve #'pcart-field y ts #'pcart-report))) - ;; Should return -;; 1.0d0 1.074911802207049d0 -0.975509986605856d0 -;; 2.0d0 -0.20563950412081608d0 -1.3992359518735706d0 +;; 1.0 0.1794874587619304 0.5317592902679298 0.46247076075331184 -1.073122136936815 +;; 2.0 0.7546873547963733 -0.6988651690131225 0.3317944848774514 -0.8667875783856668 +;; 3.0 0.8622986451947015 -1.032477629437877 -0.04382530605214124 0.1709611368887168 +;; 4.0 0.5946940935220925 -0.32928104576674966 -0.6014830191620928 1.2712646406898929 +;; 5.0 0.06363367930511842 0.8312257489444103 -0.22610367475002707 0.6709599188446416 +;; 6.0 0.015513214905789278 0.9881309198657862 0.09468174702976984 -0.3441398956696562 +;; 7.0 0.38441646507049126 0.09734614673612564 0.6971343020352996 -1.4009010238635045 +;; 8.0 0.8340714151347116 -0.9308326125800139 0.14520534336799193 -0.4863145646496502 +;; 9.0 0.8288509626732954 -0.913548305194868 -0.15954504310515494 0.5222932357142753 +;; 10.0 0.36071471240686237 0.14510464752601557 -0.6851579582975589 1.3848698186034156 ----------------------------------------------------------------------- Summary of changes: src/ffi/f77-ffi.lisp | 466 +++++++++++++++++++------------------- src/packages/odepack/dlsode.lisp | 53 +++-- 2 files changed, 268 insertions(+), 251 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-01-27 08:30:04
|
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, tensor has been updated via 23ae23fc820468a4dc5149850f7f85bd524ebd69 (commit) via d7359f74d27661c0f8674d139ff196325426e828 (commit) from dca1995e47f24342e2a8f12fdcbd4b667b76a67e (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 23ae23fc820468a4dc5149850f7f85bd524ebd69 Author: Akshay Srinivasan <aks...@gm...> Date: Sat Jan 26 19:40:55 2013 -0800 Reduced the default call bounds for L3 stuff. diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp index 60b3d8d..5113980 100644 --- a/src/base/standard-tensor.lisp +++ b/src/base/standard-tensor.lisp @@ -113,7 +113,7 @@ (assert (= (rank old) 1) nil 'tensor-not-vector :rank (rank old))) ;; -(defparameter *tensor-class-optimizations* (make-hash-table) +(defvar *tensor-class-optimizations* (make-hash-table) " Contains a either: o A property list containing: diff --git a/src/base/tweakable.lisp b/src/base/tweakable.lisp index a8cc8cf..536af42 100644 --- a/src/base/tweakable.lisp +++ b/src/base/tweakable.lisp @@ -28,7 +28,7 @@ Default set with SBCL on x86-64 linux. A reasonable value is something between 400 and 1000.") ;;Level 3--------------------------------------------------------;; -(defparameter *real-l3-fcall-lb* 100 +(defparameter *real-l3-fcall-lb* 50 " If the maximum dimension in the MM is lower than this parameter, then the lisp code is used by default, instead of @@ -37,7 +37,7 @@ Default set with SBCL on x86-64 linux. A reasonable value is something between 20 and 200.") -(defparameter *complex-l3-fcall-lb* 60 +(defparameter *complex-l3-fcall-lb* 30 " If the maximum dimension in the MM is lower than this parameter, then the lisp code is used by default, instead of commit d7359f74d27661c0f8674d139ff196325426e828 Author: Akshay Srinivasan <aks...@gm...> Date: Sat Jan 26 19:17:14 2013 -0800 Fix deprecated comments. diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp index 6c05e2c..033d36e 100644 --- a/src/level-2/gemv.lisp +++ b/src/level-2/gemv.lisp @@ -3,9 +3,6 @@ (defmacro generate-typed-gemv! (func (tensor-class blas-gemv-func fortran-call-lb)) - ;;Be very careful when using functions generated by this macro. - ;;Indexes can be tricky and this has no safety net. - ;;Use only after checking the arguments for compatibility. (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)) diff --git a/src/level-3/gemm.lisp b/src/level-3/gemm.lisp index 9ed038c..d6597e0 100644 --- a/src/level-3/gemm.lisp +++ b/src/level-3/gemm.lisp @@ -42,7 +42,6 @@ (declare (type ,(getf opt :element-type) alpha beta) (type ,matrix-class A B C) (type symbol job)) - ;;The big done-in-lisp-gemm, loop-ordering was inspired by the BLAS dgemm reference implementation. ,(let ((lisp-routine `(let-typed ((nr-C (nrows C) :type index-type) @@ -60,7 +59,7 @@ (rstp-C (row-stride C) :type index-type) (cstp-C (col-stride C) :type index-type) (hd-C (head C) :type index-type)) - ;;Replace with separate loops to maximize Row-ordered MM performance + ;; (when (eq job-A :t) (rotatef rstp-A cstp-A)) (when (eq job-B :t) ----------------------------------------------------------------------- Summary of changes: src/base/standard-tensor.lisp | 2 +- src/base/tweakable.lisp | 4 ++-- src/level-2/gemv.lisp | 3 --- src/level-3/gemm.lisp | 3 +-- 4 files changed, 4 insertions(+), 8 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-01-27 03:17:36
|
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, tensor has been updated via dca1995e47f24342e2a8f12fdcbd4b667b76a67e (commit) from 7dc82f982eff648874c67454a6f6e5a947104f3d (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 dca1995e47f24342e2a8f12fdcbd4b667b76a67e Author: Akshay Srinivasan <aks...@gm...> Date: Sat Jan 26 19:14:13 2013 -0800 Whitespace cleanup. diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp index babc897..6c05e2c 100644 --- a/src/level-2/gemv.lisp +++ b/src/level-2/gemv.lisp @@ -119,10 +119,10 @@ Purpose ======= Performs the GEneral Matrix Vector operation given by - -- - - + -- - - + + Y <- alpha * op(A) * x + beta * y - Y <- alpha * op(A) * x + beta * y - and returns y. alpha,beta are scalars, @@ -217,7 +217,7 @@ ======= Returns the GEneral Matrix Vector operation given by - alpha * op(A) * x + beta * y + alpha * op(A) * x + beta * y alpha,beta are scalars, A is a matrix, and x,y are vectors. ----------------------------------------------------------------------- Summary of changes: src/level-2/gemv.lisp | 8 ++++---- 1 files changed, 4 insertions(+), 4 deletions(-) hooks/post-receive -- matlisp |
From: Akshay S. <ak...@us...> - 2013-01-27 03:14:32
|
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, tensor has been updated via 7dc82f982eff648874c67454a6f6e5a947104f3d (commit) from 8230a46f93849cdb60ce09173ab31687f998d07d (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 7dc82f982eff648874c67454a6f6e5a947104f3d Author: Akshay Srinivasan <aks...@gm...> Date: Sat Jan 26 19:11:06 2013 -0800 Fixed bug with (:c :h) handling in GEMV. diff --git a/src/level-2/gemv.lisp b/src/level-2/gemv.lisp index c7ad77a..babc897 100644 --- a/src/level-2/gemv.lisp +++ b/src/level-2/gemv.lisp @@ -91,15 +91,16 @@ (type symbol job)) (if (member job '(:n :t)) (complex-base-typed-gemv! alpha A x beta y job) - ;;The CBLAS way. - (let-typed ((cx (let-typed ((ret (apply #'make-real-tensor (lvec->list (dimensions x))) :type complex-vector)) - (complex-typed-axpy! #c(-1d0 0d0) x ret)) - :type complex-vector)) - (complex-typed-num-scal! #c(-1d0 0d0) (tensor-realpart~ y)) + ;; + (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))) - (complex-typed-num-scal! #c(-1d0 0d0) (tensor-realpart~ y)) - y))) + (real-typed-num-scal! -1d0 y.view))) + y) ;;Symbolic #+maxima ----------------------------------------------------------------------- Summary of changes: src/level-2/gemv.lisp | 15 ++++++++------- 1 files changed, 8 insertions(+), 7 deletions(-) hooks/post-receive -- matlisp |