Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch matlisp-cffi updated. 2012-02-23-74-g8191e3e
|
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 following description gives a brief overview of how the
+C procedure is broken down into the subroutines which make up
+C the package called colnew . for further details the
+C user should refer to documentation in the various subroutines
+C and to the references cited above.
+C
+C the subroutines fall into four groups:
+C
+C part 1 - the main storage allocation and program control subr
+C
+C colnew - tests input values, does initialization and breaks up
+C the work areas, fspace and ispace, into the arrays
+C used by the program.
+C colsys - another name for colnew
+C
+C contrl - is the actual driver of the package. this routine
+C contains the strategy for nonlinear equation solving.
+C
+C skale - provides scaling for the control
+C of convergence in the nonlinear iteration.
+C
+C
+C part 2 - mesh selection and error estimation subroutines
+C
+C consts - is called once by colnew to initialize constants
+C which are used for error estimation and mesh selection.
+C
+C newmsh - generates meshes. it contains the test to decide
+C whether or not to redistribute a mesh.
+C
+C errchk - produces error estimates and checks against the
+C tolerances at each subinterval
+C
+C
+C part 3 - collocation system set-up subroutines
+C
+C lsyslv - controls the set-up and solution of the linear
+C algebraic systems of collocation equations which
+C arise at each newton iteration.
+C
+C gderiv - is used by lsyslv to set up the equation associated
+C with a side condition point.
+C
+C vwblok - is used by lsyslv to set up the equation(s) associated
+C with a collocation point.
+C
+C gblock - is used by lsyslv to construct a block of the global
+C collocation matrix or the corresponding right hand
+C side.
+C
+C
+C part 4 - service subroutines
+C
+C appsln - sets up a standard call to approx .
+C
+C approx - evaluates a piecewise polynomial solution.
+C
+C rkbas - evaluates the mesh independent runge-kutta basis
+C
+C vmonde - solves a vandermonde system for given right hand
+C side
+C
+C horder - evaluates the highest order derivatives of the
+C current collocation solution used for mesh refinement.
+C
+C
+C part 5 - linear algebra subroutines
+C
+C to solve the global linear systems of collocation equations
+C constructed in part 3, colnew uses a column oriented version
+C of the package solveblok originally due to de boor and weiss.
+C
+C to solve the linear systems for static parameter condensation
+C in each block of the collocation equations, the linpack
+C routines dgefa and dgesl are included. but these
+C may be replaced when solving problems on vector processors
+C or when solving large scale sparse jacobian problems.
+C
+C----------------------------------------------------------------------
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION M(1), ZETA(1), IPAR(1), LTOL(1), TOL(1), DUMMY(1),
+ 1 FIXPNT(1), ISPACE(1), FSPACE(1)
+C
+ COMMON /COLOUT/ PRECIS, IOUT, IPRINT
+ COMMON /COLLOC/ RHO(7), COEF(49)
+ COMMON /COLORD/ K, NC, MSTAR, KD, MMAX, MT(20)
+ COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ
+ COMMON /COLMSH/ MSHFLG, MSHNUM, MSHLMT, MSHALT
+ COMMON /COLSID/ TZETA(40), TLEFT, TRIGHT, IZETA, IDUM
+ COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS
+ COMMON /COLEST/ TTL(40), WGTMSH(40), WGTERR(40), TOLIN(40),
+ 1 ROOT(40), JTOL(40), LTTOL(40), NTOL
+C
+ EXTERNAL FSUB, DFSUB, GSUB, DGSUB, GUESS
+C
+C this subroutine can be called either COLNEW or COLSYS
+C
+ ENTRY COLSYS (NCOMP, M, ALEFT, ARIGHT, ZETA, IPAR, LTOL,
+ 1 TOL, FIXPNT, ISPACE, FSPACE, IFLAG,
+ 2 FSUB, DFSUB, GSUB, DGSUB, GUESS)
+C
+C*********************************************************************
+C
+C the actual subroutine colnew serves as an interface with
+C the package of subroutines referred to collectively as
+C colnew. the subroutine serves to test some of the input
+C parameters, rename some of the parameters (to make under-
+C standing of the coding easier), to do some initialization,
+C and to break the work areas fspace and ispace up into the
+C arrays needed by the program.
+C
+C**********************************************************************
+C
+C... specify machine dependent output unit iout and compute machine
+C... dependent constant precis = 100 * machine unit roundoff
+C
+ IF ( IPAR(7) .LE. 0 ) WRITE(6,99)
+ 99 FORMAT(//,33H VERSION *COLNEW* OF COLSYS . ,//)
+C
+ IOUT = 6
+ PRECIS = 1.D0
+ 10 PRECIS = PRECIS / 2.D0
+ PRECP1 = PRECIS + 1.D0
+ IF ( PRECP1 .GT. 1.D0 ) GO TO 10
+ PRECIS = PRECIS * 100.D0
+C
+C... in case incorrect input data is detected, the program returns
+C... immediately with iflag=-3.
+C
+ IFLAG = -3
+ IF ( NCOMP .LT. 1 .OR. NCOMP .GT. 20 ) RETURN
+ DO 20 I=1,NCOMP
+ IF ( M(I) .LT. 1 .OR. M(I) .GT. 4 ) RETURN
+ 20 CONTINUE
+C
+C... rename some of the parameters and set default values.
+C
+ NONLIN = IPAR(1)
+ K = IPAR(2)
+ N = IPAR(3)
+ IF ( N .EQ. 0 ) N = 5
+ IREAD = IPAR(8)
+ IGUESS = IPAR(9)
+ IF ( NONLIN .EQ. 0 .AND. IGUESS .EQ. 1 ) IGUESS = 0
+ IF ( IGUESS .GE. 2 .AND. IREAD .EQ. 0 ) IREAD = 1
+ ICARE = IPAR(10)
+ NTOL = IPAR(4)
+ NDIMF = IPAR(5)
+ NDIMI = IPAR(6)
+ NFXPNT = IPAR(11)
+ IPRINT = IPAR(7)
+ MSTAR = 0
+ MMAX = 0
+ DO 30 I = 1, NCOMP
+ MMAX = MAX0 ( MMAX, M(I) )
+ MSTAR = MSTAR + M(I)
+ MT(I) = M(I)
+ 30 CONTINUE
+ IF ( K .EQ. 0 ) K = MAX0( MMAX + 1 , 5 - MMAX )
+ DO 40 I = 1, MSTAR
+ 40 TZETA(I) = ZETA(I)
+ DO 50 I = 1, NTOL
+ LTTOL(I) = LTOL(I)
+ 50 TOLIN(I) = TOL(I)
+ TLEFT = ALEFT
+ TRIGHT = ARIGHT
+ NC = NCOMP
+ KD = K * NCOMP
+C
+C... print the input data for checking.
+C
+ IF ( IPRINT .GT. -1 ) GO TO 80
+ IF ( NONLIN .GT. 0 ) GO TO 60
+ WRITE (IOUT,260) NCOMP, (M(IP), IP=1,NCOMP)
+ GO TO 70
+ 60 WRITE (IOUT,270) NCOMP, (M(IP), IP=1,NCOMP)
+ 70 WRITE (IOUT,280) (ZETA(IP), IP=1,MSTAR)
+ IF ( NFXPNT .GT. 0 )
+ 1 WRITE (IOUT,340) NFXPNT, (FIXPNT(IP), IP=1,NFXPNT)
+ WRITE (IOUT,290) K
+ WRITE (IOUT,300) (LTOL(IP), IP=1,NTOL)
+ WRITE (IOUT,310) (TOL(IP), IP=1,NTOL)
+ IF (IGUESS .GE. 2) WRITE (IOUT,320)
+ IF (IREAD .EQ. 2) WRITE (IOUT,330)
+ 80 CONTINUE
+C
+C... check for correctness of data
+C
+ IF ( K .LT. 0 .OR. K .GT. 7 ) RETURN
+ IF ( N .LT. 0 ) RETURN
+ IF ( IREAD .LT. 0 .OR. IREAD .GT. 2 ) RETURN
+ IF ( IGUESS .LT. 0 .OR. IGUESS .GT. 4 ) RETURN
+ IF ( ICARE .LT. 0 .OR. ICARE .GT. 2 ) RETURN
+ IF ( NTOL .LT. 0 .OR. NTOL .GT. MSTAR ) RETURN
+ IF ( NFXPNT .LT. 0 ) RETURN
+ IF ( IPRINT .LT. (-1) .OR. IPRINT .GT. 1 ) RETURN
+ IF ( MSTAR .LT. 0 .OR. MSTAR .GT. 40 ) RETURN
+ IP = 1
+ DO 100 I = 1, MSTAR
+ IF ( DABS(ZETA(I) - ALEFT) .LT. PRECIS .OR.
+ 1 DABS(ZETA(I) - ARIGHT) .LT. PRECIS ) GO TO 100
+ 90 IF ( IP .GT. NFXPNT ) RETURN
+ IF ( ZETA(I) - PRECIS .LT. FIXPNT(IP) ) GO TO 95
+ IP = IP + 1
+ GO TO 90
+ 95 IF ( ZETA(I) + PRECIS .LT. FIXPNT(IP) ) RETURN
+ 100 CONTINUE
+C
+C... set limits on iterations and initialize counters.
+C... limit = maximum number of newton iterations per mesh.
+C... see subroutine newmsh for the roles of mshlmt , mshflg ,
+C... mshnum , and mshalt .
+C
+ MSHLMT = 3
+ MSHFLG = 0
+ MSHNUM = 1
+ MSHALT = 1
+ LIMIT = 40
+C
+C... compute the maxium possible n for the given sizes of
+C... ispace and fspace.
+C
+ NREC = 0
+ DO 110 I = 1, MSTAR
+ IB = MSTAR + 1 - I
+ IF ( ZETA(IB) .GE. ARIGHT ) NREC = I
+ 110 CONTINUE
+ NFIXI = MSTAR
+ NSIZEI = 3 + KD + MSTAR
+ NFIXF = NREC * (2*MSTAR) + 5 * MSTAR + 3
+ NSIZEF = 4 + 3 * MSTAR + (KD+5) * (KD+MSTAR) +
+ 1(2*MSTAR-NREC) * 2*MSTAR
+ NMAXF = (NDIMF - NFIXF) / NSIZEF
+ NMAXI = (NDIMI - NFIXI) / NSIZEI
+ IF ( IPRINT .LT. 1 ) WRITE(IOUT,350) NMAXF, NMAXI
+ NMAX = MIN0( NMAXF, NMAXI )
+ IF ( NMAX .LT. N ) RETURN
+ IF ( NMAX .LT. NFXPNT+1 ) RETURN
+ IF (NMAX .LT. 2*NFXPNT+2 .AND. IPRINT .LT. 1) WRITE(IOUT,360)
+C
+C... generate pointers to break up fspace and ispace .
+C
+ LXI = 1
+ LG = LXI + NMAX + 1
+ LXIOLD = LG + 2*MSTAR * (NMAX * (2*MSTAR-NREC) + NREC)
+ LW = LXIOLD + NMAX + 1
+ LV = LW + KD**2 * NMAX
+ LZ = LV + MSTAR * KD * NMAX
+ LDMZ = LZ + MSTAR * (NMAX + 1)
+ LDELZ = LDMZ + KD * NMAX
+ LDELDZ = LDELZ + MSTAR * (NMAX + 1)
+ LDQZ = LDELDZ + KD * NMAX
+ LDQDMZ = LDQZ + MSTAR * (NMAX + 1)
+ LRHS = LDQDMZ + KD * NMAX
+ LVALST = LRHS + KD * NMAX + MSTAR
+ LSLOPE = LVALST + 4 * MSTAR * NMAX
+ LACCUM = LSLOPE + NMAX
+ LSCL = LACCUM + NMAX + 1
+ LDSCL = LSCL + MSTAR * (NMAX + 1)
+ LPVTG = 1
+ LPVTW = LPVTG + MSTAR * (NMAX + 1)
+ LINTEG = LPVTW + KD * NMAX
+C
+C... if iguess .ge. 2, move xiold, z, and dmz to their proper
+C... locations in fspace.
+C
+ IF ( IGUESS .LT. 2 ) GO TO 160
+ NOLD = N
+ IF (IGUESS .EQ. 4) NOLD = ISPACE(1)
+ NZ = MSTAR * (NOLD + 1)
+ NDMZ = KD * NOLD
+ NP1 = N + 1
+ IF ( IGUESS .EQ. 4 ) NP1 = NP1 + NOLD + 1
+ DO 120 I=1,NZ
+ 120 FSPACE( LZ+I-1 ) = FSPACE( NP1+I )
+ IDMZ = NP1 + NZ
+ DO 125 I=1,NDMZ
+ 125 FSPACE( LDMZ+I-1 ) = FSPACE( IDMZ+I )
+ NP1 = NOLD + 1
+ IF ( IGUESS .EQ. 4 ) GO TO 140
+ DO 130 I=1,NP1
+ 130 FSPACE( LXIOLD+I-1 ) = FSPACE( LXI+I-1 )
+ GO TO 160
+ 140 DO 150 I=1,NP1
+ 150 FSPACE( LXIOLD+I-1 ) = FSPACE( N+1+I )
+ 160 CONTINUE
+C
+C... initialize collocation points, constants, mesh.
+C
+ CALL CONSTS ( K, RHO, COEF )
+ CALL NEWMSH (3+IREAD, FSPACE(LXI), FSPACE(LXIOLD), DUMMY,
+ 1 DUMMY, DUMMY, DUMMY, DUMMY, NFXPNT, FIXPNT)
+C
+C... determine first approximation, if the problem is nonlinear.
+C
+ IF (IGUESS .GE. 2) GO TO 230
+ NP1 = N + 1
+ DO 210 I = 1, NP1
+ 210 FSPACE( I + LXIOLD - 1 ) = FSPACE( I + LXI - 1 )
+ NOLD = N
+ IF ( NONLIN .EQ. 0 .OR. IGUESS .EQ. 1 ) GO TO 230
+C
+C... system provides first approximation of the solution.
+C... choose z(j) = 0 for j=1,...,mstar.
+C
+ DO 220 I=1, NZ
+ 220 FSPACE( LZ-1+I ) = 0.D0
+ DO 225 I=1, NDMZ
+ 225 FSPACE( LDMZ-1+I ) = 0.D0
+ 230 CONTINUE
+ IF (IGUESS .GE. 2) IGUESS = 0
+ CALL CONTRL (FSPACE(LXI),FSPACE(LXIOLD),FSPACE(LZ),FSPACE(LDMZ),
+ 1 FSPACE(LRHS),FSPACE(LDELZ),FSPACE(LDELDZ),FSPACE(LDQZ),
+ 2 FSPACE(LDQDMZ),FSPACE(LG),FSPACE(LW),FSPACE(LV),
+ 3 FSPACE(LVALST),FSPACE(LSLOPE),FSPACE(LSCL),FSPACE(LDSCL),
+ 4 FSPACE(LACCUM),ISPACE(LPVTG),ISPACE(LINTEG),ISPACE(LPVTW),
+ 5 NFXPNT,FIXPNT,IFLAG,FSUB,DFSUB,GSUB,DGSUB,GUESS )
+C
+C... prepare output
+C
+ ISPACE(1) = N
+ ISPACE(2) = K
+ ISPACE(3) = NCOMP
+ ISPACE(4) = MSTAR
+ ISPACE(5) = MMAX
+ ISPACE(6) = NZ + NDMZ + N + 2
+ K2 = K * K
+ ISPACE(7) = ISPACE(6) + K2 - 1
+ DO 240 I = 1, NCOMP
+ 240 ISPACE(7+I) = M(I)
+ DO 250 I = 1, NZ
+ 250 FSPACE( N+1+I ) = FSPACE( LZ-1+I )
+ IDMZ = N + 1 + NZ
+ DO 255 I = 1, NDMZ
+ 255 FSPACE( IDMZ+I ) = FSPACE( LDMZ-1+I )
+ IC = IDMZ + NDMZ
+ DO 258 I = 1, K2
+ 258 FSPACE( IC+I ) = COEF(I)
+ RETURN
+C----------------------------------------------------------------------
+ 260 FORMAT(/// 37H THE NUMBER OF (LINEAR) DIFF EQNS IS , I3/ 1X,
+ 1 16HTHEIR ORDERS ARE, 20I3)
+ 270 FORMAT(/// 40H THE NUMBER OF (NONLINEAR) DIFF EQNS IS , I3/ 1X,
+ 1 16HTHEIR ORDERS ARE, 20I3)
+ 280 FORMAT(27H SIDE CONDITION POINTS ZETA, 8F10.6, 4( / 27X, 8F10.6))
+ 290 FORMAT(37H NUMBER OF COLLOC PTS PER INTERVAL IS, I3)
+ 300 FORMAT(39H COMPONENTS OF Z REQUIRING TOLERANCES -,8(7X,I2,1X),
+ 1 4(/38X,8I10))
+ 310 FORMAT(33H CORRESPONDING ERROR TOLERANCES -,6X,8D10.2,
+ 1 4(/39X,8D10.2))
+ 320 FORMAT(44H INITIAL MESH(ES) AND Z,DMZ PROVIDED BY USER)
+ 330 FORMAT(27H NO ADAPTIVE MESH SELECTION)
+ 340 FORMAT(10H THERE ARE ,I5,27H FIXED POINTS IN THE MESH - ,
+ 1 10(6F10.6/))
+ 350 FORMAT(44H THE MAXIMUM NUMBER OF SUBINTERVALS IS MIN (, I4,
+ 1 23H (ALLOWED FROM FSPACE),,I4, 24H (ALLOWED FROM ISPACE) ))
+ 360 FORMAT(/53H INSUFFICIENT SPACE TO DOUBLE MESH FOR ERROR ESTIMATE)
+ END
diff --git a/lib-src/colnew/consts.f b/lib-src/colnew/consts.f
new file mode 100644
index 0000000..110df75
--- /dev/null
+++ b/lib-src/colnew/consts.f
@@ -0,0 +1,147 @@
+ SUBROUTINE CONSTS (K, RHO, COEF)
+C
+C**********************************************************************
+C
+C purpose
+C assign (once) values to various array constants.
+C
+C arrays assigned during compilation:
+C cnsts1 - weights for extrapolation error estimate
+C cnsts2 - weights for mesh selection
+C (the above weights come from the theoretical form for
+C the collocation error -- see [3])
+C
+C arrays assigned during execution:
+C wgterr - the particular values of cnsts1 used for current run
+C (depending on k, m)
+C wgtmsh - gotten from the values of cnsts2 which in turn are
+C the constants in the theoretical expression for the
+C errors. the quantities in wgtmsh are 10x the values
+C in cnsts2 so that the mesh selection algorithm
+C is aiming for errors .1x as large as the user
+C requested tolerances.
+C jtol - components of differential system to which tolerances
+C refer (viz, if ltol(i) refers to a derivative of u(j),
+C then jtol(i)=j)
+C root - reciprocals of expected rates of convergence of compo-
+C nents of z(j) for which tolerances are specified
+C rho - the k collocation points on (0,1)
+C coef -
+C acol - the runge-kutta coefficients values at collocation
+C points
+C
+C**********************************************************************
+C
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION RHO(7), COEF(K,1), CNSTS1(28), CNSTS2(28), DUMMY(1)
+C
+ COMMON /COLORD/ KDUM, NCOMP, MSTAR, KD, MMAX, M(20)
+ COMMON /COLBAS/ B(28), ACOL(28,7), ASAVE(28,4)
+ COMMON /COLEST/ TOL(40), WGTMSH(40), WGTERR(40), TOLIN(40),
+ 1 ROOT(40), JTOL(40), LTOL(40), NTOL
+C
+ DATA CNSTS1 / .25D0, .625D-1, 7.2169D-2, 1.8342D-2,
+ 1 1.9065D-2, 5.8190D-2, 5.4658D-3, 5.3370D-3, 1.8890D-2,
+ 2 2.7792D-2, 1.6095D-3, 1.4964D-3, 7.5938D-3, 5.7573D-3,
+ 3 1.8342D-2, 4.673D-3, 4.150D-4, 1.919D-3, 1.468D-3,
+ 4 6.371D-3, 4.610D-3, 1.342D-4, 1.138D-4, 4.889D-4,
+ 5 4.177D-4, 1.374D-3, 1.654D-3, 2.863D-3 /
+ DATA CNSTS2 / 1.25D-1, 2.604D-3, 8.019D-3, 2.170D-5,
+ 1 7.453D-5, 5.208D-4, 9.689D-8, 3.689D-7, 3.100D-6,
+ 2 2.451D-5, 2.691D-10, 1.120D-9, 1.076D-8, 9.405D-8,
+ 3 1.033D-6, 5.097D-13, 2.290D-12, 2.446D-11, 2.331D-10,
+ 4 2.936D-9, 3.593D-8, 7.001D-16, 3.363D-15, 3.921D-14,
+ 5 4.028D-13, 5.646D-12, 7.531D-11, 1.129D-9 /
+C
+C... assign weights for error estimate
+C
+ KOFF = K * ( K + 1 ) / 2
+ IZ = 1
+ DO 10 J = 1, NCOMP
+ MJ = M(J)
+ DO 10 L = 1, MJ
+ WGTERR(IZ) = CNSTS1(KOFF - MJ + L)
+ IZ = IZ + 1
+ 10 CONTINUE
+C
+C... assign array values for mesh selection: wgtmsh, jtol, and root
+C
+ JCOMP = 1
+ MTOT = M(1)
+ DO 40 I = 1, NTOL
+ LTOLI = LTOL(I)
+ 20 CONTINUE
+ IF ( LTOLI .LE. MTOT ) GO TO 30
+ JCOMP = JCOMP + 1
+ MTOT = MTOT + M(JCOMP)
+ GO TO 20
+ 30 CONTINUE
+ JTOL(I) = JCOMP
+ WGTMSH(I) = 1.D1 * CNSTS2(KOFF+LTOLI-MTOT) / TOLIN(I)
+ ROOT(I) = 1.D0 / DFLOAT(K+MTOT-LTOLI+1)
+ 40 CONTINUE
+C
+C... specify collocation points
+C
+ GO TO (50,60,70,80,90,100,110), K
+ 50 RHO(1) = 0.D0
+ GO TO 120
+ 60 RHO(2) = .57735026918962576451D0
+ RHO(1) = - RHO(2)
+ GO TO 120
+ 70 RHO(3) = .77459666924148337704D0
+ RHO(2) = .0D0
+ RHO(1) = - RHO(3)
+ GO TO 120
+ 80 RHO(4) = .86113631159405257523D0
+ RHO(3) = .33998104358485626480D0
+ RHO(2) = - RHO(3)
+ RHO(1) = - RHO(4)
+ GO TO 120
+ 90 RHO(5) = .90617984593866399280D0
+ RHO(4) = .53846931010568309104D0
+ RHO(3) = .0D0
+ RHO(2) = - RHO(4)
+ RHO(1) = - RHO(5)
+ GO TO 120
+ 100 RHO(6) = .93246951420315202781D0
+ RHO(5) = .66120938646626451366D0
+ RHO(4) = .23861918608319690863D0
+ RHO(3) = -RHO(4)
+ RHO(2) = -RHO(5)
+ RHO(1) = -RHO(6)
+ GO TO 120
+ 110 RHO(7) = .949107991234275852452D0
+ RHO(6) = .74153118559939443986D0
+ RHO(5) = .40584515137739716690D0
+ RHO(4) = 0.D0
+ RHO(3) = -RHO(5)
+ RHO(2) = -RHO(6)
+ RHO(1) = -RHO(7)
+ 120 CONTINUE
+C
+C... map (-1,1) to (0,1) by t = .5 * (1. + x)
+C
+ DO 130 J = 1, K
+ RHO(J) = .5D0 * (1.D0 + RHO(J))
+ 130 CONTINUE
+C
+C... now find runge-kutta coeffitients b, acol and asave
+C... the values of asave are to be used in newmsh and errchk .
+C
+ DO 140 J = 1, K
+ DO 135 I = 1, K
+ 135 COEF(I,J) = 0.D0
+ COEF(J,J) = 1.D0
+ CALL VMONDE (RHO, COEF(1,J), K)
+ 140 CONTINUE
+ CALL RKBAS ( 1.D0, COEF, K, MMAX, B, DUMMY, 0)
+ DO 150 I = 1, K
+ CALL RKBAS ( RHO(I), COEF, K, MMAX, ACOL(1,I), DUMMY, 0)
+ 150 CONTINUE
+ CALL RKBAS ( 1.D0/6.D0, COEF, K, MMAX, ASAVE(1,1), DUMMY, 0)
+ CALL RKBAS ( 1.D0/3.D0, COEF, K, MMAX, ASAVE(1,2), DUMMY, 0)
+ CALL RKBAS ( 2.D0/3.D0, COEF, K, MMAX, ASAVE(1,3), DUMMY, 0)
+ CALL RKBAS ( 5.D0/6.D0, COEF, K, MMAX, ASAVE(1,4), DUMMY, 0)
+ RETURN
+ END
diff --git a/lib-src/colnew/contrl.f b/lib-src/colnew/contrl.f
new file mode 100644
index 0000000..52a8b0a
--- /dev/null
+++ b/lib-src/colnew/contrl.f
@@ -0,0 +1,489 @@
+ SUBROUTINE CONTRL (XI, XIOLD, Z, DMZ, RHS, DELZ, DELDMZ,
+ 1 DQZ, DQDMZ, G, W, V, VALSTR, SLOPE, SCALE, DSCALE,
+ 2 ACCUM, IPVTG, INTEGS, IPVTW, NFXPNT, FIXPNT, IFLAG,
+ 3 FSUB, DFSUB, GSUB, DGSUB, GUESS )
+C
+C**********************************************************************
+C
+C purpose
+C this subroutine is the actual driver. the nonlinear iteration
+C strategy is controlled here ( see [4] ). upon convergence, errchk
+C is called to test for satisfaction of the requested tolerances.
+C
+C variables
+C
+C check - maximum tolerance value, used as part of criteria for
+C checking for nonlinear iteration convergence
+C relax - the relaxation factor for damped newton iteration
+C relmin - minimum allowable value for relax (otherwise the
+C jacobian is considered singular).
+C rlxold - previous relax
+C rstart - initial value for relax when problem is sensitive
+C ifrz - number of fixed jacobian iterations
+C lmtfrz - maximum value for ifrz before performing a reinversion
+C iter - number of iterations (counted only when jacobian
+C reinversions are performed).
+C xi - current mesh
+C xiold - previous mesh
+C ipred = 0 if relax is determined by a correction
+C = 1 if relax is determined by a prediction
+C ifreez = 0 if the jacobian is to be updated
+C = 1 if the jacobian is currently fixed (frozen)
+C iconv = 0 if no previous convergence has been obtained
+C = 1 if convergence on a previous mesh has been obtained
+C icare =-1 no convergence occurred (used for regular problems)
+C = 0 a regular problem
+C = 1 a sensitive problem
+C = 2 used for continuation (see description of ipar(10)
+C in colnew).
+C rnorm - norm of rhs (right hand side) for current iteration
+C rnold - norm of rhs for previous iteration
+C anscl - scaled norm of newton correction
+C anfix - scaled norm of newton correction at next step
+C anorm - scaled norm of a correction obtained with jacobian fixed
+C nz - number of components of z (see subroutine approx)
+C ndmz - number of components of dmz (see subroutine approx)
+C imesh - a control variable for subroutines newmsh and errchk
+C = 1 the current mesh resulted from mesh selection
+C or is the initial mesh.
+C = 2 the current mesh resulted from doubling the
+C previous mesh
+C
+C**********************************************************************
+C
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION XI(1), XIOLD(1), Z(1), DMZ(1), RHS(1)
+ DIMENSION G(1), W(1), V(1), VALSTR(1), SLOPE(1), ACCUM(1)
+ DIMENSION DELZ(1), DELDMZ(1), DQZ(1), DQDMZ(1) , FIXPNT(1)
+ DIMENSION DUMMY(1), SCALE(1), DSCALE(1)
+ DIMENSION INTEGS(1), IPVTG(1), IPVTW(1)
+C
+ COMMON /COLOUT/ PRECIS, IOUT, IPRINT
+ COMMON /COLORD/ K, NCOMP, MSTAR, KD, MMAX, M(20)
+ COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ
+ COMMON /COLMSH/ MSHFLG, MSHNUM, MSHLMT, MSHALT
+ COMMON /COLSID/ ZETA(40), ALEFT, ARIGHT, IZETA, IDUM
+ COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS
+ COMMON /COLEST/ TOL(40), WGTMSH(40), WGTERR(40), TOLIN(40),
+ 1 ROOT(40), JTOL(40), LTOL(40), NTOL
+C
+ EXTERNAL FSUB, DFSUB, GSUB, DGSUB, GUESS
+C
+C... constants for control of nonlinear iteration
+C
+ RELMIN = 1.D-3
+ RSTART = 1.D-2
+ LMTFRZ = 4
+C
+C... compute the maximum tolerance
+C
+ CHECK = 0.D0
+ DO 10 I = 1, NTOL
+ 10 CHECK = DMAX1 ( TOLIN(I), CHECK )
+ IMESH = 1
+ ICONV = 0
+ IF ( NONLIN .EQ. 0 ) ICONV = 1
+ ICOR = 0
+ NOCONV = 0
+ MSING = 0
+C
+C... the main iteration begins here .
+C... loop 20 is executed until error tolerances are satisfied or
+C... the code fails (due to a singular matrix or storage limitations)
+C
+ 20 CONTINUE
+C
+C... initialization for a new mesh
+C
+ ITER = 0
+ IF ( NONLIN .GT. 0 ) GO TO 50
+C
+C... the linear case.
+C... set up and solve equations
+C
+ CALL LSYSLV (MSING, XI, XIOLD, DUMMY, DUMMY, Z, DMZ, G,
+ 1 W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 0,
+ 2 FSUB, DFSUB, GSUB, DGSUB, GUESS )
+C
+C... check for a singular matrix
+C
+ IF ( MSING .EQ. 0 ) GO TO 400
+ 30 IF ( MSING .LT. 0 ) GO TO 40
+ IF ( IPRINT .LT. 1 ) WRITE (IOUT,495)
+ GO TO 460
+ 40 IF ( IPRINT .LT. 1 ) WRITE (IOUT,490)
+ IFLAG = 0
+ RETURN
+C
+C... iteration loop for nonlinear case
+C... define the initial relaxation parameter (= relax)
+C
+ 50 RELAX = 1.D0
+C
+C... check for previous convergence and problem sensitivity
+C
+ IF ( ICARE .EQ. 1 .OR. ICARE .EQ. (-1) ) RELAX = RSTART
+ IF ( ICONV .EQ. 0 ) GO TO 160
+C
+C... convergence on a previous mesh has been obtained. thus
+C... we have a very good initial approximation for the newton
+C... process. proceed with one full newton and then iterate
+C... with a fixed jacobian.
+C
+ IFREEZ = 0
+C
+C... evaluate right hand side and its norm and
+C... find the first newton correction
+C
+ CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,
+ 1 W, V, RHS, DQDMZ, INTEGS, IPVTG, IPVTW, RNOLD, 1,
+ 2 FSUB, DFSUB, GSUB, DGSUB, GUESS )
+C
+ IF ( IPRINT .LT. 0 ) WRITE(IOUT,530)
+ IF ( IPRINT .LT. 0 ) WRITE (IOUT,510) ITER, RNOLD
+ GO TO 70
+C
+C... solve for the next iterate .
+C... the value of ifreez determines whether this is a full
+C... newton step (=0) or a fixed jacobian iteration (=1).
+C
+ 60 IF ( IPRINT .LT. 0 ) WRITE (IOUT,510) ITER, RNORM
+ RNOLD = RNORM
+ CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,
+ 1 W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM,
+ 2 3+IFREEZ, FSUB, DFSUB, GSUB, DGSUB, GUESS )
+C
+C... check for a singular matrix
+C
+ 70 IF ( MSING .NE. 0 ) GO TO 30
+ IF ( IFREEZ .EQ. 1 ) GO TO 80
+C
+C... a full newton step
+C
+ ITER = ITER + 1
+ IFRZ = 0
+ 80 CONTINUE
+C
+C... update z and dmz , compute new rhs and its norm
+C
+ DO 90 I = 1, NZ
+ Z(I) = Z(I) + DELZ(I)
+ 90 CONTINUE
+ DO 100 I = 1, NDMZ
+ DMZ(I) = DMZ(I) + DELDMZ(I)
+ 100 CONTINUE
+ CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,
+ 1 W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 2,
+ 2 FSUB, DFSUB, GSUB, DGSUB, GUESS )
+C
+C... check monotonicity. if the norm of rhs gets smaller,
+C... proceed with a fixed jacobian; else proceed cautiously,
+C... as if convergence has not been obtained before (iconv=0).
+C
+ IF ( RNORM .LT. PRECIS ) GO TO 390
+ IF ( RNORM .GT. RNOLD ) GO TO 130
+ IF ( IFREEZ .EQ. 1 ) GO TO 110
+ IFREEZ = 1
+ GO TO 60
+C
+C... verify that the linear convergence with fixed jacobian
+C... is fast enough.
+C
+ 110 IFRZ = IFRZ + 1
+ IF ( IFRZ .GE. LMTFRZ ) IFREEZ = 0
+ IF ( RNOLD .LT. 4.D0*RNORM ) IFREEZ = 0
+C
+C... check convergence (iconv = 1).
+C
+ DO 120 IT = 1, NTOL
+ INZ = LTOL(IT)
+ DO 120 IZ = INZ, NZ, MSTAR
+ IF ( DABS(DELZ(IZ)) .GT.
+ 1 TOLIN(IT) * (DABS(Z(IZ)) + 1.D0)) GO TO 60
+ 120 CONTINUE
+C
+C... convergence obtained
+C
+ IF ( IPRINT .LT. 1 ) WRITE (IOUT,560) ITER
+ GO TO 400
+C
+C... convergence of fixed jacobian iteration failed.
+C
+ 130 IF ( IPRINT .LT. 0 ) WRITE (IOUT,510) ITER, RNORM
+ IF ( IPRINT .LT. 0 ) WRITE (IOUT,540)
+ ICONV = 0
+ RELAX = RSTART
+ DO 140 I = 1, NZ
+ Z(I) = Z(I) - DELZ(I)
+ 140 CONTINUE
+ DO 150 I = 1, NDMZ
+ DMZ(I) = DMZ(I) - DELDMZ(I)
+ 150 CONTINUE
+C
+C... update old mesh
+C
+ NP1 = N + 1
+ DO 155 I = 1, NP1
+ 155 XIOLD(I) = XI(I)
+ NOLD = N
+C
+ ITER = 0
+C
+C... no previous convergence has been obtained. proceed
+C... with the damped newton method.
+C... evaluate rhs and find the first newton correction.
+C
+ 160 IF(IPRINT .LT. 0) WRITE (IOUT,500)
+ CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,
+ 1 W, V, RHS, DQDMZ, INTEGS, IPVTG, IPVTW, RNOLD, 1,
+ 2 FSUB, DFSUB, GSUB, DGSUB, GUESS )
+C
+C... check for a singular matrix
+C
+ IF ( MSING .NE. 0 ) GO TO 30
+C
+C... bookkeeping for first mesh
+C
+ IF ( IGUESS .EQ. 1 ) IGUESS = 0
+C
+C... find initial scaling
+C
+ CALL SKALE (N, MSTAR, KD, Z, XI, SCALE, DSCALE)
+ GO TO 220
+C
+C... main iteration loop
+C
+ 170 RNOLD = RNORM
+ IF ( ITER .GE. LIMIT ) GO TO 430
+C
+C... update scaling
+C
+ CALL SKALE (N, MSTAR, KD, Z, XI, SCALE, DSCALE)
+C
+C... compute norm of newton correction with new scaling
+C
+ ANSCL = 0.D0
+ DO 180 I = 1, NZ
+ ANSCL = ANSCL + (DELZ(I) * SCALE(I))**2
+ 180 CONTINUE
+ DO 190 I = 1, NDMZ
+ ANSCL = ANSCL + (DELDMZ(I) * DSCALE(I))**2
+ 190 CONTINUE
+ ANSCL = DSQRT(ANSCL / DFLOAT(NZ+NDMZ))
+C
+C... find a newton direction
+C
+ CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,
+ 1 W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 3,
+ 2 FSUB, DFSUB, GSUB, DGSUB, GUESS )
+C
+C... check for a singular matrix
+C
+ IF ( MSING .NE. 0 ) GO TO 30
+C
+C... predict relaxation factor for newton step.
+C
+ ANDIF = 0.D0
+ DO 200 I = 1, NZ
+ ANDIF = ANDIF + ((DQZ(I) - DELZ(I)) * SCALE(I))**2
+ 200 CONTINUE
+ DO 210 I = 1, NDMZ
+ ANDIF = ANDIF + ((DQDMZ(I) - DELDMZ(I)) * DSCALE(I))**2
+ 210 CONTINUE
+ ANDIF = DSQRT(ANDIF/DFLOAT(NZ+NDMZ) + PRECIS)
+ RELAX = RELAX * ANSCL / ANDIF
+ IF ( RELAX .GT. 1.D0 ) RELAX = 1.D0
+ 220 RLXOLD = RELAX
+ IPRED = 1
+ ITER = ITER + 1
+C
+C... determine a new z and dmz and find new rhs and its norm
+C
+ DO 230 I = 1, NZ
+ Z(I) = Z(I) + RELAX * DELZ(I)
+ 230 CONTINUE
+ DO 240 I = 1, NDMZ
+ DMZ(I) = DMZ(I) + RELAX * DELDMZ(I)
+ 240 CONTINUE
+ 250 CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DQZ, DQDMZ, G,
+ 1 W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 2,
+ 2 FSUB, DFSUB, GSUB, DGSUB, GUESS )
+C
+C... compute a fixed jacobian iterate (used to...
[truncated message content] |
