1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
|
/*
* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
* Copyright (C) 2008 - INRIA - Arnaud TORSET
*
* This file must be used under the terms of the CeCILL.
* This source file is licensed as described in the file COPYING, which
* you should have received as part of this distribution. The terms
* are also available at
* http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
*
*/
#include <stdlib.h>
#include "levin.h"
#include "levinUtils.h"
#include "matrixInversion.h"
#include "matrixMultiplication.h"
#include "zeros.h"
void dlevina (int n, double* cov, int lCov, int cCov, double* la, double* sig, double* lb){
/*
[la and lb]
In Scilab, the return value la is a list of n elements. Each element is a matrix cCov*cCov,
and each element of the matrix is a polynome whose degree is n, so the polynome got n+1 elements.
The greatest size of a element of the list is : (n+1)*cCov*cCov.
Here, la is a matrix which contain all elements of the list, its size is n*(n+1)*cCov*cCov. We take the
maximum size for all elements.
The first element of the list is the first (n+1)*cCov*cCov elements of la, namely la[0] to la[(n+1)*cCov*cCov-1].
The second element of the list is the elements of la between (n+1)*cCov*cCov and 2*(n+1)*cCov*cCov,namely la[(n+1)*cCov*cCov]
to la[2*(n+1)*cCov*cCov-1],...
Enter now in a element of the list. Take the first for example.
This is, like said before, a cCov*cCov matrix. In la, it contains (n+1)*cCov*cCov. Each element of the matrix contains (n+1)
elements. As it's stocked by columns, if we represent a matrix like [a c], for example, the elements 0 to n of la represent
[b d]
a, the elements (n+1) to 2(n+1)-1 represent b,...
To finish, look at the elements of the matrix, the polynomes. The coefficients of the polynomes are stocked in increasing order.
For example, if la in Scilab is : list( [3+2x 5-2x ])
( [-5+x -2+4x])
([3+x-x^2 -1+2x ])
([1+6x+3x^2 -2x-x^2 ])
the result in dlevin will be :
-la is a table of 2*3*2*2 elements(n=2,cCov=2);
-la[]={3,2,0, -5,1,0, 5,-2,0, -2,4,0, 3,1,-1, 1,6,3 -1,2,0, 0,-2,-1}.
It's the same for lb.
[sig]
In Scilab, the return value sig is a list of n elements. Each element is a matrix cCov*cCov,
and each element of the matrix is a scalar, so 1 element.
The greatest size of a element of the list is : cCov*cCov.
Let see an example so know how it's stocked.
In Scilab, if sig is : list( [1 3])
( [2 4])
( [5 7])
( [6 8])
the result in dlevin will be :
-sig={1,2, 5,6, 3,4, 7,8}.
It's as if we put the matrix the ones under the others and we take the first column, the second,...
>>>CAREFUL<<<
la/lb and sig are stored differently
*/
int i=0;
/*version ISO C99
double tmp1[n*cCov*cCov], tmp2[n*cCov*cCov];
double sig1[cCov], gam[cCov];
double R1[n*cCov],R2[n*cCov],R3[n*cCov],R4[n*cCov];
*/
/*version pas ISO C99 */
double *tmp1, *tmp2;
double *sig1, *gam;
double *R1,*R2,*R3,*R4;
/* FIXME : malloc here */
tmp1=(double *)malloc((unsigned int)((n+1)*cCov*cCov)*sizeof(double));
tmp2=(double *)malloc((unsigned int)((n+1)*cCov*cCov)*sizeof(double));
sig1=(double *)malloc((unsigned int)(cCov*cCov)*sizeof(double));
gam=(double *)malloc((unsigned int)(cCov*cCov)*sizeof(double));
R1=(double *)malloc((unsigned int)(n*cCov*cCov)*sizeof(double));
R2=(double *)malloc((unsigned int)(n*cCov*cCov)*sizeof(double));
R3=(double *)malloc((unsigned int)(n*cCov*cCov)*sizeof(double));
R4=(double *)malloc((unsigned int)(n*cCov*cCov)*sizeof(double));
/*
* Initializations
* */
dzerosa(sig,n*cCov*cCov,1);
dzerosa(la,n*(n+1)*cCov*cCov,1);
dzerosa(lb,n*(n+1)*cCov*cCov,1);
/*equal to eye(la) and eye(lb)
but we can't use eye cause to the indexation*/
for (i=0;i<cCov;i++){
la[i*((n+1)*(cCov+1))]=1;
lb[i*((n+1)*(cCov+1))]=1;
}
dr1(cov,lCov,cCov,n,R1);
dr2(cov,lCov,cCov,n,R2);
dr3(cov,lCov,cCov,n,R3);
dr4(cov,lCov,cCov,n,R4);
/* case i=0 */
/*computation of sig */
dlevinmul(la,R4,n,cCov,0,sig,n*cCov,0,'d');
/*computation of gam1 */
dlevinmul(lb,R2,n,cCov,0,gam,cCov,0,'u');
/*computation of c1*r1 */
dlevinmul(la,R1,n,cCov,0,tmp1,cCov,0,'u');
/*computation of inv(gam1) */
dinverma(gam,sig1,cCov);
/*computation of k1 = c1*r1*inv(gam1) */
dmulma(tmp1,cCov,cCov,sig1,cCov,cCov,tmp2);
/*computation of k1*lb */
dlevinmul2(tmp2,lb,0,n,cCov,tmp1);
/*computation of k1*lb*z */
ddecalage(tmp1,0,n,cCov,tmp1);
/*computation of la */
dlevinsub(la,tmp1,n,cCov,0,0,la);
/*computation of sig1 (we extract the value if sig at time 0)*/
dlevinsig(sig,0,cCov,n*cCov,sig1);
/*computation of c2*r3 */
dlevinmul(lb,R3,n,cCov,0,tmp1,cCov,0,'d');
/*computation of inv(sig1)*/
dinverma(sig1,gam,cCov);
/*computation of k2 = c2*r3*inv(sig1) */
dmulma(tmp1,cCov,cCov,gam,cCov,cCov,tmp2);
/*computation of k2*la (here it's lb cause la have been modified
and the precedent values hadn't been saved)*/
dlevinmul2(tmp2,lb,0,n,cCov,tmp1);
/*computation of lb*z */
ddecalage(lb,0,n,cCov,lb);
/*computation of lb */
dlevinsub(lb,tmp1,n,cCov,0,0,lb);
for (i=1;i<n;i++){
dlevinmul(la,R4,n,cCov,i,sig,n*cCov,1,'d');
dlevinmul(lb,R2,n,cCov,i,gam,cCov,0,'u');
dlevinmul(la,R1,n,cCov,i,tmp1,cCov,0,'u');
dinverma(gam,sig1,cCov);
dmulma(tmp1,cCov,cCov,sig1,cCov,cCov,tmp2);
dlevinmul2(tmp2,lb,i-1,n,cCov,tmp1);
ddecalage(tmp1,0,n,cCov,tmp1);
dlevinsub(la,tmp1,n,cCov,i,i,la);/*a*/
/*computation of sig1 (we extract the value if sig at time i)*/
dlevinsig(sig,i,cCov,n*cCov,sig1);
dlevinmul(lb,R3,n,cCov,i,tmp1,cCov,0,'d');
dinverma(sig1,gam,cCov);
dmulma(tmp1,cCov,cCov,gam,cCov,cCov,tmp2);
/*computation of k2*la (now it's la at time (i-1))*/
dlevinmul2(tmp2,la,i-1,n,cCov,tmp1);
ddecalage(lb,(i-1)*(n+1)*cCov*cCov,n,cCov,tmp2);
dlevinsub(tmp2,tmp1,n,cCov,0,i,lb);
}
free(R4);
free(R3);
free(R2);
free(R1);
free(gam);
free(sig1);
free(tmp2);
free(tmp1);
}
|
