share/lapack/eigensys.lisp

   1 (in-package :maxima)
   2
   3 (defun maxima-matrix-dims (a)
   4   (let ((row (second a)))
   5     ;; Should we check that all rows have the same length?
   6     (values (length (cdr a)) (length (cdr row)))))
   7
   8 (defun complex-maxima-matrix-p (a)
   9   (dolist (row (cdr a))
  10     (dolist (col (cdr row))
  11       (unless (eql ($imagpart col) 0)
  12         (return-from complex-maxima-matrix-p t))))
  13   nil)
  14
  15 (defun lapack-lispify-matrix (a nrow ncol &optional (assume-complex-maxima-matrix-p nil))
  16   "Convert a Maxima matrix A of dimension NROW and NCOL to Lisp matrix
  17   suitable for use with LAPACK"
  18   (setq a ($float a))
  19   (let* ((array-type (if (or assume-complex-maxima-matrix-p (complex-maxima-matrix-p a))
  20                          '(complex flonum)
  21                          'flonum))
  22          (mat (make-array (* nrow ncol)
  23                           :element-type array-type))
  24          (mat-2d (make-array (list ncol nrow)
  25                              :element-type array-type
  26                              :displaced-to mat))
  27          (r 0))
  28     (dolist (row (cdr a))
  29       (let ((c 0))
  30         (dolist (col (cdr row))
  31           ;; Fortran matrices are in column-major order!
  32           (setf (aref mat-2d c r) (if (eql array-type 'flonum)
  33                                       (coerce col 'flonum)
  34                                       (coerce (complex ($realpart col) ($imagpart col))
  35                                               '(complex flonum))
  36                                       ))
  37           (incf c)))
  38       (incf r))
  39     mat))
  40
  41 (defun lapack-maxify-matrix (nrow ncol a)
  42   "Convert an LAPACK matrix of dimensions NROW and NCOL into a Maxima
  43 matrix (list of lists)"
  44   (let ((2d (make-array (list ncol nrow) :element-type (array-element-type a)
  45                         :displaced-to a)))
  46     (let (res)
  47       (dotimes (r nrow)
  48         (let (row)
  49           (dotimes (c ncol)
  50             ;; Fortran arrays are column-major order!
  51             (let ((v (aref 2d c r)))
  52               (push (add (realpart v) (mul '$%i (imagpart v))) row)))
  53           (push `((mlist) ,@(nreverse row)) res)))
  54       `(($matrix) ,@(nreverse res)))))
  55
  56 (defun $dgeev (a &optional right-vec-p left-vec-p)
  57   "
  58 DGEEV computes for an N-by-N real nonsymmetric matrix A, the
  59 eigenvalues and, optionally, the left and/or right eigenvectors.
  60
  61 The right eigenvector v(j) of A satisfies
  62                  A * v(j) = lambda(j) * v(j)
  63 where lambda(j) is its eigenvalue.
  64 The left eigenvector u(j) of A satisfies
  65               u(j)**H * A = lambda(j) * u(j)**H
  66 where u(j)**H denotes the conjugate transpose of u(j).
  67
  68 The computed eigenvectors are normalized to have Euclidean norm
  69 equal to 1 and largest component real.
  70
  71 A list of three items is returned.  The first item is a list of the
  72 eigenvectors.  The second item is false or the matrix of right
  73 eigenvectors.  The last itme is false or the matrix of left
  74 eigenvectors."
  75   (flet ((make-eigval (wr wi)
  76            `((mlist) ,@(map 'list #'(lambda (r i)
  77                                       (add r (mul '$%i i)))
  78                             wr wi)))
  79          (make-eigvec (n vr wi)
  80            ;; dgeev stores the eigen vectors in a special way.  Undo
  81            ;; that.  For simplicity, we create a 2-D matrix and store
  82            ;; the eigenvectors there.  Then we convert that matrix
  83            ;; into a form that maxima wants.  Somewhat inefficient.
  84            (let ((evec (make-array (list n n))))
  85              (do ((col 0 (incf col))
  86                   (posn 0))
  87                  ((>= col n))
  88                (cond ((zerop (aref wi col))
  89                       (dotimes (row n)
  90                         (setf (aref evec row col) (aref vr posn))
  91                         (incf posn)))
  92                      (t
  93                       (dotimes (row n)
  94                         (let* ((next-posn (+ posn n))
  95                                (val+ (add (aref vr posn)
  96                                           (mul '$%i (aref vr next-posn))))
  97                                (val- (sub (aref vr posn)
  98                                           (mul '$%i (aref vr next-posn)))))
  99                           (setf (aref evec row col) val+)
 100                           (setf (aref evec row (1+ col)) val-)
 101                           (incf posn)))
 102                       ;; Skip over the next column, which we've already used
 103                       (incf col)
 104                       (incf posn n))))
 105              ;; Now convert this 2-D Lisp matrix into a maxima matrix
 106              (let (res)
 107                (dotimes (r n)
 108                  (let (row)
 109                    (dotimes (c n)
 110                      (push (aref evec r c) row))
 111                    (push `((mlist) ,@(nreverse row)) res)))
 112                `(($matrix) ,@(nreverse res)))
 113              )))
 114
 115     (let* ((n (maxima-matrix-dims a))
 116            (a-mat (lapack-lispify-matrix a n n))
 117            (wr (make-array n :element-type 'flonum))
 118            (wi (make-array n :element-type 'flonum))
 119            (vl (make-array (if left-vec-p (* n n) 0)
 120                            :element-type 'flonum))
 121            (vr (make-array (if right-vec-p (* n n) 0)
 122                            :element-type 'flonum)))
 123       ;; XXX: FIXME: We need to do more error checking in the calls to
 124       ;; dgeev!
 125       (multiple-value-bind (z-jobvl z-jobvr z-n z-a z-lda z-wr z-wi z-vl
 126                                     z-ldvl z-vr z-ldvr z-work z-lwork info)
 127           ;; Figure out how much space we need in the work array.
 128           (lapack:dgeev (if left-vec-p "V" "N")
 129                         (if right-vec-p "V" "N")
 130                         n a-mat n wr wi vl n vr n wr -1 0)
 131         (declare (ignore z-jobvl z-jobvr z-n z-a z-lda z-wr z-wi z-vl
 132                          z-ldvl z-vr z-ldvr z-work z-lwork info))
 133         (let* ((opt-lwork (truncate (aref wr 0)))
 134                (work (make-array opt-lwork :element-type 'flonum)))
 135           ;; Now do the work with the optimum size of the work space.
 136           (multiple-value-bind (z-jobvl z-jobvr z-n z-a z-lda z-wr z-wi z-vl
 137                                         z-ldvl z-vr z-ldvr z-work z-lwork info)
 138               (lapack:dgeev (if left-vec-p "V" "N")
 139                             (if right-vec-p "V" "N")
 140                             n a-mat n wr wi vl n vr n work opt-lwork 0)
 141             (declare (ignore z-jobvl z-jobvr z-n z-a z-lda z-wr z-wi z-vl
 142                              z-ldvl z-vr z-ldvr z-work z-lwork))
 143             (cond ((< info 0)
 144                    (merror "DGEEV: invalid arguments: ~D" info))
 145                   ((> info 0)
 146                    (merror "DGEEV: failed to converge: ~D" info)))
 147             ;; Convert wr+%i*wi to maxima form
 148             #+nil
 149             (progn
 150               (format t "info = ~A~%" info)
 151               (format t "lwork = ~A~%" (aref work 0))
 152               (format t "vr = ~A~%" vr))
 153             (let ((e-val (make-eigval wr wi))
 154                   (e-vec-right (if right-vec-p
 155                                    (make-eigvec n vr wi)
 156                                    nil))
 157                   (e-vec-left (if left-vec-p
 158                                   (make-eigvec n vl wi)
 159                                   nil)))
 160               `((mlist) ,e-val ,e-vec-right ,e-vec-left))))))))
 161
 162 (defun $dgesvd (a &optional jobu jobvt)
 163   "
 164 DGESVD computes the singular value decomposition (SVD) of a real
 165 M-by-N matrix A, optionally computing the left and/or right singular
 166 vectors. The SVD is written
 167
 168      A = U * SIGMA * transpose(V)
 169
 170 where SIGMA is an M-by-N matrix which is zero except for its
 171 min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
 172 V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
 173 are the singular values of A; they are real and non-negative, and
 174 are returned in descending order.  The first min(m,n) columns of
 175 U and V are the left and right singular vectors of A.
 176
 177 Note that the routine returns V**T, not V.
 178
 179 A list of three items is returned.  The first is a list containing the
 180 non-zero elements of SIGMA.  If jobu is not false, The second element
 181 is the matrix U.  Otherwise it is false.  Similarly, the third element
 182 is V**T or false, depending on jobvt."
 183
 184   (flet ((maxify-vector (v)
 185            `((mlist) ,@(coerce v 'list)))
 186          (fixup-jobu (arg)
 187            (if arg "All columns of U" "No columns of U"))
 188          (fixup-jobvt (arg)
 189            (if arg "All columns of V^T" "No columns of V^T")))
 190     (multiple-value-bind (nrow ncol)
 191         (maxima-matrix-dims a)
 192       (let* ((a-mat (lapack-lispify-matrix a nrow ncol))
 193              (s (make-array (min nrow ncol) :element-type 'flonum))
 194              (u (make-array (* nrow nrow) :element-type 'flonum))
 195              (u1 (make-array (list nrow nrow) :element-type 'flonum
 196                              :displaced-to u))
 197              (vt (make-array (* ncol ncol)
 198                              :element-type 'flonum))
 199              (vt1 (make-array (list ncol ncol) :element-type 'flonum
 200                               :displaced-to vt))
 201              (wr (make-array 1 :element-type 'flonum)))
 202         ;; XXX: FIXME: We need to do more error checking in the calls to
 203         ;; dgesvd!
 204         (multiple-value-bind (z-jobu z-jobvt z-m z-n z-a z-lda z-s z-u z-ldu
 205                                      z-vt z-ldvt z-work z-lwork info)
 206             ;; Figure out the optimum size for the work array
 207             (lapack::dgesvd (fixup-jobu jobu)
 208                             (fixup-jobvt jobvt)
 209                             nrow ncol a-mat nrow
 210                             s u nrow
 211                             vt ncol
 212                             wr -1
 213                             0)
 214           (declare (ignore z-jobu z-jobvt z-m z-n z-a z-lda z-s z-u z-ldu
 215                            z-vt z-ldvt z-work z-lwork info))
 216           (let* ((opt-lwork (truncate (aref wr 0)))
 217                  (work (make-array opt-lwork :element-type 'flonum)))
 218             ;; Allocate the optimum work array and do the requested
 219             ;; computation.
 220             (multiple-value-bind (z-jobu z-jobvt z-m z-n z-a z-lda z-s z-u
 221                                          z-ldu z-vt z-ldvt z-work z-lwork info)
 222                 (lapack::dgesvd (fixup-jobu jobu)
 223                                 (fixup-jobvt jobvt)
 224                                 nrow ncol a-mat nrow
 225                                 s u nrow
 226                                 vt ncol
 227                                 work opt-lwork
 228                                 0)
 229               (declare (ignore z-jobu z-jobvt z-m z-n z-a z-lda z-s z-u z-ldu
 230                                z-vt z-ldvt z-work z-lwork))
 231               (cond ((< info 0)
 232                      (merror "DGESVD: invalid arguments: ~D" info))
 233                     ((> info 0)
 234                      (merror "DGESVD: failed to converge: ~D" info)))
 235               (let ((u-max (if jobu
 236                                (lapack-maxify-matrix nrow nrow u1)
 237                                nil))
 238                     (vt-max (if jobvt
 239                                 (lapack-maxify-matrix ncol ncol vt1)
 240                                 nil))
 241                     (s-max (maxify-vector s)))
 242                 `((mlist) ,s-max ,u-max ,vt-max)))))))))
 243
 244
 245 \f
 246 (defun $dlange (norm a)
 247   "
 248 DLANGE returns the value
 249
 250    DLANGE = ( max(abs(A(i,j))), NORM = '$max
 251             (
 252             ( norm1(A),         NORM = '$one_norm
 253             (
 254             ( normI(A),         NORM = '$inf_norm
 255             (
 256             ( normF(A),         NORM = '$frobenius
 257
 258 where  norm1  denotes the  one norm of a matrix (maximum column sum),
 259 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
 260 normF  denotes the  Frobenius norm of a matrix (square root of sum of
 261 squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm."
 262
 263   ;; Norm should be '$max, '$one_norm, '$inf_norm, '$frobenius
 264   (multiple-value-bind (nrows ncols)
 265       (maxima-matrix-dims a)
 266     (let* ((a-mat (lapack-lispify-matrix a nrows ncols))
 267            (norm-type (ecase norm
 268                         ($max "M")
 269                         ($one_norm "O")
 270                         ($inf_norm "I")
 271                         ($frobenius "F")))
 272            (work (make-array (if (equal norm-type "I") nrows 0)
 273                              :element-type 'flonum)))
 274       (lapack::dlange norm-type nrows ncols a-mat nrows work))))
 275
 276
 277 (defun $zlange (norm a)
 278   "
 279 DLANGE returns the value
 280
 281    DLANGE = ( max(abs(A(i,j))), NORM = '$max
 282             (
 283             ( norm1(A),         NORM = '$one_norm
 284             (
 285             ( normI(A),         NORM = '$inf_norm
 286             (
 287             ( normF(A),         NORM = '$frobenius
 288
 289 where  norm1  denotes the  one norm of a matrix (maximum column sum),
 290 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
 291 normF  denotes the  Frobenius norm of a matrix (square root of sum of
 292 squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm."
 293
 294   ;; Norm should be '$max, '$one_norm, '$inf_norm, '$frobenius
 295   (multiple-value-bind (nrows ncols)
 296       (maxima-matrix-dims a)
 297     (let* ((a-mat (lapack-lispify-matrix a nrows ncols))
 298            (norm-type (ecase norm
 299                         ($max "M")
 300                         ($one_norm "O")
 301                         ($inf_norm "I")
 302                         ($frobenius "F")))
 303            (work (make-array (if (equal norm-type "I") nrows 0)
 304                              :element-type 'flonum)))
 305       (lapack::zlange norm-type nrows ncols a-mat nrows work))))
 306
 307 (defun $zgeev (a &optional right-vec-p left-vec-p)
 308   "
 309 ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the
 310 eigenvalues and, optionally, the left and/or right eigenvectors.
 311
 312 The right eigenvector v(j) of A satisfies
 313                  A * v(j) = lambda(j) * v(j)
 314 where lambda(j) is its eigenvalue.
 315 The left eigenvector u(j) of A satisfies
 316               u(j)**H * A = lambda(j) * u(j)**H
 317 where u(j)**H denotes the conjugate transpose of u(j).
 318
 319 The computed eigenvectors are normalized to have Euclidean norm
 320 equal to 1 and largest component real.
 321
 322 A list of three items is returned.  The first item is a list of the
 323 eigenvectors.  The second item is false or the matrix of right
 324 eigenvectors.  The last itme is false or the matrix of left
 325 eigenvectors."
 326   (flet ((make-eigval (w)
 327            `((mlist) ,@(map 'list #'(lambda (z)
 328                                       (add (realpart z) (mul '$%i (imagpart z))))
 329                             w))))
 330
 331     (let* ((n (maxima-matrix-dims a))
 332            (a-mat (lapack-lispify-matrix a n n t))
 333            (w (make-array n :element-type '(complex flonum)))
 334            (vl (make-array (if left-vec-p (* n n) 0)
 335                            :element-type '(complex flonum)))
 336            (vr (make-array (if right-vec-p (* n n) 0)
 337                            :element-type '(complex flonum)))
 338            (rwork (make-array 1 :element-type 'flonum)))
 339       ;; XXX: FIXME: We need to do more error checking in the calls to
 340       ;; zgeev!
 341       (multiple-value-bind (z-jobvl z-jobvr z-n z-a z-lda z-wr z-wi z-vl
 342                                     z-ldvl z-vr z-ldvr z-work z-lwork info)
 343           ;; Figure out how much space we need in the work array.
 344           (lapack:zgeev (if left-vec-p "V" "N")
 345                         (if right-vec-p "V" "N")
 346                         n a-mat n w vl n vr n w -1 rwork 0)
 347         (declare (ignore z-jobvl z-jobvr z-n z-a z-lda z-wr z-wi z-vl
 348                          z-ldvl z-vr z-ldvr z-work z-lwork info))
 349         (let* ((opt-lwork (truncate (realpart (aref w 0))))
 350                (work (make-array opt-lwork :element-type '(complex flonum)))
 351                (rwork (make-array opt-lwork :element-type 'flonum)))
 352           ;; Now do the work with the optimum size of the work space.
 353           (multiple-value-bind (z-jobvl z-jobvr z-n z-a z-lda z-wr z-wi z-vl
 354                                         z-ldvl z-vr z-ldvr z-work z-lwork info)
 355               (lapack:zgeev (if left-vec-p "V" "N")
 356                             (if right-vec-p "V" "N")
 357                             n a-mat n w vl n vr n work opt-lwork rwork 0)
 358             (declare (ignore z-jobvl z-jobvr z-n z-a z-lda z-wr z-wi z-vl
 359                              z-ldvl z-vr z-ldvr z-work z-lwork))
 360             (cond ((< info 0)
 361                    (merror "ZGEEV: invalid arguments: ~D" info))
 362                   ((> info 0)
 363                    (merror "ZGEEV: failed to converge: ~D" info)))
 364             (let ((e-val (make-eigval w))
 365                   (e-vec-right (if right-vec-p
 366                                    (lapack-maxify-matrix n n vr)
 367                                    nil))
 368                   (e-vec-left (if left-vec-p
 369                                   (lapack-maxify-matrix n n vl)
 370                                   nil)))
 371               `((mlist) ,e-val ,e-vec-right ,e-vec-left))))))))
 372
 373 (defun $zheev (a &optional eigen-vector-p)
 374   (flet ((make-eigval (w)
 375            `((mlist) ,@(map 'list #'(lambda (z)
 376                                       (add (realpart z) (mul '$%i (imagpart z))))
 377                             w))))
 378
 379     (let* ((n (maxima-matrix-dims a))
 380            (a-mat (lapack-lispify-matrix a n n t))
 381            (w (make-array n :element-type 'flonum))
 382            (work (make-array 1 :element-type '(complex flonum)))
 383            (rwork (make-array (max 1 (- (* 3 n) 2))
 384                               :element-type 'flonum)))
 385       ;; XXX: FIXME: We need to do more error checking in the calls to
 386       ;; zgeev!
 387       (multiple-value-bind (z-jobz z-uplo z-n z-a z-lda z-w z-work
 388                             z-lwork z-rwork info)
 389           ;; Figure out how much space we need in the work array.
 390           (lapack:zheev (if eigen-vector-p "V" "N")
 391                         "U"
 392                         n
 393                         a-mat
 394                         n
 395                         w
 396                         work
 397                         -1
 398                         rwork
 399                         0)
 400         (declare (ignore z-jobz z-uplo z-n z-a z-lda z-w z-work
 401                          z-lwork z-rwork))
 402         (let* ((opt-lwork (truncate (realpart (aref work 0))))
 403                (work (make-array opt-lwork :element-type '(complex flonum))))
 404           ;; Now do the work with the optimum size of the work space.
 405           (multiple-value-bind (z-jobz z-uplo z-n z-a z-lda z-w z-work
 406                                 z-lwork z-rwork info)
 407               (lapack:zheev (if eigen-vector-p "V" "N")
 408                             "U"
 409                             n
 410                             a-mat
 411                             n
 412                             w
 413                             work
 414                             opt-lwork
 415                             rwork
 416                             0)
 417             (declare (ignore z-jobz z-uplo z-n z-a z-lda z-w z-work
 418                              z-lwork z-rwork))
 419             (cond ((< info 0)
 420                    (merror "ZHEEV: invalid arguments: ~D" info))
 421                   ((> info 0)
 422                    (merror "ZHEEV: failed to converge: ~D" info)))
 423             (let ((e-val (make-eigval w))
 424                   (e-vec (if eigen-vector-p
 425                                    (lapack-maxify-matrix n n a-mat)
 426                                    nil)))
 427               `((mlist) ,e-val ,e-vec))))))))