RiDMC/src/lexp_aux.c

   1 /*
   2 iDMC C library
   3
   4 Adapted from iDMC, Copyright (C) 2004-2006 Marji Lines and Alfredo Medio
   5 Copyright (C) 2006,2007 Marji Lines and Alfredo Medio.
   6
   7 This program is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU General Public License as published by
   9 the Free Software Foundation; either version 2 of the License, or any
  10 later version.
  11
  12 This program is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15 General Public License for more details.
  16
  17 Last modified: $Date$
  18 */
  19 #include <string.h>
  20 #include <math.h>
  21 #include <gsl/gsl_block.h>
  22 #include <gsl/gsl_vector.h>
  23 #include <gsl/gsl_matrix.h>
  24 #include <gsl/gsl_blas.h>
  25 #include <gsl/gsl_linalg.h>
  26 #include <gsl/gsl_odeiv.h>
  27 #include "lexp_aux.h"
  28 #include "matrices.h"
  29 #include "model.h"
  30
  31 struct combined{
  32         idmc_model* model;
  33         int dim;
  34         double* pars;
  35         double* data;
  36         double* draft1;
  37         double* draft2;
  38         double* draft3;
  39 };
  40
  41 typedef struct combined Combined;
  42
  43 int func(double t,const double y[],double f[],void* params){
  44         int status;
  45         Combined* c=(Combined*) params;
  46         status=idmc_model_f((*c).model,(*c).pars,y,f);
  47         if (status==IDMC_OK)
  48                 return GSL_SUCCESS;
  49         else
  50                 return GSL_SUCCESS-1;
  51 }
  52
  53 inline int func_jac(double t,double* y,double* f,Combined* params){
  54         int status;
  55         Combined* par=(Combined*) params;
  56         status=idmc_model_Jf((*par).model,(*par).pars,y,f);
  57         if (status==IDMC_OK)
  58                 return GSL_SUCCESS;
  59         else
  60                 return GSL_SUCCESS-1;
  61 }
  62
  63 //(*c).draft should hold dim*dim
  64 inline int interpolated_jac(Combined* c, double t0, double* y0, double t1, double* y1, double t, double f[])
  65 {
  66         int status1,status2;
  67         int dim=(*c).dim;
  68         double* draft=(*c).draft1;//(double*) malloc(dim*dim*sizeof(double));
  69         //assumes that y0 y1 hold dim elements each
  70
  71         double a;
  72         if (t1==t0)
  73                 a=0;
  74         else
  75                 a=(t-t0)/(t1-t0);
  76         status1=func_jac(t0,y0,f,c);
  77         status2=func_jac(t1,y1,draft,c);
  78
  79         gsl_matrix_view m1v = gsl_matrix_view_array(draft,dim,dim);
  80         gsl_matrix_view m0v = gsl_matrix_view_array(f,dim,dim);
  81         gsl_matrix_scale(&m1v.matrix,a);
  82         gsl_matrix_scale(&m0v.matrix,1-a);
  83         gsl_matrix_add(&m0v.matrix,&m1v.matrix);
  84         if (status1==GSL_SUCCESS && status2==GSL_SUCCESS)
  85                 return GSL_SUCCESS;
  86         else
  87                 return GSL_SUCCESS-1;
  88 }
  89
  90 // (*c).data is 2dim+2, (*c).draft should hold 5*dim*dim (sizeof(double))
  91 int Qsystem (double t, const double Q[], double RHS[], void *params)
  92 {
  93         int status;
  94         Combined* c=(Combined*) params;
  95         int dim=(*c).dim;
  96         int dd=dim*dim;
  97
  98         double* ptr=(*c).data;
  99         double tt0=*ptr;
 100         double* yy0=ptr+1;
 101         double tt1=*(ptr+dim+1);
 102         double* yy1=ptr+(dim+2);
 103
 104
 105         double* draft=(*c).draft2;//(double*) malloc(4*dim*dim*sizeof(double));
 106         double* Jac=draft;
 107
 108         status=interpolated_jac(c,tt0, yy0, tt1, yy1, t, Jac);
 109
 110         double* Qt=Jac+dd;
 111         transpose(Q,Qt,dim,dim);
 112
 113         double* tmp1=Qt+dd;
 114         mmproduct(Jac,Q,tmp1,dim,dim,dim);
 115         double* tmp2=tmp1+dd;
 116         mmproduct(Qt,tmp1,tmp2,dim,dim,dim);
 117         asym(tmp2,dim);
 118         mmproduct(Q,tmp2,RHS,dim,dim,dim);
 119
 120         if (status==GSL_SUCCESS)
 121                 return GSL_SUCCESS;
 122         else
 123                 return GSL_SUCCESS-1;
 124 }
 125
 126 inline int elementary_step(Combined* c, double step, gsl_odeiv_step* ode_step, gsl_odeiv_step* Q_step, gsl_odeiv_system* sys,
 127                                                    double* t, double* y, double* Q, double* l, double* y_err, double* Q_err ){
 128         int status1,status2,status3;
 129         int dim=(*c).dim;
 130         int dd=dim*dim;
 131
 132         double* Jac;
 133         double* Qt;
 134         double* QtJ;
 135         double* QtJQ;
 136         double* R;
 137         double* Qcopy;
 138         double* tmp1;
 139
 140         *((*c).data)=*t;
 141         memcpy((*c).data+1,y,dim*sizeof(double));
 142
 143
 144         status1 = gsl_odeiv_step_apply (ode_step, *t, step, y, y_err, NULL, NULL, sys);
 145
 146         *((*c).data+dim+1)=*t+step;
 147         memcpy((*c).data+dim+2,y,dim*sizeof(double));
 148
 149         gsl_odeiv_system Q_sys={Qsystem, NULL,dim*dim, c};
 150
 151         status2 = gsl_odeiv_step_apply (Q_step, *t, step,  Q, Q_err, NULL, NULL, &Q_sys);
 152         *t=*t+step;
 153
 154         double* draft=(*c).draft3;//(double*) malloc((6*dd+dim)*sizeof(double));
 155         R=draft;
 156         tmp1=R+dd;
 157         Qcopy=tmp1+dim;
 158         Qt=Qcopy+dd;
 159         Jac=Qt+dd;
 160         QtJ=Jac+dd;
 161         QtJQ=QtJ+dd;
 162
 163         memcpy(Qcopy,Q,dd*sizeof(double));
 164         qr_coded( Qcopy, tmp1, dim,dim);
 165         qr_decomp(Qcopy, tmp1, Q, R,dim,dim);
 166
 167         status3=idmc_model_Jf((*c).model,(*c).pars,y,Jac);
 168         transpose(Q,Qt,dim,dim);
 169         mmproduct(Qt,Jac,QtJ,dim,dim,dim);
 170         mmproduct(QtJ,Q,QtJQ,dim,dim,dim);
 171         gsl_vector_view lv=gsl_vector_view_array(l,dim);
 172         gsl_matrix_view QtJQv=gsl_matrix_view_array(QtJQ,dim,dim);
 173         gsl_vector_view diagv=gsl_matrix_diagonal(&QtJQv.matrix);
 174         gsl_vector_scale(&diagv.vector,step);
 175         gsl_vector_add(&lv.vector,&diagv.vector);
 176
 177         if (status1==GSL_SUCCESS && status2==GSL_SUCCESS && status3==IDMC_OK)
 178                 return GSL_SUCCESS;
 179         else
 180                 return GSL_SUCCESS-1;
 181 }
 182
 183
 184 inline int time_plot_step(idmc_model* model, int dim, double step, double* t, double* pars,
 185                                                   double* y, double* Q,  double* l, double* alloc_memory ){
 186         int s,status;
 187         status=IDMC_OK;
 188
 189         int dd=dim*dim;
 190
 191         const  gsl_odeiv_step_type* step_type= gsl_odeiv_step_rk2;
 192
 193         gsl_odeiv_step* ode_step=gsl_odeiv_step_alloc (step_type, dim);
 194         gsl_odeiv_step* Q_step=gsl_odeiv_step_alloc (step_type, dim*dim);
 195
 196         gsl_vector_view lv=gsl_vector_view_array(l,dim);
 197         if (t==0)
 198                 gsl_vector_set_zero(&lv.vector);
 199
 200         double* work_space=alloc_memory;
 201         double* y_err=work_space;
 202
 203         double* Q_err=y_err+dim;
 204         double* point_pair=Q_err+dd;
 205
 206         double* draft1=alloc_memory+3*dim+2+2*dd;
 207         double* draft2=draft1+dd;
 208         double* draft3=draft2+4*dd;
 209         gsl_matrix_view Qv=gsl_matrix_view_array(Q,dim,dim);
 210         if (t==0)
 211                 gsl_matrix_set_identity(&Qv.matrix);
 212
 213         Combined c;
 214         c.model=model;
 215         c.dim=dim;
 216         c.pars=pars;
 217         c.data=point_pair;
 218         c.draft1=draft1;
 219         c.draft2=draft2;
 220         c.draft3=draft3;
 221         gsl_odeiv_system sys = {func, NULL, dim, &c};
 222
 223         s=elementary_step(&c, step, ode_step, Q_step, &sys, t, y, Q, l, y_err, Q_err);
 224
 225         if (s!=GSL_SUCCESS)
 226                 status=IDMC_EMATH;
 227
 228         gsl_odeiv_step_free(ode_step);
 229         gsl_odeiv_step_free(Q_step);
 230         return(status);
 231 }
 232
 233 //y0=y, t0=0, t1 is the final time, l is the array containing the Lyap.exps.
 234 inline int compute_lexp(idmc_model* model, double* par, int dim,
 235                                                 double step, double* y, double t1, double* l, double* alloc_memory){
 236         int s,status;
 237         status=IDMC_OK;
 238
 239         int dd=dim*dim;
 240
 241         const  gsl_odeiv_step_type* step_type= gsl_odeiv_step_rk2;
 242
 243         gsl_odeiv_step* ode_step=gsl_odeiv_step_alloc (step_type, dim);
 244         gsl_odeiv_step* Q_step=gsl_odeiv_step_alloc (step_type, dim*dim);
 245
 246         double t = 0;
 247
 248         gsl_vector_view lv=gsl_vector_view_array(l,dim);
 249         gsl_vector_set_zero(&lv.vector);
 250
 251         double* work_space=alloc_memory;
 252         double* y_err=work_space;
 253
 254         double* Q=y_err+dim;
 255         double* Q_err=Q+dd;
 256         double* point_pair=Q_err+dd;
 257
 258         double* draft1=alloc_memory+3*dim+2+2*dd;
 259         double* draft2=draft1+dd;
 260         double* draft3=draft2+4*dd;
 261
 262         gsl_matrix_view Qv=gsl_matrix_view_array(Q,dim,dim);
 263         gsl_matrix_set_identity(&Qv.matrix);
 264
 265         Combined c;
 266         c.model=model;
 267         c.dim=dim;
 268         c.pars=par;
 269         c.data=point_pair;
 270         c.draft1=draft1;
 271         c.draft2=draft2;
 272         c.draft3=draft3;
 273         gsl_odeiv_system sys = {func, NULL, dim, &c};
 274
 275         while (t < t1){
 276                 s = elementary_step(&c, step, ode_step, Q_step, &sys, &t,
 277                                                          y, Q, l, y_err, Q_err);
 278                 if (s!=GSL_SUCCESS)
 279                         status=IDMC_OK-1;
 280         }
 281
 282         gsl_vector_scale(&lv.vector,1/(t1+step));//if t1 is set to zero, does not give an error
 283         gsl_odeiv_step_free (ode_step);
 284         gsl_odeiv_step_free (Q_step);
 285         return(status);
 286 }
 287