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