Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
refactor out matrix decomp/factorization methods into smaller chunks.
Signed-off-by: AJ Rossini <blindglobe@gmail.com>
- Loading branch information
1 parent
b6a26f5
commit 445cfc2
Showing
5 changed files
with
175 additions
and
156 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,74 @@ | ||
|
|
||
| (in-package :lisp-matrix) | ||
|
|
||
| ;;; CHOLESKY | ||
|
|
||
| ;;; POTRF - compute the Cholesky Factorization of a real sym pos-def | ||
| ;;; matrix A. | ||
| ;;; Returns Matrix, upper/lower triang char, info | ||
| (def-lapack-method potrf ((a !matrix-type)) | ||
| (let ((info (make-fnv-int32 1 :initial-value 0))) | ||
| (assert (<= (ncols a) (nrows a))) ; make sure A supports options | ||
| (with-copies ((a (or (not unit-strides-p) | ||
| transposed-p) | ||
| t)) | ||
| (list a | ||
| "U" | ||
| (check-info (fnv-int32-ref info 0) "POTRF")) | ||
| (!function "U" ; store in Upper section | ||
| (ncols a) ; N | ||
| a ; matrix (in/out) | ||
| (real-nrows a) ; LDA | ||
| info)))) ; info | ||
|
|
||
| ;;; POTRI - compute the inverse of a real symmetric positive definite | ||
| ;;; matrix A using the Cholesky factorization A = U**T*U or A = L*L**T | ||
| (def-lapack-method potri ((a !matrix-type)) | ||
| (assert (= (ncols a) (nrows a))) ;; only works with square matrices | ||
| (let ((info (make-fnv-int32 1 :initial-value 0))) | ||
| (with-copies ((a (or (not unit-strides-p) | ||
| transposed-p) | ||
| t)) | ||
| ;; Returning: | ||
| ;; - inverse, | ||
| ;; - "U" since upper format trangular, | ||
| ;; - info, for correctness of results. | ||
| ;; Should we put INFO first?! | ||
| (list a | ||
| "U" ; not useful until we add option for lowercase. | ||
| (check-info (fnv-int32-ref info 0) "POTRI")) | ||
| (!function "U" ; "L" (in) is lower an option? | ||
| (ncols a) ; N (in) (order of matrix, columns, 2nd index | ||
| a ; a (in/out) matrix | ||
| (nrows a) ; LDA (in) leading dimension, LDA >= max(1,N) | ||
| ; above was "(real-nrows a)" ? | ||
| info)))) ; info (out) | ||
|
|
||
|
|
||
|
|
||
| ;;; CHOLESKY | ||
| ;; | ||
| ;; POTRI - compute the inverse of a real symmetric positive definite | ||
| ;; matrix A using the Cholesky factorization A = U**T*U or A = L*L**T | ||
| (def-lapack-method potri ((a !matrix-type)) | ||
| (assert (= (ncols a) (nrows a))) ; only square matrices | ||
| (let ((info (make-fnv-int32 1 :initial-value 0))) | ||
| (with-copies ((a (or (not unit-strides-p) | ||
| transposed-p) | ||
| t)) | ||
| ;; Returning: | ||
| ;; - inverse, | ||
| ;; - "U" since upper format trangular, | ||
| ;; - info, for correctness of results. | ||
| ;; Should we put INFO first?! | ||
| (list a | ||
| "U" ; not useful until we add option for lowercase. | ||
| (check-info (fnv-int32-ref info 0) "POTRI")) | ||
| (!function "U" ; "L" (in) is lower an option? | ||
| (ncols a) ; N (in) (order of matrix, columns, 2nd index | ||
| a ; a (in/out) matrix | ||
| (nrows a) ; LDA (in) leading dimension, LDA >= max(1,N) | ||
| ; above was "(real-nrows a)" ? | ||
| info)))) ; info (out) | ||
|
|
||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,26 @@ | ||
|
|
||
|
|
||
| ;;; LU | ||
|
|
||
| ;;; GETRF - compute the LU Factorization of a matrix. | ||
| ;;; Returns Matrix, upper/lower triang char, info | ||
| (def-lapack-method getrf ((a !matrix-type)) | ||
| (let ((info (make-fnv-int32 1 :initial-value 0))) | ||
| (assert (<= (ncols a) (nrows a))) ; make sure A supports options | ||
| (unless ipvt ; make it the bigger of # cols/ # rows | ||
| (setf ipvt (make-fnv-int32 (nrows a) :initial-value 0))) | ||
| (with-copies ((a (or (not unit-strides-p) | ||
| transposed-p) | ||
| t)) | ||
| (list a | ||
| (check-info (fnv-int32-ref info 0) "GETRF")) | ||
| (call-with-work (lwork work !data-type) | ||
| (!function (nrows a) ; N | ||
| (ncols a) ; M | ||
| a ; A | ||
| (max 1 (ncols a)); LDA | ||
| ipiv | ||
| ;; IPIV (output) INTEGER array, dimension (min(M,N)) | ||
| ;; The pivot indices; for 1 <= i <= min(M,N), row i of | ||
| ;; the matrix was interchanged with row IPIV(i). | ||
| info))))); info |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,64 @@ | ||
|
|
||
| (in-package :lisp-matrix) | ||
|
|
||
| ;;; QR | ||
|
|
||
|
|
||
|
|
||
| ;;; QR decomposition. Need one more front end to provide appropriate | ||
| ;;; processing. A and TAU will have different values at the end, more | ||
| ;;; appropriate to the transformation (i.e. will be the QR, but stored | ||
| ;;; in common compact form). | ||
| ;;; | ||
| ;; M (input) INTEGER | ||
| ;; The number of rows of the matrix A. M >= 0. | ||
| ;; N (input) INTEGER | ||
| ;; The number of columns of the matrix A. N >= 0. | ||
| ;; A (input/output) DOUBLE PRECISION array, dimension (LDA,N) | ||
| ;; On entry, the M-by-N matrix A. On exit, the elements on and | ||
| ;; above the diagonal of the array contain the min(M,N)-by-N upper | ||
| ;; trapezoidal matrix R (R is upper triangular if m >= n); the | ||
| ;; elements below the diagonal, with the array TAU, represent the | ||
| ;; orthogonal matrix Q as a product of min(m,n) elementary reflec‐ | ||
| ;; tors (see Further Details). | ||
| ;; LDA (input) INTEGER | ||
| ;; The leading dimension of the array A. LDA >= max(1,M). | ||
| ;; TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) | ||
| ;; The scalar factors of the elementary reflectors (see Further | ||
| ;; Details). | ||
| ;; WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) | ||
| ;; On exit, if INFO = 0, WORK(1) returns the optimal LWORK. | ||
| ;; LWORK (input) INTEGER | ||
| ;; The dimension of the array WORK. LWORK >= max(1,N). For opti‐ | ||
| ;; mum performance LWORK >= N*NB, where NB is the optimal block‐ | ||
| ;; size. | ||
| ;; If LWORK = -1, then a workspace query is assumed; the routine | ||
| ;; only calculates the optimal size of the WORK array, returns | ||
| ;; this value as the first entry of the WORK array, and no error | ||
| ;; message related to LWORK is issued by XERBLA. | ||
| ;; INFO (output) INTEGER | ||
| ;; = 0: successful exit | ||
| ;; < 0: if INFO = -i, the i-th argument had an illegal value | ||
|
|
||
| (def-lapack-method geqrf ((a !matrix-type) | ||
| (tau !matrix-type)) ; tau is for output | ||
| (assert (<= (ncols a) (nrows a))) ; make sure A supports options | ||
| (let ((info (make-fnv-int32 1 :initial-value 0))) | ||
| (with-copies ((a (or (not unit-strides-p) | ||
| transposed-p) | ||
| t) | ||
| (tau (or (not unit-strides-p) | ||
| transposed-p) | ||
| t)) | ||
| (list a | ||
| tau | ||
| (check-info (fnv-int32-ref info 0) "GEQRF")) | ||
| (call-with-work (lwork work !data-type) | ||
| (!function (nrows a) | ||
| (ncols a) | ||
| a | ||
| (max 1 (nrows a)) | ||
| (data tau) | ||
| (data work) | ||
| lwork | ||
| info))))) |