提出 #24569066
ソースコード 拡げる
(in-package :cl-user)
(eval-when (:compile-toplevel :load-toplevel :execute)
(defparameter *opt*
#+swank '(optimize (speed 3) (safety 2))
#-swank '(optimize (speed 3) (safety 0) (debug 0)))
#+swank (ql:quickload '(:cl-debug-print :fiveam :cp/util) :silent t)
#+swank (use-package :cp/util :cl-user)
#-swank (set-dispatch-macro-character
#\# #\> (lambda (s c p) (declare (ignore c p)) `(values ,(read s nil nil t))))
#+sbcl (dolist (f '(:popcnt :sse4)) (pushnew f sb-c:*backend-subfeatures*))
(setq *random-state* (make-random-state t)))
#-swank (eval-when (:compile-toplevel)
(setq *break-on-signals* '(and warning (not style-warning))))
#+swank (set-dispatch-macro-character #\# #\> #'cl-debug-print:debug-print-reader)
(macrolet ((def (b)
`(progn (deftype ,(intern (format nil "UINT~A" b)) () '(unsigned-byte ,b))
(deftype ,(intern (format nil "INT~A" b)) () '(signed-byte ,b))))
(define-int-types (&rest bits) `(progn ,@(mapcar (lambda (b) `(def ,b)) bits))))
(define-int-types 2 4 7 8 15 16 31 32 62 63 64))
(defconstant +mod+ 1000000007)
(defmacro dbg (&rest forms)
(declare (ignorable forms))
#+swank (if (= (length forms) 1)
`(format *error-output* "~A => ~A~%" ',(car forms) ,(car forms))
`(format *error-output* "~A => ~A~%" ',forms `(,,@forms))))
(declaim (inline println))
(defun println (obj &optional (stream *standard-output*))
(let ((*read-default-float-format*
(if (typep obj 'double-float) 'double-float *read-default-float-format*)))
(prog1 (princ obj stream) (terpri stream))))
;; BEGIN_INSERTED_CONTENTS
(defpackage :cp/csc
(:use :cl)
(:export
#:coo #:coo-m #:coo-n #:make-coo #:coo-insert!
#:coo-rows #:coo-cols #:coo-values #:coo-nz
#:csc #:make-csc #:csc-to-array #:make-csc-from-array #:make-csc-from-coo
#:csc-gemv #:csc-gemv-with-basis #:csc-transpose #:csc-m #:csc-n #:csc-nz
#:csc-colstarts #:csc-rows #:csc-values
#:sparse-vector #:make-sparse-vector #:make-sparse-vector-from
#:make-array-from-sparse-vector
#:sparse-vector-nz #:sparse-vector-values #:sparse-vector-indices)
(:documentation "Provides some representations of a sparse matrix."))
(in-package :cp/csc)
(deftype csc-float () 'rational)
(defconstant +zero+ (coerce 0 'csc-float))
(defconstant +one+ (coerce 1 'csc-float))
(defstruct (coo (:constructor %make-coo (m n nz rows cols values))
(:copier nil)
(:predicate nil))
"Stores a sparse matrix with coordinate list representation (aka COO
format). Note that you can increase M and N after construction, but cannot
decrease them."
(m nil :type (mod #.array-dimension-limit))
(n nil :type (mod #.array-dimension-limit))
(nz nil :type (mod #.array-dimension-limit))
(rows nil :type (simple-array fixnum (*)))
(cols nil :type (simple-array fixnum (*)))
(values nil :type (simple-array csc-float (*))))
(defun make-coo (m n &optional (initial-size 4))
(declare ((integer 1 (#.array-dimension-limit)) initial-size))
(let ((initial-size (max 1 initial-size)))
(%make-coo m n 0
(make-array initial-size :element-type 'fixnum)
(make-array initial-size :element-type 'fixnum)
(make-array initial-size :element-type 'csc-float))))
(defun coo-insert! (coo row col value)
"Destructively executes an assignment operation: COO[ROW][COL] := VALUE."
(declare (optimize (speed 3))
((mod #.array-dimension-limit) row col))
(assert (and (< row (coo-m coo))
(< col (coo-n coo))))
(symbol-macrolet ((rows (coo-rows coo))
(cols (coo-cols coo))
(values (coo-values coo))
(nz (coo-nz coo)))
(when (= nz (length rows))
(let ((new-size (max 1 (the (mod #.array-dimension-limit) (* 2 nz)))))
(setq rows (adjust-array rows new-size)
cols (adjust-array cols new-size)
values (adjust-array values new-size))))
(setf (aref rows nz) row
(aref cols nz) col
(aref values nz) value
nz (+ nz 1))
coo))
(defstruct (csc (:constructor make-csc (m n nz colstarts rows values))
(:copier nil)
(:predicate nil))
"Stores a sparse matrix with compressed sparse column representation (aka CSC
format). Note that you can increase M after construction, but cannot decrease it
or change N."
(m nil :type (mod #.array-dimension-limit))
(n nil :type (mod #.array-dimension-limit))
(nz nil :type (mod #.array-dimension-limit))
(colstarts nil :type (simple-array fixnum (*)))
(rows nil :type (simple-array fixnum (*)))
(values nil :type (simple-array csc-float (*))))
(defun csc-to-array (csc &optional rowperm colperm)
"Makes a 2-dimensional array from a CSC."
(declare (optimize (speed 3))
((or null vector) rowperm colperm))
(let* ((m (csc-m csc))
(n (csc-n csc))
(res (make-array (list m n) :element-type 'csc-float :initial-element +zero+))
(colstarts (csc-colstarts csc))
(rows (csc-rows csc))
(values (csc-values csc)))
(dotimes (col n)
(let ((dest-col (if colperm (aref colperm col) col)))
(loop for i from (aref colstarts col) below (aref colstarts (+ col 1))
for row = (aref rows i)
for dest-row = (if rowperm (aref rowperm row) row)
for value = (aref values i)
do (setf (aref res dest-row dest-col) value))))
res))
(defun make-csc-from-array (array)
"Makes a CSC from a 2-dimensional array."
(declare (optimize (speed 3))
((array * (* *)) array))
(destructuring-bind (m n) (array-dimensions array)
(declare ((mod #.array-dimension-limit) m n))
(let* ((colstarts (make-array (+ n 1) :element-type 'fixnum))
(nz (count +zero+ (sb-ext:array-storage-vector array) :test-not #'=))
(rows (make-array nz :element-type 'fixnum))
(values (make-array nz :element-type 'csc-float))
(end 0))
(declare ((mod #.array-dimension-limit) end))
(dotimes (col n)
(setf (aref colstarts col) end)
(dotimes (row m)
(let ((value (aref array row col)))
(unless (zerop value)
(setf (aref rows end) row
(aref values end) value)
(incf end)))))
(setf (aref colstarts n) end)
(make-csc m n nz colstarts rows values))))
(declaim (inline make-csc-from-coo))
(defun make-csc-from-coo (coo)
"Makes a CSC from a COO. Note that the returned CSC contains zero when COO
contains it."
(declare (optimize (speed 3))
(inline stable-sort))
(let* ((m (coo-m coo))
(n (coo-n coo))
(rows (coo-rows coo))
(cols (coo-cols coo))
(values (coo-values coo))
(indices (let ((tmp (make-array (coo-nz coo) :element-type 'fixnum)))
(dotimes (i (length tmp))
(setf (aref tmp i) i))
(stable-sort tmp (lambda (i1 i2)
(if (= (aref cols i1) (aref cols i2))
(< (aref rows i1) (aref rows i2))
(< (aref cols i1) (aref cols i2)))))))
(nz 0))
(declare ((simple-array fixnum (*)) indices)
((mod #.array-dimension-limit) nz))
;; drop duplicate elements
(dotimes (i* (length indices))
(let ((i (aref indices i*)))
(if (and (> i* 0)
(let ((prev-i (aref indices (- i* 1))))
(and (= (aref rows i) (aref rows prev-i))
(= (aref cols i) (aref cols prev-i)))))
(setf (aref indices (- nz 1)) i)
(setf (aref indices nz) i
nz (+ nz 1)))))
(let ((colstarts (make-array (+ n 1) :element-type 'fixnum))
(csc-rows (make-array nz :element-type 'fixnum))
(csc-values (make-array nz :element-type 'csc-float))
(end 0)
(prev-col -1))
(declare ((mod #.array-dimension-limit) end)
((integer -1 (#.array-dimension-limit)) prev-col))
(dotimes (i* nz)
(let* ((i (aref indices i*))
(row (aref rows i))
(col (aref cols i))
(value (aref values i)))
(loop for j from col above prev-col
do (setf (aref colstarts j) end))
(setf (aref csc-rows end) row
(aref csc-values end) value)
(incf end)
(setq prev-col col)))
(loop for j from n above prev-col
do (setf (aref colstarts j) end))
(make-csc m n nz colstarts csc-rows csc-values))))
(declaim (inline csc-gemv-with-basis))
(defun csc-gemv-with-basis (csc vector basis)
"Multiplies CSC and VECTOR and returns a resultant vector."
(declare (vector vector basis))
(let* ((m (csc-m csc))
(colstarts (csc-colstarts csc))
(rows (csc-rows csc))
(values (csc-values csc))
(res (make-array m :element-type 'csc-float :initial-element +zero+)))
(dotimes (i m)
(let ((bi (aref basis i)))
(loop for k from (aref colstarts bi) below (aref colstarts (+ bi 1))
do (incf (aref res (aref rows k))
(* (aref values k) (coerce (aref vector i) 'csc-float))))))
res))
(declaim (inline csc-gemv))
(defun csc-gemv (csc vector)
"Multiplies CSC and VECTOR and returns a resultant vector."
(declare (vector vector))
(let* ((m (csc-m csc))
(n (length vector))
(colstarts (csc-colstarts csc))
(rows (csc-rows csc))
(values (csc-values csc))
(res (make-array m :element-type 'csc-float :initial-element +zero+)))
(dotimes (j n)
(loop for k from (aref colstarts j) below (aref colstarts (+ j 1))
do (incf (aref res (aref rows k))
(* (aref values k) (coerce (aref vector j) 'csc-float)))))
res))
(defun csc-transpose (csc)
"Returns the transposed matrix of CSC. This function is non-destructive."
(declare (optimize (speed 3)))
(let* ((m (csc-m csc))
(n (csc-n csc))
(nz (csc-nz csc))
(colstarts (csc-colstarts csc))
(rows (csc-rows csc))
(values (csc-values csc))
(tmp (make-array m :element-type 'fixnum :initial-element 0))
(new-colstarts (make-array (+ m 1) :element-type 'fixnum :initial-element 0))
(new-rows (make-array nz :element-type 'fixnum :initial-element 0))
(new-values (make-array nz :element-type 'csc-float :initial-element +zero+)))
(dotimes (k nz)
(let ((row (aref rows k)))
(incf (aref tmp row))))
;; TMP[ROW] == the number of the elements in ROW
(dotimes (i m)
(setf (aref new-colstarts (+ i 1))
(+ (aref new-colstarts i) (aref tmp i))))
(fill tmp 0)
(dotimes (j n)
(loop for k from (aref colstarts j) below (aref colstarts (+ j 1))
for row = (aref rows k)
for new-pos = (+ (aref new-colstarts row) (aref tmp row))
do (incf (aref tmp row))
(setf (aref new-rows new-pos) j
(aref new-values new-pos) (aref values k))))
(make-csc n m nz new-colstarts new-rows new-values)))
;; TODO: Should we store the intended dimension of vector?
(defstruct (sparse-vector (:constructor %make-sparse-vector (nz values indices)))
(nz nil :type (mod #.array-dimension-limit))
(values nil :type (simple-array csc-float (*)))
(indices nil :type (simple-array fixnum (*))))
(defun make-sparse-vector (size)
(%make-sparse-vector
0
(make-array size :element-type 'csc-float)
(make-array size :element-type 'fixnum)))
(defun make-sparse-vector-from (vector)
(declare (vector vector))
(let* ((nz (count +zero+ vector :test-not #'=))
(values (make-array nz :element-type 'csc-float))
(indices (make-array nz :element-type 'fixnum))
(end 0))
(declare ((mod #.array-dimension-limit) end))
(dotimes (i (length vector))
(unless (zerop (aref vector i))
(setf (aref values end) (aref vector i)
(aref indices end) i)
(incf end)))
(%make-sparse-vector nz values indices)))
(defun make-array-from-sparse-vector (sparse-vector)
(let* ((nz (sparse-vector-nz sparse-vector))
(indices (sparse-vector-indices sparse-vector))
(values (sparse-vector-values sparse-vector))
(m (if (zerop (length indices))
0
(+ 1 (reduce #'max indices))))
(res (make-array m :element-type 'csc-float :initial-element +zero+)))
(dotimes (i nz)
(let ((index (aref indices i))
(value (aref values i)))
(setf (aref res index) value)))
res))
(defpackage :cp/extend-vector
(:use :cl)
(:export #:extend-vector #:extend-vectorf))
(in-package :cp/extend-vector)
(declaim (inline %power-of-two-ceiling))
(defun %power-of-two-ceiling (x)
(ash 1 (integer-length (- x 1))))
(defun extend-vector (vector size &optional (initial-element nil supplied-p))
(declare (optimize (speed 3))
((mod #.array-dimension-limit) size)
(vector vector))
(let ((new-size (%power-of-two-ceiling (max size 1))))
(declare ((mod #.array-dimension-limit) new-size))
(if (< (length vector) new-size)
(if supplied-p
(adjust-array vector new-size :initial-element initial-element)
(adjust-array vector new-size))
vector)))
(defconstant +not-supplied+ (if (boundp '+not-supplied+)
(symbol-value '+not-supplied+)
(make-symbol "NOT-SUPPLIED")))
(define-modify-macro extend-vectorf (new-size &optional (initial-element '+not-supplied+))
(lambda (vector new-size initial-element)
(declare ((mod #.array-dimension-limit) new-size))
(if (< (length vector) new-size)
(if (eq initial-element +not-supplied+)
(extend-vector vector new-size)
(extend-vector vector new-size initial-element))
vector)))
(defpackage :cp/rdtscp
(:use :cl :sb-assem)
(:import-from :sb-c #:defknown #:define-vop #:move)
(:import-from :sb-vm #:unsigned-reg #:rax-offset #:rdx-offset #:rcx-offset #:unsigned-num)
(:export #:read-tsc)
(:documentation "Provides a reader for time-stamp counter using RDTSCP
instruction."))
(in-package :cp/rdtscp)
(eval-when (:compile-toplevel :load-toplevel :execute)
(when (find-package :sb-x86-64-asm)
(sb-ext:unlock-package :sb-x86-64-asm)))
(defmacro emit-bytes (segment &rest bytes)
`(progn ,@(mapcar (lambda (x) `(emit-byte ,segment ,x)) bytes)))
(eval-when (:compile-toplevel :load-toplevel :execute)
(define-instruction rdtscp (segment)
(:emitter (emit-bytes segment #x0f #x01 #xf9))))
(eval-when (:compile-toplevel :load-toplevel :execute)
(defknown %read-tsc () (values (unsigned-byte 32) (unsigned-byte 32)) ()
:overwrite-fndb-silently t)
(define-vop (%read-tsc)
(:policy :fast-safe)
(:translate %read-tsc)
(:temporary (:sc unsigned-reg :offset rax-offset :target lo) eax)
(:temporary (:sc unsigned-reg :offset rdx-offset :target hi) edx)
;; RDTSCP instruction reads IA32_TSC_AUX value into ECX
(:temporary (:sc unsigned-reg :offset rcx-offset) ecx)
(:ignore ecx)
(:results (hi :scs (unsigned-reg))
(lo :scs (unsigned-reg)))
(:result-types unsigned-num unsigned-num)
(:generator 3
(inst rdtscp)
(move lo eax)
(move hi edx))))
(declaim (inline read-tsc))
(defun read-tsc ()
(multiple-value-bind (hi lo) (%read-tsc)
(dpb hi (byte 32 32) lo)))
(eval-when (:compile-toplevel :load-toplevel :execute)
(when (find-package :sb-x86-64-asm)
(sb-ext:lock-package :sb-x86-64-asm)))
(defpackage :cp/iterset
(:use :cl)
(:export #:node-key #:iterset-empty-p #:iterset-first #:iterset-last #:iterset-insert
#:make-iterset #:iterset-map #:node-next #:node-prev))
(in-package :cp/iterset)
(deftype key-type () 'fixnum)
(declaim (inline order))
(defun order (x y)
(declare (key-type x y))
(< x y))
(defstruct (node (:constructor %make-node (key priority &optional left right parent))
(:copier nil)
(:predicate nil)
(:conc-name %node-))
(key nil :type key-type)
(priority nil :type (integer 0 #.most-positive-fixnum))
(left nil :type (or null node))
(right nil :type (or null node))
(parent nil :type (or null node)))
(declaim (inline node-key))
(defun node-key (node)
(%node-key node))
(defstruct (iterset (:constructor make-iterset ())
(:copier nil)
(:predicate nil)
(:conc-name %iterset-))
(root nil :type (or null node)))
(declaim (inline iterset-empty-p))
(defun iterset-empty-p (iterset)
(null (%iterset-root iterset)))
(defun node-split (node key)
"Destructively splits NODE with reference to KEY and returns two trees, the
smaller tree (< KEY) and the larger one (>= KEY)."
(declare (optimize (speed 3)))
(labels ((recur (node)
(cond ((null node)
(values nil nil))
((order (%node-key node) key)
(multiple-value-bind (left right)
(recur (%node-right node))
(when left
(setf (%node-parent left) node))
(setf (%node-right node) left)
(values node right)))
(t
(multiple-value-bind (left right)
(recur (%node-left node))
(when right
(setf (%node-parent right) node))
(setf (%node-left node) right)
(values left node))))))
(recur node)))
(defun node-insert (node key)
(declare (optimize (speed 3))
((or null node) node))
"Destructively inserts KEY into NODE and returns the resultant tree. You
**cannot** rely on the side effect. Use the returned value. The behavior is
undefined when duplicate keys are inserted."
(let ((new-node (%make-node key (random #.(+ 1 most-positive-fixnum)))))
(labels ((recur (node)
(cond ((null node) new-node)
((> (%node-priority new-node) (%node-priority node))
(multiple-value-bind (left right)
(node-split node key)
(when left
(setf (%node-parent left) new-node))
(when right
(setf (%node-parent right) new-node))
(setf (%node-left new-node) left
(%node-right new-node) right)
new-node))
((order key (%node-key node))
(let ((left (recur (%node-left node))))
(when left
(setf (%node-parent left) node))
(setf (%node-left node) left))
node)
(t
(let ((right (recur (%node-right node))))
(when right
(setf (%node-parent right) node))
(setf (%node-right node) right))
node))))
(recur node))))
(defun iterset-insert (iterset key)
"Destructively inserts KEY into ITERSET. You **can** rely on the side
effect. The behavior is undefined when duplicate keys are inserted."
(setf (%iterset-root iterset)
(node-insert (%iterset-root iterset) key))
iterset)
(declaim (inline node-first))
(defun node-first (node)
(declare (node node))
(loop while (%node-left node)
do (setq node (%node-left node)))
node)
(declaim (inline iterset-first))
(defun iterset-first (iterset)
"Returns the first **node** of ITERSET."
(when (%iterset-root iterset)
(node-first (%iterset-root iterset))))
(declaim (inline node-last))
(defun node-last (node)
(declare (node node))
(loop while (%node-right node)
do (setq node (%node-right node)))
node)
(declaim (inline iterset-last))
(defun iterset-last (iterset)
"Returns the last **node** of ITERSET."
(when (%iterset-root iterset)
(node-last (%iterset-root iterset))))
(declaim (inline node-map))
(defun node-map (function node)
"Successively applies FUNCTION to each key of NODE. FUNCTION must take one
argument."
(declare (function function))
(labels ((recur (node)
(when node
(recur (%node-left node))
(funcall function (%node-key node))
(recur (%node-right node)))))
(recur node)))
(declaim (inline iterset-map))
(defun iterset-map (function iterset)
"Successively applies FUNCTION to each key of ITERSET. FUNCTION must take one
argument."
(node-map function (%iterset-root iterset)))
(defmethod print-object ((object node) stream)
(print-unreadable-object (object stream :type t)
(let ((init t))
(node-map (lambda (key)
(if init
(setf init nil)
(write-char #\ stream))
(write key :stream stream))
object))))
(defun node-next (node)
"Returns the next node of NODE. Returns NIL instead if NODE is the last node."
(declare (optimize (speed 3))
(node node))
(let ((key (%node-key node)))
(if (%node-right node)
(node-first (%node-right node))
(loop while (%node-parent node)
do (setq node (%node-parent node))
when (order key (%node-key node))
do (return node)
finally (return nil)))))
(defun node-prev (node)
"Returns the previous node of NODE. Returns NIL instead if NODE is the first
node."
(declare (optimize (speed 3))
(node node))
(let ((key (%node-key node)))
(if (%node-left node)
(node-last (%node-left node))
(loop while (%node-parent node)
do (setq node (%node-parent node))
when (order (%node-key node) key)
do (return node)
finally (return nil)))))
(defpackage :cp/movable-binary-heap
(:use :cl)
(:export #:make-heap #:heap-empty-error #:heap-pop
#:heap-ensure-key #:heap-clear #:heap-peek #:heap-empty-p #:heap-count)
(:documentation "Provides binary heap implementation with decrease-key and
increase-key operations."))
(in-package :cp/movable-binary-heap)
;; TODO: abstraction
;; NOTE: not tested
(declaim (inline compare))
(defun compare (x y)
(declare (fixnum x y))
(< x y))
(define-condition heap-empty-error (error)
((heap :initarg :heap :reader heap-empty-error-heap))
(:report
(lambda (condition stream)
(format stream "Attempted to get an element from empty heap ~W"
(heap-empty-error-heap condition)))))
(defstruct (heap
(:constructor make-heap
(size
&aux
(keys (make-array (+ 1 size) :element-type 'fixnum))
(priorities (make-array (+ 1 size) :element-type 'fixnum))
(positions (make-array (+ 1 size)
:element-type 'fixnum
:initial-element -1)))))
(keys nil :type (simple-array fixnum (*)))
(priorities nil :type (simple-array fixnum (*)))
(positions nil :type (simple-array fixnum (*)))
(end 1 :type (integer 1 (#.array-dimension-limit))))
(defun %heapify-up (heap pos)
(declare (optimize (speed 3))
((mod #.array-dimension-limit) pos))
(let* ((keys (heap-keys heap))
(priorities (heap-priorities heap))
(positions (heap-positions heap))
(key (aref keys pos)))
(labels ((recur (pos)
(unless (= pos 1)
(let ((parent-pos (ash pos -1)))
(when (compare (aref priorities pos) (aref priorities parent-pos))
(rotatef (aref positions key) (aref positions (aref keys parent-pos)))
(rotatef (aref keys pos) (aref keys parent-pos))
(rotatef (aref priorities pos) (aref priorities parent-pos))
(recur parent-pos))))))
(recur pos))))
(defun %heapify-down (heap pos)
(declare (optimize (speed 3))
((mod #.array-dimension-limit) pos))
(let* ((keys (heap-keys heap))
(priorities (heap-priorities heap))
(positions (heap-positions heap))
(end (heap-end heap))
(key (aref keys pos)))
(labels
((recur (pos)
(let* ((cpos1 (+ pos pos))
(cpos2 (+ 1 cpos1)))
(declare ((mod #.array-dimension-limit) cpos1 cpos2))
(when (< cpos1 end)
(if (< cpos2 end)
(if (compare (aref priorities cpos1) (aref priorities cpos2))
(when (compare (aref priorities cpos1) (aref priorities pos))
(rotatef (aref positions key) (aref positions (aref keys cpos1)))
(rotatef (aref keys pos) (aref keys cpos1))
(rotatef (aref priorities pos) (aref priorities cpos1))
(recur cpos1))
(when (compare (aref priorities cpos2) (aref priorities pos))
(rotatef (aref positions key) (aref positions (aref keys cpos2)))
(rotatef (aref keys pos) (aref keys cpos2))
(rotatef (aref priorities pos) (aref priorities cpos2))
(recur cpos2)))
(when (compare (aref priorities cpos1) (aref priorities pos))
(rotatef (aref positions key) (aref positions (aref keys cpos1)))
(rotatef (aref keys pos) (aref keys cpos1))
(rotatef (aref priorities pos) (aref priorities cpos1))))))))
(recur pos))))
(declaim #+sbcl (sb-ext:maybe-inline heap-ensure-key))
(defun heap-ensure-key (heap key priority &optional if-exists)
"Adds KEY to HEAP if it doesn't exist and sets its priority to
PRIORITY. Otherwise updates its priority by IF-EXISTS."
(declare (optimize (speed 3))
((mod #.array-dimension-limit) key))
(symbol-macrolet ((end (heap-end heap)))
(when (>= end (length (heap-keys heap)))
(let ((new-size (min (* end 2) #.(- array-dimension-limit 1))))
(setf (heap-keys heap) (adjust-array (heap-keys heap) new-size)
(heap-priorities heap) (adjust-array (heap-priorities heap) new-size))))
(let ((keys (heap-keys heap))
(priorities (heap-priorities heap))
(positions (heap-positions heap)))
(let ((pos (aref positions key)))
(if (= pos -1)
(progn
(setf (aref keys end) key
(aref priorities end) priority
(aref positions key) end)
(%heapify-up heap end)
(incf end))
(let ((prev-priority (aref priorities pos)))
(setf (aref priorities pos)
(if if-exists
(funcall if-exists prev-priority)
priority))
(if (compare prev-priority (aref priorities pos))
(%heapify-down heap pos)
(%heapify-up heap pos)))))
heap)))
(declaim #+sbcl (sb-ext:maybe-inline heap-pop))
(defun heap-pop (heap)
"Deletes the topmost element of HEAP and returns its key and priority."
(declare (optimize (speed 3)))
(symbol-macrolet ((end (heap-end heap)))
(let ((keys (heap-keys heap))
(priorities (heap-priorities heap))
(positions (heap-positions heap)))
(when (= end 1)
(error 'heap-empty-error :heap heap))
(multiple-value-prog1 (values (aref keys 1) (aref priorities 1))
(decf end)
(setf (aref positions (aref keys 1)) -1)
(unless (= end 1)
(setf (aref keys 1) (aref keys end)
(aref priorities 1) (aref priorities end)
(aref positions (aref keys 1)) 1))
(%heapify-down heap 1)))))
(declaim (inline heap-clear))
(defun heap-clear (heap)
"Makes heap empty. Note that it takes linear time."
(setf (heap-end heap) 1)
(fill (heap-positions heap) -1)
heap)
(declaim (inline heap-empty-p))
(defun heap-empty-p (heap)
"Returns true iff HEAP is empty."
(= 1 (heap-end heap)))
(declaim (inline heap-count))
(defun heap-count (heap)
"Returns the current number of the elements in HEAP."
(- (heap-end heap) 1))
(declaim (inline heap-peek))
(defun heap-peek (heap)
"Returns the topmost element of HEAP, i.e. its key and priority."
(if (= 1 (heap-end heap))
(error 'heap-empty-error :heap heap)
(values (aref (heap-keys heap) 1)
(aref (heap-priorities heap) 1))))
(defpackage :cp/lud
(:use :cl :cp/csc #:cp/movable-binary-heap #:cp/iterset #:cp/rdtscp #:cp/extend-vector)
(:import-from :cp/csc #:csc-float #:+zero+ #:+one+)
(:export #:lud #:lud-eta #:lu-factor
#:lud-lower #:lud-upper #:lud-tlower #:lud-tupper #:lud-diagu
#:lud-rank #:lud-colperm #:lud-icolperm #:lud-rowperm #:lud-irowperm #:lud-m
#:make-lud-eta #:dense-solve! #:add-eta! #:lud-eta-lud #:lud-eta-count
#:refactor #:refactor-p #:*refactor-threshold* #:*refactor-by-time*
#:gauss-eta! #:gauss-eta-transposed!
#:sparse-solve! #:sparse-solve-transposed!)
(:documentation "Provides LU decomposition for sparse matrix. This module is
served for simplex method.
Reference:
Robert J. Vanderbei. Linear Programming: Foundations and Extensions. 5th edition."))
(in-package :cp/lud)
(defconstant +eps+ (coerce 0 'csc-float))
(defconstant +epsnum+ (coerce 0 'csc-float))
(defstruct (lud (:conc-name lud-))
(m nil :type (mod #.array-dimension-limit))
;; lower triangular matrix
(lower nil :type csc)
;; transposition of LOWER
(tlower nil :type csc)
;; upper triangular matrix
(upper nil :type csc)
;; transposition of UPPER
(tupper nil :type csc)
;; diagonal elements of upper triagular matrix
(diagu nil :type (simple-array csc-float (*)))
(rank nil :type (mod #.array-dimension-limit))
(colperm nil :type (simple-array fixnum (*)))
(icolperm nil :type (simple-array fixnum (*)))
(rowperm nil :type (simple-array fixnum (*)))
(irowperm nil :type (simple-array fixnum (*))))
(defmacro vector-set* (vector index new-element)
(let ((i (gensym))
(elm (gensym)))
`(let ((,i ,index)
(,elm ,new-element))
(extend-vectorf ,vector (+ 1 ,i))
(setf (aref ,vector ,i) ,elm))))
(deftype ivec () '(simple-array fixnum (*)))
(deftype fvec () '(simple-array csc-float (*)))
(defconstant +nan+ -1)
(declaim (inline %power-of-two-ceiling))
(defun %power-of-two-ceiling (x)
(ash 1 (integer-length (- x 1))))
(defun lu-factor (matrix basis)
(declare (optimize (speed 3))
(vector basis))
;; B: submatrix of MATRIX w.r.t. BASIS
;; BT: transposition of B
(let* ((m (csc-m matrix))
(colstarts (csc-colstarts matrix))
(rows (csc-rows matrix))
(values (csc-values matrix))
(colperm (make-array m :element-type 'fixnum :initial-element 0))
(icolperm (make-array m :element-type 'fixnum :initial-element +nan+))
(rowperm (make-array m :element-type 'fixnum :initial-element 0))
(irowperm (make-array m :element-type 'fixnum :initial-element +nan+))
;; the number of non-zero elements in each column
(b-degs (make-array m :element-type 'fixnum :initial-element 0))
(bt-degs (make-array m :element-type 'fixnum :initial-element 0))
;; HEAP keeps the column numbers in ascending order of degree.
;; Priority of col is M+1 iff B-DEGS[col] is zero;
;; Priority of col is M+2 iff col has no elements which are larger than +EPSNUM+.
(heap (make-heap m))
(tmp (make-array m :element-type 'fixnum :initial-element 0))
(tmp2 (make-array m :element-type 'fixnum :initial-element +nan+)))
(assert (= m (length basis)))
(dotimes (i m)
(let ((start (aref colstarts (aref basis i)))
(end (aref colstarts (+ 1 (the fixnum (aref basis i))))))
(declare ((mod #.array-dimension-limit) start end))
(setf (aref b-degs i) (- end start))
(loop for k from start below end
do (incf (aref bt-degs (aref rows k))))))
(let* ((estimated-nz (%power-of-two-ceiling
(ash (the fixnum (reduce #'+ b-degs)) -1)))
(lower-colstarts (make-array (+ m 1) :element-type 'fixnum :initial-element 0))
(lower-rows (make-array estimated-nz :element-type 'fixnum))
(lower-values (make-array estimated-nz :element-type 'csc-float))
(tupper-colstarts (make-array (+ m 1) :element-type 'fixnum :initial-element 0))
(tupper-rows (make-array estimated-nz :element-type 'fixnum))
(tupper-values (make-array estimated-nz :element-type 'csc-float))
(diagu (make-array m :element-type 'csc-float))
(b-indices (make-array m :element-type t))
(b-values (make-array m :element-type t))
(bt-indices (make-array m :element-type t))
(bt-values (make-array m :element-type t))
(rank m)
(tag 1)
(lower-nz 0)
(upper-nz 0))
(declare (ivec lower-rows tupper-rows)
(fvec lower-values tupper-values)
((mod #.array-dimension-limit) rank tag lower-nz upper-nz))
;; initialize B and BT
(dotimes (i m)
(setf (aref b-indices i) (make-array (aref b-degs i) :element-type 'fixnum)
(aref b-values i) (make-array (aref b-degs i) :element-type 'csc-float)
(aref bt-indices i) (make-array (aref bt-degs i) :element-type 'fixnum)
(aref bt-values i) (make-array (aref bt-degs i) :element-type 'csc-float)))
(let ((bt-poses tmp))
(dotimes (i m)
(let ((bpos 0)
(start (aref colstarts (aref basis i)))
(end (aref colstarts (+ 1 (the fixnum (aref basis i))))))
(declare ((mod #.array-dimension-limit) start end))
(loop for k from start below end
for row = (aref rows k)
for bt-pos = (aref bt-poses row)
do (setf (aref (the ivec (aref b-indices i)) bpos) row
(aref (the fvec (aref b-values i)) bpos) (aref values k)
(aref (the ivec (aref bt-indices row)) bt-pos) i
(aref (the fvec (aref bt-values row)) bt-pos) (aref values k))
(incf (aref bt-poses row))
(incf bpos))))
(fill bt-poses 0))
;; initialize heap
(dotimes (j m)
(heap-ensure-key heap j
(if (zerop (aref b-degs j)) (+ m 1) (aref b-degs j))))
(block decomposition
(dotimes (i m)
(multiple-value-bind (pivot-col coldeg pivot-row rowdeg)
(loop
(multiple-value-bind (col priority) (heap-pop heap)
(when (> priority m)
;; singular matrix
(when (= priority (+ m 1))
(assert (zerop (aref b-degs col))))
(setq rank i)
(return-from decomposition))
(let ((coldeg (aref b-degs col))
(rowdeg (+ m 1))
row)
;; select row of minimal degree
(dotimes (k coldeg)
(when (and (< (aref bt-degs (aref (the ivec (aref b-indices col)) k))
rowdeg)
(> (abs (aref (the fvec (aref b-values col)) k))
+epsnum+))
(setq row (aref (the ivec (aref b-indices col)) k)
rowdeg (aref bt-degs row))))
(when (< rowdeg (+ m 1))
(return (values col coldeg row rowdeg)))
;; col has no elements which are sufficiently distant from zero
(heap-ensure-key heap col (+ m 2)))))
(declare ((mod #.array-dimension-limit) pivot-col coldeg pivot-row rowdeg))
(setf (aref colperm i) pivot-col
(aref icolperm pivot-col) i
(aref rowperm i) pivot-row
(aref irowperm pivot-row) i)
;; copy pivot column into LOWER and pivot row into TUPPER
(setf (aref lower-colstarts (+ i 1))
(+ (aref lower-colstarts i) coldeg -1))
(dotimes (k coldeg)
(let ((row (aref (the ivec (aref b-indices pivot-col)) k))
(val (aref (the fvec (aref b-values pivot-col)) k)))
(unless (= row pivot-row)
(vector-set* lower-rows lower-nz row)
(vector-set* lower-values lower-nz val)
(incf lower-nz))))
(setf (aref tupper-colstarts (+ i 1))
(+ (aref tupper-colstarts i) rowdeg -1))
(dotimes (k rowdeg)
(let ((col (aref (the ivec (aref bt-indices pivot-row)) k))
(val (aref (the fvec (aref bt-values pivot-row)) k)))
(if (= col pivot-col)
(setf (aref diagu i) val)
(progn
(vector-set* tupper-rows upper-nz col)
(vector-set* tupper-values upper-nz val)
(incf upper-nz)))))
;; remove pivot rows and element from B and BT
(dotimes (k coldeg)
(let ((row (aref (the ivec (aref b-indices pivot-col)) k)))
(decf (aref bt-degs row))
(let ((bt-pos 0))
(loop until (= pivot-col (aref (the ivec (aref bt-indices row)) bt-pos))
do (incf bt-pos))
(rotatef (aref (the ivec (aref bt-indices row)) bt-pos)
(aref (the ivec (aref bt-indices row)) (aref bt-degs row)))
(rotatef (aref (the fvec (aref bt-values row)) bt-pos)
(aref (the fvec (aref bt-values row)) (aref bt-degs row))))))
(dotimes (k rowdeg)
(let ((col (aref (the ivec (aref bt-indices pivot-row)) k)))
(decf (aref b-degs col))
(let ((b-pos 0))
(loop until (= pivot-row (aref (the ivec (aref b-indices col)) b-pos))
do (incf b-pos))
(rotatef (aref (the ivec (aref b-indices col)) b-pos)
(aref (the ivec (aref b-indices col)) (aref b-degs col)))
(rotatef (aref (the fvec (aref b-values col)) b-pos)
(aref (the fvec (aref b-values col)) (aref b-degs col))))))
(setf (aref b-degs pivot-col) 0
(aref bt-degs pivot-row) 0)
;; update B and BT
(loop
with tags = tmp
with bt-poses = tmp2
for k from (aref lower-colstarts i) below (aref lower-colstarts (+ i 1))
for row = (aref lower-rows k)
do (dotimes (bt-pos (aref bt-degs row))
(setf (aref tags (aref (the ivec (aref bt-indices row)) bt-pos)) tag
(aref bt-poses (aref (the ivec (aref bt-indices row)) bt-pos)) bt-pos))
(loop
for kk from (aref tupper-colstarts i) below (aref tupper-colstarts (+ i 1))
for col = (aref tupper-rows kk)
for decr of-type csc-float = (/ (* (aref lower-values k)
(aref tupper-values kk))
(aref diagu i))
when (= (aref tags col) tag)
do (decf (aref (the fvec (aref bt-values row)) (aref bt-poses col))
decr)
else
do (let ((deg (aref bt-degs row)))
(vector-set* (the ivec (aref bt-indices row)) deg col)
(vector-set* (the fvec (aref bt-values row)) deg (- decr))
(incf (aref bt-degs row))))
(incf tag))
(loop
with tags = tmp
with b-poses = tmp2
for k from (aref tupper-colstarts i) below (aref tupper-colstarts (+ i 1))
for col = (aref tupper-rows k)
do (dotimes (b-pos (aref b-degs col))
(setf (aref tags (aref (the ivec (aref b-indices col)) b-pos)) tag
(aref b-poses (aref (the ivec (aref b-indices col)) b-pos)) b-pos))
(loop
for kk from (aref lower-colstarts i) below (aref lower-colstarts (+ i 1))
for row = (aref lower-rows kk)
for decr of-type csc-float = (/ (* (aref lower-values kk)
(aref tupper-values k))
(aref diagu i))
when (= (aref tags row) tag)
do (decf (aref (the fvec (aref b-values col)) (aref b-poses row))
decr)
else
do (let ((deg (aref b-degs col)))
(vector-set* (the ivec (aref b-indices col)) deg row)
(vector-set* (the fvec (aref b-values col)) deg (- decr))
(incf (aref b-degs col))))
(incf tag))
(loop for k from (aref tupper-colstarts i) below (aref tupper-colstarts (+ i 1))
for col = (aref tupper-rows k)
do (heap-ensure-key heap col
(if (zerop (aref b-degs col))
(+ m 1)
(aref b-degs col)))))))
;; fill dependent rows and cols
(let ((end rank))
(declare ((mod #.array-dimension-limit) end))
(dotimes (col m)
(when (= +nan+ (aref icolperm col))
(setf (aref colperm end) col
(aref icolperm col) end
end (+ end 1)))))
(let ((end rank))
(declare ((mod #.array-dimension-limit) end))
(dotimes (row m)
(when (= +nan+ (aref irowperm row))
(setf (aref rowperm end) row
(aref irowperm row) end
end (+ end 1)))))
(loop for i from rank below m
do (setf (aref lower-colstarts (+ i 1)) (aref lower-colstarts i)
(aref tupper-colstarts (+ i 1)) (aref tupper-colstarts i)
(aref diagu i) +zero+))
(dotimes (k (aref lower-colstarts m))
(setf (aref lower-rows k)
(aref irowperm (aref lower-rows k))))
(dotimes (k (aref tupper-colstarts m))
(setf (aref tupper-rows k)
(aref icolperm (aref tupper-rows k))))
(dotimes (i m)
(loop for k from (aref lower-colstarts i) below (aref lower-colstarts (+ i 1))
do (setf (aref lower-values k)
(/ (aref lower-values k) (aref diagu i)))))
(let* ((lower (make-csc m m lower-nz lower-colstarts lower-rows lower-values))
(tupper (make-csc m m upper-nz tupper-colstarts tupper-rows tupper-values))
(tlower (csc-transpose lower))
(upper (csc-transpose tupper)))
(make-lud :m m
:lower lower
:tlower tlower
:upper upper
:tupper tupper
:diagu diagu
:rank rank
:colperm colperm
:icolperm icolperm
:rowperm rowperm
:irowperm irowperm)))))
(defun dense-solve! (lud y)
"Solves LUx = y and stores the solution to Y. The consequence is undefined
when it is infeasible."
(declare (optimize (speed 3))
(fvec y))
(let* ((m (lud-m lud))
(irowperm (lud-irowperm lud))
(colperm (lud-colperm lud))
(lower (lud-lower lud))
(lower-colstarts (csc-colstarts lower))
(lower-rows (csc-rows lower))
(lower-values (csc-values lower))
(tupper (lud-tupper lud))
(tupper-colstarts (csc-colstarts tupper))
(tupper-rows (csc-rows tupper))
(tupper-values (csc-values tupper))
(diagu (lud-diagu lud))
(rank (lud-rank lud))
(tmp (make-array m :element-type 'csc-float)))
(dotimes (i m)
(setf (aref tmp i) (aref y i)))
(dotimes (i m)
(setf (aref y (aref irowperm i)) (aref tmp i)))
;; y := L^(-1)y
(dotimes (i rank)
(let ((val (aref y i)))
(loop for k from (aref lower-colstarts i) below (aref lower-colstarts (+ i 1))
for row = (aref lower-rows k)
do (decf (aref y row) (* (aref lower-values k) val)))))
(loop for i from (- m 1) downto rank
do (setf (aref y i) +zero+))
(loop for i from (- rank 1) downto 0
for val of-type csc-float = (aref y i)
do (loop for k from (aref tupper-colstarts i) below (aref tupper-colstarts (+ i 1))
do (decf val (* (aref tupper-values k) (aref y (aref tupper-rows k)))))
(setf (aref y i) (/ val (aref diagu i))))
(dotimes (i m)
(setf (aref tmp i) (aref y i)))
(dotimes (i m)
(setf (aref y (aref colperm i)) (aref tmp i)))
y))
(defstruct (lud-eta (:constructor %make-lud-eta))
(lud nil :type lud)
(count 0 :type (mod #.array-dimension-limit))
(nz 0 :type (mod #.array-dimension-limit))
;; This vector stores the leaving columns of each pivotting. Note that this
;; `column' doesn't mean the index of some variable but the column number of
;; the matrix B.
(leaving-cols nil :type (simple-array fixnum (*)))
(colstarts nil :type (simple-array fixnum (*)))
(rows nil :type (simple-array fixnum (*)))
(values nil :type (simple-array csc-float (*)))
(cumtime 0 :type (unsigned-byte 64))
(prev-cumtime 0 :type (unsigned-byte 64)))
(defun make-lud-eta (lud &key (initial-size 0) cumtime prev-cumtime)
(declare ((or null (unsigned-byte 64)) cumtime prev-cumtime))
(%make-lud-eta :lud lud
:leaving-cols (make-array 0 :element-type 'fixnum)
:colstarts (make-array 1 :element-type 'fixnum :initial-element 0)
:rows (make-array initial-size :element-type 'fixnum)
:values (make-array initial-size :element-type 'csc-float)
:cumtime (or cumtime 0)
:prev-cumtime (or prev-cumtime 0)))
(defun refactor (matrix basis)
(let* ((start-time (read-tsc))
(lud (lu-factor matrix basis))
(timedelta (ldb (byte 64 0) (- (read-tsc) start-time))))
(make-lud-eta lud :cumtime timedelta)))
(defun add-eta! (lud-eta leaving-col sparse-vector)
(declare (optimize (speed 3)))
(symbol-macrolet ((enz (lud-eta-nz lud-eta))
(ecount (lud-eta-count lud-eta))
(vector-nz (sparse-vector-nz sparse-vector))
(leaving-cols (lud-eta-leaving-cols lud-eta))
(colstarts (lud-eta-colstarts lud-eta))
(rows (lud-eta-rows lud-eta))
(values (lud-eta-values lud-eta)))
(vector-set* leaving-cols ecount leaving-col)
(loop for k from enz below (+ enz vector-nz)
for pos across (sparse-vector-indices sparse-vector)
for value of-type csc-float across (sparse-vector-values sparse-vector)
do (vector-set* rows k pos)
(vector-set* values k value))
(incf ecount)
(incf enz vector-nz)
(vector-set* colstarts ecount enz)
lud-eta))
(defmacro with-updating-time (lud-eta &body body)
(let ((start-time (gensym "START"))
(timedelta (gensym "DELTA"))
(lude (gensym "LUDE")))
`(let ((,start-time (read-tsc))
(,lude ,lud-eta))
(unwind-protect (progn ,@body)
(let ((,timedelta (ldb (byte 64 0) (- (read-tsc) ,start-time))))
(setf (lud-eta-cumtime ,lude)
(ldb (byte 64 0) (+ (lud-eta-cumtime ,lude) ,timedelta))))))))
(defconstant +initial-size+ 16
"Is an initial size of temporary storage.")
(declaim (fvec *gauss-eta-values*)
(ivec *gauss-eta-tags* *gauss-eta-rows*)
((integer 0 #.most-positive-fixnum) *gauss-eta-tag*))
(defparameter *gauss-eta-values* (make-array +initial-size+ :element-type 'csc-float))
(defparameter *gauss-eta-tags* (make-array +initial-size+ :element-type 'fixnum))
(defparameter *gauss-eta-tag* 1)
(defparameter *gauss-eta-rows* (make-array +initial-size+ :element-type 'fixnum))
(defun gauss-eta! (lud-eta xb)
(declare (optimize (speed 3))
(sparse-vector xb))
(let* ((lud (lud-eta-lud lud-eta))
(m (lud-m lud)))
(extend-vectorf *gauss-eta-values* m)
(extend-vectorf *gauss-eta-tags* m)
(extend-vectorf *gauss-eta-rows* m))
(with-updating-time lud-eta
(let* ((vector-nz (sparse-vector-nz xb))
(vector-values (sparse-vector-values xb))
(vector-indices (sparse-vector-indices xb))
(tmp-values *gauss-eta-values*)
(tmp-tags *gauss-eta-tags*)
(tmp-rows *gauss-eta-rows*)
(tmp-end 0)
(leaving-cols (lud-eta-leaving-cols lud-eta))
(colstarts (lud-eta-colstarts lud-eta))
(values (lud-eta-values lud-eta))
(rows (lud-eta-rows lud-eta)))
(declare (fvec tmp-values vector-values)
(ivec tmp-rows tmp-tags vector-indices)
((mod #.array-dimension-limit) tmp-end))
(when (zerop (lud-eta-count lud-eta))
(return-from gauss-eta!))
(labels ((add-to-tmp (nzpos value)
(setf (aref tmp-values nzpos) value
(aref tmp-tags nzpos) *gauss-eta-tag*
(aref tmp-rows tmp-end) nzpos)
(incf tmp-end)))
(dotimes (k vector-nz)
(add-to-tmp (aref vector-indices k) (aref vector-values k)))
(dotimes (j (lud-eta-count lud-eta))
(let ((leaving-col (aref leaving-cols j))
leaving-col-k)
(loop for k from (aref colstarts j) below (aref colstarts (+ j 1))
for row = (aref rows k)
unless (= (aref tmp-tags row) *gauss-eta-tag*)
do (add-to-tmp row +zero+)
when (= row leaving-col)
do (setq leaving-col-k k))
(let ((coef (/ (aref tmp-values leaving-col) (aref values leaving-col-k))))
(unless (zerop coef)
(loop for k from (aref colstarts j) below leaving-col-k
for row = (aref rows k)
do (decf (aref tmp-values row) (* (aref values k) coef)))
(setf (aref tmp-values leaving-col) coef)
(loop for k from (+ leaving-col-k 1) below (aref colstarts (+ j 1))
for row = (aref rows k)
do (decf (aref tmp-values row) (* (aref values k) coef)))))))
(let ((nz 0))
(dotimes (k tmp-end)
(let ((row (aref tmp-rows k)))
(when (> (abs (aref tmp-values row)) +eps+)
(vector-set* vector-values nz (aref tmp-values row))
(vector-set* vector-indices nz row)
(incf nz))))
(setf (sparse-vector-values xb) vector-values
(sparse-vector-indices xb) vector-indices
(sparse-vector-nz xb) nz))
(incf *gauss-eta-tag*)
xb))))
(defun gauss-eta-transposed! (lud-eta xb)
(declare (optimize (speed 3)))
(let* ((lud (lud-eta-lud lud-eta))
(m (lud-m lud)))
(extend-vectorf *gauss-eta-values* m)
(extend-vectorf *gauss-eta-tags* m))
(with-updating-time lud-eta
(let* ((vector-nz (sparse-vector-nz xb))
(vector-values (sparse-vector-values xb))
(vector-indices (sparse-vector-indices xb))
(tmp-values *gauss-eta-values*)
(tmp-tags *gauss-eta-tags*)
(leaving-cols (lud-eta-leaving-cols lud-eta))
(colstarts (lud-eta-colstarts lud-eta))
(values (lud-eta-values lud-eta))
(rows (lud-eta-rows lud-eta)))
(declare (fvec tmp-values vector-values)
(ivec tmp-tags vector-indices))
(loop for j from (- (lud-eta-count lud-eta) 1) downto 0
for leaving-col = (aref leaving-cols j)
for leaving-col-k = nil
do (dotimes (k vector-nz)
(let ((row (aref vector-indices k)))
(when (= row leaving-col)
(setq leaving-col-k k))
(setf (aref tmp-values row) (aref vector-values k)
(aref tmp-tags row) *gauss-eta-tag*)))
(unless (= *gauss-eta-tag* (aref tmp-tags leaving-col))
(vector-set* vector-values vector-nz +zero+)
(vector-set* vector-indices vector-nz leaving-col)
(setq leaving-col-k vector-nz)
(incf vector-nz)
(setf (aref tmp-values leaving-col) +zero+
(aref tmp-tags leaving-col) *gauss-eta-tag*))
(let ((tmp (aref vector-values leaving-col-k))
leaving-col-eta-k)
(declare (csc-float tmp))
(loop for k from (aref colstarts j) below (aref colstarts (+ j 1))
for row = (aref rows k)
when (= row leaving-col)
do (setq leaving-col-eta-k k)
when (and (= (aref tmp-tags row) *gauss-eta-tag*)
(/= row leaving-col))
do (decf tmp (* (aref values k) (aref tmp-values row))))
(incf *gauss-eta-tag*)
(setf (aref vector-values leaving-col-k)
(/ tmp (aref values leaving-col-eta-k)))))
(setf (sparse-vector-values xb) vector-values
(sparse-vector-indices xb) vector-indices
(sparse-vector-nz xb) vector-nz)
xb)))
(declaim (fvec *sparse-solve-values*)
(ivec *sparse-solve-tags* *sparse-solve-rows*)
((integer 0 #.most-positive-fixnum) *sparse-solve-tag*))
(defparameter *sparse-solve-values* (make-array +initial-size+ :element-type 'csc-float))
(defparameter *sparse-solve-tags* (make-array +initial-size+ :element-type 'fixnum))
(defparameter *sparse-solve-tag* 1)
(defun sparse-solve! (lud-eta y)
(declare (optimize (speed 3)))
(let* ((lud (lud-eta-lud lud-eta))
(m (lud-m lud)))
(extend-vectorf *sparse-solve-values* m)
(extend-vectorf *sparse-solve-tags* m))
(with-updating-time lud-eta
(let* ((lud (lud-eta-lud lud-eta))
(rank (lud-rank lud))
(vector-nz (sparse-vector-nz y))
(vector-indices (sparse-vector-indices y))
(vector-values (sparse-vector-values y))
(tree (make-iterset))
;; vector values sorted w.r.t. the order of LU decomposition
(tmp-values *sparse-solve-values*)
(tmp-tags *sparse-solve-tags*)
(tag *sparse-solve-tag*))
(declare (fvec tmp-values vector-values)
(ivec tmp-tags vector-indices))
(labels ((add-to-tmp (lu-pos value)
(setf (aref tmp-values lu-pos) value
(aref tmp-tags lu-pos) tag)
(iterset-insert tree lu-pos)))
(let ((irowperm (lud-irowperm lud)))
(dotimes (k vector-nz)
(let ((lu-pos (aref irowperm (aref vector-indices k))))
(add-to-tmp lu-pos (aref vector-values k)))))
;; y := L^(-1)y
(let* ((lower (lud-lower lud))
(lower-colstarts (csc-colstarts lower))
(lower-rows (csc-rows lower))
(lower-values (csc-values lower))
(node (iterset-first tree)))
(loop while (and node (< (node-key node) rank))
for lu-pos = (node-key node)
for beta of-type csc-float = (aref tmp-values lu-pos)
for start = (aref lower-colstarts lu-pos)
for end = (aref lower-colstarts (+ lu-pos 1))
do (loop for k from start below end
for row = (aref lower-rows k)
unless (= (aref tmp-tags row) tag)
do (add-to-tmp row +zero+)
do (decf (aref tmp-values row) (* (aref lower-values k) beta)))
(setq node (node-next node))))
;; fill free variable
;; y := U^(-1)y
(let* ((upper (lud-upper lud))
(upper-colstarts (csc-colstarts upper))
(upper-rows (csc-rows upper))
(upper-values (csc-values upper))
(diagu (lud-diagu lud))
(node (iterset-last tree)))
(loop while (and node (>= (node-key node) rank))
for lu-pos = (node-key node)
do (setf (aref tmp-values lu-pos) +zero+)
(setq node (node-prev node)))
(loop while node
for lu-pos = (node-key node)
for beta of-type csc-float = (/ (aref tmp-values lu-pos) (aref diagu lu-pos))
for start = (aref upper-colstarts lu-pos)
for end = (aref upper-colstarts (+ lu-pos 1))
do (loop for k from start below end
for row = (aref upper-rows k)
unless (= (aref tmp-tags row) tag)
do (add-to-tmp row +zero+)
do (decf (aref tmp-values row) (* (aref upper-values k) beta)))
(setf (aref tmp-values lu-pos) beta)
(setq node (node-prev node)))))
;; update y
(let ((nz 0)
(node (iterset-first tree))
(colperm (lud-colperm lud)))
(declare ((mod #.array-dimension-limit) nz))
(loop while node
for lu-pos = (node-key node)
when (> (abs (aref tmp-values lu-pos)) +eps+)
do (vector-set* vector-values nz (aref tmp-values lu-pos))
(vector-set* vector-indices nz (aref colperm lu-pos))
(incf nz)
do (setq node (node-next node)))
(setf (sparse-vector-nz y) nz
(sparse-vector-indices y) vector-indices
(sparse-vector-values y) vector-values))
(gauss-eta! lud-eta y)
(incf *sparse-solve-tag*)
y)))
(defun sparse-solve-transposed! (lud-eta y)
(declare (optimize (speed 3)))
(let* ((lud (lud-eta-lud lud-eta))
(m (lud-m lud)))
(extend-vectorf *sparse-solve-values* m)
(extend-vectorf *sparse-solve-tags* m))
(with-updating-time lud-eta
(gauss-eta-transposed! lud-eta y)
(let* ((lud (lud-eta-lud lud-eta))
(rank (lud-rank lud))
(vector-nz (sparse-vector-nz y))
(vector-indices (sparse-vector-indices y))
(vector-values (sparse-vector-values y))
(tree (make-iterset))
;; vector values sorted w.r.t. the order of LU decomposition
(tmp-values *sparse-solve-values*)
(tmp-tags *sparse-solve-tags*)
(tag *sparse-solve-tag*))
(declare (fvec tmp-values vector-values)
(ivec tmp-tags vector-indices))
(labels ((add-to-tmp (lu-pos value)
(setf (aref tmp-values lu-pos) value
(aref tmp-tags lu-pos) tag)
(iterset-insert tree lu-pos)))
(let ((icolperm (lud-icolperm lud)))
(dotimes (k vector-nz)
(let ((lu-pos (aref icolperm (aref vector-indices k))))
(add-to-tmp lu-pos (aref vector-values k)))))
;; y := U^(-T)y
(let* ((tupper (lud-tupper lud))
(tupper-colstarts (csc-colstarts tupper))
(tupper-rows (csc-rows tupper))
(tupper-values (csc-values tupper))
(diagu (lud-diagu lud))
(node (iterset-first tree)))
(loop while (and node (< (node-key node) rank))
for lu-pos = (node-key node)
for beta of-type csc-float = (/ (aref tmp-values lu-pos) (aref diagu lu-pos))
for start = (aref tupper-colstarts lu-pos)
for end = (aref tupper-colstarts (+ 1 lu-pos))
do (loop for k from start below end
for row = (aref tupper-rows k)
unless (= (aref tmp-tags row) tag)
do (add-to-tmp row +zero+)
do (decf (aref tmp-values row) (* (aref tupper-values k) beta)))
(setf (aref tmp-values lu-pos) beta)
(setq node (node-next node))))
;; fill free varaiable
;; y := L^(-T)y
(let* ((tlower (lud-tlower lud))
(tlower-colstarts (csc-colstarts tlower))
(tlower-rows (csc-rows tlower))
(tlower-values (csc-values tlower))
(node (iterset-last tree)))
(loop while (and node (>= (node-key node) rank))
for lu-pos = (node-key node)
do (setf (aref tmp-values lu-pos) +zero+)
(setq node (node-prev node)))
(loop while node
for lu-pos = (node-key node)
for beta of-type csc-float = (aref tmp-values lu-pos)
for start = (aref tlower-colstarts lu-pos)
for end = (aref tlower-colstarts (+ lu-pos 1))
do (loop for k from start below end
for row = (aref tlower-rows k)
unless (= (aref tmp-tags row) tag)
do (add-to-tmp row +zero+)
do (decf (aref tmp-values row) (* (aref tlower-values k) beta)))
(setq node (node-prev node)))))
;; update y
(let ((nz 0)
(node (iterset-first tree))
(rowperm (lud-rowperm lud)))
(declare ((mod #.array-dimension-limit) nz))
(loop while node
for lu-pos = (node-key node)
when (> (abs (aref tmp-values lu-pos)) +eps+)
do (vector-set* vector-values nz (aref tmp-values lu-pos))
(vector-set* vector-indices nz (aref rowperm lu-pos))
(incf nz)
do (setq node (node-next node)))
(setf (sparse-vector-nz y) nz
(sparse-vector-indices y) vector-indices
(sparse-vector-values y) vector-values))
(incf *sparse-solve-tag*)
;; #>y
y)))
(declaim ((integer 0 #.most-positive-fixnum) *refactor-threshold*))
(defparameter *refactor-threshold* 100)
(defparameter *refactor-by-time* t)
(defun refactor-p (lud-eta leaving-col)
(declare (optimize (speed 3))
((mod #.array-dimension-limit) leaving-col))
(let ((count (lud-eta-count lud-eta))
(colstarts (lud-eta-colstarts lud-eta))
(rows (lud-eta-rows lud-eta)))
(assert (> count 0))
(prog1 (or (>= count *refactor-threshold*)
;; leaving column vanishes at the last E
(loop for k from (aref colstarts (- count 1)) below (aref colstarts count)
never (= (aref rows k) leaving-col))
(and *refactor-by-time*
(>= count 2)
(>= (floor (lud-eta-cumtime lud-eta) count)
(floor (lud-eta-prev-cumtime lud-eta) (- count 1)))))
(setf (lud-eta-prev-cumtime lud-eta)
(lud-eta-cumtime lud-eta)))))
(defpackage :cp/sparse-simplex
(:use :cl :cp/csc :cp/lud)
(:import-from :cp/csc #:+zero+ #:+one+ #:csc-float)
(:import-from :cp/lud #:vector-set* #:extend-vectorf)
(:use :cl)
(:export #:make-dictionary #:make-sparse-lp #:slp-restore
#:slp-lude #:slp-m #:slp-n #:slp-dictionary
#:slp-mat #:slp-tmat #:slp-b #:slp-c #:slp-x-basic #:slp-y-nonbasic
#:correct-x-basic! #:correct-y-nonbasic!
#:dictionary-basics #:dictionary-nonbasics #:dictionary-basic-flag
#:slp-primal! #:slp-dual! #:slp-dual-primal! #:slp-self-dual!)
(:documentation
"Provides two kinds of simplex method for sparse LP:
- two-phase (dual-then-primal) simplex method using Dantzig's pivot rule;
- TODO: parametric self-dual simplex method.
Usage procedure:
1. MAKE-SPARSE-LP
2. SPARSE-DUAL-PRIMAL!
3. SRARSE-LP-RESTORE
Reference:
Robert J. Vanderbei. Linear Programming: Foundations and Extensions. 5th edition."))
(in-package :cp/sparse-simplex)
(defconstant +eps-large+ 0)
(defconstant +eps-middle+ 0)
(defconstant +eps-small+ 0)
(defconstant +inf+ most-positive-fixnum)
(defconstant +neg-inf+ most-negative-fixnum)
(defun add-slack! (a)
(declare (optimize (speed 3))
(csc a))
(symbol-macrolet ((m (csc-m a))
(n (csc-n a))
(colstarts (csc-colstarts a))
(rows (csc-rows a))
(values (csc-values a))
(nz (csc-nz a)))
;; Add slack variable
(loop for row below m
for col of-type (mod #.array-dimension-limit) from n
for k from (aref colstarts n)
do (vector-set* values k +one+)
(vector-set* rows k row)
(vector-set* colstarts (+ col 1) (+ k 1))
(incf nz)
finally (setf n (+ n m)))
a))
(defstruct (dictionary (:constructor %make-dictionary))
(m nil :type (mod #.array-dimension-limit))
(n nil :type (mod #.array-dimension-limit))
(basics nil :type (simple-array fixnum (*)))
(nonbasics nil :type (simple-array fixnum (*)))
(basic-flag nil :type (simple-array fixnum (*))))
(defconstant +nan+ most-positive-fixnum)
(defun make-dictionary (m n basics)
"BASICS[i] = column number of the constraints matrix (of size m * (n + m))
which is currently basic.
0, 1, 2, ....., n-1, n, n+1, ..., n+m-1
|- non-slack vars -| |-- slack vars --|
"
(declare (optimize (speed 3))
(vector basics)
((mod #.array-dimension-limit) m n))
(assert (= m (length basics)))
(let* ((basics (coerce basics '(simple-array fixnum (*))))
(basic-flag (make-array (+ n m) :element-type 'fixnum :initial-element +nan+))
(nonbasics (make-array n :element-type 'fixnum)))
(dotimes (i m)
(let ((col (aref basics i)))
(setf (aref basic-flag col) i)))
(let ((j 0))
(dotimes (col (length basic-flag))
(when (= (aref basic-flag col) +nan+)
(setf (aref nonbasics j) col
(aref basic-flag col) (lognot j))
(incf j))))
(%make-dictionary :m m :n n :basics basics :nonbasics nonbasics :basic-flag basic-flag)))
(defun dictionary-swap! (dictionary col-out col-in)
(declare (optimize (speed 3))
((mod #.array-dimension-limit) col-out col-in))
(let* ((basics (dictionary-basics dictionary))
(nonbasics (dictionary-nonbasics dictionary))
(basic-flag (dictionary-basic-flag dictionary))
(i (aref basics col-out))
(j (aref nonbasics col-in)))
(setf (aref basics col-out) j
(aref nonbasics col-in) i
(aref basic-flag i) (lognot col-in)
(aref basic-flag j) col-out))
dictionary)
(declaim (ftype (function * (values csc-float &optional)) dot*))
(defun dot* (coefs x-basic basics)
(declare (optimize (speed 3))
((simple-array csc-float (*)) coefs x-basic)
((simple-array fixnum (*)) basics))
(let ((res +zero+))
(declare (csc-float res))
(dotimes (i (length x-basic))
(incf res (* (aref coefs (aref basics i)) (aref x-basic i))))
res))
(defconstant +initial-size+ 16)
(declaim ((simple-array csc-float (*)) *tmp-values*)
((simple-array fixnum (*)) *tmp-tags* *tmp-rows*)
((integer 0 #.most-positive-fixnum) *tmp-tag*))
(defparameter *tmp-values*
(make-array +initial-size+ :element-type 'csc-float))
(defparameter *tmp-tags*
(make-array +initial-size+ :element-type 'fixnum :initial-element 0))
(defparameter *tmp-tag* 1)
(defparameter *tmp-rows* (make-array +initial-size+ :element-type 'fixnum))
(defun tmat-times-vec! (tmat vec basic-flag &optional res)
(declare (optimize (speed 3))
(csc tmat)
(sparse-vector vec)
((or null sparse-vector) res)
((simple-array fixnum (*)) basic-flag))
(let ((m (csc-m tmat)))
(extend-vectorf *tmp-values* m)
(extend-vectorf *tmp-tags* m)
(extend-vectorf *tmp-rows* m))
(let ((tmp-values *tmp-values*)
(tmp-tags *tmp-tags*)
(tag *tmp-tag*)
(tmp-rows *tmp-rows*)
(vector-indices (sparse-vector-indices vec))
(vector-values (sparse-vector-values vec))
(tmat-values (csc-values tmat))
(tmat-colstarts (csc-colstarts tmat))
(tmat-rows (csc-rows tmat))
(end 0))
(declare ((mod #.array-dimension-limit)))
(dotimes (k1 (sparse-vector-nz vec))
(let ((col (aref vector-indices k1)))
(loop for k2 from (aref tmat-colstarts col) below (aref tmat-colstarts (+ col 1))
for row = (aref tmat-rows k2)
when (< (aref basic-flag row) 0)
do (unless (eql tag (aref tmp-tags row))
(setf (aref tmp-values row) +zero+
(aref tmp-tags row) tag
(aref tmp-rows end) row
end (+ end 1)))
(incf (aref tmp-values row)
(* (aref vector-values k1) (aref tmat-values k2))))))
(let* ((res (or res (make-sparse-vector end)))
(res-values (sparse-vector-values res))
(res-indices (sparse-vector-indices res))
(nz 0))
(declare ((simple-array csc-float (*)) res-values)
((simple-array fixnum (*)) res-indices)
((mod #.array-dimension-limit) nz))
(extend-vectorf res-values end)
(extend-vectorf res-indices end)
(dotimes (k end)
(let ((row (aref tmp-rows k)))
(when (> (abs (aref tmp-values row)) +eps-large+)
(setf (aref res-values nz) (aref tmp-values row)
(aref res-indices nz) (lognot (aref basic-flag row))
nz (+ nz 1)))))
(setf (sparse-vector-values res) res-values
(sparse-vector-indices res) res-indices
(sparse-vector-nz res) nz)
(incf *tmp-tag*)
res)))
(defstruct (sparse-lp (:constructor %make-sparse-lp)
(:conc-name slp-))
(m nil :type (mod #.array-dimension-limit))
(n nil :type (mod #.array-dimension-limit))
(mat nil :type csc)
(tmat nil :type csc)
(b nil :type (simple-array csc-float (*)))
(c nil :type (simple-array csc-float (*)))
(x-basic nil :type (simple-array csc-float (*)))
(y-nonbasic nil :type (simple-array csc-float (*)))
(dictionary nil :type dictionary)
(lude nil :type lud-eta))
(defun correct-x-basic! (lude x-basic)
(assert (zerop (lud-eta-count lude)))
(dense-solve! (lud-eta-lud lude) x-basic))
(defun correct-y-nonbasic! (sparse-lp)
(declare (optimize (speed 3)))
(let* ((lude (slp-lude sparse-lp))
(m (slp-m sparse-lp))
(n (slp-n sparse-lp))
(tmat (slp-tmat sparse-lp))
(c (slp-c sparse-lp))
(tmp (make-sparse-vector m))
(dictionary (slp-dictionary sparse-lp))
(basics (dictionary-basics dictionary))
(nonbasics (dictionary-nonbasics dictionary))
(basic-flag (dictionary-basic-flag dictionary))
(y-nonbasic (slp-y-nonbasic sparse-lp))
(tmp-values (sparse-vector-values tmp))
(tmp-indices (sparse-vector-indices tmp)))
(symbol-macrolet ((tmp-nz (sparse-vector-nz tmp)))
(dotimes (i m)
(let ((coef (aref c (aref basics i))))
(when (> (abs coef) +eps-large+)
(setf (aref tmp-values tmp-nz) coef
(aref tmp-indices tmp-nz) i)
(incf tmp-nz))))
(sparse-solve-transposed! lude tmp)
(let* ((tmp (tmat-times-vec! tmat tmp basic-flag))
(tmp-values (sparse-vector-values tmp))
(tmp-indices (sparse-vector-indices tmp)))
(dotimes (j n)
(setf (aref y-nonbasic j) (- (aref c (aref nonbasics j)))))
(dotimes (k (sparse-vector-nz tmp))
(incf (aref y-nonbasic (aref tmp-indices k)) (aref tmp-values k)))))
sparse-lp))
(defun make-sparse-lp (a b c &key (add-slack t) (dictionary nil supplied-p))
"Creates SPARSE-LP from a sparse matrix, which has the standard form: maximize
c'x subject to Ax <= b, x >= 0.
This function translates a given LP to an equality form Ax + w = b by adding
slack variables and changes A to (A E). If you want to give an equality form
directly, just disable ADD-SLACK.
You can set DICTIONARY to an arbitrary initial dictionary, but please note that
the consequence is undefined when it is rank-deficient."
(declare (optimize (speed 3))
(csc a)
((simple-array csc-float (*)) b c))
(let* ((m (csc-m a))
(n (if add-slack (csc-n a) (- (csc-n a) (csc-m a)))))
(assert (= m (length b)))
(when add-slack
(setq a (add-slack! a)))
;; Add coefficients for basic variables
(unless (= (length c) (+ m n))
(assert (= n (length c)))
(setq c (adjust-array c (the (mod #.array-dimension-limit) (+ n m))
:initial-element +zero+)))
(let* ((x-basic (make-array m :element-type 'csc-float))
(y-nonbasic (make-array n :element-type 'csc-float))
(dictionary (or dictionary
(let ((basics (make-array m :element-type 'fixnum)))
(dotimes (i m)
(setf (aref basics i) (+ n i)))
(make-dictionary m n basics))))
(basics (dictionary-basics dictionary))
(a-transposed (csc-transpose a)))
(unless supplied-p
(dotimes (j n)
(setf (aref y-nonbasic j) (- (aref c j)))))
(dotimes (i m)
(setf (aref x-basic i) (aref b i)))
(let* ((lude (refactor a basics))
(slp (%make-sparse-lp :m m :n n
:mat a :tmat a-transposed
:b b :c c
:x-basic x-basic
:y-nonbasic y-nonbasic
:dictionary dictionary
:lude lude)))
(when supplied-p
(correct-x-basic! lude x-basic)
(correct-y-nonbasic! slp))
slp))))
(defun slp-restore (sparse-lp)
"Restores the current solution of LP and returns five values: objective value,
primal solution, dual solution, values of primal slack variables, and values of
dual slack variables. (Note that they are not necessarily feasible solutions if
the current dictionary is not feasible.)"
(declare (optimize (speed 3)))
(let* ((m (slp-m sparse-lp))
(n (slp-n sparse-lp))
(c (slp-c sparse-lp))
(x-basic (slp-x-basic sparse-lp))
(y-nonbasic (slp-y-nonbasic sparse-lp))
(dictionary (slp-dictionary sparse-lp))
(basics (dictionary-basics dictionary))
(nonbasics (dictionary-nonbasics dictionary))
(x (make-array (+ n m) :element-type 'csc-float :initial-element +zero+))
(y (make-array (+ n m) :element-type 'csc-float :initial-element +zero+)))
(dotimes (i m)
(setf (aref x (aref basics i)) (aref x-basic i)))
(dotimes (i n)
(setf (aref y (aref nonbasics i)) (aref y-nonbasic i)))
(values (dot* c x-basic basics)
(subseq x 0 n)
(subseq y n)
(subseq x n)
(subseq y 0 n))))
(defun pick-negative (vector)
(declare (optimize (speed 3))
((simple-array csc-float (*)) vector))
(let ((min (- +eps-small+))
res)
(dotimes (i (length vector))
(when (< (aref vector i) min)
(setq min (aref vector i)
res i)))
res))
(defun ratio-test (x dx)
(declare (optimize (speed 3))
((simple-array csc-float (*)) x)
(sparse-vector dx))
(let ((min +inf+)
(dx-indices (sparse-vector-indices dx))
(dx-values (sparse-vector-values dx))
res)
(dotimes (k (sparse-vector-nz dx))
(when (> (aref dx-values k) +eps-large+)
(let* ((index (aref dx-indices k))
(rate (/ (aref x index) (aref dx-values k))))
(when (< rate min)
(setq min rate
res index)))))
res))
(defun slp-primal! (sparse-lp)
"Applies primal simplex method to SPARSE-LP and returns the terminal state:
:optimal or :unbounded."
(declare (optimize (speed 3)))
(let* ((m (slp-m sparse-lp))
(n (slp-n sparse-lp))
(x-basic (slp-x-basic sparse-lp))
(y-nonbasic (slp-y-nonbasic sparse-lp))
(dictionary (slp-dictionary sparse-lp))
(basics (dictionary-basics dictionary))
(nonbasics (dictionary-nonbasics dictionary))
(basic-flag (dictionary-basic-flag dictionary))
(mat (slp-mat sparse-lp))
(tmat (slp-tmat sparse-lp))
(dx (make-sparse-vector m))
(dy (make-sparse-vector n))
(tmp (make-sparse-vector m)))
(symbol-macrolet ((lude (slp-lude sparse-lp))
(dx-values (sparse-vector-values dx))
(dx-indices (sparse-vector-indices dx))
(dx-nz (sparse-vector-nz dx))
(dy-values (sparse-vector-values dy))
(dy-indices (sparse-vector-indices dy))
(dy-nz (sparse-vector-nz dy))
(tmp-values (sparse-vector-values tmp))
(tmp-indices (sparse-vector-indices tmp))
(tmp-nz (sparse-vector-nz tmp)))
(loop
;; find entering column
(let* ((col-in (pick-negative y-nonbasic)))
(unless col-in
(return :optimal))
;; dx_B := B^(-1)Ne_j (j = col-in)
(let ((acolstarts (csc-colstarts mat))
(arows (csc-rows mat))
(avalues (csc-values mat))
(j (aref nonbasics col-in)))
(loop for i from 0
for k from (aref acolstarts j) below (aref acolstarts (+ j 1))
do (setf (aref dx-values i) (aref avalues k)
(aref dx-indices i) (aref arows k))
finally (setq dx-nz i)))
(sparse-solve! lude dx)
;; find leaving column
(let ((col-out (ratio-test x-basic dx)))
(unless col-out
(return :unbounded))
;; dy_N := -(B^(-1)N)^Te_i (i = col-out)
(setf (aref tmp-values 0) (- +one+)
(aref tmp-indices 0) col-out
tmp-nz 1)
(sparse-solve-transposed! lude tmp)
(tmat-times-vec! tmat tmp basic-flag dy)
;; t := x_i/dx_i
;; s := y_j/dy_j
(let ((rate-t (loop for k below dx-nz
when (= (aref dx-indices k) col-out)
do (return (/ (aref x-basic col-out)
(aref dx-values k)))
finally (error "Huh?")))
(rate-s (loop for k below dy-nz
when (= (aref dy-indices k) col-in)
do (return (/ (aref y-nonbasic col-in)
(aref dy-values k)))
finally (error "Huh?"))))
;; y_N := y_N - s dy_N
;; y_i := s
;; x_B := x_B - t dx_B
;; x_j := t
(dotimes (k dy-nz)
(let ((j (aref dy-indices k)))
(decf (aref y-nonbasic j) (* rate-s (aref dy-values k)))))
(setf (aref y-nonbasic col-in) rate-s)
(dotimes (k dx-nz)
(let ((i (aref dx-indices k)))
(decf (aref x-basic i) (* rate-t (aref dx-values k)))))
(setf (aref x-basic col-out) rate-t)
;; Update basis
(dictionary-swap! dictionary col-out col-in)
(add-eta! lude col-out dx)
(when (refactor-p lude col-out)
(setq lude (refactor mat basics))))))))))
(defun slp-dual! (sparse-lp)
"Applies dual simplex method to SPARSE-LP and returns the terminal state:
:optimal or :infeasible."
(declare (optimize (speed 3)))
(let* ((m (slp-m sparse-lp))
(n (slp-n sparse-lp))
(x-basic (slp-x-basic sparse-lp))
(y-nonbasic (slp-y-nonbasic sparse-lp))
(dictionary (slp-dictionary sparse-lp))
(basics (dictionary-basics dictionary))
(nonbasics (dictionary-nonbasics dictionary))
(basic-flag (dictionary-basic-flag dictionary))
(mat (slp-mat sparse-lp))
(tmat (slp-tmat sparse-lp))
(dx (make-sparse-vector m))
(dy (make-sparse-vector n))
(tmp (make-sparse-vector m)))
(symbol-macrolet ((lude (slp-lude sparse-lp))
(dx-values (sparse-vector-values dx))
(dx-indices (sparse-vector-indices dx))
(dx-nz (sparse-vector-nz dx))
(dy-values (sparse-vector-values dy))
(dy-indices (sparse-vector-indices dy))
(dy-nz (sparse-vector-nz dy))
(tmp-values (sparse-vector-values tmp))
(tmp-indices (sparse-vector-indices tmp))
(tmp-nz (sparse-vector-nz tmp)))
(loop
;; find leaving column
(let ((col-out (pick-negative x-basic)))
(unless col-out
(return :optimal))
;; dy_N := -(B^(-1)N)^Te_i (i = col-out)
(setf (aref tmp-values 0) (- +one+)
(aref tmp-indices 0) col-out
tmp-nz 1)
(sparse-solve-transposed! lude tmp)
(tmat-times-vec! tmat tmp basic-flag dy)
;; find entering column
(let ((col-in (ratio-test y-nonbasic dy)))
(unless col-in
(return :infeasible))
;; dx_B := B^(-1)Ne_j (j = col-in)
(let ((acolstarts (csc-colstarts mat))
(arows (csc-rows mat))
(avalues (csc-values mat))
(j (aref nonbasics col-in)))
(loop for i from 0
for k from (aref acolstarts j) below (aref acolstarts (+ j 1))
do (setf (aref dx-values i) (aref avalues k)
(aref dx-indices i) (aref arows k))
finally (setq dx-nz i)))
(sparse-solve! lude dx)
;; t := x_i/dx_i
;; s := y_j/dy_j
(let ((rate-t (loop for k below dx-nz
when (= (aref dx-indices k) col-out)
do (return (/ (aref x-basic col-out)
(aref dx-values k)))
finally (error "Huh?")))
(rate-s (loop for k below dy-nz
when (= (aref dy-indices k) col-in)
do (return (/ (aref y-nonbasic col-in)
(aref dy-values k)))
finally (error "Huh?"))))
;; y_N := y_N - s dy_N
;; y_i := s
;; x_B := x_B - t dx_B
;; x_j := t
(dotimes (k dy-nz)
(let ((j (aref dy-indices k)))
(decf (aref y-nonbasic j) (* rate-s (aref dy-values k)))))
(setf (aref y-nonbasic col-in) rate-s)
(dotimes (k dx-nz)
(let ((i (aref dx-indices k)))
(decf (aref x-basic i) (* rate-t (aref dx-values k)))))
(setf (aref x-basic col-out) rate-t)
;; Update basis
(dictionary-swap! dictionary col-out col-in)
(add-eta! lude col-out dx)
(when (refactor-p lude col-out)
(setq lude (refactor mat basics))))))))))
(defun slp-dual-primal! (sparse-lp)
"Applies two-phase simplex method to SPARSE-LP and returns the terminal state:
:optimal, :unbounded, or :infeasible. "
(declare (optimize (speed 3)))
(let* ((n (slp-n sparse-lp))
(dictionary (slp-dictionary sparse-lp))
(nonbasics (dictionary-nonbasics dictionary))
(y-nonbasic (slp-y-nonbasic sparse-lp))
(c (slp-c sparse-lp)))
;; Set all the coefficiets of objective to negative values.
(dotimes (j n)
(let ((col (aref nonbasics j)))
(setf (aref y-nonbasic j)
(+ (max (if (< col n) (aref c col) +zero+) +one+)
(random +one+)))))
(let ((state-dual (slp-dual! sparse-lp)))
(correct-y-nonbasic! sparse-lp)
(unless (eql state-dual :optimal)
(return-from slp-dual-primal! state-dual))
(slp-primal! sparse-lp))))
;;;
;;; self-dual simplex method
;;;
(defun self-dual-ratio-test (x dx mu x-params)
(declare (optimize (speed 3))
((simple-array csc-float (*)) x x-params)
(sparse-vector dx)
(csc-float mu))
(let ((min +inf+)
(dx-indices (sparse-vector-indices dx))
(dx-values (sparse-vector-values dx))
res)
(dotimes (k (sparse-vector-nz dx))
(when (> (aref dx-values k) +eps-large+)
(let* ((index (aref dx-indices k))
(rate (/ (+ (aref x index) (* mu (aref x-params index)))
(aref dx-values k))))
(when (< rate min)
(setq min rate
res index)))))
res))
;; not tested
(defun slp-self-dual! (sparse-lp)
"Applies self-dual simplex method to SPARSE-LP and returns the terminal state:
:optimal, :infeasible, or :dual-infeasible."
(declare (optimize (speed 3)))
(let* ((m (slp-m sparse-lp))
(n (slp-n sparse-lp))
(x-basic (slp-x-basic sparse-lp))
(y-nonbasic (slp-y-nonbasic sparse-lp))
(dictionary (slp-dictionary sparse-lp))
(basics (dictionary-basics dictionary))
(nonbasics (dictionary-nonbasics dictionary))
(basic-flag (dictionary-basic-flag dictionary))
(mat (slp-mat sparse-lp))
(tmat (slp-tmat sparse-lp))
(dx (make-sparse-vector m))
(dy (make-sparse-vector n))
(tmp (make-sparse-vector m))
(x-params (make-array n :element-type 'csc-float :initial-element +zero+))
(y-params (make-array m :element-type 'csc-float :initial-element +zero+)))
(symbol-macrolet ((lude (slp-lude sparse-lp))
(dx-values (sparse-vector-values dx))
(dx-indices (sparse-vector-indices dx))
(dx-nz (sparse-vector-nz dx))
(dy-values (sparse-vector-values dy))
(dy-indices (sparse-vector-indices dy))
(dy-nz (sparse-vector-nz dy))
(tmp-values (sparse-vector-values tmp))
(tmp-indices (sparse-vector-indices tmp))
(tmp-nz (sparse-vector-nz tmp)))
;; initialize parameters for x and y
(dotimes (j n)
(let ((colstarts (csc-colstarts mat))
(rows (csc-rows mat))
(values (csc-values mat)))
(loop for k from (aref colstarts j) below (aref colstarts (+ j 1))
for a2 of-type csc-float = (expt (aref values k) 2)
do (incf (aref x-params (aref rows k)) a2)
(incf (aref y-params j) a2))))
(map-into x-params
(lambda (x) (+ (random 1d0) (sqrt (the (double-float 0d0) x))))
x-params)
(map-into y-params
(lambda (x) (+ (random 1d0) (sqrt (the (double-float 0d0) x))))
y-params)
(loop
(let ((mu +neg-inf+)
col-in
col-out)
(dotimes (j n)
(when (and (> (aref y-params j) +eps-small+)
(< mu (/ (- (aref y-nonbasic j)) (aref y-params j))))
(setq mu (/ (- (aref y-nonbasic j)) (aref y-params j))
col-in j)))
(dotimes (i m)
(when (and (> (aref x-params i) +eps-small+)
(< mu (/ (- (aref x-basic i)) (aref x-params i))))
(setq mu (/ (- (aref x-basic i)) (aref x-params i))
col-out i)))
(when (<= mu +eps-middle+)
(return :optimal))
(assert (or (and col-in (not col-out))
(and (not col-in) col-out)))
(when col-out
;; dy_N := -(B^(-1)N)^T e_i where i is leaving column
(setf (aref tmp-values 0) (- +one+)
(aref tmp-indices 0) col-out
tmp-nz 1)
(sparse-solve-transposed! lude tmp)
(tmat-times-vec! tmat tmp basic-flag dy)
(let ((col-in (self-dual-ratio-test y-nonbasic dy mu y-params)))
(unless col-in
(return :infeasible))
;; dx_B := B^(-1)Ne_j where j is entering column
(let* ((j (aref nonbasics col-in))
(colstarts (csc-colstarts mat))
(rows (csc-rows mat))
(values (csc-values mat))
(start (aref colstarts j))
(end (aref colstarts (+ j 1))))
(loop for k from start below end
for i of-type (mod #.array-dimension-limit) = (- k start)
do (setf (aref dx-values i) (aref values k)
(aref dx-indices i) (aref rows k)))
(setq dx-nz (- end start))
(sparse-solve! lude dx))))
(when col-in
;; dx_B := B^(-1)Ne_j where j is entering column
(let* ((j (aref nonbasics col-in))
(colstarts (csc-colstarts mat))
(rows (csc-rows mat))
(values (csc-values mat))
(start (aref colstarts j))
(end (aref colstarts (+ j 1))))
(loop for k from start below end
for i of-type (mod #.array-dimension-limit) = (- k start)
do (setf (aref dx-values i) (aref values k)
(aref dx-indices i) (aref rows k)))
(setq dx-nz (- end start))
(sparse-solve! lude dx)
(let ((col-out (self-dual-ratio-test x-basic dx mu x-params)))
(unless col-out
(return :dual-infeasible))
;; dy_N := -(B^(-1)N)^T e_i where i is leaving column
(setf (aref tmp-values 0) (- +one+)
(aref tmp-indices 0) col-out
tmp-nz 1)
(sparse-solve-transposed! lude tmp)
(tmat-times-vec! tmat tmp basic-flag dy))))
;; t := x_i/dx_i
;; tparam := xparam_i/dx_i
;; s := y_j/dy_j
;; sparam := y_param_j/dy_j
(multiple-value-bind (rate-t rate-tparam)
(dotimes (k dx-nz (error "Huh?"))
(when (= (aref dx-indices k) col-out)
(return (values (/ (aref x-basic col-out)
(aref dx-values k))
(/ (aref x-params col-out)
(aref dx-values k))))))
(multiple-value-bind (rate-s rate-sparam)
(dotimes (k dy-nz (error "Huh?"))
(when (= (aref dy-indices k) col-in)
(return (values (/ (aref y-nonbasic col-in)
(aref dy-values k))
(/ (aref y-params col-in)
(aref dy-values k))))))
;; y_N := y_N - s dy_N
;; yparam := yparam - s dy_N
;; y_i := s
;; x_B := x_B - t dx_B
;; xparam_B := xparam - t dx_B
;; x_j := t
(dotimes (k dy-nz)
(let ((j (aref dy-indices k)))
(decf (aref y-nonbasic j) (* rate-s (aref dy-values k)))
(decf (aref y-params j) (* rate-sparam (aref dy-values k)))))
(setf (aref y-nonbasic col-in) rate-s
(aref y-params col-in) rate-sparam)
(dotimes (k dx-nz)
(let ((i (aref dx-indices k)))
(decf (aref x-basic i) (* rate-t (aref dx-values k)))
(decf (aref x-params i) (* rate-tparam (aref dx-values k)))))
(setf (aref x-basic col-out) rate-t
(aref x-params col-out) rate-tparam)
;; Update basis
(dictionary-swap! dictionary col-out col-in)
(add-eta! lude col-out dx)
(when (refactor-p lude col-out)
(setq lude (refactor mat basics))))))))))
(defpackage :cp/read-fixnum
(:use :cl)
(:export #:read-fixnum))
(in-package :cp/read-fixnum)
(declaim (ftype (function * (values fixnum &optional)) read-fixnum))
(defun read-fixnum (&optional (in *standard-input*))
"NOTE: cannot read -2^62"
(macrolet ((%read-byte ()
`(the (unsigned-byte 8)
#+swank (char-code (read-char in nil #\Nul))
#-swank (sb-impl::ansi-stream-read-byte in nil #.(char-code #\Nul) nil))))
(let* ((minus nil)
(result (loop (let ((byte (%read-byte)))
(cond ((<= 48 byte 57)
(return (- byte 48)))
((zerop byte) ; #\Nul
(error "Read EOF or #\Nul."))
((= byte #.(char-code #\-))
(setq minus t)))))))
(declare ((integer 0 #.most-positive-fixnum) result))
(loop
(let* ((byte (%read-byte)))
(if (<= 48 byte 57)
(setq result (+ (- byte 48)
(* 10 (the (integer 0 #.(floor most-positive-fixnum 10))
result))))
(return (if minus (- result) result))))))))
;; BEGIN_USE_PACKAGE
(eval-when (:compile-toplevel :load-toplevel :execute)
(use-package :cp/read-fixnum :cl-user))
(eval-when (:compile-toplevel :load-toplevel :execute)
(use-package :cp/sparse-simplex :cl-user))
(eval-when (:compile-toplevel :load-toplevel :execute)
(use-package :cp/csc :cl-user))
(in-package :cl-user)
;;;
;;; Body
;;;
(defun main ()
(let* ((n (read))
(coo (make-coo (+ 1 (* 3 (- n 1))) (* 3 (- n 1))))
(as (make-array n :element-type 'uint62 :initial-element 0))
(bs (make-array n :element-type 'uint62 :initial-element 0))
(coefs (make-array (* 3 (- n 1)) :element-type 'rational :initial-element 0))
(rhs (make-array (+ 1 (* 3 (- n 1))) :element-type 'rational :initial-element 0)))
(dotimes (i n)
(setf (aref as i) (read-fixnum)))
(dotimes (i n)
(setf (aref bs i) (read-fixnum)))
(labels ((%getvar (index op)
(ecase op
(:y+ index)
(:y- (+ index n -1))
(:z (+ index (* 2 (- n 1)))))))
(let ((row 0))
(dotimes (i n)
(when (< i (- n 1))
(coo-insert! coo row (%getvar i :y+) 1)
(coo-insert! coo row (%getvar i :y-) -1))
(when (>= (- i 1) 0)
(coo-insert! coo row (%getvar (- i 1) :y+) -1)
(coo-insert! coo row (%getvar (- i 1) :y-) 1))
(setf (aref rhs row) (- (aref as i) (aref bs i)))
;; #>(csc-to-array (make-csc-from-coo coo))
(incf row))
;; #>row
(dotimes (i (- n 1))
(coo-insert! coo row (%getvar i :y+) 1)
(coo-insert! coo row (%getvar i :y-) -1)
(coo-insert! coo row (%getvar i :z) -1)
(incf row)
(coo-insert! coo row (%getvar i :y+) -1)
(coo-insert! coo row (%getvar i :y-) 1)
(coo-insert! coo row (%getvar i :z) -1)
(incf row))
(dotimes (i (- n 1))
(setf (aref coefs (%getvar i :z)) -1))
;; #>row
)
;; #>rhs
(let* ((csc (make-csc-from-coo coo))
;; (mat (csc-to-array csc))
(lp (make-sparse-lp csc rhs coefs)))
;; #>mat
(slp-dual-primal! lp)
(println (- (slp-restore lp)))))))
#-swank (main)
;;;
;;; Test
;;;
#+swank
(progn
(defparameter *lisp-file-pathname* (uiop:current-lisp-file-pathname))
(setq *default-pathname-defaults* (uiop:pathname-directory-pathname *lisp-file-pathname*))
(uiop:chdir *default-pathname-defaults*)
(defparameter *dat-pathname* (uiop:merge-pathnames* "test.dat" *lisp-file-pathname*))
(defparameter *problem-url* "https://atcoder.jp/contests/kupc2016/tasks/kupc2016_h"))
#+swank
(defun gen-dat ()
(uiop:with-output-file (out *dat-pathname* :if-exists :supersede)
(format out "")))
#+swank
(defun bench (&optional (out (make-broadcast-stream)))
(time (run *dat-pathname* out)))
#+(and sbcl (not swank))
(eval-when (:compile-toplevel)
(when sb-c::*undefined-warnings*
(error "undefined warnings: ~{~A~^ ~}" sb-c::*undefined-warnings*)))
;; To run: (5am:run! :sample)
#+swank
(5am:test :sample
(5am:is
(equal "2
"
(run "2
1 5
3 1
" nil)))
(5am:is
(equal "6
"
(run "5
1 2 3 4 5
3 3 1 1 1
" nil)))
(5am:is
(equal "48
"
(run "27
46 3 4 2 10 2 5 2 6 7 20 13 9 49 3 8 4 3 19 9 3 5 4 13 9 5 7
10 2 5 6 2 6 3 2 2 5 3 11 13 2 2 7 7 3 9 5 13 4 17 2 2 2 4
" nil)))
(5am:is
(equal "6302172
"
(run "18
3878348 423911 8031742 1035156 24256 10344593 19379 3867285 4481365 1475384 1959412 1383457 164869 4633165 6674637 9732852 10459147 2810788
1236501 770807 4003004 131688 1965412 266841 3980782 565060 816313 192940 541896 250801 217586 3806049 1220252 1161079 31168 2008961
" nil)))
(5am:is
(equal "1234567890
"
(run "2
1 99999999999
1234567891 1
" nil))))
提出情報
| 提出日時 | |
|---|---|
| 問題 | H - 壁壁壁壁壁壁壁 |
| ユーザ | sansaqua |
| 言語 | Common Lisp (SBCL 2.0.3) |
| 得点 | 60 |
| コード長 | 99250 Byte |
| 結果 | TLE |
| 実行時間 | 2218 ms |
| メモリ | 518712 KiB |
コンパイルエラー
; file: /imojudge/sandbox/Main.lisp ; in: DEFUN CSC-TO-ARRAY ; (AREF CP/CSC::COLPERM CP/CSC::COL) ; ==> ; (SB-KERNEL:HAIRY-DATA-VECTOR-REF/CHECK-BOUNDS ARRAY SB-INT:INDEX) ; ; note: unable to ; optimize ; because: ; Upgraded element type of array is not known at compile time. ; (AREF CP/CSC::ROWPERM CP/CSC::ROW) ; ==> ; (SB-KERNEL:HAIRY-DATA-VECTOR-REF/CHECK-BOUNDS ARRAY SB-INT:INDEX) ; ; note: unable to ; optimize ; because: ; Upgraded element type of array is not known at compile time. ; (MAKE-ARRAY (LIST CP/CSC::M CP/CSC::N) :ELEMENT-TYPE 'CP/CSC::CSC-FLOAT ; :INITIAL-ELEMENT CP/CSC::+ZERO+) ; --> LET* THE ; ==> ; (* #:D0 #:D1) ; ; note: forced to do GENERIC-* (cost 30) ; unable to do inline fixnum arithmetic (cost 2) because: ; The result is a (VALUES (MOD 21267647932558653929567424817066410001) ; &OPTIONAL), not a (VALUES FIXNUM &REST T). ; unable to do inline (signed-byte 64) arithmetic (cost 4) because: ; The result is a (VALUES (MOD 21267647932558653929567424817066410001) ; &OPTIONAL), not a (VALUES (SIGNED-BYTE 64) &REST ; T). ; etc. ; file: /imojudge/sandbox/Main.lisp ; in: DEFUN MAKE-CSC-FROM-ARRAY ; (AREF ARRAY CP/CSC::ROW CP/CSC::COL) ; --> LET* ; ==> ; (SB-KERNEL:HAIRY-DATA-VECTOR-REF ARRAY SB-INT:INDEX) ; ; note: unable to ; avoid runtime dispatch on array element type ; due to type uncertainty: ; The first argument is a (ARRAY * (* *)), not a SIMPLE-ARRAY. ; (ZEROP CP/CSC::VALUE) ; ==> ; (= CP/CSC::VALUE 0) ; ; note: unable to ; open-code FLOAT to RATIONAL comparison ; due to type uncertainty: ; The first argument is a NUMBER, not a FLOAT. ; ; note: unable to ; open-code FLOAT to RATIONAL comparison ; due to type uncertainty: ; The first argument is a NUMBER, not a (OR (COMPLEX SINGLE-FLOAT) ; (COMPLEX DOUBLE-FLOAT)). ; ; note: unab...
ジャッジ結果
| セット名 | Subtask1 | Subtask2 | All | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| 得点 / 配点 | 30 / 30 | 30 / 30 | 0 / 140 | ||||||||
| 結果 |
|
|
|
| セット名 | テストケース |
|---|---|
| Subtask1 | 00_00_sample.txt, 00_01_sample.txt, 00_02_sample.txt, 00_03_random.txt, 00_04_random.txt, 00_05_random.txt, 00_06_random.txt, 00_07_random.txt, 00_08_random.txt, 00_09_random.txt, 00_10_random.txt, 00_11_random.txt, 00_12_random.txt, 00_13_random.txt, 00_14_random.txt, 00_15_random.txt, 00_16_random.txt, 00_17_random.txt, 00_18_random.txt, 00_19_random.txt, 00_20_random.txt, 00_21_random.txt, 00_22_random.txt, 00_23_freedom.txt, 00_24_freedom.txt, 00_25_freedom.txt, 00_26_full.txt, 00_27_full.txt, 00_28_full.txt, 00_29_min.txt, 00_30_min.txt, 00_31_min.txt, 00_32_max.txt, 00_33_max.txt, 00_34_max.txt, 00_35_max.txt, 00_36_max.txt |
| Subtask2 | 01_37_sample.txt, 01_38_sample.txt, 01_39_random.txt, 01_40_random.txt, 01_41_random.txt, 01_42_random.txt, 01_43_random.txt, 01_44_random.txt, 01_45_random.txt, 01_46_random.txt, 01_47_random.txt, 01_48_random.txt, 01_49_random.txt, 01_50_random.txt, 01_51_random.txt, 01_52_random.txt, 01_53_random.txt, 01_54_random.txt, 01_55_random.txt, 01_56_random.txt, 01_57_random.txt, 01_58_random.txt, 01_59_freedom.txt, 01_60_freedom.txt, 01_61_freedom.txt, 01_62_full.txt, 01_63_full.txt, 01_64_full.txt, 01_65_min.txt, 01_66_min.txt, 01_67_min.txt, 01_68_max.txt, 01_69_max.txt, 01_70_max.txt, 01_71_max.txt, 01_72_max.txt |
| All | 00_00_sample.txt, 00_01_sample.txt, 00_02_sample.txt, 00_03_random.txt, 00_04_random.txt, 00_05_random.txt, 00_06_random.txt, 00_07_random.txt, 00_08_random.txt, 00_09_random.txt, 00_10_random.txt, 00_11_random.txt, 00_12_random.txt, 00_13_random.txt, 00_14_random.txt, 00_15_random.txt, 00_16_random.txt, 00_17_random.txt, 00_18_random.txt, 00_19_random.txt, 00_20_random.txt, 00_21_random.txt, 00_22_random.txt, 00_23_freedom.txt, 00_24_freedom.txt, 00_25_freedom.txt, 00_26_full.txt, 00_27_full.txt, 00_28_full.txt, 00_29_min.txt, 00_30_min.txt, 00_31_min.txt, 00_32_max.txt, 00_33_max.txt, 00_34_max.txt, 00_35_max.txt, 00_36_max.txt, 01_37_sample.txt, 01_38_sample.txt, 01_39_random.txt, 01_40_random.txt, 01_41_random.txt, 01_42_random.txt, 01_43_random.txt, 01_44_random.txt, 01_45_random.txt, 01_46_random.txt, 01_47_random.txt, 01_48_random.txt, 01_49_random.txt, 01_50_random.txt, 01_51_random.txt, 01_52_random.txt, 01_53_random.txt, 01_54_random.txt, 01_55_random.txt, 01_56_random.txt, 01_57_random.txt, 01_58_random.txt, 01_59_freedom.txt, 01_60_freedom.txt, 01_61_freedom.txt, 01_62_full.txt, 01_63_full.txt, 01_64_full.txt, 01_65_min.txt, 01_66_min.txt, 01_67_min.txt, 01_68_max.txt, 01_69_max.txt, 01_70_max.txt, 01_71_max.txt, 01_72_max.txt, 02_100_min.txt, 02_101_min.txt, 02_102_max.txt, 02_103_max.txt, 02_104_max.txt, 02_105_max.txt, 02_106_max.txt, 02_73_random.txt, 02_74_random.txt, 02_75_random.txt, 02_76_random.txt, 02_77_random.txt, 02_78_random.txt, 02_79_random.txt, 02_80_random.txt, 02_81_random.txt, 02_82_random.txt, 02_83_random.txt, 02_84_random.txt, 02_85_random.txt, 02_86_random.txt, 02_87_random.txt, 02_88_random.txt, 02_89_random.txt, 02_90_random.txt, 02_91_random.txt, 02_92_random.txt, 02_93_freedom.txt, 02_94_freedom.txt, 02_95_freedom.txt, 02_96_full.txt, 02_97_full.txt, 02_98_full.txt, 02_99_min.txt |
| ケース名 | 結果 | 実行時間 | メモリ |
|---|---|---|---|
| 00_00_sample.txt | AC | 27 ms | 28172 KiB |
| 00_01_sample.txt | AC | 24 ms | 28228 KiB |
| 00_02_sample.txt | AC | 26 ms | 28452 KiB |
| 00_03_random.txt | AC | 18 ms | 28324 KiB |
| 00_04_random.txt | AC | 21 ms | 29044 KiB |
| 00_05_random.txt | AC | 22 ms | 28712 KiB |
| 00_06_random.txt | AC | 25 ms | 28488 KiB |
| 00_07_random.txt | AC | 27 ms | 28212 KiB |
| 00_08_random.txt | AC | 23 ms | 28396 KiB |
| 00_09_random.txt | AC | 20 ms | 28460 KiB |
| 00_10_random.txt | AC | 20 ms | 28544 KiB |
| 00_11_random.txt | AC | 19 ms | 28376 KiB |
| 00_12_random.txt | AC | 23 ms | 30740 KiB |
| 00_13_random.txt | AC | 20 ms | 28168 KiB |
| 00_14_random.txt | AC | 22 ms | 28972 KiB |
| 00_15_random.txt | AC | 19 ms | 28236 KiB |
| 00_16_random.txt | AC | 27 ms | 28432 KiB |
| 00_17_random.txt | AC | 20 ms | 28368 KiB |
| 00_18_random.txt | AC | 20 ms | 28572 KiB |
| 00_19_random.txt | AC | 21 ms | 29300 KiB |
| 00_20_random.txt | AC | 20 ms | 28424 KiB |
| 00_21_random.txt | AC | 20 ms | 28812 KiB |
| 00_22_random.txt | AC | 19 ms | 28312 KiB |
| 00_23_freedom.txt | AC | 25 ms | 28484 KiB |
| 00_24_freedom.txt | AC | 20 ms | 28416 KiB |
| 00_25_freedom.txt | AC | 18 ms | 28648 KiB |
| 00_26_full.txt | AC | 19 ms | 28776 KiB |
| 00_27_full.txt | AC | 24 ms | 29776 KiB |
| 00_28_full.txt | AC | 28 ms | 29152 KiB |
| 00_29_min.txt | AC | 20 ms | 28240 KiB |
| 00_30_min.txt | AC | 18 ms | 28184 KiB |
| 00_31_min.txt | AC | 18 ms | 28188 KiB |
| 00_32_max.txt | AC | 26 ms | 28612 KiB |
| 00_33_max.txt | AC | 28 ms | 31584 KiB |
| 00_34_max.txt | AC | 19 ms | 28588 KiB |
| 00_35_max.txt | AC | 26 ms | 29868 KiB |
| 00_36_max.txt | AC | 27 ms | 29864 KiB |
| 01_37_sample.txt | AC | 21 ms | 28224 KiB |
| 01_38_sample.txt | AC | 18 ms | 28240 KiB |
| 01_39_random.txt | AC | 272 ms | 89160 KiB |
| 01_40_random.txt | AC | 688 ms | 96828 KiB |
| 01_41_random.txt | AC | 54 ms | 41712 KiB |
| 01_42_random.txt | AC | 38 ms | 34796 KiB |
| 01_43_random.txt | AC | 19 ms | 28140 KiB |
| 01_44_random.txt | AC | 235 ms | 90048 KiB |
| 01_45_random.txt | AC | 34 ms | 31992 KiB |
| 01_46_random.txt | AC | 91 ms | 58964 KiB |
| 01_47_random.txt | AC | 26 ms | 28588 KiB |
| 01_48_random.txt | AC | 34 ms | 34288 KiB |
| 01_49_random.txt | AC | 103 ms | 64760 KiB |
| 01_50_random.txt | AC | 71 ms | 48784 KiB |
| 01_51_random.txt | AC | 19 ms | 29184 KiB |
| 01_52_random.txt | AC | 107 ms | 65808 KiB |
| 01_53_random.txt | AC | 65 ms | 48028 KiB |
| 01_54_random.txt | AC | 38 ms | 34456 KiB |
| 01_55_random.txt | AC | 86 ms | 56148 KiB |
| 01_56_random.txt | AC | 152 ms | 83956 KiB |
| 01_57_random.txt | AC | 113 ms | 68696 KiB |
| 01_58_random.txt | AC | 118 ms | 71876 KiB |
| 01_59_freedom.txt | AC | 29 ms | 29684 KiB |
| 01_60_freedom.txt | AC | 32 ms | 31080 KiB |
| 01_61_freedom.txt | AC | 22 ms | 31112 KiB |
| 01_62_full.txt | AC | 197 ms | 81424 KiB |
| 01_63_full.txt | AC | 320 ms | 89200 KiB |
| 01_64_full.txt | AC | 197 ms | 82060 KiB |
| 01_65_min.txt | AC | 18 ms | 28236 KiB |
| 01_66_min.txt | AC | 17 ms | 28248 KiB |
| 01_67_min.txt | AC | 20 ms | 28184 KiB |
| 01_68_max.txt | AC | 87 ms | 56444 KiB |
| 01_69_max.txt | AC | 100 ms | 61516 KiB |
| 01_70_max.txt | AC | 83 ms | 54704 KiB |
| 01_71_max.txt | AC | 60 ms | 45608 KiB |
| 01_72_max.txt | AC | 206 ms | 88136 KiB |
| 02_100_min.txt | AC | 26 ms | 28164 KiB |
| 02_101_min.txt | AC | 17 ms | 28188 KiB |
| 02_102_max.txt | TLE | 2218 ms | 464428 KiB |
| 02_103_max.txt | TLE | 2217 ms | 462636 KiB |
| 02_104_max.txt | TLE | 2213 ms | 435500 KiB |
| 02_105_max.txt | TLE | 2218 ms | 470100 KiB |
| 02_106_max.txt | TLE | 2217 ms | 466108 KiB |
| 02_73_random.txt | TLE | 2216 ms | 415820 KiB |
| 02_74_random.txt | TLE | 2217 ms | 464076 KiB |
| 02_75_random.txt | TLE | 2212 ms | 271900 KiB |
| 02_76_random.txt | TLE | 2213 ms | 464716 KiB |
| 02_77_random.txt | TLE | 2210 ms | 302216 KiB |
| 02_78_random.txt | TLE | 2214 ms | 353756 KiB |
| 02_79_random.txt | TLE | 2210 ms | 201872 KiB |
| 02_80_random.txt | TLE | 2212 ms | 446280 KiB |
| 02_81_random.txt | TLE | 2212 ms | 449056 KiB |
| 02_82_random.txt | TLE | 2214 ms | 374396 KiB |
| 02_83_random.txt | TLE | 2213 ms | 278024 KiB |
| 02_84_random.txt | TLE | 2213 ms | 347460 KiB |
| 02_85_random.txt | TLE | 2213 ms | 436228 KiB |
| 02_86_random.txt | TLE | 2214 ms | 346228 KiB |
| 02_87_random.txt | TLE | 2210 ms | 207640 KiB |
| 02_88_random.txt | TLE | 2212 ms | 465748 KiB |
| 02_89_random.txt | TLE | 2215 ms | 518712 KiB |
| 02_90_random.txt | TLE | 2210 ms | 237908 KiB |
| 02_91_random.txt | TLE | 2209 ms | 150804 KiB |
| 02_92_random.txt | TLE | 2212 ms | 311488 KiB |
| 02_93_freedom.txt | AC | 439 ms | 226408 KiB |
| 02_94_freedom.txt | AC | 235 ms | 127612 KiB |
| 02_95_freedom.txt | AC | 334 ms | 175620 KiB |
| 02_96_full.txt | TLE | 2208 ms | 139356 KiB |
| 02_97_full.txt | TLE | 2216 ms | 472816 KiB |
| 02_98_full.txt | TLE | 2212 ms | 478604 KiB |
| 02_99_min.txt | AC | 28 ms | 28200 KiB |