diff --git a/lisp/hermitian-matrix.lisp b/lisp/hermitian-matrix.lisp index 2839f34..85795cd 100644 --- a/lisp/hermitian-matrix.lisp +++ b/lisp/hermitian-matrix.lisp @@ -384,3 +384,19 @@ matrix2." "Permutation matrix~A and dense matrix~A sizes incompatible." (matrix-dimensions matrix) (matrix-dimensions permutation)))) + +(defmethod solve ((matrix hermitian-matrix) (vector column-vector)) + "Return the solution to the system of equations." + (make-instance + 'column-vector + :contents + (hermitian-cholesky-solver + (copy-array (contents matrix)) (copy-array (contents vector))))) + +(defmethod nsolve ((matrix hermitian-matrix) (vector column-vector)) + "Return the solution to the system of equations." + (setf + (contents vector) + (hermitian-cholesky-solver (contents matrix) (contents vector))) + ;; Return the solution vector + vector) diff --git a/lisp/symmetric-matrix.lisp b/lisp/symmetric-matrix.lisp index c34238a..081b8a4 100644 --- a/lisp/symmetric-matrix.lisp +++ b/lisp/symmetric-matrix.lisp @@ -357,3 +357,19 @@ subtracted to a symmetric matrix." (nsubtract-array (contents matrix1) (contents matrix2) scalar1 scalar2) matrix1) + +(defmethod solve ((matrix symmetric-matrix) (vector column-vector)) + "Return the solution to the system of equations." + (make-instance + 'column-vector + :contents + (symmetric-cholesky-solver + (copy-array (contents matrix)) (copy-array (contents vector))))) + +(defmethod nsolve ((matrix symmetric-matrix) (vector column-vector)) + "Return the solution to the system of equations." + (setf + (contents vector) + (symmetric-cholesky-solver (contents matrix) (contents vector))) + ;; Return the solution vector + vector) diff --git a/test/hermitian-matrix.lisp b/test/hermitian-matrix.lisp index e088744..beffb4e 100644 --- a/test/hermitian-matrix.lisp +++ b/test/hermitian-matrix.lisp @@ -1332,3 +1332,88 @@ 'error (linear-algebra:product (unit-hermitian-matrix 3) (unit-hermitian-matrix 4)))) + +(define-test solve-hermitian-matrix + (:tag :hermitian-matrix :solve) + (let ((vector2 (linear-algebra:column-vector 2.0 1.0)) + (vector3 (linear-algebra:column-vector 2.3 1.2 2.2)) + (matrix2 + (linear-algebra:make-matrix + 2 2 :matrix-type 'linear-algebra:hermitian-matrix + :initial-contents + #2A((#C(2.0 0.0) #C(1.0 -2.0)) + (#C(1.0 2.0) #C(3.0 0.0))))) + (matrix3 + (linear-algebra:make-matrix + 3 3 :matrix-type 'linear-algebra:hermitian-matrix + :initial-contents + #2A((#C(3.31 0.0) #C(1.26 -2.0) #C(1.37 -3.0)) + (#C(1.26 2.0) #C(2.23 0.0) #C(2.31 -1.5)) + (#C(1.37 3.0) #C(2.31 1.5) #C(8.15 0.0)))))) + ;; 2x2 + (assert-float-equal + #(#C(5.0 2.0) #C(0.0 -4.0)) + (linear-algebra:solve matrix2 vector2)) + (assert-float-equal + #2A((#C(2.0 0.0) #C(1.0 -2.0)) + (#C(1.0 2.0) #C(3.0 0.0))) + matrix2) + (assert-float-equal + #(2.0 1.0) vector2 (linear-algebra::contents vector2)) + ;; 3x3 + (assert-float-equal + #(#C( 3.5175734 3.4673646) + #C( 3.3198433 -4.3366637) + #C(-0.78414906 -1.2595192)) + (linear-algebra:solve matrix3 vector3)) + (assert-float-equal + #2A((#C(3.31 0.0) #C(1.26 -2.0) #C(1.37 -3.0)) + (#C(1.26 2.0) #C(2.23 0.0) #C(2.31 -1.5)) + (#C(1.37 3.0) #C(2.31 1.5) #C(8.15 0.0))) + matrix3) + (assert-float-equal + #(2.3 1.2 2.2) vector3 (linear-algebra::contents vector3)))) + +(define-test nsolve-hermitian-matrix + (:tag :hermitian-matrix :nsolve) + (let ((vector2 (linear-algebra:column-vector 2.0 1.0)) + (vector3 (linear-algebra:column-vector 2.3 1.2 2.2)) + (matrix2 + (linear-algebra:make-matrix + 2 2 :matrix-type 'linear-algebra:hermitian-matrix + :initial-contents + #2A((#C(2.0 0.0) #C(1.0 -2.0)) + (#C(1.0 2.0) #C(3.0 0.0))))) + (matrix3 + (linear-algebra:make-matrix + 3 3 :matrix-type 'linear-algebra:hermitian-matrix + :initial-contents + #2A((#C(3.31 0.0) #C(1.26 -2.0) #C(1.37 -3.0)) + (#C(1.26 2.0) #C(2.23 0.0) #C(2.31 -1.5)) + (#C(1.37 3.0) #C(2.31 1.5) #C(8.15 0.0)))))) + ;; 2x2 + (assert-float-equal + #(#C(5.0 2.0) #C(0.0 -4.0)) + (linear-algebra:nsolve matrix2 vector2)) + (assert-float-equal + #2A((#C(2.0 0.0) #C(0.5 -1.0)) + (#C(0.5 1.0) #C(0.5 0.0))) + matrix2) + (assert-float-equal + #(#C(5.0 2.0) #C(0.0 -4.0)) vector2) + ;; 3x3 + (assert-float-equal + #(#C( 3.5175734 3.4673646) + #C( 3.3198433 -4.3366637) + #C(-0.78414906 -1.2595192)) + (linear-algebra:nsolve matrix3 vector3)) + (assert-float-equal + #2A((#C(3.31 0.0) #C( 0.38066468 -0.6042296) #C( 0.4138973 -0.9063445)) + (#C(0.38066468 0.6042296) #C( 0.54190326 0.0) #C(-0.044656467 -2.1882146)) + (#C(0.4138973 0.9063445) #C(-0.044656467 2.1882146) #C( 2.2680602 -3.7252904E-9))) + matrix3) + (assert-float-equal + #(#C( 3.5175734 3.4673646) + #C( 3.3198433 -4.3366637) + #C(-0.78414906 -1.2595192)) + vector3))) diff --git a/test/symmetric-matrix.lisp b/test/symmetric-matrix.lisp index 9cd3d6e..2660003 100644 --- a/test/symmetric-matrix.lisp +++ b/test/symmetric-matrix.lisp @@ -1244,3 +1244,69 @@ (1.2 2.2 2.3) (1.3 2.3 3.3))) 2.1))) + +(define-test solve-symmetric-matrix + (:tag :symmetric-matrix :solve) + (let ((vector2 (linear-algebra:column-vector 2.0 1.0)) + (vector3 (linear-algebra:column-vector 2.3 1.2 2.2)) + (matrix2 + (linear-algebra:make-matrix + 2 2 :matrix-type 'linear-algebra:symmetric-matrix + :initial-contents #2A((1.1 1.2) (1.2 2.2)))) + (matrix3 + (linear-algebra:make-matrix + 3 3 :matrix-type 'linear-algebra:symmetric-matrix + :initial-contents + #2A((1.15 1.26 1.37) + (1.26 2.23 2.31) + (1.37 2.31 3.31))))) + ;; 2x2 + (assert-float-equal + #(3.2653065 -1.3265308) + (linear-algebra:solve matrix2 vector2)) + (assert-float-equal + #2A((1.1 1.2) (1.2 2.2)) matrix2) + (assert-float-equal #(2.0 1.0)vector2) + ;; 3x3 + (assert-float-equal + #(3.5856622 -2.306286 0.79007966) + (linear-algebra:solve matrix3 vector3)) + (assert-float-equal + #2A((1.15 1.26 1.37) + (1.26 2.23 2.31) + (1.37 2.31 3.31)) + matrix3) + (assert-float-equal #(2.3 1.2 2.2) vector3))) + +(define-test nsolve-symmetric-matrix + (:tag :symmetric-matrix :nsolve) + (let ((vector2 (linear-algebra:column-vector 2.0 1.0)) + (vector3 (linear-algebra:column-vector 2.3 1.2 2.2)) + (matrix2 + (linear-algebra:make-matrix + 2 2 :matrix-type 'linear-algebra:symmetric-matrix + :initial-contents #2A((1.1 1.2) (1.2 2.2)))) + (matrix3 + (linear-algebra:make-matrix + 3 3 :matrix-type 'linear-algebra:symmetric-matrix + :initial-contents + #2A((1.15 1.26 1.37) + (1.26 2.23 2.31) + (1.37 2.31 3.31))))) + ;; 2x2 + (assert-float-equal + #(3.2653065 -1.3265308) + (linear-algebra:nsolve matrix2 vector2)) + (assert-float-equal + #2A((1.1 1.0909091) (1.0909091 0.8909091)) matrix2) + (assert-float-equal #(3.2653065 -1.3265308) vector2) + ;; 3x3 + (assert-float-equal + #(3.5856622 -2.306286 0.79007966) + (linear-algebra:nsolve matrix3 vector3)) + (assert-float-equal + #2A((1.15 1.0956522 1.1913043) + (1.0956522 0.84947825 0.9522979) + (1.1913043 0.9522979 0.90754557)) + matrix3) + (assert-float-equal #(3.5856622 -2.306286 0.79007966) vector3)))