Function: math-nlfit-givens

math-nlfit-givens is a byte-compiled function defined in calc-nlfit.el.gz.

Signature

(math-nlfit-givens C D)

Source Code

;; Defined in /usr/src/emacs/lisp/calc/calc-nlfit.el.gz
;;; For solving the linearized system.
;;; (The Givens method, from Schwarz.)

(defun math-nlfit-givens (C d)
  (let* ((C (copy-tree C))
         (d (copy-tree d))
         (n (length (car C)))
         (N (length C))
         (j 1)
         (r (make-list N 0))
         (x (make-list N 0))
         w
         gamma
         sigma
         rho)
    (while (<= j n)
      (let ((i (1+ j)))
        (while (<= i N)
          (let ((cij (math-nlfit-get-matx-elt C i j))
                (cjj (math-nlfit-get-matx-elt C j j)))
            (when (not (math-equal 0 cij))
                (if (math-lessp (calcFunc-abs cjj)
                                (math-mul math-nlfit-delta (calcFunc-abs cij)))
                    (setq w (math-neg cij)
                          gamma 0
                          sigma 1
                          rho 1)
                  (setq w (math-mul
                           (calcFunc-sign cjj)
                           (calcFunc-sqrt
                            (math-add
                             (math-mul cjj cjj)
                             (math-mul cij cij))))
                        gamma (math-div cjj w)
                        sigma (math-neg (math-div cij w)))
                  (if (math-lessp (calcFunc-abs sigma) gamma)
                      (setq rho sigma)
                    (setq rho (math-div (calcFunc-sign sigma) gamma))))
              (setq cjj w
                    cij rho)
              (math-nlfit-set-matx-elt C j j w)
              (math-nlfit-set-matx-elt C i j rho)
              (let ((k (1+ j)))
                (while (<= k n)
                  (let* ((cjk (math-nlfit-get-matx-elt C j k))
                         (cik (math-nlfit-get-matx-elt C i k))
                         (h (math-sub
                             (math-mul gamma cjk) (math-mul sigma cik))))
                    (setq cik (math-add
                               (math-mul sigma cjk)
                               (math-mul gamma cik)))
                    (setq cjk h)
                    (math-nlfit-set-matx-elt C i k cik)
                    (math-nlfit-set-matx-elt C j k cjk)
                    (setq k (1+ k)))))
              (let* ((di (math-nlfit-get-elt d i))
                     (dj (math-nlfit-get-elt d j))
                     (h (math-sub
                         (math-mul gamma dj)
                         (math-mul sigma di))))
                (setq di (math-add
                          (math-mul sigma dj)
                          (math-mul gamma di)))
                (setq dj h)
                (math-nlfit-set-elt d i di)
                (math-nlfit-set-elt d j dj))))
          (setq i (1+ i))))
      (setq j (1+ j)))
    (let ((i n)
          s)
      (while (>= i 1)
        (math-nlfit-set-elt r i 0)
        (setq s (math-nlfit-get-elt d i))
        (let ((k (1+ i)))
          (while (<= k n)
            (setq s (math-add s (math-mul (math-nlfit-get-matx-elt C i k)
                            (math-nlfit-get-elt x k))))
            (setq k (1+ k))))
        (math-nlfit-set-elt x i
                            (math-neg
                             (math-div s
                                       (math-nlfit-get-matx-elt C i i))))
        (setq i (1- i))))
    (let ((i (1+ n)))
      (while (<= i N)
        (math-nlfit-set-elt r i (math-nlfit-get-elt d i))
        (setq i (1+ i))))
    (let ((j n))
      (while (>= j 1)
        (let ((i N))
          (while (>= i (1+ j))
            (setq rho (math-nlfit-get-matx-elt C i j))
            (if (math-equal rho 1)
                (setq gamma 0
                      sigma 1)
              (if (math-lessp (calcFunc-abs rho) 1)
                  (setq sigma rho
                        gamma (calcFunc-sqrt
                               (math-sub 1 (math-mul sigma sigma))))
                (setq gamma (math-div 1 (calcFunc-abs rho))
                      sigma (math-mul (calcFunc-sign rho)
                                       (calcFunc-sqrt
                                        (math-sub 1 (math-mul gamma gamma)))))))
            (let ((ri (math-nlfit-get-elt r i))
                  (rj (math-nlfit-get-elt r j))
                  h)
              (setq h (math-add (math-mul gamma rj)
                         (math-mul sigma ri)))
              (setq ri (math-sub
                        (math-mul gamma ri)
                        (math-mul sigma rj)))
              (setq rj h)
              (math-nlfit-set-elt r i ri)
              (math-nlfit-set-elt r j rj))
            (setq i (1- i))))
        (setq j (1- j))))

    x))