share/lapack/dgesv.lisp

   1 ;; dgesv.lisp -- Maxima interface to lapack::dgesv
   2 ;; copyright 2010 by Robert Dodier
   3 ;; I release this work under terms of the GNU General Public License.
   4 (in-package :maxima)
   5
   6 ;; dgesv(a, b) returns solution x of linear equations a . x = b
   7 ;; as computed by the LU decomposition.
   8 ;; a is a n-by-n Maxima matrix, b is a n-by-m Maxima matrix,
   9 ;; where m maybe be greater than or equal to 1.
  10 ;; a and b are not modified.
  11
  12 (defun $dgesv (a b)
  13
  14   (multiple-value-bind (a-nrow a-ncol) (maxima-matrix-dims a)
  15     (multiple-value-bind (b-nrow b-ncol) (maxima-matrix-dims b)
  16
  17       (let
  18         ((a-mat (lapack-lispify-matrix a a-nrow a-ncol))
  19          (b-mat (lapack-lispify-matrix b b-nrow b-ncol))
  20          (ipiv (make-array a-nrow :element-type 'f2cl-lib:integer4)))
  21
  22         (multiple-value-bind (z-n z-nrhs z-a z-lda z-ipiv z-b z-ldb$ z-info)
  23           (lapack::dgesv a-nrow b-ncol a-mat a-nrow ipiv b-mat b-nrow 0)
  24
  25           (declare (ignore z-n z-nrhs z-a z-lda z-ipiv z-b z-ldb$))
  26           (cond
  27             ((< z-info 0)
  28              (merror "dgesv: ~M-th argument has an illegal value." (- z-info)))
  29             ((> z-info 0)
  30              (merror "dgesv: U(~M, ~M) is exactly zero; cannot compute a solution." z-info z-info))
  31             (t
  32               (lapack-maxify-matrix b-nrow b-ncol b-mat))))))))