(DEFINE-FILE-INFO READTABLE "XCL" PACKAGE (DEFPACKAGE "BLAS" ( USE "LISP" "PT")))
(il:filecreated " 5-Dec-88 19:03:20" il:{qv}<idl>next>blas-fns.\;13 12925
il:|changes| il:|to:| (il:functions givens-rotate) (il:vars il:blas-fnscoms)
il:|previous| il:|date:| " 1-Dec-88 18:46:17" il:{qv}<idl>next>blas-fns.\;12)
; Copyright (c) 1988 by Xerox Corporation. All rights reserved.
(il:prettycomprint il:blas-fnscoms)
(il:rpaqq il:blas-fnscoms ((il:functions check->=-0 check->-0 compute-count compute-limit limit-inside-p) (il:functions a+x ax x+y ax+y) (il:functions sum-x max-index max-abs-index <xy> l2-norm) (il:functions compute-givens-coefficients givens-rotate) (il:functions copy-x-to-y swap-x-and-y fill-x) (eval-when (load) (il:p (export (quote (a+x ax x+y ax+y sum-x max-index max-abs-index <xy> l2-norm compute-givens-coefficients givens-rotate copy-x-to-y swap-x-and-y fill-x)) (find-package "BLAS")))) (xcl:file-environments "BLAS-FNS")))
(defmacro check->=-0 (x) (il:bquote (if (not (typep (il:\\\, x) (quote (integer 0)))) (error "Not a non-negative integer: ~s" (il:\\\, x)))))
(defmacro check->-0 (x) (il:bquote (if (not (typep (il:\\\, x) (quote (integer 1)))) (error "Not a positive integer: ~s" (il:\\\, x)))))
(defmacro compute-count (offset skip size) (xcl:once-only (offset skip size) (il:bquote (if (eq (il:\\\, skip) 1) (- (il:\\\, size) (il:\\\, offset)) (ceiling (- (il:\\\, size) (il:\\\, offset)) (il:\\\, skip))))))
(defmacro compute-limit (offset count skip) (il:bquote (+ (il:\\\, offset) (* (il:\\\, count) (il:\\\, skip)))))
(defmacro limit-inside-p (limit skip count) (il:bquote (<= (- (il:\\\, limit) (il:\\\, skip)) (il:\\\, count))))
(defun a+x (a x &key (offset 0) (skip 1) count) (il:* il:|;;;| "Implements BLAS scalar plus vector. COUNT defaults to largest vector commensurate with OFFSET and SKIP. Returns X which is side-effected.") (let ((size (length x))) (il:* il:|;;| "Arg checking") (check->=-0 offset) (check->-0 skip) (when (null count) (setq count (compute-count offset skip size)) (check->=-0 count)) (let ((limit (compute-limit offset count skip))) (if (not (limit-inside-p limit skip size)) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((i offset (+ i skip))) ((eq i limit)) (setf (aref x i) (+ a (aref x i))))) x))
(defun ax (a x &key (offset 0) (skip 1) count) (il:* il:|;;;| "Implements BLAS scalar times vector. COUNT defaults to largest vector commensurate with OFFSET and SKIP. Returns X, which is side effected.") (let ((size (length x))) (il:* il:|;;| "Arg checking") (check->=-0 offset) (check->-0 skip) (when (null count) (setq count (compute-count offset skip size)) (check->=-0 count)) (let ((limit (compute-limit offset count skip))) (if (not (limit-inside-p limit skip size)) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((i offset (+ i skip))) ((eq i limit)) (setf (aref x i) (* a (aref x i))))) x))
(defun x+y (x y &key (x-offset 0) (x-skip 1) (y-offset 0) (y-skip 1) count) (il:* il:|;;;| "Implements BLAS elementary operation, X + Y. ") (il:* il:|;;;| "Count defaults to largest count commensurate with Offsets and skips. ") (il:* il:|;;;| "Returns Y, which is side-effected.") (let ((x-size (length x)) (y-size (length y))) (il:* il:|;;| "Arg checking") (check->=-0 x-offset) (check->=-0 y-offset) (check->-0 x-skip) (check->-0 y-skip) (when (null count) (setq count (min (compute-count x-offset x-skip x-size) (compute-count y-offset y-skip y-size))) (check->=-0 count)) (let ((limit-x (compute-limit x-offset count x-skip)) (limit-y (compute-limit y-offset count y-skip))) (if (not (and (limit-inside-p limit-x x-skip x-size) (limit-inside-p limit-y y-skip y-size))) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((i x-offset (+ i x-skip)) (j y-offset (+ j y-skip))) ((eq i limit-x)) (setf (aref y j) (+ (aref x i) (aref y j))))) y))
(defun ax+y (a x y &key (x-offset 0) (x-skip 1) (y-offset 0) (y-skip 1) count) (il:* il:|;;;| "Implements BLAS elementary operation, a * X + Y. ") (il:* il:|;;;| "Count defaults to largest count commensurate with Offsets and skips. ") (il:* il:|;;;| "Returns Result Y, which is side-effectd.") (let ((x-size (length x)) (y-size (length y))) (il:* il:|;;| "Arg checking") (check->=-0 x-offset) (check->=-0 y-offset) (check->-0 x-skip) (check->-0 y-skip) (when (null count) (setq count (min (compute-count x-offset x-skip x-size) (compute-count y-offset y-skip y-size))) (check->=-0 count)) (let ((limit-x (compute-limit x-offset count x-skip)) (limit-y (compute-limit y-offset count y-skip))) (if (not (and (limit-inside-p limit-x x-skip x-size) (limit-inside-p limit-y y-skip y-size))) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((i x-offset (+ i x-skip)) (j y-offset (+ j y-skip))) ((eq i limit-x)) (setf (aref y j) (+ (* a (aref x i)) (aref y j))))) y))
(defun sum-x (x &key (offset 0) (skip 1) count) (il:* il:|;;;| "Implements BLAS sum reduce X. COUNT defaults to largest vector commensurate with OFFSET and SKIP.") (let ((size (length x))) (il:* il:|;;| "Arg checking") (check->=-0 offset) (check->-0 skip) (when (null count) (setq count (compute-count offset skip size)) (check->=-0 count)) (let ((limit (compute-limit offset count skip))) (if (not (limit-inside-p limit skip size)) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((result 0) (i offset (+ i skip))) ((eq i limit) result) (incf result (aref x i))))))
(defun max-index (x &key (offset 0) (skip 1) count) (il:* il:|;;;| "Implements BLAS index of entry with maximum value. COUNT defaults to largest vector commensurate with OFFSET and SKIP.") (let ((size (length x))) (il:* il:|;;| "Arg checking") (check->=-0 offset) (check->-0 skip) (when (null count) (setq count (compute-count offset skip size))) (check->-0 count) (let ((limit (compute-limit offset count skip))) (if (not (limit-inside-p limit skip size)) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((i offset (+ i skip)) (index 0 (1+ index)) (max-elt (aref x 0)) (max-index 0) elt) ((eq index count) max-index) (setq elt (aref x i)) (when (> elt max-elt) (setq max-elt elt) (setq max-index index))))))
(defun max-abs-index (x &key (offset 0) (skip 1) count) (il:* il:|;;;| "Implements BLAS index of entry with maximum absolute value. COUNT defaults to largest vector commensurate with OFFSET and SKIP.") (let ((size (length x))) (il:* il:|;;| "Arg checking") (check->=-0 offset) (check->-0 skip) (when (null count) (setq count (compute-count offset skip size))) (check->-0 count) (let ((limit (compute-limit offset count skip))) (if (not (limit-inside-p limit skip size)) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((i offset (+ i skip)) (index 0 (1+ index)) (max-elt (aref x 0)) (max-index 0) elt) ((eq index count) max-index) (setq elt (abs (aref x i))) (when (> elt max-elt) (setq max-elt elt) (setq max-index index))))))
(defun <xy> (x y &key (x-offset 0) (x-skip 1) (y-offset 0) (y-skip 1) count) (il:* il:|;;;| "Implements BLAS elementary operation, <X ,Y> (inner product). ") (il:* il:|;;;| "Count defaults to largest count commensurate with Offsets and skips. ") (let ((x-size (length x)) (y-size (length y))) (il:* il:|;;| "Arg checking") (check->=-0 x-offset) (check->=-0 y-offset) (check->-0 x-skip) (check->-0 y-skip) (when (null count) (setq count (min (compute-count x-offset x-skip x-size) (compute-count y-offset y-skip y-size))) (check->=-0 count)) (let ((limit-x (compute-limit x-offset count x-skip)) (limit-y (compute-limit y-offset count y-skip))) (if (not (and (limit-inside-p limit-x x-skip x-size) (limit-inside-p limit-y y-skip y-size))) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((i x-offset (+ i x-skip)) (j y-offset (+ j y-skip)) (sum 0)) ((eq i limit-x) sum) (incf sum (* (aref x i) (aref y j)))))))
(defun l2-norm (x &key (offset 0) (skip 1) count) (il:* il:|;;;| "Implements BLAS vector euclidean norm. COUNT defaults to largest vector commensurate with OFFSET and SKIP.") (il:* il:|;;| "Should be implemented more carefully to avoid overflow on intermediate results.") (sqrt (<xy> x x :x-offset offset :x-skip skip :y-offset offset :y-skip skip :count count)))
(defun compute-givens-coefficients (a b c s) (il:* il:|;;| "Compute the Givens rotation that will eliminate a and b. Returns the coefficients C and S as multiple values. ") (let* ((sigma (if (> (abs a) (abs b)) (if (< a 0) -1.0 1.0) (if (< b 0) -1.0 1.0))) (radius (* sigma (sqrt (* a a) (* b b))))) (if (= radius 0.0) (values 1.0 0.0) (values (/ a radius) (/ b radius)))))
(defun givens-rotate (c s x y &key (x-offset 0) (x-skip 1) (y-offset 0) (y-skip 1) count) (il:* il:|;;;| "Implements BLAS elementary operation, givens rotation defined by c and s. ") (il:* il:|;;;| "Count defaults to largest count commensurate with Offsets and skips. ") (il:* il:|;;;| "Returns Y, but side effects both X and Y.") (let ((x-size (length x)) (y-size (length y))) (il:* il:|;;| "Arg checking") (check->=-0 x-offset) (check->=-0 y-offset) (check->-0 x-skip) (check->-0 y-skip) (when (null count) (setq count (min (compute-count x-offset x-skip x-size) (compute-count y-offset y-skip y-size))) (check->=-0 count)) (let ((limit-x (compute-limit x-offset count x-skip)) (limit-y (compute-limit y-offset count y-skip))) (if (not (and (limit-inside-p limit-x x-skip x-size) (limit-inside-p limit-y y-skip y-size))) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((i x-offset (+ i x-skip)) (j y-offset (+ j y-skip)) x-elt y-elt) ((eq i limit-x)) (setq x-elt (aref x i)) (setq y-elt (aref y j)) (setf (aref x i) (+ (* c x-elt) (* s y-elt))) (setf (aref y j) (- (* c y-elt) (* s x-elt))))) y))
(defun copy-x-to-y (x y &key (x-offset 0) (x-skip 1) (y-offset 0) (y-skip 1) count) (il:* il:|;;;| "Implements BLAS elementary operation, copy x to y. ") (il:* il:|;;;| "Count defaults to largest count commensurate with Offsets and skips. ") (il:* il:|;;;| "Returns Y. Note that result is ill-determined if (eq x y).") (let ((x-size (length x)) (y-size (length y))) (il:* il:|;;| "Arg checking") (check->=-0 x-offset) (check->=-0 y-offset) (check->-0 x-skip) (check->-0 y-skip) (when (null count) (setq count (min (compute-count x-offset x-skip x-size) (compute-count y-offset y-skip y-size))) (check->=-0 count)) (let ((limit-x (compute-limit x-offset count x-skip)) (limit-y (compute-limit y-offset count y-skip))) (if (not (and (limit-inside-p limit-x x-skip x-size) (limit-inside-p limit-y y-skip y-size))) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((i x-offset (+ i x-skip)) (j y-offset (+ j y-skip))) ((eq i limit-x)) (setf (aref y j) (aref x i)))) y))
(defun swap-x-and-y (x y &key (x-offset 0) (x-skip 1) (y-offset 0) (y-skip 1) count) (il:* il:|;;;| "Implements BLAS elementary operation, swap x and y. ") (il:* il:|;;;| "Count defaults to largest count commensurate with Offsets and skips. ") (il:* il:|;;;| "Returns Y, but side effects both X and Y. Note that result is ill-determined if (eq x y).") (let ((x-size (length x)) (y-size (length y))) (il:* il:|;;| "Arg checking") (check->=-0 x-offset) (check->=-0 y-offset) (check->-0 x-skip) (check->-0 y-skip) (when (null count) (setq count (min (compute-count x-offset x-skip x-size) (compute-count y-offset y-skip y-size))) (check->=-0 count)) (let ((limit-x (compute-limit x-offset count x-skip)) (limit-y (compute-limit y-offset count y-skip))) (if (not (and (limit-inside-p limit-x x-skip x-size) (limit-inside-p limit-y y-skip y-size))) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((i x-offset (+ i x-skip)) (j y-offset (+ j y-skip))) ((eq i limit-x)) (rotatef (aref y j) (aref x i)))) y))
(defun fill-x (a x &key (offset 0) (skip 1) count) (il:* il:|;;;| "Implements BLAS fill x with a. COUNT defaults to largest vector commensurate with OFFSET and SKIP. Returns X.") (let ((size (length x))) (il:* il:|;;| "Arg checking") (check->=-0 offset) (check->-0 skip) (when (null count) (setq count (compute-count offset skip size)) (check->=-0 count)) (let ((limit (compute-limit offset count skip))) (if (not (limit-inside-p limit skip size)) (error "Count out of bounds: ~s" count)) (il:* il:|;;| "Compute the result") (do ((i offset (+ i skip))) ((eq i limit)) (setf (aref x i) a))) x))
(eval-when (load)
(export (quote (a+x ax x+y ax+y sum-x max-index max-abs-index <xy> l2-norm compute-givens-coefficients givens-rotate copy-x-to-y swap-x-and-y fill-x)) (find-package "BLAS"))
)
(xcl:define-file-environment "BLAS-FNS" :package (xcl:defpackage "BLAS" (:use "LISP" "PT")) :readtable "XCL" :compiler :compile-file)
(il:putprops il:blas-fns il:copyright ("Xerox Corporation" 1988))
(il:declare\: il:dontcopy
(il:filemap (nil)))
il:stop