matlisp-commit

[Raymond T. <rt...@us...>]

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".

The branch, matlisp-cffi has been updated
       via  45e57ae1f888c8c271a91d42b9231b99bd55691e (commit)
      from  d54f0adc84bb64c23dbbc4bbb1c7885e8cd610e6 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 45e57ae1f888c8c271a91d42b9231b99bd55691e
Author: Raymond Toy <toy...@gm...>
Date:   Mon Mar 26 22:39:54 2012 -0700

    Add demo2 for colnew
    
    src/colnew-demo2.lisp:
    o Second example program for colnew.
    
    src/colnew.lisp:
    o Compute the dimensions of the vectors in the callbacks correctly.
    
    matlisp.asd:
    o Add colnew-demo2.lisp.

diff --git a/matlisp.asd b/matlisp.asd
index b72039e..adf4dd6 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -320,8 +320,7 @@
     :components
     ((:file "colnew")
      (:file "colnew-demo1" :depends-on ("colnew"))
-     #+nil
-     (:file "colnew-demo4" :depends-on ("colnew"))))))
+     (:file "colnew-demo2" :depends-on ("colnew"))))))
 
 (defmethod perform ((op asdf:test-op) (c (eql (asdf:find-system :matlisp))))
   (oos 'asdf:test-op 'matlisp-tests))
diff --git a/src/colnew-demo2.lisp b/src/colnew-demo2.lisp
new file mode 100644
index 0000000..f358b9c
--- /dev/null
+++ b/src/colnew-demo2.lisp
@@ -0,0 +1,145 @@
+;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*-
+
+(in-package #:matlisp)
+
+(defvar *gamma* 1.1d0)
+(defvar *eps* .01d0)
+(defvar *dmu* *eps*)
+(defvar *eps4mu* (/ (expt *eps* 4) *dmu*))
+(defvar *xt* (sqrt (/ (* 2 (- *gamma* 1)) *gamma*)))
+
+
+(defun fsub (x z f)
+  (setf (fv-ref f 0)
+	(+ (/ (fv-ref z 0) x x)
+	   (- (/ (fv-ref z 1) x))
+	   (/ (- (fv-ref z 0)
+		 (* (fv-ref z 2)
+		    (- 1 (/ (fv-ref z 0)
+			    x)))
+		 (* *gamma* x (- 1
+				 (* x x 0.5d0))))
+	      *eps4mu*)))
+  (setf (fv-ref f 1)
+	(+ (/ (fv-ref z 2) x x)
+	   (/ (- (fv-ref z 3)) x)
+	   (* (fv-ref z 0)
+	      (/ (- 1 (/ (fv-ref z 0) 2 x))
+		 *dmu*)))))
+
+(defun dfsub (x z df)
+  (let ((nrows 2))
+    (flet ((column-major-index (r c)
+	     (+ (- r 1)
+		(* (- c 1) nrows))))
+      (setf (fv-ref df (column-major-index 1 1))
+	    (+ (/ 1 x x)
+	       (/ (+ 1
+		     (/ (fv-ref z 2)
+			x))
+		  *eps4mu*)))
+      (setf (fv-ref df (column-major-index 1 2))
+	    (/ -1 x))
+      (setf (fv-ref df (column-major-index 1 3))
+	    (- (/ (- 1 (/ (fv-ref z 0)
+			  x))
+		  *eps4mu*)))
+      (setf (fv-ref df (column-major-index 1 4))
+	    0d0)
+      (setf (fv-ref df (column-major-index 2 1))
+	    (/ (- 1
+		  (/ (fv-ref z 0)
+		     x))
+	       *dmu*))
+      (setf (fv-ref df (column-major-index 2 2))
+	    0d0)
+      (setf (fv-ref df (column-major-index 2 3))
+	    (/ 1 x x))
+      (setf (fv-ref df (column-major-index 2 4))
+	    (/ -1 x)))))
+
+(defun gsub (i z g)
+  (case i
+    ((or 1 3)
+     (setf (fv-ref g 0) (fv-ref z 0)))
+    (2
+     (setf (fv-ref g 0) (fv-ref z 2)))
+    (4
+     (setf (fv-ref g 0) (+ (fv-ref z 3)
+			   (* -0.3d0 (fv-ref z 2))
+			   0.7d0)))))
+
+(defun dgsub (i z dg)
+  (dotimes (k 4)
+    (setf (fv-ref dg k) 0d0))
+  (case i
+    ((or 1 3)
+     (setf (fv-ref dg 0) 1d0))
+    (2
+     (setf (fv-ref dg 2) 1d0))
+    (4
+     (setf (fv-ref dg 3) 1d0)
+     (setf (fv-ref dg 2) -0.3d0))))
+
+(defun guess (x z dmval)
+  (let ((con (* *gamma* x (- 1 (* 0.5d0 x x))))
+	(dcon (* *gamma* (- 1 (* 1.5d0 x x))))
+	(d2con (* -3 *gamma* x)))
+    (cond ((<= x *xt*)
+	   (setf (fv-ref z 0) (* 2 x))
+	   (setf (fv-ref z 1) 2d0)
+	   (setf (fv-ref z 2) (+ (* -2 x) con))
+	   (setf (fv-ref z 3) (+ dcon -2d0))
+	   (setf (fv-ref dmval 1) d2con))
+	  (t
+	   (setf (fv-ref z 0) 0d0)
+	   (setf (fv-ref z 1) 0d0)
+	   (setf (fv-ref z 2) (- con))
+	   (setf (fv-ref z 3) (- dcon))
+	   (setf (fv-ref dmval 1) (- d2con))))
+    (setf (fv-ref dmval 0) 0d0)))
+
+
+(defun colnew-prob2 ()
+  (let* (
+	 (aleft 0d0)
+	 (aright 1d0)
+	 ;; Two differential equations
+	 (ncomp 2)
+	 ;; Orders of each equation
+	 (m (make-array 2 :element-type '(signed-byte 32)
+			  :initial-contents '(2 2)))
+	 ;; Locations of side conditions
+	 (zeta (make-array 4 :element-type 'double-float
+			     :initial-contents '(0d0 0d0 1d0 1d0)))
+	 (ipar (make-array 11 :element-type '(signed-byte 32)
+			      :initial-element 0))
+	 ;; Error tolerances on u and its second derivative
+	 (ltol (make-array 4 :element-type '(signed-byte 32)
+			     :initial-contents '(1 2 3 4)))
+	 (tol (make-array 4 :element-type 'double-float
+			    :initial-contents '(1d-5 1d-5 1d-5 1d-5)))
+	 (fspace (make-array 40000 :element-type 'double-float))
+	 (ispace (make-array 2500 :element-type '(signed-byte 32)))
+	 (fixpnt (make-array 1 :element-type 'double-float)))
+    ;; Set up parameters of the problem.
+    (setf (aref ipar 0) 1)		; nonlinear problem
+    (setf (aref ipar 1) 4)		; 4 collocation points per subinterval
+    (setf (aref ipar 2) 10)		; Initial uniform mesh of 10 subintervals
+    (setf (aref ipar 7) 0)
+    (setf (aref ipar 4) 40000)		; Size of fspace
+    (setf (aref ipar 5) 2500)		; Size of ispace
+    (setf (aref ipar 6) -1)		; Full output
+    (setf (aref ipar 8) 1)		; Initial approx provided
+    (setf (aref ipar 9) 0)		; Regular problem
+    (setf (aref ipar 10) 0)		; No fixed points in mesh
+    (setf (aref ipar 3) 4)		; Tolerances on all components
+    
+    (colnew ncomp m aleft aright zeta ipar ltol tol fixpnt ispace fspace 0
+	    #'fsub #'dfsub #'gsub #'dgsub #'guess)
+    (let ((x 0d0)
+	  (z (make-array 4 :element-type 'double-float)))
+      (dotimes (j 21)
+	(appsln x z fspace ispace)
+	(format t "~5,2f ~{~15,5e~}~%" x (coerce z 'list))
+	(incf x 0.05d0)))))
diff --git a/src/colnew.lisp b/src/colnew.lisp
index d47f76c..adefa60 100644
--- a/src/colnew.lisp
+++ b/src/colnew.lisp
@@ -4,6 +4,10 @@
 
 (cffi:use-foreign-library colnew)
 
+(defun m* (ncomp m)
+  (loop for k from 0 below ncomp
+	sum (aref m k)))
+
 (def-fortran-routine colnew :void
   "COLNEW"
   (ncomp :integer :input)
@@ -20,24 +24,24 @@
   (iflag :integer :output)
   (fsub (:callback :void
 		   (x :double-float :input)
-		   (z (* :double-float :size (aref m 0)) :input)
-		   (f (* :double-float :size (aref m 0)) :output)))
+		   (z (* :double-float :size (m* ncomp m)) :input)
+		   (f (* :double-float :size ncomp) :output)))
   (dfsub (:callback :void
 		    (x :double-float :input)
-		    (z (* :double-float :size (aref m 0)) :input)
-		    (df (* :double-float :size (aref m 0)) :output)))
+		    (z (* :double-float :size (m* ncomp m)) :input)
+		    (df (* :double-float :size (* ncomp (m* ncomp m))) :output)))
   (gsub (:callback :void
 		   (i :integer :input)
-		   (z (* :double-float :size (aref m 0)) :input)
-		   (g (* :double-float :size (aref m 0)) :output)))
+		   (z (* :double-float :size (m* ncomp m)) :input)
+		   (g (* :double-float :size (m* ncomp m)) :output)))
   (dgsub (:callback :void
 		    (i :integer :input)
-		    (z (* :double-float :size (aref m 0)) :input)
-		    (dg (* :double-float :size (aref m 0)) :output)))
+		    (z (* :double-float :size (m* ncomp m)) :input)
+		    (dg (* :double-float :size (expt (m* ncomp m) 2)) :output)))
   (guess (:callback :void
 		    (x :double-float :input)
-		    (z (* :double-float) :output)
-		    (dmval (* :double-float) :output))))
+		    (z (* :double-float :size (m* ncomp m)) :output)
+		    (dmval (* :double-float :size ncomp) :output))))
 
 
 (def-fortran-routine appsln :void

-----------------------------------------------------------------------

Summary of changes:
 matlisp.asd           |    3 +-
 src/colnew-demo2.lisp |  145 +++++++++++++++++++++++++++++++++++++++++++++++++
 src/colnew.lisp       |   24 +++++----
 3 files changed, 160 insertions(+), 12 deletions(-)
 create mode 100644 src/colnew-demo2.lisp


hooks/post-receive
-- 
matlisp