Menu
▾
▴
matlisp-commit
|
From: Raymond T. <rt...@us...> - 2012-03-24 22:08:31
|
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 8191e3e96c966fc09587e44f98a09e42cd5985e9 (commit)
via b53c930a62d5eadcb565e1c77b13a33ac3f24297 (commit)
from 7f20064540e1c4bbb9ba535c37fb1533831cb217 (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 8191e3e96c966fc09587e44f98a09e42cd5985e9
Author: Raymond Toy <toy...@gm...>
Date: Sat Mar 24 15:08:07 2012 -0700
Add support for colnew.
Makefile.am:
o Build and install the colnew library.
configure.ac:
o Add lib-src/colnew/Makefile to list of output files.
lib/lazy-loader.lisp.in:
o Define the colnew library.
matlisp.asd:
o Add defsystem for colnew
src/colnew.lisp:
o New file defining interface to main routines in the colnew package.
src/colnew-demo1.lisp:
o New file giving a simple demo of colnew.
diff --git a/Makefile.am b/Makefile.am
index b073e06..35fb5c5 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -16,6 +16,7 @@ all :
(cd lib-src/toms715; $(MAKE) install)
(cd lib-src/odepack; $(MAKE) install)
(cd lib-src/compat; $(MAKE) install)
+ (cd lib-src/colnew; $(MAKE) install)
if !ATLAS
(cd LAPACK/BLAS/SRC; $(MAKE) install)
(cd LAPACK/SRC; $(MAKE) install)
diff --git a/configure.ac b/configure.ac
index 9bbb440..89f0f1b 100644
--- a/configure.ac
+++ b/configure.ac
@@ -421,13 +421,14 @@ AC_CONFIG_FILES([
start.lisp
config.lisp
lib/lazy-loader.lisp
+ src/f77-mangling.lisp
LAPACK/SRC/Makefile
LAPACK/BLAS/SRC/Makefile
dfftpack/Makefile
lib-src/toms715/Makefile
- lib-src/odepack/Makefile
lib-src/compat/Makefile
- src/f77-mangling.lisp
+ lib-src/odepack/Makefile
+ lib-src/colnew/Makefile
])
echo FLIBS = $FLIBS
diff --git a/lib/lazy-loader.lisp.in b/lib/lazy-loader.lisp.in
index ba1296c..46afb45 100644
--- a/lib/lazy-loader.lisp.in
+++ b/lib/lazy-loader.lisp.in
@@ -137,6 +137,10 @@
(:darwin "libodepack.dylib")
(t (:default "@libdir@/libodepack")))
+(cffi:define-foreign-library colnew
+ (:darwin "libcolnew.dylib")
+ (t (:default "@libdir@/libcolnew")))
+
(cffi:define-foreign-library matlisp
(:darwin "libmatlisp.dylib")
(t (:default "@libdir@/libmatlisp")))
diff --git a/matlisp.asd b/matlisp.asd
index 2c8019c..b72039e 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -311,7 +311,17 @@
:components
((:module "src"
:components
- ((:file "dlsode")))))
+ ((:file "dlsode")))))
+
+(asdf:defsystem matlisp-colnew
+ :pathname #.(translate-logical-pathname "matlisp:srcdir;")
+ :components
+ ((:module "src"
+ :components
+ ((:file "colnew")
+ (:file "colnew-demo1" :depends-on ("colnew"))
+ #+nil
+ (:file "colnew-demo4" :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-demo1.lisp b/src/colnew-demo1.lisp
new file mode 100644
index 0000000..ae5d53c
--- /dev/null
+++ b/src/colnew-demo1.lisp
@@ -0,0 +1,91 @@
+;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*-
+
+(in-package #:matlisp)
+
+(defun fsub (x z f)
+ (setf (fv-ref f 0) (/ (- 1
+ (* 6 x x (fv-ref z 3))
+ (* 6 x (fv-ref z 2)))
+ (* x x x))))
+(defun dfsub (x z df)
+ (setf (fv-ref df 0) 0d0)
+ (setf (fv-ref df 1) 0d0)
+ (setf (fv-ref df 2) (/ -6 x x))
+ (setf (fv-ref df 3) (/ -6 x)))
+
+(defun gsub (i z g)
+ (setf (fv-ref g 0)
+ (if (or (= i 1) (= i 3))
+ (fv-ref z 0)
+ (fv-ref z 2))))
+
+(defun dgsub (i z dg)
+ (dotimes (k 4)
+ (setf (fv-ref dg k) 0d0))
+ (if (or (= i 1) (= i 3))
+ (setf (fv-ref dg 0) 1d0)
+ (setf (fv-ref dg 2) 1d0)))
+
+(defun guess (x z dmval)
+ )
+
+(defun exact (x)
+ (declare (type double-float x))
+ (let ((result (make-array 4 :element-type 'double-float)))
+ (setf (aref result 0)
+ (+ (* .25d0 (- (* 10 (log 2d0)) 3)
+ (- 1 x))
+ (* 0.5d0 (+ (/ x)
+ (* (+ 3 x) (log x))
+ (- x)))))
+ (setf (aref result 1)
+ (+ (* -0.25d0 (- (* 10 (log 2d0)) 3))
+ (* .5d0
+ (+ (/ -1 x x)
+ (log x)
+ (/ (+ 3 x) x)
+ -1))))
+ (setf (aref result 2)
+ (* 0.5d0
+ (+ (/ 2 (expt x 3))
+ (/ x)
+ (/ -3 x x))))
+ (setf (aref result 3)
+ (* 0.5d0
+ (+ (/ -6 (expt x 4))
+ (/ -1 x x)
+ (/ 6 (expt x 3)))))
+ result))
+
+(defun colnew-prob1 ()
+ (let ((m (make-array 1 :element-type '(signed-byte 32)
+ :initial-element 4))
+ (zeta (make-array 4 :element-type 'double-float
+ :initial-contents '(1d0 1d0 2d0 2d0)))
+ (ipar (make-array 11 :element-type '(signed-byte 32)
+ :initial-element 0))
+ (ltol (make-array 2 :element-type '(signed-byte 32)
+ :initial-contents '(1 3)))
+ (tol (make-array 2 :element-type 'double-float
+ :initial-contents '(1d-7 1d-7)))
+ (fspace (make-array 2000 :element-type 'double-float))
+ (ispace (make-array 200 :element-type '(signed-byte 32)))
+ (fixpnt (make-array 1 :element-type 'double-float))
+ (errors (make-array 4 :element-type 'double-float)))
+ (setf (aref ipar 2) 1)
+ (setf (aref ipar 3) 2)
+ (setf (aref ipar 4) 2000)
+ (setf (aref ipar 5) 200)
+ (colnew 1 m 1d0 2d0 zeta ipar ltol tol fixpnt ispace fspace 0
+ #'fsub #'dfsub #'gsub #'dgsub #'guess)
+ (let ((x 1d0)
+ (z (make-array 4 :element-type 'double-float)))
+ (dotimes (j 100)
+ (appsln x z fspace ispace)
+ (map-into errors #'(lambda (a b c)
+ (max a (abs (- b c))))
+ errors
+ (exact x)
+ z))
+ (format t "The exact errors are: ~{ ~11,4e~}~%" (coerce errors 'list)))))
+
\ No newline at end of file
diff --git a/src/colnew.lisp b/src/colnew.lisp
new file mode 100644
index 0000000..d47f76c
--- /dev/null
+++ b/src/colnew.lisp
@@ -0,0 +1,48 @@
+;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*-
+
+(in-package #:matlisp)
+
+(cffi:use-foreign-library colnew)
+
+(def-fortran-routine colnew :void
+ "COLNEW"
+ (ncomp :integer :input)
+ (m (* :integer) :input)
+ (aleft :double-float :input)
+ (aright :double-float :input)
+ (zeta (* :double-float) :input)
+ (ipar (* :integer) :input)
+ (ltol (* :integer) :input)
+ (tol (* :double-float) :input)
+ (fixpnt (* :double-float) :input)
+ (ispace (* :integer) :input-output)
+ (fspace (* :double-float) :input-output)
+ (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)))
+ (dfsub (:callback :void
+ (x :double-float :input)
+ (z (* :double-float :size (aref m 0)) :input)
+ (df (* :double-float :size (aref m 0)) :output)))
+ (gsub (:callback :void
+ (i :integer :input)
+ (z (* :double-float :size (aref m 0)) :input)
+ (g (* :double-float :size (aref m 0)) :output)))
+ (dgsub (:callback :void
+ (i :integer :input)
+ (z (* :double-float :size (aref m 0)) :input)
+ (dg (* :double-float :size (aref m 0)) :output)))
+ (guess (:callback :void
+ (x :double-float :input)
+ (z (* :double-float) :output)
+ (dmval (* :double-float) :output))))
+
+
+(def-fortran-routine appsln :void
+ (x :double-float :input)
+ (z (* :double-float) :output)
+ (fspace (* :double-float) :input)
+ (ispace (* :double-float) :input))
+
commit b53c930a62d5eadcb565e1c77b13a33ac3f24297
Author: Raymond Toy <toy...@gm...>
Date: Sat Mar 24 15:04:56 2012 -0700
Add source files for colnew.
diff --git a/lib-src/colnew/Makefile.am b/lib-src/colnew/Makefile.am
new file mode 100644
index 0000000..9562b0b
--- /dev/null
+++ b/lib-src/colnew/Makefile.am
@@ -0,0 +1,35 @@
+lib_LTLIBRARIES = libcolnew.la
+
+AM_FFLAGS = $(F2C)
+if LIB32
+AM_FFLAGS += -m32
+endif
+
+libcolnew_la_LDFLAGS = $(FLIBS_OPTS)
+libcolnew_la_LIBADD = $(FLIBS)
+
+libcolnew_la_SOURCES = \
+approx.f \
+appsln.f \
+colnew.f \
+consts.f \
+contrl.f \
+dgefa.f \
+dgesl.f \
+dmzsol.f \
+errchk.f \
+factrb.f \
+fcblok.f \
+gblock.f \
+gderiv.f \
+horder.f \
+lsyslv.f \
+newmsh.f \
+rkbas.f \
+sbblok.f \
+shiftb.f \
+skale.f \
+subbak.f \
+subfor.f \
+vmonde.f \
+vwblok.f
\ No newline at end of file
diff --git a/lib-src/colnew/approx.f b/lib-src/colnew/approx.f
new file mode 100644
index 0000000..1f9b93c
--- /dev/null
+++ b/lib-src/colnew/approx.f
@@ -0,0 +1,120 @@
+ SUBROUTINE APPROX (I, X, ZVAL, A, COEF, XI, N, Z, DMZ, K,
+ 1 NCOMP, MMAX, M, MSTAR, MODE, DMVAL, MODM )
+C
+C**********************************************************************
+C
+C purpose
+C (1) (m1-1) (mncomp-1)
+C evaluate z(u(x))=(u (x),u (x),...,u (x),...,u (x) )
+C 1 1 1 mncomp
+C at one point x.
+C
+C variables
+C a - array of mesh independent rk-basis coefficients
+C basm - array of mesh dependent monomial coefficients
+C xi - the current mesh (having n subintervals)
+C z - the current solution vector
+C dmz - the array of mj-th derivatives of the current solution
+C mode - determines the amount of initialization needed
+C = 4 forms z(u(x)) using z, dmz and ha
+C = 3 as in =4, but computes local rk-basis
+C = 2 as in =3, but determines i such that
+C xi(i) .le. x .lt. xi(i+1) (unless x=xi(n+1))
+C = 1 retrieve z=z(u(x(i))) directly
+C
+C**********************************************************************
+C
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION ZVAL(1), DMVAL(1), XI(1), M(1), A(7,1), DM(7)
+ DIMENSION Z(1), DMZ(1), BM(4), COEF(1)
+C
+ COMMON /COLOUT/ PRECIS, IOUT, IPRINT
+C
+ GO TO (10, 30, 80, 90), MODE
+C
+C... mode = 1 , retrieve z( u(x) ) directly for x = xi(i).
+C
+ 10 X = XI(I)
+ IZ = (I-1) * MSTAR
+ DO 20 J = 1, MSTAR
+ IZ = IZ + 1
+ ZVAL(J) = Z(IZ)
+ 20 CONTINUE
+ RETURN
+C
+C... mode = 2 , locate i so xi(i) .le. x .lt. xi(i+1)
+C
+ 30 CONTINUE
+ IF ( X .GE. XI(1)-PRECIS .AND. X .LE. XI(N+1)+PRECIS )
+ 1 GO TO 40
+ IF (IPRINT .LT. 1) WRITE(IOUT,900) X, XI(1), XI(N+1)
+ IF ( X .LT. XI(1) ) X = XI(1)
+ IF ( X .GT. XI(N+1) ) X = XI(N+1)
+ 40 IF ( I .GT. N .OR. I .LT. 1 ) I = (N+1) / 2
+ ILEFT = I
+ IF ( X .LT. XI(ILEFT) ) GO TO 60
+ DO 50 L = ILEFT, N
+ I = L
+ IF ( X .LT. XI(L+1) ) GO TO 80
+ 50 CONTINUE
+ GO TO 80
+ 60 IRIGHT = ILEFT - 1
+ DO 70 L = 1, IRIGHT
+ I = IRIGHT + 1 - L
+ IF ( X .GE. XI(I) ) GO TO 80
+ 70 CONTINUE
+C
+C... mode = 2 or 3 , compute mesh independent rk-basis.
+C
+ 80 CONTINUE
+ S = (X - XI(I)) / (XI(I+1) - XI(I))
+ CALL RKBAS ( S, COEF, K, MMAX, A, DM, MODM )
+C
+C... mode = 2, 3, or 4 , compute mesh dependent rk-basis.
+C
+ 90 CONTINUE
+ BM(1) = X - XI(I)
+ DO 95 L = 2, MMAX
+ BM(L) = BM(1) / DFLOAT(L)
+ 95 CONTINUE
+C
+C... evaluate z( u(x) ).
+C
+ 100 IR = 1
+ IZ = (I-1) * MSTAR + 1
+ IDMZ = (I-1) * K * NCOMP
+ DO 140 JCOMP = 1, NCOMP
+ MJ = M(JCOMP)
+ IR = IR + MJ
+ IZ = IZ + MJ
+ DO 130 L = 1, MJ
+ IND = IDMZ + JCOMP
+ ZSUM = 0.D0
+ DO 110 J = 1, K
+ ZSUM = ZSUM + A(J,L) * DMZ(IND)
+ 110 IND = IND + NCOMP
+ DO 120 LL = 1, L
+ LB = L + 1 - LL
+ 120 ZSUM = ZSUM * BM(LB) + Z(IZ-LL)
+ 130 ZVAL(IR-L) = ZSUM
+ 140 CONTINUE
+ IF ( MODM .EQ. 0 ) RETURN
+C
+C... for modm = 1 evaluate dmval(j) = mj-th derivative of uj.
+C
+ DO 150 JCOMP = 1, NCOMP
+ 150 DMVAL(JCOMP) = 0.D0
+ IDMZ = IDMZ + 1
+ DO 170 J = 1, K
+ FACT = DM(J)
+ DO 160 JCOMP = 1, NCOMP
+ DMVAL(JCOMP) = DMVAL(JCOMP) + FACT * DMZ(IDMZ)
+ IDMZ = IDMZ + 1
+ 160 CONTINUE
+ 170 CONTINUE
+ RETURN
+C--------------------------------------------------------------------
+ 900 FORMAT(37H ****** DOMAIN ERROR IN APPROX ******
+ 1 /4H X =,D20.10, 10H ALEFT =,D20.10,
+ 2 11H ARIGHT =,D20.10)
+ END
diff --git a/lib-src/colnew/appsln.f b/lib-src/colnew/appsln.f
new file mode 100644
index 0000000..7edb726
--- /dev/null
+++ b/lib-src/colnew/appsln.f
@@ -0,0 +1,31 @@
+C
+C----------------------------------------------------------------------
+C p a r t 4
+C polynomial and service routines
+C----------------------------------------------------------------------
+C
+ SUBROUTINE APPSLN (X, Z, FSPACE, ISPACE)
+C
+C*****************************************************************
+C
+C purpose
+C
+C set up a standard call to approx to evaluate the
+C approximate solution z = z( u(x) ) at a point x
+C (it has been computed by a call to colnew ).
+C the parameters needed for approx are retrieved
+C from the work arrays ispace and fspace .
+C
+C*****************************************************************
+C
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION Z(1), FSPACE(1), ISPACE(1), A(28), DUMMY(1)
+ IS6 = ISPACE(6)
+ IS5 = ISPACE(1) + 2
+ IS4 = IS5 + ISPACE(4) * (ISPACE(1) + 1)
+ I = 1
+ CALL APPROX (I, X, Z, A, FSPACE(IS6), FSPACE(1), ISPACE(1),
+ 1 FSPACE(IS5), FSPACE(IS4), ISPACE(2), ISPACE(3),
+ 2 ISPACE(5), ISPACE(8), ISPACE(4), 2, DUMMY, 0)
+ RETURN
+ END
diff --git a/lib-src/colnew/colnew.f b/lib-src/colnew/colnew.f
new file mode 100644
index 0000000..06aa822
--- /dev/null
+++ b/lib-src/colnew/colnew.f
@@ -0,0 +1,740 @@
+c From research!csnet!CSNET-RELAY!mit-multics.arpa!UBC.mailnet!USER=NBAF Tue, 3 Feb 87 15:25:36 PST
+C**********************************************************************
+C this package solves boundary value problems for
+C ordinary differential equations, as described below.
+C
+C COLNEW is a modification of the package COLSYS by ascher,
+C christiansen and russell [1]. It incorporates a new basis
+C representation replacing b-splines, and improvements for
+C the linear and nonlinear algebraic equation solvers.
+C the package can be referenced as either COLNEW or COLSYS.
+C**********************************************************************
+C----------------------------------------------------------------------
+C p a r t 1
+C main storage allocation and program control subroutines
+C----------------------------------------------------------------------
+C
+ SUBROUTINE COLNEW (NCOMP, M, ALEFT, ARIGHT, ZETA, IPAR, LTOL,
+ 1 TOL, FIXPNT, ISPACE, FSPACE, IFLAG,
+ 2 FSUB, DFSUB, GSUB, DGSUB, GUESS)
+C
+C
+C**********************************************************************
+C
+C written by
+C u. ascher,
+C department of computer science,
+C university of british columbia,
+C vancouver, b. c., canada v6t 1w5
+C g. bader,
+C institut f. angewandte mathematik
+C university of heidelberg
+C im neuenheimer feld 294
+C d-6900 heidelberg 1
+C
+C**********************************************************************
+C
+C purpose
+C
+C this package solves a multi-point boundary value
+C problem for a mixed order system of ode-s given by
+C
+C (m(i))
+C u = f ( x; z(u(x)) ) i = 1, ... ,ncomp
+C i i
+C
+C aleft .lt. x .lt. aright,
+C
+C
+C g ( zeta(j); z(u(zeta(j))) ) = 0 j = 1, ... ,mstar
+C j
+C mstar = m(1)+m(2)+...+m(ncomp),
+C
+C
+C where t
+C u = (u , u , ... ,u ) is the exact solution vector
+C 1 2 ncomp
+C
+C (mi)
+C u is the mi=m(i) th derivative of u
+C i i
+C
+C (1) (m1-1) (mncomp-1)
+C z(u(x)) = ( u (x),u (x),...,u (x),...,u (x) )
+C 1 1 1 ncomp
+C
+C f (x,z(u)) is a (generally) nonlinear function of
+C i
+C z(u)=z(u(x)).
+C
+C g (zeta(j);z(u)) is a (generally) nonlinear function
+C j
+C used to represent a boundary condition.
+C
+C the boundary points satisfy
+C aleft .le. zeta(1) .le. .. .le. zeta(mstar) .le. aright
+C
+C the orders mi of the differential equations satisfy
+C 1 .le. m(i) .le. 4.
+C
+C
+C**********************************************************************
+C
+C method
+C
+C the method used to approximate the solution u is
+C collocation at gaussian points, requiring m(i)-1 continuous
+C derivatives in the i-th component, i = 1, ..., ncomp.
+C here, k is the number of collocation points (stages) per
+C subinterval and is chosen such that k .ge. max m(i).
+C a runge-kutta-monomial solution representation is utilized.
+C
+C references
+C
+C [1] u. ascher, j. christiansen and r.d. russell,
+C collocation software for boundary-value odes,
+C acm trans. math software 7 (1981), 209-222.
+C this paper contains EXAMPLES where use of the code
+C is demonstrated.
+C
+C [2] g. bader and u. ascher,
+C a new basis implementation for a mixed order
+C boundary value ode solver,
+C siam j. scient. stat. comput. (1987).
+C
+C [3] u. ascher, j. christiansen and r.d. russell,
+C a collocation solver for mixed order
+C systems of boundary value problems,
+C math. comp. 33 (1979), 659-679.
+C
+C [4] u. ascher, j. christiansen and r.d. russell,
+C colsys - a collocation code for boundary
+C value problems,
+C lecture notes comp.sc. 76, springer verlag,
+C b. childs et. al. (eds.) (1979), 164-185.
+C
+C [5] c. deboor and r. weiss,
+C solveblok: a package for solving almost block diagonal
+C linear systems,
+C acm trans. math. software 6 (1980), 80-87.
+C
+C**********************************************************************
+C
+C *************** input to colnew ***************
+C
+C variables
+C
+C ncomp - no. of differential equations (ncomp .le. 20)
+C
+C m(j) - order of the j-th differential equation
+C ( mstar = m(1) + ... + m(ncomp) .le. 40 )
+C
+C aleft - left end of interval
+C
+C aright - right end of interval
+C
+C zeta(j) - j-th side condition point (boundary point). must
+C have zeta(j) .le. zeta(j+1). all side condition
+C points must be mesh points in all meshes used,
+C see description of ipar(11) and fixpnt below.
+C
+C ipar - an integer array dimensioned at least 11.
+C a list of the parameters in ipar and their meaning follows
+C some parameters are renamed in colnew; these new names are
+C given in parentheses.
+C
+C ipar(1) ( = nonlin )
+C = 0 if the problem is linear
+C = 1 if the problem is nonlinear
+C
+C ipar(2) = no. of collocation points per subinterval (= k )
+C where max m(i) .le. k .le. 7 . if ipar(2)=0 then
+C colnew sets k = max ( max m(i)+1, 5-max m(i) )
+C
+C ipar(3) = no. of subintervals in the initial mesh ( = n ).
+C if ipar(3) = 0 then colnew arbitrarily sets n = 5.
+C
+C ipar(4) = no. of solution and derivative tolerances. ( = ntol )
+C we require 0 .lt. ntol .le. mstar.
+C
+C ipar(5) = dimension of fspace. ( = ndimf )
+C
+C ipar(6) = dimension of ispace. ( = ndimi )
+C
+C ipar(7) - output control ( = iprint )
+C = -1 for full diagnostic printout
+C = 0 for selected printout
+C = 1 for no printout
+C
+C ipar(8) ( = iread )
+C = 0 causes colnew to generate a uniform initial mesh.
+C = 1 if the initial mesh is provided by the user. it
+C is defined in fspace as follows: the mesh
+C aleft=x(1).lt.x(2).lt. ... .lt.x(n).lt.x(n+1)=aright
+C will occupy fspace(1), ..., fspace(n+1). the
+C user needs to supply only the interior mesh
+C points fspace(j) = x(j), j = 2, ..., n.
+C = 2 if the initial mesh is supplied by the user
+C as with ipar(8)=1, and in addition no adaptive
+C mesh selection is to be done.
+C
+C ipar(9) ( = iguess )
+C = 0 if no initial guess for the solution is
+C provided.
+C = 1 if an initial guess is provided by the user
+C in subroutine guess.
+C = 2 if an initial mesh and approximate solution
+C coefficients are provided by the user in fspace.
+C (the former and new mesh are the same).
+C = 3 if a former mesh and approximate solution
+C coefficients are provided by the user in fspace,
+C and the new mesh is to be taken twice as coarse;
+C i.e.,every second point from the former mesh.
+C = 4 if in addition to a former initial mesh and
+C approximate solution coefficients, a new mesh
+C is provided in fspace as well.
+C (see description of output for further details
+C on iguess = 2, 3, and 4.)
+C
+C ipar(10)= 0 if the problem is regular
+C = 1 if the first relax factor is =rstart, and the
+C nonlinear iteration does not rely on past covergence
+C (use for an extra sensitive nonlinear problem only).
+C = 2 if we are to return immediately upon (a) two
+C successive nonconvergences, or (b) after obtaining
+C error estimate for the first time.
+C
+C ipar(11)= no. of fixed points in the mesh other than aleft
+C and aright. ( = nfxpnt , the dimension of fixpnt)
+C the code requires that all side condition points
+C other than aleft and aright (see description of
+C zeta ) be included as fixed points in fixpnt.
+C
+C ltol - an array of dimension ipar(4). ltol(j) = l specifies
+C that the j-th tolerance in tol controls the error
+C in the l-th component of z(u). also require that
+C 1.le.ltol(1).lt.ltol(2).lt. ... .lt.ltol(ntol).le.mstar
+C
+C tol - an array of dimension ipar(4). tol(j) is the
+C error tolerance on the ltol(j) -th component
+C of z(u). thus, the code attempts to satisfy
+C for j=1,...,ntol on each subinterval
+C abs(z(v)-z(u)) .le. tol(j)*abs(z(u)) +tol(j)
+C ltol(j) ltol(j)
+C
+C if v(x) is the approximate solution vector.
+C
+C fixpnt - an array of dimension ipar(11). it contains
+C the points, other than aleft and aright, which
+C are to be included in every mesh.
+C
+C ispace - an integer work array of dimension ipar(6).
+C its size provides a constraint on nmax,
+C the maximum number of subintervals. choose
+C ipar(6) according to the formula
+C ipar(6) .ge. nmax*nsizei
+C where
+C nsizei = 3 + kdm
+C with
+C kdm = kd + mstar ; kd = k * ncomp ;
+C nrec = no. of right end boundary conditions.
+C
+C
+C fspace - a real work array of dimension ipar(5).
+C its size provides a constraint on nmax.
+C choose ipar(5) according to the formula
+C ipar(5) .ge. nmax*nsizef
+C where
+C nsizef = 4 + 3 * mstar + (5+kd) * kdm +
+C (2*mstar-nrec) * 2*mstar.
+C
+C
+C iflag - the mode of return from colnew.
+C = 1 for normal return
+C = 0 if the collocation matrix is singular.
+C =-1 if the expected no. of subintervals exceeds storage
+C specifications.
+C =-2 if the nonlinear iteration has not converged.
+C =-3 if there is an input data error.
+C
+C
+C**********************************************************************
+C
+C ************* user supplied subroutines *************
+C
+C
+C the following subroutines must be declared external in the
+C main program which calls colnew.
+C
+C
+C fsub - name of subroutine for evaluating f(x,z(u(x))) =
+C t
+C (f ,...,f ) at a point x in (aleft,aright). it
+C 1 ncomp
+C should have the heading
+C
+C subroutine fsub (x , z , f)
+C
+C where f is the vector containing the value of fi(x,z(u))
+C in the i-th component and t
+C z(u(x))=(z(1),...,z(mstar))
+C is defined as above under purpose .
+C
+C
+C dfsub - name of subroutine for evaluating the jacobian of
+C f(x,z(u)) at a point x. it should have the heading
+C
+C subroutine dfsub (x , z , df)
+C
+C where z(u(x)) is defined as for fsub and the (ncomp) by
+C (mstar) array df should be filled by the partial deriv-
+C atives of f, viz, for a particular call one calculates
+C df(i,j) = dfi / dzj, i=1,...,ncomp
+C j=1,...,mstar.
+C
+C
+C gsub - name of subroutine for evaluating the i-th component of
+C g(x,z(u(x))) = g (zeta(i),z(u(zeta(i)))) at a point x =
+C i
+C zeta(i) where 1.le.i.le.mstar. it should have the heading
+C
+C subroutine gsub (i , z , g)
+C
+C where z(u) is as for fsub, and i and g=g are as above.
+C i
+C note that in contrast to f in fsub , here
+C only one value per call is returned in g.
+C
+C
+C dgsub - name of subroutine for evaluating the i-th row of
+C the jacobian of g(x,u(x)). it should have the heading
+C
+C subroutine dgsub (i , z , dg)
+C
+C where z(u) is as for fsub, i as for gsub and the mstar-
+C vector dg should be filled with the partial derivatives
+C of g, viz, for a particular call one calculates
+C dg(i,j) = dgi / dzj j=1,...,mstar.
+C
+C
+C guess - name of subroutine to evaluate the initial
+C approximation for z(u(x)) and for dmval(u(x))= vector
+C of the mj-th derivatives of u(x). it should have the
+C heading
+C
+C subroutine guess (x , z , dmval)
+C
+C note that this subroutine is needed only if using
+C ipar(9) = 1, and then all mstar components of z
+C and ncomp components of dmval should be specified
+C for any x, aleft .le. x .le. aright .
+C
+C
+C**********************************************************************
+C
+C ************ use of output from colnew ************
+C
+C *** solution evaluation ***
+C
+C on return from colnew, the arrays fspace and ispace
+C contain information specifying the approximate solution.
+C the user can produce the solution vector z( u(x) ) at
+C any point x, aleft .le. x .le. aright, by the statement,
+C
+C call appsln (x, z, fspace, ispace)
+C
+C when saving the coefficients for later reference, only
+C ispace(1),...,ispace(7+ncomp) and
+C fspace(1),...,fspace(ispace(7)) need to be saved as
+C these are the quantities used by appsln.
+C
+C
+C *** simple continuation ***
+C
+C
+C a formerly obtained solution can easily be used as the
+C first approximation for the nonlinear iteration for a
+C new problem by setting (iguess =) ipar(9) = 2, 3 or 4.
+C
+C if the former solution has just been obtained then the
+C values needed to define the first approximation are
+C already in ispace and fspace.
+C alternatively, if the former solution was obtained in a
+C previous run and its coefficients were saved then those
+C coefficients must be put back into
+C ispace(1),..., ispace(7+ncomp) and
+C fspace(1),..., fspace(ispace(7)).
+C
+C for ipar(9) = 2 or 3 set ipar(3) = ispace(1) ( = the
+C size of the previous mesh ).
+C
+C for ipar(9) = 4 the user specifies a new mesh of n subintervals
+C as follows.
+C the values in fspace(1),...,fspace(ispace(7)) have to be
+C shifted by n+1 locations to fspace(n+2),..,fspace(ispace(7)+n+1)
+C and the new mesh is then specified in fspace(1),..., fspace(n+1).
+C also set ipar(3) = n.
+C
+C
+C**********************************************************************
+C
+C *************** package subroutines ***************
+C
+C the...
[truncated message content] |
Thanks for helping keep SourceForge clean.
X