Permalink
Browse files

Some tweaks

  • Loading branch information...
pkhuong committed Dec 18, 2013
1 parent 2fcda3b commit 3a489beaa77dc3f0d95ae4822dd0544819d61754
Showing with 42 additions and 35 deletions.
  1. +32 −27 affine-scaling.lisp
  2. +10 −8 sparse-cholesky.lisp
@@ -102,14 +102,17 @@
;; Let B = [constraints][scale], N = BB', c = -[scale]c
;;
;; x = c - B'N^-1 B c
(let* ((c (matlisp:m.*! scale
(matlisp:scal -1 c)))
(B (scale-sparse! constraints scale))
(Bc (sparse-m* B c))
(N^-1Bc (solve-symmetric! B bc sparse-state)))
;; c - Bt N^-1 B c
(sparse-m* B N^-1Bc :transpose t :y c
:alpha -1d0 :beta 1d0)))
(let* ((sc (matlisp:m.*! scale
(matlisp:scal -1 c)))
(AD (scale-sparse! constraints scale))
(AD2c (sparse-m* AD sc))
(N^-1AD2c (solve-symmetric! AD AD2c
sparse-state)))
;; Dc - (AD)t N^-1 AD Dc
(values (sparse-m* AD N^-1AD2c :transpose t :y sc
:alpha -1d0 :beta 1d0
:output sc)
N^-1AD2c)))

(defvar *max-slack* 1d8)

@@ -148,25 +151,27 @@
(let* ((x (affine-x state))
(l (affine-l state))
(u (affine-u state))
(slack (slack l x u *max-slack*))
(dg (project slack
(affine-c state)
(affine-A-copy state)
(affine-sparse-state state)))
(g (matlisp:m.* dg slack)))
(let ((step (* *gamma* (max-step l x u g)))
(norm-g (matlisp:norm g)))
(when (> step 1d10)
(error "Unbounded problem"))
(format t "~12,5g " norm-g)
(when (or (< (min (* step norm-g) norm-g)
(min (* 1d-10 (matlisp:nrows g))
1d-6))
(plusp (matlisp:dot g (affine-c state))))
(return-from one-affine-scaling-iteration
(values state nil)))
(setf (affine-x state)
(matlisp:axpy! step g x)))
(slack (slack l x u *max-slack*)))
(multiple-value-bind (dg ye)
(project slack (affine-c state)
(affine-A-copy state)
(affine-sparse-state state))
(declare (ignore ye))
(let* ((g (matlisp:m.* dg slack))
(step (* *gamma* (max-step l x u g)))
(norm-g (matlisp:norm g))
(norm-dg (matlisp:norm dg)))
(when (> step 1d10)
(error "Unbounded problem"))
(format t "~12,5g ~12,5g" norm-g norm-dg)
(when (or (< (min (* step norm-g) norm-dg)
(max (* 1d-10 (matlisp:nrows g))
1d-6))
(plusp (matlisp:dot g (affine-c state))))
(return-from one-affine-scaling-iteration
(values state nil)))
(setf (affine-x state)
(matlisp:axpy! step g x))))
(values state t)))

(defun residual (state)
@@ -521,8 +521,9 @@
(cholmod-free-factor (addr factor) common)
(cholmod-free-dense (addr b) common)))))

(defun solve-sparse-recycle (As b state &aux (common *cholmod-common*)
(matlisp-b b))
(defun solve-sparse-recycle (As b state factorized
&aux (common *cholmod-common*)
(matlisp-b b))
(with-alien ((As (* cholmod-sparse) :local As)
(b (* cholmod-dense) :local (make-dense-from-matlisp
b
@@ -537,10 +538,11 @@
(solve-sparse-state-workspace-y state))
(E (* cholmod-dense) :local
(solve-sparse-state-workspace-e state)))
(cholmod-set-status common 0)
(cholmod-factorize As factor common)
(when (/= (cholmod-get-status common) 0)
(return-from solve-sparse-recycle))
(unless factorized
(cholmod-set-status common 0)
(cholmod-factorize As factor common)
(when (/= (cholmod-get-status common) 0)
(return-from solve-sparse-recycle)))
(unless (plusp (cholmod-solve2
0
factor
@@ -557,9 +559,9 @@
(solve-sparse-state-workspace-e state) E)
(dense-to-matlisp x matlisp-b)))

(defun solve-sparse (As b &optional state)
(defun solve-sparse (As b &optional state factorized)
(if state
(solve-sparse-recycle As b state)
(solve-sparse-recycle As b state factorized)
(solve-sparse-one-shot As b)))

(defun sparse-m* (sparse x &key transpose

0 comments on commit 3a489be

Please sign in to comment.