Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch tensor updated. 2012-02-23-199-g1367160
|
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
|
