summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
Diffstat (limited to 'packages')
-rw-r--r--packages/base/src/Data/Packed/Matrix.hs17
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Data.hs5
-rw-r--r--packages/gsl/src/Numeric/GSL/ODE.hs126
-rw-r--r--packages/gsl/src/Numeric/GSL/gsl-ode.c31
4 files changed, 121 insertions, 58 deletions
diff --git a/packages/base/src/Data/Packed/Matrix.hs b/packages/base/src/Data/Packed/Matrix.hs
index 70b9232..3690293 100644
--- a/packages/base/src/Data/Packed/Matrix.hs
+++ b/packages/base/src/Data/Packed/Matrix.hs
@@ -33,8 +33,9 @@ module Data.Packed.Matrix (
33 fromRows, toRows, fromColumns, toColumns, 33 fromRows, toRows, fromColumns, toColumns,
34 fromBlocks, diagBlock, toBlocks, toBlocksEvery, 34 fromBlocks, diagBlock, toBlocks, toBlocksEvery,
35 repmat, 35 repmat,
36 flipud, fliprl, 36 flipud, fliprl, subMatrix,
37 subMatrix, takeRows, dropRows, takeColumns, dropColumns, 37 takeRows, takeLastRows, dropRows, dropLastRows,
38 takeColumns, takeLastColumns, dropColumns, dropLastColumns,
38 extractRows, extractColumns, 39 extractRows, extractColumns,
39 diagRect, takeDiag, 40 diagRect, takeDiag,
40 mapMatrix, mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_, 41 mapMatrix, mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_,
@@ -253,15 +254,27 @@ r >< c = f where
253-- | Creates a matrix with the first n rows of another matrix 254-- | Creates a matrix with the first n rows of another matrix
254takeRows :: Element t => Int -> Matrix t -> Matrix t 255takeRows :: Element t => Int -> Matrix t -> Matrix t
255takeRows n mt = subMatrix (0,0) (n, cols mt) mt 256takeRows n mt = subMatrix (0,0) (n, cols mt) mt
257-- | Creates a matrix with the last n rows of another matrix
258takeLastRows :: Element t => Int -> Matrix t -> Matrix t
259takeLastRows n mt = subMatrix (rows mt - n, 0) (n, cols mt) mt
256-- | Creates a copy of a matrix without the first n rows 260-- | Creates a copy of a matrix without the first n rows
257dropRows :: Element t => Int -> Matrix t -> Matrix t 261dropRows :: Element t => Int -> Matrix t -> Matrix t
258dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt 262dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt
263-- | Creates a copy of a matrix without the last n rows
264dropLastRows :: Element t => Int -> Matrix t -> Matrix t
265dropLastRows n mt = subMatrix (0,0) (rows mt - n, cols mt) mt
259-- |Creates a matrix with the first n columns of another matrix 266-- |Creates a matrix with the first n columns of another matrix
260takeColumns :: Element t => Int -> Matrix t -> Matrix t 267takeColumns :: Element t => Int -> Matrix t -> Matrix t
261takeColumns n mt = subMatrix (0,0) (rows mt, n) mt 268takeColumns n mt = subMatrix (0,0) (rows mt, n) mt
269-- |Creates a matrix with the last n columns of another matrix
270takeLastColumns :: Element t => Int -> Matrix t -> Matrix t
271takeLastColumns n mt = subMatrix (0, cols mt - n) (rows mt, n) mt
262-- | Creates a copy of a matrix without the first n columns 272-- | Creates a copy of a matrix without the first n columns
263dropColumns :: Element t => Int -> Matrix t -> Matrix t 273dropColumns :: Element t => Int -> Matrix t -> Matrix t
264dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt 274dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt
275-- | Creates a copy of a matrix without the last n columns
276dropLastColumns :: Element t => Int -> Matrix t -> Matrix t
277dropLastColumns n mt = subMatrix (0,0) (rows mt, cols mt - n) mt
265 278
266---------------------------------------------------------------- 279----------------------------------------------------------------
267 280
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index 9935c15..2161e75 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -53,8 +53,9 @@ module Numeric.LinearAlgebra.Data(
53 53
54 -- * Matrix extraction 54 -- * Matrix extraction
55 Extractor(..), (??), 55 Extractor(..), (??),
56 takeRows, dropRows, takeColumns, dropColumns, subMatrix, (?), (¿), fliprl, flipud, 56 takeRows, takeLastRows, dropRows, dropLastRows,
57 remap, 57 takeColumns, takeLastColumns, dropColumns, dropLastColumns,
58 subMatrix, (?), (¿), fliprl, flipud, remap,
58 59
59 -- * Block matrix 60 -- * Block matrix
60 fromBlocks, (|||), (===), diagBlock, repmat, toBlocks, toBlocksEvery, 61 fromBlocks, (|||), (===), diagBlock, repmat, toBlocks, toBlocksEvery,
diff --git a/packages/gsl/src/Numeric/GSL/ODE.hs b/packages/gsl/src/Numeric/GSL/ODE.hs
index 3258b83..9e52873 100644
--- a/packages/gsl/src/Numeric/GSL/ODE.hs
+++ b/packages/gsl/src/Numeric/GSL/ODE.hs
@@ -32,7 +32,7 @@ main = mplot (ts : toColumns sol)
32----------------------------------------------------------------------------- 32-----------------------------------------------------------------------------
33 33
34module Numeric.GSL.ODE ( 34module Numeric.GSL.ODE (
35 odeSolve, odeSolveV, ODEMethod(..), Jacobian 35 odeSolve, odeSolveV, odeSolveVWith, ODEMethod(..), Jacobian, StepControl(..)
36) where 36) where
37 37
38import Numeric.LinearAlgebra.HMatrix 38import Numeric.LinearAlgebra.HMatrix
@@ -44,9 +44,10 @@ import System.IO.Unsafe(unsafePerformIO)
44 44
45------------------------------------------------------------------------- 45-------------------------------------------------------------------------
46 46
47type TVV = TV (TV Res) 47type TVV = TV (TV Res)
48type TVM = TV (TM Res) 48type TVM = TV (TM Res)
49type TVVM = TV (TV (TM Res)) 49type TVVM = TV (TV (TM Res))
50type TVVVM = TV (TV (TV (TM Res)))
50 51
51type Jacobian = Double -> Vector Double -> Matrix Double 52type Jacobian = Double -> Vector Double -> Matrix Double
52 53
@@ -63,68 +64,100 @@ data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
63 | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12. 64 | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12.
64 | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems. 65 | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems.
65 66
67-- | Adaptive step-size control functions
68data StepControl = X Double Double -- ^ abs. and rel. tolerance for x(t)
69 | X' Double Double -- ^ abs. and rel. tolerance for x'(t)
70 | XX' Double Double Double Double -- ^ include both via rel. tolerance scaling factors a_x, a_x'
71 | ScXX' Double Double Double Double (Vector Double) -- ^ scale abs. tolerance of x(t) components
66 72
67-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists. 73-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
68odeSolve 74odeSolve
69 :: (Double -> [Double] -> [Double]) -- ^ xdot(t,x) 75 :: (Double -> [Double] -> [Double]) -- ^ x'(t,x)
70 -> [Double] -- ^ initial conditions 76 -> [Double] -- ^ initial conditions
71 -> Vector Double -- ^ desired solution times 77 -> Vector Double -- ^ desired solution times
72 -> Matrix Double -- ^ solution 78 -> Matrix Double -- ^ solution
73odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts 79odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
74 where hi = (ts!1 - ts!0)/100 80 where hi = (ts!1 - ts!0)/100
75 epsAbs = 1.49012e-08 81 epsAbs = 1.49012e-08
76 epsRel = 1.49012e-08 82 epsRel = epsAbs
77 l2v f = \t -> fromList . f t . toList 83 l2v f = \t -> fromList . f t . toList
78 84
79-- | Evolution of the system with adaptive step-size control. 85-- | A version of 'odeSolveVWith' with reasonable default step control.
80odeSolveV 86odeSolveV
81 :: ODEMethod 87 :: ODEMethod
82 -> Double -- ^ initial step size 88 -> Double -- ^ initial step size
83 -> Double -- ^ absolute tolerance for the state vector 89 -> Double -- ^ absolute tolerance for the state vector
84 -> Double -- ^ relative tolerance for the state vector 90 -> Double -- ^ relative tolerance for the state vector
85 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) 91 -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x)
86 -> Vector Double -- ^ initial conditions 92 -> Vector Double -- ^ initial conditions
87 -> Vector Double -- ^ desired solution times 93 -> Vector Double -- ^ desired solution times
88 -> Matrix Double -- ^ solution 94 -> Matrix Double -- ^ solution
89odeSolveV RK2 = odeSolveV' 0 Nothing 95odeSolveV meth hi epsAbs epsRel = odeSolveVWith meth (XX' epsAbs epsRel 1 1) hi
90odeSolveV RK4 = odeSolveV' 1 Nothing 96
91odeSolveV RKf45 = odeSolveV' 2 Nothing 97-- | Evolution of the system with adaptive step-size control.
92odeSolveV RKck = odeSolveV' 3 Nothing 98odeSolveVWith
93odeSolveV RK8pd = odeSolveV' 4 Nothing 99 :: ODEMethod
94odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac) 100 -> StepControl
95odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac) 101 -> Double -- ^ initial step size
96odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac) 102 -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x)
97odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac)
98odeSolveV MSAdams = odeSolveV' 9 Nothing
99odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac)
100
101
102odeSolveV'
103 :: CInt
104 -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian
105 -> Double -- ^ initial step size
106 -> Double -- ^ absolute tolerance for the state vector
107 -> Double -- ^ relative tolerance for the state vector
108 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x)
109 -> Vector Double -- ^ initial conditions 103 -> Vector Double -- ^ initial conditions
110 -> Vector Double -- ^ desired solution times 104 -> Vector Double -- ^ desired solution times
111 -> Matrix Double -- ^ solution 105 -> Matrix Double -- ^ solution
112odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do 106odeSolveVWith method control = odeSolveVWith' m mbj c epsAbs epsRel aX aX' mbsc
113 let n = size xiv 107 where (m, mbj) = case method of
114 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) 108 RK2 -> (0 , Nothing )
115 jp <- case mbjac of 109 RK4 -> (1 , Nothing )
116 Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) 110 RKf45 -> (2 , Nothing )
117 Nothing -> return nullFunPtr 111 RKck -> (3 , Nothing )
118 sol <- vec xiv $ \xiv' -> 112 RK8pd -> (4 , Nothing )
119 vec (checkTimes ts) $ \ts' -> 113 RK2imp jac -> (5 , Just jac)
120 createMIO (size ts) n 114 RK4imp jac -> (6 , Just jac)
121 (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' ) 115 BSimp jac -> (7 , Just jac)
122 "ode" 116 RK1imp jac -> (8 , Just jac)
123 freeHaskellFunPtr fp 117 MSAdams -> (9 , Nothing )
124 return sol 118 MSBDF jac -> (10, Just jac)
119 (c, epsAbs, epsRel, aX, aX', mbsc) = case control of
120 X ea er -> (0, ea, er, 1 , 0 , Nothing)
121 X' ea er -> (0, ea, er, 0 , 1 , Nothing)
122 XX' ea er ax ax' -> (0, ea, er, ax, ax', Nothing)
123 ScXX' ea er ax ax' sc -> (1, ea, er, ax, ax', Just sc)
124
125odeSolveVWith'
126 :: CInt -- ^ stepping function
127 -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian
128 -> CInt -- ^ step-size control function
129 -> Double -- ^ absolute tolerance for step-size control
130 -> Double -- ^ relative tolerance for step-size control
131 -> Double -- ^ scaling factor for relative tolerance of x(t)
132 -> Double -- ^ scaling factor for relative tolerance of x'(t)
133 -> Maybe (Vector Double) -- ^ optional scaling for absolute error
134 -> Double -- ^ initial step size
135 -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x)
136 -> Vector Double -- ^ initial conditions
137 -> Vector Double -- ^ desired solution times
138 -> Matrix Double -- ^ solution
139odeSolveVWith' method mbjac control epsAbs epsRel aX aX' mbsc h f xiv ts =
140 unsafePerformIO $ do
141 let n = size xiv
142 sc = case mbsc of
143 Just scv -> checkdim1 n scv
144 Nothing -> xiv
145 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
146 jp <- case mbjac of
147 Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
148 Nothing -> return nullFunPtr
149 sol <- vec sc $ \sc' -> vec xiv $ \xiv' ->
150 vec (checkTimes ts) $ \ts' -> createMIO (size ts) n
151 (ode_c method control h epsAbs epsRel aX aX' fp jp
152 // sc' // xiv' // ts' )
153 "ode"
154 freeHaskellFunPtr fp
155 return sol
125 156
126foreign import ccall safe "ode" 157foreign import ccall safe "ode"
127 ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM 158 ode_c :: CInt -> CInt -> Double
159 -> Double -> Double -> Double -> Double
160 -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVVM
128 161
129------------------------------------------------------- 162-------------------------------------------------------
130 163
@@ -141,4 +174,3 @@ checkdim2 n m
141checkTimes ts | size ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts 174checkTimes ts | size ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
142 | otherwise = error "odeSolve requires increasing times" 175 | otherwise = error "odeSolve requires increasing times"
143 where ts' = toList ts 176 where ts' = toList ts
144
diff --git a/packages/gsl/src/Numeric/GSL/gsl-ode.c b/packages/gsl/src/Numeric/GSL/gsl-ode.c
index 3f2771b..a6bdb55 100644
--- a/packages/gsl/src/Numeric/GSL/gsl-ode.c
+++ b/packages/gsl/src/Numeric/GSL/gsl-ode.c
@@ -23,10 +23,11 @@ int odejac (double t, const double y[], double *dfdy, double dfdt[], void *param
23} 23}
24 24
25 25
26int ode(int method, double h, double eps_abs, double eps_rel, 26int ode(int method, int control, double h,
27 double eps_abs, double eps_rel, double a_y, double a_dydt,
27 int f(double, int, const double*, int, double*), 28 int f(double, int, const double*, int, double*),
28 int jac(double, int, const double*, int, int, double*), 29 int jac(double, int, const double*, int, int, double*),
29 KRVEC(xi), KRVEC(ts), RMAT(sol)) { 30 KRVEC(sc), KRVEC(xi), KRVEC(ts), RMAT(sol)) {
30 31
31 const gsl_odeiv_step_type * T; 32 const gsl_odeiv_step_type * T;
32 33
@@ -46,8 +47,16 @@ int ode(int method, double h, double eps_abs, double eps_rel,
46 } 47 }
47 48
48 gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); 49 gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);
49 gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);
50 gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); 50 gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);
51 gsl_odeiv_control * c;
52
53 switch(control) {
54 case 0: { c = gsl_odeiv_control_standard_new
55 (eps_abs, eps_rel, a_y, a_dydt); break; }
56 case 1: { c = gsl_odeiv_control_scaled_new
57 (eps_abs, eps_rel, a_y, a_dydt, scp, scn); break; }
58 default: ERROR(BAD_CODE);
59 }
51 60
52 Tode P; 61 Tode P;
53 P.f = f; 62 P.f = f;
@@ -112,10 +121,11 @@ int odejac (double t, const double y[], double *dfdy, double dfdt[], void *param
112} 121}
113 122
114 123
115int ode(int method, double h, double eps_abs, double eps_rel, 124int ode(int method, int control, double h,
125 double eps_abs, double eps_rel, double a_y, double a_dydt,
116 int f(double, int, const double*, int, double*), 126 int f(double, int, const double*, int, double*),
117 int jac(double, int, const double*, int, int, double*), 127 int jac(double, int, const double*, int, int, double*),
118 KRVEC(xi), KRVEC(ts), RMAT(sol)) { 128 KRVEC(sc), KRVEC(xi), KRVEC(ts), RMAT(sol)) {
119 129
120 const gsl_odeiv2_step_type * T; 130 const gsl_odeiv2_step_type * T;
121 131
@@ -141,8 +151,15 @@ int ode(int method, double h, double eps_abs, double eps_rel,
141 151
142 gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; 152 gsl_odeiv2_system sys = {odefunc, odejac, xin, &P};
143 153
144 gsl_odeiv2_driver * d = 154 gsl_odeiv2_driver * d;
145 gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); 155
156 switch(control) {
157 case 0: { d = gsl_odeiv2_driver_alloc_standard_new
158 (&sys, T, h, eps_abs, eps_rel, a_y, a_dydt); break; }
159 case 1: { d = gsl_odeiv2_driver_alloc_scaled_new
160 (&sys, T, h, eps_abs, eps_rel, a_y, a_dydt, scp); break; }
161 default: ERROR(BAD_CODE);
162 }
146 163
147 double t = tsp[0]; 164 double t = tsp[0];
148 165