Index: mttroot/mtt/bin/mtt
==================================================================
--- mttroot/mtt/bin/mtt
+++ mttroot/mtt/bin/mtt
@@ -13,10 +13,13 @@
###############################################################
## Version control history
###############################################################
## $Header$
## $Log$
+## Revision 1.292 2001/02/05 17:27:40 gawthrop
+## Make sure _def.r exists before creating _state.txt
+##
## Revision 1.291 2000/12/27 14:50:40 peterg
## This is the first CVS version (4.9).
## Commented out code now deleted
##
## Revision 1.290 2000/12/05 09:59:37 peterg
@@ -1747,13 +1750,13 @@
if [ -f "$filename" ]; then
echo $filename exists
else
if [ -n "$Verbose" ]; then
- echo make -s -f $2_rep.make "SYS=$1" "LANG=$3" "ARG=$4"
+ echo make -s -f $2_rep.make "SYS=$1" "LANG=$3" "ARG=$4" "OPTS=$mtt_switches"
fi
- make -s -f $2_rep.make "SYS=$1" "LANG=$3" "ARG=$4"
+ make -s -f $2_rep.make "SYS=$1" "LANG=$3" "ARG=$4" "OPTS=$mtt_switches"
if [ -n "$4" ]; then
echo Copying $1_$2$__ARGS.$ps
cp $1_$2$__ARGS.$ps ..
fi
fi
ADDED mttroot/mtt/lib/cc/mtt_euler.cc
Index: mttroot/mtt/lib/cc/mtt_euler.cc
==================================================================
--- /dev/null
+++ mttroot/mtt/lib/cc/mtt_euler.cc
@@ -0,0 +1,36 @@
+#include <octave/oct.h>
+
+DEFUN_DLD (mtt_euler, args, ,
+ "euler integration method")
+{
+#ifdef OCTAVE_DEV
+ ColumnVector x = args(0).column_vector_value ();
+ const ColumnVector dx = args(1).column_vector_value ();
+ const double ddt = args(2).double_value ();
+ const int Nx = static_cast<int> (args(3).double_value ());
+ const ColumnVector openx = args(4).column_vector_value ();
+#else
+ const ColumnVector x = args(0).vector_value ();
+ const ColumnVector dx = args(1).vector_value ();
+ const double ddt = args(2).double_value ();
+ const int Nx = static_cast<int> (args(3).double_value ());
+ const ColumnVector openx = args(4).vector_value ();
+#endif
+
+ register int i, n;
+
+ n = Nx;
+ for (i = 0; i < Nx; i++)
+ {
+ if (0 != openx (i))
+ {
+ x (i) = 0.0;
+ }
+ else
+ {
+ x (i) += dx (i) * ddt;
+ }
+ }
+
+ return octave_value (x);
+}
ADDED mttroot/mtt/lib/cc/mtt_implicit.cc
Index: mttroot/mtt/lib/cc/mtt_implicit.cc
==================================================================
--- /dev/null
+++ mttroot/mtt/lib/cc/mtt_implicit.cc
@@ -0,0 +1,139 @@
+#include <octave/oct.h>
+#include "useful-functions.hh"
+
+#include <octave/xdiv.h>
+static inline int
+result_ok (int info, double rcond, int warn = 1)
+{
+ assert (info != -1);
+
+ if (info == -2)
+ {
+ if (warn)
+ warning ("matrix singular to machine precision, rcond = %g", rcond);
+ else
+ error ("matrix singular to machine precision, rcond = %g", rcond);
+
+ return 0;
+ }
+ else
+ return 1;
+}
+bool
+mx_leftdiv_conform (const Matrix& a, const ColumnVector& b)
+{
+ int a_nr = a.rows ();
+ int b_nr = b.length ();
+
+ if (a_nr != b_nr)
+ {
+ int a_nc = a.cols ();
+ int b_nc = 1;
+
+ gripe_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc);
+ return false;
+ }
+
+ return true;
+}
+
+// need to update xdiv.cc ?
+inline ColumnVector
+ldiv (const Matrix &a, const ColumnVector &b)
+{
+ if (! mx_leftdiv_conform (a, b))
+ return ColumnVector ();
+
+ int info;
+ if (a.rows () == a.columns ())
+ {
+ double rcond = 0.0;
+ ColumnVector result = a.solve (b, info, rcond);
+ if (result_ok (info, rcond))
+ return result;
+ }
+
+ int rank;
+ return a.lssolve (b, info, rank);
+}
+
+DEFUN_DLD (mtt_implicit, args, ,
+ "implicit integration method")
+{
+#ifdef OCTAVE_DEV
+ ColumnVector x = args(0).column_vector_value ();
+ const ColumnVector dx = args(1).column_vector_value ();
+ const Matrix AA = args(2).matrix_value ();
+ const ColumnVector AAx = args(3).column_vector_value ();
+ const double t = args(4).double_value ();
+ const int Nx = (int) (args(5).double_value ());
+ const ColumnVector openx = args(6).column_vector_value ();
+#else
+ ColumnVector x = args(0).vector_value ();
+ const ColumnVector dx = args(1).vector_value ();
+ const Matrix AA = args(2).matrix_value ();
+ const ColumnVector AAx = args(3).vector_value ();
+ const double t = args(4).double_value ();
+ const int Nx = (int) (args(5).double_value ());
+ const ColumnVector openx = args(6).vector_value ();
+#endif
+
+ register int i, n;
+ register int col_old, col_new;
+ register int row_old, row_new;
+
+ n = Nx;
+ for (i = 0; i < Nx; i++)
+ {
+ if (0 != openx (i))
+ {
+ n--;
+ }
+ }
+
+ ColumnVector tmp_dx (n, 0.0);
+ ColumnVector tmp_x (n, 0.0);
+ ColumnVector tmp_AAx (n, 0.0);
+ Matrix tmp_AA (n, n, 0.0);
+
+ for (row_new = row_old = 0; row_old < Nx; row_old++)
+ {
+ if (0 == openx (row_old))
+ {
+ tmp_dx (row_new) = dx (row_old);
+ tmp_AAx (row_new) = AAx (row_old);
+ for (col_new = col_old = 0; col_old < Nx; col_old++)
+ {
+ if (0 == openx (col_old))
+ {
+ // xxx: this can be improved by symmetry
+ tmp_AA (row_new,col_new) = AA (row_old,col_old);
+ col_new++;
+ }
+ }
+ row_new++;
+ }
+ }
+
+ // can't get ldiv to work - doesn't like ColVector
+ // tmp_x = tmp_AA.pseudo_inverse () * tmp_AAx + t * tmp_dx;
+ tmp_x = ldiv (tmp_AA, (tmp_AAx + t * tmp_dx));
+
+ row_new = 0;
+ for (row_old = 0; row_old < Nx; row_old++)
+ {
+ if (0 == openx (row_old))
+ {
+ x (row_old) = tmp_x (row_new);
+ row_new++;
+ }
+ else
+ {
+ x (row_old) = 0.0;
+ }
+ }
+
+ return octave_value (x);
+}
+
+
ADDED mttroot/mtt/lib/cc/mtt_write.cc
Index: mttroot/mtt/lib/cc/mtt_write.cc
==================================================================
--- /dev/null
+++ mttroot/mtt/lib/cc/mtt_write.cc
@@ -0,0 +1,54 @@
+#include <octave/oct.h>
+#include <octave/variables.h>
+
+#ifdef STANDALONE
+static Matrix MTT_data;
+#endif
+
+DEFUN_DLD (mtt_write, args, ,
+ "append current data to output")
+{
+ const double t = args(0).double_value ();
+#ifdef OCTAVE_DEV
+ const ColumnVector x = args(1).column_vector_value ();
+ const ColumnVector y = args(2).column_vector_value ();
+#else
+ const ColumnVector x = args(1).vvector_value ();
+ const ColumnVector y = args(2).vector_value ();
+#endif
+
+ const int nx = x.length ();
+ const int ny = y.length ();
+
+
+ ColumnVector Output (2+nx+ny, 0.0);
+
+ Output (0) = Output (1+nx) = t;
+ Output.insert (x.transpose (), 1);
+ Output.insert (y.transpose (), 2+nx);
+
+ Matrix data;
+
+ if (0.0 == t)
+ {
+ data = static_cast<Matrix> (Output.transpose ());
+ }
+ else
+ {
+#ifdef STANDALONE
+ data = MTT_data.transpose ();
+#else
+ data = get_global_value ("MTT_data").matrix_value ().transpose ();
+#endif
+ data = data.append (Output);
+ }
+ data = data.transpose ();
+#ifdef STANDALONE
+ MTT_data = data;
+ cout << Output.transpose () << endl;
+#else
+ set_global_value ("MTT_data", data);
+#endif
+ return data;
+}
+
ADDED mttroot/mtt/lib/rep/standalone_rep.make
Index: mttroot/mtt/lib/rep/standalone_rep.make
==================================================================
--- /dev/null
+++ mttroot/mtt/lib/rep/standalone_rep.make
@@ -0,0 +1,57 @@
+# -*-makefile-*-
+
+.POSIX:
+
+MTTFLAGS = -q -u -oct $(OPTS)
+
+# Adapt according to local set-up and mkoctfile
+CC = g++
+CXXFLAGS = $(DEBUG) $(OPTIM) $(DEFINES) $(ARCHFLAGS) -fno-rtti -fno-exceptions -fno-implicit-templates
+
+DEBUG = -v -pg -g
+OPTIM = -O3
+
+PREFIX = /usr/local
+
+INCLUDES = -I$(PREFIX)/include/octave
+
+OCTAVE_SRC_PATH = $(PREFIX)/src/octave
+
+LIBOCTAVE = -L$(PREFIX)/lib/octave -loctave -lcruft -loctinterp
+LIBKPATHSEA = -L$(OCTAVE_SRC_PATH)/kpathsea -lkpathsea
+LIBREADLINE = -L$(OCTAVE_SRC_PATH)/readline -lreadline
+LIBBLAS = -L/usr/local/src/ATLAS/lib/Linux_PIII -lcblas -lf77blas -llapack -latlas -ltstatlas
+LIBF2C = -lg2c
+LIBRARIES = -ldl -lm -lncurses
+
+ARCHFLAGS = $(i386FLAGS)
+i386FLAGS = -mieee-fp
+
+# Define -DOCTAVE_DEV for octave 2.1.x
+ifeq (0, $(shell octave --version | awk -F\. '{print $2}'))
+DEFINES = -DSTANDALONE
+else
+DEFINES = -DSTANDALONE -DOCTAVE_DEV
+endif
+
+all: $(SYS)_standalone.$(LANG)
+
+$(SYS)_standalone.exe: $(SYS)_ode2odes.cc $(SYS)_def.h $(SYS)_sympar.h
+ cp $(MTT_LIB)/cc/*.cc .
+ $(CC) *.cc -o $@ $(CXXFLAGS) $(INCLUDES) $(LIBOCTAVE) $(LIBKPATHSEA) $(LIBREADLINE) $(LIBBLAS) $(LIBF2C) $(LIBRARIES)
+
+.PHONY: $(SYS)_standalone.clean
+
+$(SYS)_ode2odes.cc:
+ mtt $(MTTFLAGS) $(SYS) ode2odes m
+
+$(SYS)_def.h:
+ mtt $(MTTFLAGS) $(SYS) def h
+
+$(SYS)_sympar.h:
+ mtt $(MTTFLAGS) $(SYS) sympar h
+
+
+$(SYS)_standalone.clean:
+ cd .. ; mtt Clean
+ rm -f ../$(SYS)_standalone.exe