From mboxrd@z Thu Jan 1 00:00:00 1970 Path: news.gmane.org!.POSTED.blaine.gmane.org!not-for-mail From: Matt Wette Newsgroups: gmane.lisp.guile.user Subject: Re: Matrix or array operations library Date: Sun, 27 Jan 2019 10:41:57 -0800 Message-ID: References: <1C6C55C9-24B1-4952-BC1A-18CF80749565@bluewin.ch> <2D1247B3-4B56-45C3-9AC4-5EA5E1B0858B@bluewin.ch> Mime-Version: 1.0 Content-Type: text/plain; charset=utf-8; format=flowed Content-Transfer-Encoding: 7bit Injection-Info: blaine.gmane.org; posting-host="blaine.gmane.org:195.159.176.226"; logging-data="179526"; mail-complaints-to="usenet@blaine.gmane.org" User-Agent: Mozilla/5.0 (X11; Linux x86_64; rv:60.0) Gecko/20100101 Thunderbird/60.4.0 To: guile-user@gnu.org Original-X-From: guile-user-bounces+guile-user=m.gmane.org@gnu.org Sun Jan 27 19:50:30 2019 Return-path: Envelope-to: guile-user@m.gmane.org Original-Received: from lists.gnu.org ([209.51.188.17]) by blaine.gmane.org with esmtps (TLS1.0:RSA_AES_256_CBC_SHA1:256) (Exim 4.89) (envelope-from ) id 1gnpVl-000kb6-NT for guile-user@m.gmane.org; Sun, 27 Jan 2019 19:50:29 +0100 Original-Received: from localhost ([127.0.0.1]:49197 helo=lists.gnu.org) by lists.gnu.org with esmtp (Exim 4.71) (envelope-from ) id 1gnpVk-0001hu-MW for guile-user@m.gmane.org; Sun, 27 Jan 2019 13:50:28 -0500 Original-Received: from eggs.gnu.org ([209.51.188.92]:34534) by lists.gnu.org with esmtp (Exim 4.71) (envelope-from ) id 1gnpSQ-0007nq-Mt for guile-user@gnu.org; Sun, 27 Jan 2019 13:47:03 -0500 Original-Received: from Debian-exim by eggs.gnu.org with spam-scanned (Exim 4.71) (envelope-from ) id 1gnpNZ-0003c1-M3 for guile-user@gnu.org; Sun, 27 Jan 2019 13:42:02 -0500 Original-Received: from mail-pl1-x630.google.com ([2607:f8b0:4864:20::630]:44029) by eggs.gnu.org with esmtps (TLS1.0:RSA_AES_128_CBC_SHA1:16) (Exim 4.71) (envelope-from ) id 1gnpNZ-0003Zs-Dp for guile-user@gnu.org; Sun, 27 Jan 2019 13:42:01 -0500 Original-Received: by mail-pl1-x630.google.com with SMTP id gn14so6730425plb.10 for ; Sun, 27 Jan 2019 10:42:00 -0800 (PST) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20161025; h=subject:to:references:from:message-id:date:user-agent:mime-version :in-reply-to:content-transfer-encoding:content-language; bh=oSnXVw+q9FhY1uZzJqOJN0BPMcKFV36kNmDAIpt0IF8=; b=AGbe6Sc6sDxVYo2Oxv0g5oJr7svRxRVj7s+Gi6MEu7xPbEAQCQzNg5iM5QIje5o6hj 0xEInocyjUXkhEECva9gibhrI0/IVuVQND7Ij9FUxgsoeUCnU+3xn0JWZeq3swNgIPbf Oyhzi6Rlq73zW07EMz6uw7DgY+hxUzfX6SYvevOO2l8hcrp00vTxXf7cRDD4WOwr35DR MysZpokqrKocSJsd0YaL0x9Tc5XDpFSo+BNhX8RzC5z/S68szfH+du28t0eQXB36emFX ujnngNMEuoy2WmxeaKvFiB9cIJfYdDBnTrTbdI2tlnEK/UB0c0LcTyaKvB3WENWJtsuR 77/w== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20161025; h=x-gm-message-state:subject:to:references:from:message-id:date :user-agent:mime-version:in-reply-to:content-transfer-encoding :content-language; bh=oSnXVw+q9FhY1uZzJqOJN0BPMcKFV36kNmDAIpt0IF8=; b=WHF+nhN+kL6X7p09bIF9MVbmzZ2RInvZeWI/Kq9nsG6ZGCjEwHZP4lrX0LW+FO5AeG 7l5gGpBGyhGlN/+9z7vNgru+HRxhQspwVg0E7D1FOV/E6Pj1j2pSH2eHvFRg5ndWKhT9 BZCu3MK4W7wpKvrJs3Ht/AhcvfIGQaozVMp1aqjGzFPvwp8gTyLoIYYugQ4nRAVuZOxQ XfTR6WmbLWi9EGjDXYB854J7RZF9w/waHsugvwiRIRllXPcPxTK7gzmcjpT0WFo8YIqg Xi/EMQYctJON2B323R+rYz4uTQhrXmrovmldf6FEDFMbTX4x5wD5+Wy7QrOT5FsQPDRM 7L+g== X-Gm-Message-State: AJcUuke70tZnlVvBi9RJ/7etfUhnndnyBVkj+qAFRM2ZaXKcjk2n5p6x jvq4e5beKlMn2SrCJBoSgaMbvnik X-Google-Smtp-Source: ALg8bN64ej5ljmuftp5bXVSTkM/N9zV8vNvUXrSnUjWKRTxnMEfV4hWn9bhr3P0m5w1wKaEpWhPlSA== X-Received: by 2002:a17:902:1102:: with SMTP id d2mr19030052pla.138.1548614519319; Sun, 27 Jan 2019 10:41:59 -0800 (PST) Original-Received: from [192.168.2.183] (216-165-228-38.championbroadband.com. [216.165.228.38]) by smtp.gmail.com with ESMTPSA id q187sm43828869pfq.128.2019.01.27.10.41.58 for (version=TLS1_2 cipher=ECDHE-RSA-AES128-GCM-SHA256 bits=128/128); Sun, 27 Jan 2019 10:41:58 -0800 (PST) In-Reply-To: <2D1247B3-4B56-45C3-9AC4-5EA5E1B0858B@bluewin.ch> Content-Language: en-US X-detected-operating-system: by eggs.gnu.org: Genre and OS details not recognized. X-Received-From: 2607:f8b0:4864:20::630 X-BeenThere: guile-user@gnu.org X-Mailman-Version: 2.1.21 Precedence: list List-Id: General Guile related discussions List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Errors-To: guile-user-bounces+guile-user=m.gmane.org@gnu.org Original-Sender: "guile-user" Xref: news.gmane.org gmane.lisp.guile.user:15265 Archived-At: So I tried an example using Guile to perform Gaussian elimination (in situ, no column pivoting so unstable) using shared-arrays and using explicit computation of offsets. The explicit approach is much more grungy but I'm guessing it does a lot less memory allocation. Also, is it possible to get a pointer to the first element of an array (i.e., it's #f64() root vector), for the purpose of using with FFI? (use-modules ((srfi srfi-1) #:select (fold))) (use-modules (srfi srfi-11)) ;; === method 1: use share arrays (define (axpy1 a x y) (let ((n (car (array-dimensions x)))) (if (not (= n (car (array-dimensions y)))) (throw 'error)) (let loop ((i 0)) (when (< i n) (array-set! y (+ (* a (array-ref x i)) (array-ref y i)) i) (loop (1+ i)))))) ;; modify A, b in-situ so that A is upper triangular (define (reduce1 A b n) (let ((n-1 (1- n))) (let loop ((kl (iota n))) (when (pair? kl) (let* ((k (car kl)) (alpha (- (/ 1.0 (array-ref A k k)))) (x (make-shared-array A (lambda (j) (list k j)) n))) (let loop1 ((i (1+ k))) (when (< i n) (let ((a (* alpha (array-ref A i k))) (y (make-shared-array A (lambda (j) (list i j)) n))) (axpy1 a x y) (array-set! b (+ (* a (array-ref b k)) (array-ref b i)) i)) (loop1 (1+ i))))) (loop (cdr kl)))))) ;; solve A*x = b for x where A is upper triangular (define (bksolv1 A b x n) (let loop ((k (1- n))) (unless (negative? k) (let loop1 ((r 0.0) (j (1- n))) (cond ((< k j) (loop1 (+ r (* (array-ref A k j) (array-ref x j))) (1- j))) (else (array-set! x (/ (- (array-ref b k) r) (array-ref A k k)) k)))) (loop (1- k))))) ;; === method 2: use root vector and computing offsets explicitly ;; based on root vectors ;; x, y are root vectors (from shared-array-root) ;; x0, y0 are starting indices ;; xi, yi are increments (define (axpy2 n a x x0 xi y y0 yi) (let loop ((ix x0) (iy y0) (i 0)) (when (< i n) (array-set! y (+ (* a (array-ref x ix)) (array-ref y iy)) iy) (loop (+ ix xi) (+ iy yi) (1+ i))))) ;; @deffn {Procedure} array-vec-model a vec-idx . idxs ;; Generate params for accessing a vector in an array. ;; Returned expression is ;; (values rv off inc) ;; where @code{rv} is the root-vector, @code{off} is the offset of the ;; first element and @code{inc} is the increment. ;; @end deffn (define (array-vec-model a vec-idx . idxs) (let ((ic (shared-array-increments a))) (values (shared-array-root a) (fold (lambda (bd ic ix of) (+ of (* ic (- ix (car bd))))) (shared-array-offset a) (array-shape a) ic idxs) (list-ref ic vec-idx)))) ;; @deffn {Procedure} array-mat-model a row-idx col-idx . idxs ;; Generate params for accessing a vector in an array. ;; Returned expression is ;; (values rv off inc) ;; where @code{rv} is the root-vector, @code{off} is the offset of the ;; first element and @code{inc} is the increment. ;; @end deffn (define (array-mat-model a row-idx col-idx . idxs) (let ((ic (shared-array-increments a))) (values (shared-array-root a) (fold (lambda (bd ic ix of) (+ of (* ic (- ix (car bd))))) (shared-array-offset a) (array-shape a) ic idxs) (list-ref ic row-idx) (list-ref ic col-idx)))) ;; modify A, b in-situ so that A is upper triangular (define (reduce2 A b n) (let-values (((av a0 arinc acinc) (array-mat-model A 0 1 0 0)) ((bv b0 binc) (array-vec-model b 0 0))) (let loop ((kl (iota n)) (akk a0) (bk b0)) (when (pair? kl) (let* ((k (car kl)) (alpha (- (/ 1.0 (array-ref av akk))))) (let loop1 ((i (1+ k)) (aik (+ akk arinc)) (bi (+ bk binc))) (when (< i n) (let ((a (* alpha (array-ref av aik)))) (axpy2 (- n k) a av akk acinc av aik acinc) (array-set! bv (+ (* a (array-ref bv bk)) (array-ref bv bi)) bi)) (loop1 (1+ i) (+ aik arinc) (+ bi binc))))) (loop (cdr kl) (+ akk arinc acinc) (+ bk binc)))))) ;; solve A*x = b for x where A is upper triangular (define (bksolv2 A b x n) (let*-values (((av a00 arinc acinc) (array-mat-model A 0 1 0 0)) ((bv b0 binc) (array-vec-model b 0 0)) ((xv x0 xinc) (array-vec-model x 0 0)) ((bn) (+ b0 (* (1- n) binc)))) (let loop ((k (1- n)) (akn (+ a00 (* (1- n) arinc) (* (1- n) acinc))) (xk (+ x0 (* (1- n) xinc)))) (unless (negative? k) (let loop1 ((r 0.0) (j (1- n)) (akj akn) (bj bn)) (cond ((< k j) (loop1 (+ r (* (array-ref av akj) (array-ref x j))) (1- j) (- akj acinc) (- bj binc))) (else (array-set! xv (/ (- (array-ref b j) r) (array-ref av akj)) xk)))) (loop (1- k) (- akn arinc) (- xk xinc)))))) ;; === example problem (define A (list->typed-array 'f64 2 '((11.0 1.2 1.3 1.4) (2.1 22.0 2.3 2.4) (3.1 3.2 33.0 3.4) (4.1 4.2 4.3 44.0)))) (define b (list->typed-array 'f64 1 '(4.0 3.0 2.0 1.0))) (define x (make-typed-array 'f64 0.0 4)) ;; expected result (define x-soln (list->typed-array 'f64 1 '(0.352856 0.103010 0.019728 -0.021913))) (display "expect:\n") (vec-disp x-soln) (newline) (when #f (reduce1 A b 4) ;;(display "A:\n") (mat-disp A) ;;(display "b:\n") (vec-disp b) (bksolv1 A b x 4) (display "method1:\n") (vec-disp x) ) (when #t (reduce2 A b 4) ;;(display "A:\n") (mat-disp A) ;;(display "b:\n") (vec-disp b) (bksolv2 A b x 4) (display "method2:\n") (vec-disp x) )