🚀 go-pugleaf

Archive-name: Lsoda.mex
Submitted-by: fitz@rpi.edu (Brian P. Fitzgerald)
Version-number: 1.0

Pramit Sarma writes:
>In article <1993Apr3.184136.13054@vax.oxford.ac.uk>, reese@vax.oxford.ac.uk writes:
>|> 
...
>An alternative - harass the hell of out of (a) yourself, or (b)
>one of the local systems gurus to convert one of easily available
>(and far, far better performing than ODE**) stiff integrators, such 
>as LSODE,LSODA,DDASSL,DDASAC, which are FORTRANed, into a MEX

Good idea.  (Maybe I'll finish this project someday.
It works now, though.)

Distribution: na
              ^^ hmm. Original thread posted from Oxford, but distribution of
followup was "na".  "na" means "North America", not "Not Applicable"

(You compile the matlab gateway objects and link to the real odepack,
which you can get from research.att.com.  ok?)

Brian


#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  The tool that generated this
# shell archive is called "shar", and is available by anonymous ftp
# from ftp.uu.net in subdirectory /usenet/comp.sources.unix, and from many
# other places. Check 'archie' for the latest locations.  If this archive
# is complete, you will see the following message at the end:
#		"End of shell archive."
# Contents:  Makefile Lsoda.c Odepack.h f.c jac.c
# Wrapped by fitzgb@mml0.meche.rpi.edu on Fri Apr  9 20:53:44 1993
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'Makefile' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'Makefile'\"
else
echo shar: Extracting \"'Makefile'\" \(798 characters\)
sed "s/^X//" >'Makefile' <<'END_OF_FILE'
X# Makefile by Brian P. Fitzgerald 10/15/92
X# Mechanics of Materials Laboratory
X# RPI, Troy, NY
X
X# support matlab 4.0a 3/3/93
X
XML_ROOT		= /campus/mathworks/matlab-4.0a/@sys/src
X
XINCL		= -I$(ML_ROOT)/extern/include
XCCOPTS		= -O
XCFLAGS		= $(DEFINES) $(INCL) $(CCOPTS)
XLINTFLAGS	= $(DEFINES) $(INCL)
X
XLIBROOT		= Odepack
XLIBNAME		= lib$(LIBROOT).a
X
XSYS		= @sys
XDIR		= /dept/meche/netlib/odepack/0.2/$(SYS)/matlab/lib
XLIB		= $(DIR)/$(LIBNAME)
X
XLNS		= Lsoda.ln f.ln jac.ln
XMEMS		= \
X		$(LIB)(f.o) \
X		$(LIB)(jac.o)
X
XOBJS		= Lsoda.o
X
Xall: $(LIB) $(OBJS)
X
X$(LIB):	$(MEMS)
X	ranlib $(LIB)
X
Xinstall: $(OBJS)
X	cp $(OBJS) $(DIR)
X
Xlint: $(LNS)
X	lint $(LNS)
X
XDIST		= Makefile Lsoda.c Odepack.h f.c jac.c
X
Xshar:	$(DIST)
X	shar $(DIST) > Odepack.shar
X
Xclean:
X	rm -f *.o *.ln core typescript *.BAK
X
X$(MEMS):	Odepack.h
END_OF_FILE
if test 798 -ne `wc -c <'Makefile'`; then
    echo shar: \"'Makefile'\" unpacked with wrong size!
fi
# end of 'Makefile'
fi
if test -f 'Lsoda.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'Lsoda.c'\"
else
echo shar: Extracting \"'Lsoda.c'\" \(5222 characters\)
sed "s/^X//" >'Lsoda.c' <<'END_OF_FILE'
X/*-
X * Brian P. Fitzgerald
X * Mechanics of Materials Laboratory
X * Rensselaer Polytechnic Institute
X * Troy, NY
X * October 16, 1992
X *
X * Lsoda.c	MEX file implementation of the fortran subroutine:
X *
X * subroutine lsoda (f, neq, y, t, tout, itol, rtol, atol, itask,
X *	istate, iopt, rwork, lrw, iwork, liw, jac, jt)
X *
X * Call from matlab as:
X *
X * [yret tret]=Lsoda(f, neq, y, t, tout, itol, rtol, atol, itask,
X *	istate, iopt, rwork, iwork, jac, jt);
X *
X * If tout is a vector of desired output points, then yret will contain the
X * solution at each of these points.
X *
X * omit lrw, liw.  Only first 20 elements of rwork and iwork used.
X *
X * f and jac should be strings -- the names of matlab functions (.m or .mex
X * files).
X *
X * Return:
X *
X * tret  - integration time points (row-vector).
X *
X * yret  - solution, one column vector per time point
X *
X * assume a fixed length rwork (not dynamic)
X *
X */
X
X#include "Odepack.h"
X
X/* Input Arguments */
X
X#define	F		prhs[0]
X#define	NEQ		prhs[1]
X#define	Y		prhs[2]
X#define	T		prhs[3]
X#define	TOUT		prhs[4]
X#define	ITOL		prhs[5]
X#define	RTOL		prhs[6]
X#define	ATOL		prhs[7]
X#define	ITASK		prhs[8]
X#define	ISTATE		prhs[9]
X#define	IOPT		prhs[10]
X#define	RWORK		prhs[11]
X#define	IWORK		prhs[12]
X#define	JAC		prhs[13]
X#define	JT		prhs[14]
X
X#define NEQN		*neq
X
X/* Output Arguments */
X
X#define	YRET		plhs[0]
X#define	TRET		plhs[1]
X
Xchar            f_name[64];
Xchar            jac_name[64];
X
Xint
Xget_lrw(neq, iopt, iwork, jt)
X    int             neq, iopt, *iwork, jt;
X{
X    int             lmat, lrn, lrs;
X    int             nyh = neq;
X
X    int             mxordn =
X    iopt == 1 && *(iwork + 7) && *(iwork + 7) < 12 ? *(iwork + 7) : 12;
X    int             mxords =
X    iopt == 1 && *(iwork + 8) && *(iwork + 8) < 5 ? *(iwork + 8) : 5;
X
X    if (jt == 1 || jt == 2)
X	lmat = neq * neq + 2;
X    else if (jt == 4 || jt == 5) {
X	int             ml = *iwork;
X	int             mu = *(iwork + 1);
X	lmat = (2 * ml + mu + 1) * neq + 2;
X    } else
X	mexErrMsgTxt("Invalid jt for Lsoda");
X
X    lrn = 20 + nyh * (mxordn + 1) + 3 * neq;
X    lrs = 20 + nyh * (mxords + 1) + 3 * neq + lmat;
X    return max(lrn, lrs);
X}
X
XmexFunction(nlhs, plhs, nrhs, prhs)
X    int             nlhs, nrhs;
X    Matrix         *plhs[], *prhs[];
X{
X    double         *t, *y, *tout, *rtol, *atol, *rwork;
X    double         *neq_in, *rwork_in, *iwork_in;	/* input "integers" */
X    double         *yret, *tret;
X    double         *dp;
X
X    int             itol, itask, istate, iopt, lrw, liw, jt;
X    int			l_neq, l_iwork, l_rwork, l_tout;
X    int            *neq, *iwork;
X    register unsigned int l, m;
X
X    if (nrhs != 15)
X	mexErrMsgTxt("Wrong number of input arguments for Lsoda");
X    else if (nlhs != 2)
X	mexErrMsgTxt("Wrong number of output arguments for Lsoda");
X
X    if (!IS_VECTOR(NEQ))
X	mexErrMsgTxt("Lsoda: neq: must be a vector or a scalar");
X    if (!IS_VECTOR(Y))
X	mexErrMsgTxt("Lsoda: y: must be a vector");
X    if (!IS_VECTOR(TOUT))
X	mexErrMsgTxt("Lsoda: tout: must be a vector or a scalar");
X    if (!IS_VECTOR(ATOL))
X	mexErrMsgTxt("Lsoda: atol: must be a vector or a scalar");
X    if (!IS_VECTOR(RWORK))
X	mexErrMsgTxt("Lsoda: rwork: must be a vector");
X    if (!IS_VECTOR(IWORK))
X	mexErrMsgTxt("Lsoda: iwork: must be a vector");
X
X    for (m = 0; m < mxGetN(F); m++)	/* Get user derivative function name */
X	f_name[m] = (char) *(mxGetPr(F) + m);
X    f_name[m] = '\0';
X    neq_in = mxGetPr(NEQ);
X    l_neq = LENGTH(NEQ);
X    neq = create_int_array(l_neq);	/* (could be a vector) */
X    for (m = 0; m < l_neq; m++)
X	*(neq + m) = *(neq_in + m);
X    dp = mxGetPr(ITOL);
X    itol = (int) *dp;
X    rtol = mxGetPr(RTOL);
X    atol = mxGetPr(ATOL);
X    dp = mxGetPr(ITASK);
X    itask = (int) *dp;
X    dp = mxGetPr(ISTATE);
X    istate = (int) *dp;
X    dp = mxGetPr(IOPT);
X    iopt = (int) *dp;
X    iwork_in = mxGetPr(IWORK);
X    liw = 20 + NEQN;
X    iwork = create_int_array(liw);
X    l_iwork = LENGTH(IWORK);
X    for (m = 0; m < min(20, l_iwork); m++)
X	*(iwork + m) = *(iwork_in + m);
X    dp = mxGetPr(JT);
X    jt = (int) *dp;
X    rwork_in = mxGetPr(RWORK);
X    lrw = get_lrw(NEQN, iopt, iwork, jt);
X    rwork = create_real_array(lrw);
X    l_rwork = LENGTH(RWORK);
X    for (m = 0; m < min(20, l_rwork); m++)
X	*(rwork + m) = *(rwork_in + m);
X    for (m = 0; m < mxGetN(JAC); m++)/* Get user jacobian function name */
X	jac_name[m] = (char) *(mxGetPr(JAC) + m);
X    jac_name[m] = '\0';
X
X    l_tout = LENGTH(TOUT);
X    YRET = mxCreateFull(NEQN, l_tout, REAL);
X    TRET = mxCreateFull(l_tout, 1, REAL);
X
X    yret = mxGetPr(YRET);
X    tret = mxGetPr(TRET);
X    y = mxGetPr(Y);
X    t = mxGetPr(T);
X    tout = mxGetPr(TOUT);
X
X    l = l_tout;
X    while (l--) {		/* do a loop over the tout vector */
X	LSODA(f_, neq, y, t, tout++, &itol, rtol, atol, &itask, &istate,
X	       &iopt, rwork, &lrw, iwork, &liw, jac_, &jt);
X
X	if (istate != 2) {
X	    mexPrintf("return from lsoda istate = %d\n", istate);
X	    mexErrMsgTxt("error in lsoda");
X	}
X	dp = y;
X	for (m = 0; m < NEQN; m++)
X	    *yret++ = *dp++;
X	*tret++ = *t;
X    }
X    /* pass back some interesting information */
X    for (m = 0; m < min(20, l_rwork); m++)
X	*(rwork_in + m) = *(rwork + m);
X    for (m = 0; m < min(20, l_iwork); m++)
X	*(iwork_in + m) = *(iwork + m);
X}
END_OF_FILE
if test 5222 -ne `wc -c <'Lsoda.c'`; then
    echo shar: \"'Lsoda.c'\" unpacked with wrong size!
fi
# end of 'Lsoda.c'
fi
if test -f 'Odepack.h' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'Odepack.h'\"
else
echo shar: Extracting \"'Odepack.h'\" \(419 characters\)
sed "s/^X//" >'Odepack.h' <<'END_OF_FILE'
X#include <math.h>
X#include "mex.h"
X
X#define	max(A, B)	((A) > (B) ? (A) : (B))
X#define	min(A, B)	((A) < (B) ? (A) : (B))
X
X#define IS_VECTOR(p)	(min( mxGetM(p) , mxGetN(p)) == 1)
X#define LENGTH(p)	(max( mxGetM(p) , mxGetN(p)))
X
Xextern char     f_name[];
Xextern char     jac_name[];
X
Xvoid            f_();
Xvoid            jac_();
X
X#ifdef sun
X#define LSODA lsoda_
X#else
X#define LSODA lsoda
X#endif
X
Xvoid            LSODA();
END_OF_FILE
if test 419 -ne `wc -c <'Odepack.h'`; then
    echo shar: \"'Odepack.h'\" unpacked with wrong size!
fi
# end of 'Odepack.h'
fi
if test -f 'f.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'f.c'\"
else
echo shar: Extracting \"'f.c'\" \(975 characters\)
sed "s/^X//" >'f.c' <<'END_OF_FILE'
X/*-
X * Brian P. Fitzgerald
X * Mechanics of Materials Laboratory Rensselaer
X * Polytechnic Institute Troy, NY
X * October 16, 1992
X * 
X * user needs to write a matlab function
X * 
X * function ydot=fex(neq,t,y)
X * 
X */
X
X#include "Odepack.h"
X
X#define	YDOT		f_plhs[0]
X
X#define	NEQ		f_prhs[0]
X#define	T		f_prhs[1]
X#define	Y		f_prhs[2]
X
Xvoid
Xf_(neq, t, y, ydot)
X    int            *neq;
X    double         *t, *y, *ydot;
X{
X    Matrix         *f_plhs[1], *f_prhs[3];
X    double         *dp;
X    register unsigned int m;
X
X    NEQ = mxCreateFull(1, 1, REAL);
X    T = mxCreateFull(1, 1, REAL);
X    Y = mxCreateFull(*neq, 1, REAL);
X
X    dp = mxGetPr(NEQ);
X    *dp = (double) *neq;
X    dp = mxGetPr(T);
X    *dp = *t;
X    m = LENGTH(Y);
X    dp = mxGetPr(Y);
X    while (m--)
X	*dp++ = *y++;
X
X    matlab_fcn(1, f_plhs, 3, f_prhs, f_name);
X
X    if (!IS_VECTOR(YDOT))
X	mexErrMsgTxt("f: ydot: must be a vector");
X    m = LENGTH(YDOT);
X    dp = mxGetPr(YDOT);
X    while (m--)
X	*ydot++ = *dp++;
X}
END_OF_FILE
if test 975 -ne `wc -c <'f.c'`; then
    echo shar: \"'f.c'\" unpacked with wrong size!
fi
# end of 'f.c'
fi
if test -f 'jac.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'jac.c'\"
else
echo shar: Extracting \"'jac.c'\" \(1506 characters\)
sed "s/^X//" >'jac.c' <<'END_OF_FILE'
X/*-
X * Brian P. Fitzgerald
X * Mechanics of Materials Laboratory
X * Rensselaer Polytechnic Institute
X * Troy, NY
X * October 16, 1992
X *
X */
X
X#include "Odepack.h"
X
X#define	PD		f_plhs[0]
X
X#define NEQ		f_prhs[0]
X#define T		f_prhs[1]
X#define Y		f_prhs[2]
X#define ML		f_prhs[3]
X#define MU		f_prhs[4]
X#define NROWPD		f_prhs[5]
X
Xvoid
Xjac_(neq, t, y, ml, mu, pd, nrowpd)
X    int            *neq, *ml, *mu, *nrowpd;
X    double         *t, *y, *pd;
X{
X    Matrix         *f_plhs[1], *f_prhs[6];
X    double         *dp;
X    register unsigned int m;
X
X    NEQ = mxCreateFull(1, 1, REAL);
X    T = mxCreateFull(1, 1, REAL);
X    Y = mxCreateFull(*neq, 1, REAL);
X    ML = mxCreateFull(1, 1, REAL);
X    MU = mxCreateFull(1, 1, REAL);
X    NROWPD = mxCreateFull(1, 1, REAL);
X
X    dp = mxGetPr(NEQ);
X    *dp = (double) *neq;
X    dp = mxGetPr(T);
X    *dp = *t;
X    m = LENGTH(Y);
X    dp = mxGetPr(Y);
X    while (m--)
X	*dp++ = *y++;
X    dp = mxGetPr(ML);
X    *dp = (double) *ml;
X    dp = mxGetPr(MU);
X    *dp = (double) *mu;
X    /*-
X	nrowpd = neq if jt = 1
X	nrowpd = meband if jt = 4
X	where,
X	mband = ml + mu + 1
X	meband = mband + ml
X     */
X    dp = mxGetPr(NROWPD);
X    *dp = (double) *nrowpd;
X
X    mexCallMATLAB(1, f_plhs, 6, f_prhs, jac_name);
X
X    m = mxGetM(PD) * mxGetN(PD);
X    if (m != *nrowpd * *neq) {
X	mexPrintf("%s returned %d x %d array pd\n",
X		jac_name, mxGetM(PD), mxGetN(PD));
X	mexPrintf(" needed %d x %d array pd\n", *nrowpd, *neq);
X	mexErrMsgTxt();
X    }
X    dp = mxGetPr(PD);
X    while (m--)
X	*pd++ = *dp++;
X}
END_OF_FILE
if test 1506 -ne `wc -c <'jac.c'`; then
    echo shar: \"'jac.c'\" unpacked with wrong size!
fi
# end of 'jac.c'
fi
echo shar: End of shell archive.
exit 0