share/lapack/dgeqrf.lisp

   1 ;; dgeqrf.lisp -- Maxima interface to lapack::dgeqrf
   2 ;; copyright 2011 by Robert Dodier
   3 ;; I release this work under terms of the GNU General Public License.
   4 (in-package :maxima)
   5
   6 ;; dgeqrf(a) computes QR factorization of m-by-n Maxima matrix,
   7 ;; where m and n might be equal or unequal.
   8 ;; a is not modified.
   9
  10 (defun $dgeqrf (a)
  11
  12   (multiple-value-bind (a-nrow a-ncol)
  13       (maxima-matrix-dims a)
  14     (let
  15       ((a-mat (lapack-lispify-matrix a a-nrow a-ncol))
  16        (work (make-array 1 :element-type 'flonum))
  17        (tau (make-array (min a-nrow a-ncol) :element-type 'flonum)))
  18
  19        ;; First call DGEQRF with LWORK = -1 to find optimal block size NB.
  20        ;; Use that to compute LWORK.
  21
  22        (multiple-value-bind (foo-1 foo-2 foo-3 foo-4 foo-5 foo-6 foo-7 foo-8)
  23          (lapack::dgeqrf a-nrow a-ncol a-mat a-nrow tau work -1 0)
  24          (declare (ignore foo-1 foo-2 foo-3 foo-4 foo-5 foo-6 foo-7 foo-8))
  25
  26          (let*
  27            ((lwork (floor (aref work 0)))
  28             (work (make-array lwork :element-type 'flonum)))
  29
  30            (multiple-value-bind (foo-1 foo-2 foo-3 foo-4 foo-5 foo-6 foo-7 z-info)
  31              (lapack::dgeqrf a-nrow a-ncol a-mat a-nrow tau work lwork 0)
  32              (declare (ignore foo-1 foo-2 foo-3 foo-4 foo-5 foo-6 foo-7))
  33
  34              (cond
  35                ((< z-info 0)
  36                 (merror "dgeqrf: ~M-th argument has an illegal value." (- z-info)))
  37                (t
  38                  (dgeqrf-unpack a-nrow a-ncol a-mat tau)))))))))
  39
  40 (defun dgeqrf-unpack (a-nrow a-ncol a-mat tau)
  41   (let
  42     ((r (dgeqrf-unpack-r a-nrow a-ncol a-mat))
  43      (q (dgeqrf-unpack-q a-nrow a-ncol a-mat tau)))
  44     (list '(mlist) q r)))
  45
  46 (defun dgeqrf-unpack-r (a-nrow a-ncol a-mat)
  47   (let*
  48     ((r-nrow a-nrow)
  49      (r-ncol a-ncol)
  50      (r-mat (make-array (* r-ncol r-nrow) :element-type 'flonum
  51                                           :initial-element 0e0)))
  52     (dotimes (j a-ncol)
  53       (dotimes (i (1+ j))
  54         (if (< i r-nrow)
  55           (let
  56             ((a-ij (+ (* j a-nrow) i))
  57              (r-ij (+ (* j r-nrow) i)))
  58             (setf (aref r-mat r-ij) (aref a-mat a-ij))))))
  59     (lapack-maxify-matrix r-nrow r-ncol r-mat)))
  60
  61 (defun dgeqrf-unpack-q (a-nrow a-ncol a-mat tau)
  62   (let ((q-mat (make-array (* a-nrow a-nrow) :element-type 'flonum
  63                                              :initial-element 0e0)))
  64     (dotimes (i a-nrow)
  65       (setf (aref q-mat (+ (* i a-nrow) i)) 1e0))
  66     (dotimes (i (min a-nrow a-ncol))
  67       (let ((h-mat-i (dgeqrf-h i tau a-nrow a-mat)))
  68         (dgeqrf-multiply-into a-nrow q-mat h-mat-i)))
  69     (lapack-maxify-matrix a-nrow a-nrow q-mat)))
  70
  71 (defun dgeqrf-h (i tau a-nrow a-mat)
  72   (let ((h-mat (make-array (* a-nrow a-nrow) :element-type 'flonum)))
  73     (dotimes (j a-nrow)
  74       (dotimes (k a-nrow)
  75         (let*
  76           ((v-j (dgeqrf-v i j a-nrow a-mat))
  77            (v-k (dgeqrf-v i k a-nrow a-mat))
  78            (foo (* (aref tau i) v-j v-k))
  79            (bar (if (= j k) (1+ (- foo)) (- foo))))
  80           (setf (aref h-mat (+ (* k a-nrow) j)) bar))))
  81     h-mat))
  82
  83 (defun dgeqrf-v (i j a-nrow a-mat)
  84   (cond
  85     ((= j i) 1)
  86     ((> j i) (aref a-mat (+ (* i a-nrow) j)))
  87     (t 0)))
  88
  89 (defun dgeqrf-multiply-into (n mat-1 mat-2)
  90   (let ((row (make-array n :element-type 'flonum :initial-element 0e0)))
  91     (dotimes (i n)
  92       (dotimes (j n)
  93         (setf (aref row j) (dgeqrf-inner-product n mat-1 mat-2 i j)))
  94       (dotimes (j n)
  95         (setf (aref mat-1 (+ (* j n) i)) (aref row j))))))
  96
  97 (defun dgeqrf-inner-product (n mat-1 mat-2 i j)
  98   (let ((sum 0))
  99     (dotimes (k n)
 100       (setq sum (+ sum (* (aref mat-1 (+ (* k n) i)) (aref mat-2 (+ (* j n) k))))))
 101     sum))