SOLVE/NSOLVE implementations for symmetric and hermitian matrices.

Uses the Cholesky decomposition, so the matrices have to be positive
definite. May need to add logic to handle other types of symmetric/
hermitian matrices.

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)

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)

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)))

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)))