| 1 |
// $Header: $ |
| 2 |
|
| 3 |
#include "ctypes_def.h" |
| 4 |
|
| 5 |
#if TYPE == FLOAT |
| 6 |
#define GAPSTRUCTURE sgapstructure |
| 7 |
#define BURG smulti_burg |
| 8 |
#define FAHLMAN_ULRYCH sfahlman_ulrych |
| 9 |
#define FILL_GAPS sfill_gaps |
| 10 |
#define TOEPLITZ_MULT stoeplitz_mult |
| 11 |
#define PRECONDITION sprecondition |
| 12 |
#define PCG spcg |
| 13 |
#define LEVINSON slevinson |
| 14 |
#define EPSILON FLT_EPSILON |
| 15 |
#elif TYPE == DOUBLE |
| 16 |
#define GAPSTRUCTURE dgapstructure |
| 17 |
#define BURG dmulti_burg |
| 18 |
#define FAHLMAN_ULRYCH dfahlman_ulrych |
| 19 |
#define FILL_GAPS dfill_gaps |
| 20 |
#define TOEPLITZ_MULT dtoeplitz_mult |
| 21 |
#define PRECONDITION dprecondition |
| 22 |
#define PCG dpcg |
| 23 |
#define LEVINSON dlevinson |
| 24 |
#define EPSILON DBL_EPSILON |
| 25 |
#elif TYPE == COMPLEXFLOAT |
| 26 |
#define GAPSTRUCTURE cgapstructure |
| 27 |
#define BURG cmulti_burg |
| 28 |
#define FAHLMAN_ULRYCH cfahlman_ulrych |
| 29 |
#define FILL_GAPS cfill_gaps |
| 30 |
#define TOEPLITZ_MULT ctoeplitz_mult |
| 31 |
#define PRECONDITION cprecondition |
| 32 |
#define PCG cpcg |
| 33 |
#define LEVINSON clevinson |
| 34 |
#define EPSILON FLT_EPSILON |
| 35 |
#elif TYPE == COMPLEXDOUBLE |
| 36 |
#define GAPSTRUCTURE zgapstructure |
| 37 |
#define BURG zmulti_burg |
| 38 |
#define FAHLMAN_ULRYCH zfahlman_ulrych |
| 39 |
#define FILL_GAPS zfill_gaps |
| 40 |
#define TOEPLITZ_MULT ztoeplitz_mult |
| 41 |
#define PRECONDITION zprecondition |
| 42 |
#define PCG zpcg |
| 43 |
#define LEVINSON zlevinson |
| 44 |
#define EPSILON DBL_EPSILON |
| 45 |
#endif |
| 46 |
|
| 47 |
#define PCG_MAXIT n_gap |
| 48 |
#define PCG_TOL 1e-6 |
| 49 |
|
| 50 |
static int GAPSTRUCTURE(int n, CTYPE *data, int *isgood, |
| 51 |
gapped_timeseries *ts, int **data_length, |
| 52 |
CTYPE ***good_data, int minpercentage); |
| 53 |
static void FILL_GAPS(int n, CTYPE *data, int *isgood, int order, CTYPE *ar_coeff); |
| 54 |
static void TOEPLITZ_MULT(int n, CTYPE *x, CTYPE *y, void **data); |
| 55 |
static void PRECONDITION(int n, CTYPE *x, CTYPE *y, void **data); |
| 56 |
|
| 57 |
//static void identity(int n, CTYPE *x, CTYPE *y, void **data); |
| 58 |
|
| 59 |
|
| 60 |
|
| 61 |
int FAHLMAN_ULRYCH(int n, CTYPE *data_in, int *isgood_in, int minpercentage, |
| 62 |
int maxorder, int iterations, int padends, |
| 63 |
int *order, CTYPE *ar_coeff_in) |
| 64 |
{ |
| 65 |
int i,j, first, last, iter, effective_order; |
| 66 |
RTYPE *err ; |
| 67 |
gapped_timeseries ts,tsm; |
| 68 |
CTYPE **good_data; |
| 69 |
int *isgood; |
| 70 |
CTYPE *data, nold; |
| 71 |
int *data_length; |
| 72 |
int *model_gaps, oldorder; |
| 73 |
int model_order; |
| 74 |
CTYPE *ar_coeff; |
| 75 |
|
| 76 |
|
| 77 |
if (ar_coeff_in) |
| 78 |
ar_coeff = ar_coeff_in; |
| 79 |
else |
| 80 |
ar_coeff = malloc(maxorder*sizeof(CTYPE)); |
| 81 |
|
| 82 |
// the following check will be skipped in production executables since they will be built with NDEBUG defined. |
| 83 |
// it doesn't matter because the check is already done in jtsfiddle.c which calls this function. |
| 84 |
assert(minpercentage>0 && minpercentage<=100); |
| 85 |
// Execute the Fahlman-Ulrych algorithm: |
| 86 |
// FOR i=1 TO iterations DO |
| 87 |
// Step 1. Compute AR model using known sample values. |
| 88 |
// Step 2. For the unknown samples, compute predicted values that |
| 89 |
// minimize the forward and backward prediction errors in |
| 90 |
// the least squares sense. |
| 91 |
// END FOR |
| 92 |
|
| 93 |
if (padends) |
| 94 |
{ |
| 95 |
nold = n; |
| 96 |
n = nold + 2*maxorder; |
| 97 |
data = calloc(n,sizeof(CTYPE)); |
| 98 |
isgood = calloc(n,sizeof(int)); |
| 99 |
memcpy(&data[maxorder],data_in,nold*sizeof(CTYPE)); |
| 100 |
memcpy(&isgood[maxorder],isgood_in,nold*sizeof(int)); |
| 101 |
} |
| 102 |
else |
| 103 |
{ |
| 104 |
isgood = isgood_in; |
| 105 |
data = data_in; |
| 106 |
} |
| 107 |
|
| 108 |
oldorder = -1; |
| 109 |
model_order = GAPSTRUCTURE(n, data, isgood, &ts, &data_length, |
| 110 |
&good_data, minpercentage); |
| 111 |
tsm = ts; |
| 112 |
model_gaps = malloc(n*sizeof(int)); |
| 113 |
memcpy(model_gaps, isgood, n*sizeof(int)); |
| 114 |
for (iter=0; iter<iterations; iter++) |
| 115 |
{ |
| 116 |
// Analyze gapstructure. |
| 117 |
if (iter>0) |
| 118 |
model_order = GAPSTRUCTURE(n, data, model_gaps, &tsm, &data_length, |
| 119 |
&good_data, minpercentage); |
| 120 |
|
| 121 |
|
| 122 |
model_order = min( model_order,maxorder); |
| 123 |
|
| 124 |
// Sort gaps by length and determine the largest AR model order |
| 125 |
// such that "minpercentage" percent of the data will be included |
| 126 |
// in the estimation of the AR coefficients. |
| 127 |
// model_order = maxorder(&ts,data_length,minpercentage); |
| 128 |
err = (RTYPE *) malloc((model_order+1)*sizeof(RTYPE)); |
| 129 |
|
| 130 |
#ifdef DEBUGOUT |
| 131 |
{ |
| 132 |
printf("Determining AR(%d) model...",model_order); |
| 133 |
fflush(stdout); |
| 134 |
} |
| 135 |
#endif |
| 136 |
|
| 137 |
effective_order = BURG(tsm.m_data, data_length, good_data, model_order, |
| 138 |
ar_coeff, err); |
| 139 |
|
| 140 |
#ifdef DEBUGOUT |
| 141 |
{ |
| 142 |
printf("done\n"); |
| 143 |
printf("Specified order = %d, Effective order = %d.\n", model_order, effective_order); |
| 144 |
printf("Filling gaps...\n"); |
| 145 |
fflush(stdout); |
| 146 |
} |
| 147 |
#endif |
| 148 |
|
| 149 |
model_order = effective_order; |
| 150 |
if (ts.first_is_good) |
| 151 |
{ |
| 152 |
i = 1; |
| 153 |
first = max(0,ts.data_int[0].last-(model_order)+1); |
| 154 |
last=-1; |
| 155 |
} |
| 156 |
else |
| 157 |
{ |
| 158 |
i = 0; |
| 159 |
first = 0; |
| 160 |
last=-1; |
| 161 |
} |
| 162 |
memcpy(model_gaps, isgood, n*sizeof(int)); |
| 163 |
while ( i<ts.m_data && last<n-1) |
| 164 |
{ |
| 165 |
// move the right-hand end of the interval until we find a |
| 166 |
// data section longer than the model order, i.e. which decouples |
| 167 |
// the unknowns in the gaps on either side. |
| 168 |
while (i<ts.m_data && |
| 169 |
(ts.data_int[i].last-ts.data_int[i].first+1)<(model_order) && |
| 170 |
(i==0 || (ts.data_int[i].first-ts.data_int[i-1].last-1)<model_order) |
| 171 |
) |
| 172 |
{ i++; } |
| 173 |
|
| 174 |
|
| 175 |
if (i!=0 && (ts.data_int[i].first-ts.data_int[i-1].last)>=model_order) |
| 176 |
last = ts.data_int[i-1].last+(model_order)-1; |
| 177 |
else if ( i>=ts.m_data-1 ) |
| 178 |
last = n-1; |
| 179 |
else |
| 180 |
last = ts.data_int[i].first+(model_order)-1; |
| 181 |
|
| 182 |
#ifdef DEBUGOUT |
| 183 |
printf("Filling missing data in [ %d : %d ].\n",first,last); |
| 184 |
#endif |
| 185 |
|
| 186 |
FILL_GAPS(last-first+1, &data[first], &isgood[first], |
| 187 |
(model_order), ar_coeff); |
| 188 |
|
| 189 |
first = ts.data_int[i].last-(model_order)+1; |
| 190 |
i++; |
| 191 |
} |
| 192 |
// Mark filled gaps shorter than effective order as good. |
| 193 |
for(i=0; i<ts.m_gap; i++) |
| 194 |
{ |
| 195 |
if ( ts.gap_int[i].last-ts.gap_int[i].first+1 <= model_order) |
| 196 |
for (j = ts.gap_int[i].first; j<=ts.gap_int[i].last; j++) |
| 197 |
model_gaps[j] = 1; |
| 198 |
} |
| 199 |
#ifdef DEBUGOUT |
| 200 |
printf("done\n"); |
| 201 |
#endif |
| 202 |
if (iter>0) |
| 203 |
{ |
| 204 |
free(tsm.gap_int); free(tsm.data_int); |
| 205 |
} |
| 206 |
free(data_length); free(good_data); |
| 207 |
free(err); |
| 208 |
|
| 209 |
if (oldorder==model_order || tsm.m_gap==0) |
| 210 |
break; |
| 211 |
|
| 212 |
oldorder = model_order; |
| 213 |
} |
| 214 |
free(ts.gap_int); free(ts.data_int); |
| 215 |
memcpy(isgood,model_gaps,n*sizeof(int)); |
| 216 |
for (i=0;i<n;i++) |
| 217 |
if (model_gaps[i]==0) |
| 218 |
data[i] = 0; |
| 219 |
|
| 220 |
free(model_gaps); |
| 221 |
|
| 222 |
|
| 223 |
if (padends) |
| 224 |
{ |
| 225 |
memcpy(data_in, &data[maxorder],nold*sizeof(CTYPE)); |
| 226 |
memcpy(isgood_in,&isgood[maxorder],nold*sizeof(int)); |
| 227 |
free(data); |
| 228 |
free(isgood); |
| 229 |
} |
| 230 |
|
| 231 |
if (ar_coeff_in==NULL) |
| 232 |
free(ar_coeff); |
| 233 |
|
| 234 |
if (order) |
| 235 |
*order = model_order; |
| 236 |
|
| 237 |
return 0; |
| 238 |
} |
| 239 |
|
| 240 |
|
| 241 |
static void FILL_GAPS(int n, CTYPE *data, int *isgood, int order, |
| 242 |
CTYPE *ar_coeff) |
| 243 |
{ |
| 244 |
int iter, idx, i, shift, *gap_idx, *gap_idx_inv, *block; |
| 245 |
int n_gap, n_data, n_blocks, maxblock; |
| 246 |
CTYPE *rhs, *rhs2, *data_gap, *gamma2, *scratch; |
| 247 |
RTYPE rnorm, tol; |
| 248 |
void *aarg[3], *marg[3]; |
| 249 |
CTYPE *gamma; |
| 250 |
|
| 251 |
|
| 252 |
// for (i=0; i<10; i++) |
| 253 |
// printf("data[%d] = %f%+fi.\n",i,creal(data[i]),cimag(data[i])); |
| 254 |
|
| 255 |
// Compute gamma containing the first row of T^H*T. |
| 256 |
gamma = (CTYPE *)malloc((order+1)*sizeof(CTYPE)); |
| 257 |
rhs = (CTYPE *)malloc(2*(n-order)*sizeof(CTYPE)); |
| 258 |
for (shift = 0; shift<=order; shift++) |
| 259 |
{ |
| 260 |
gamma[shift] = 0; |
| 261 |
for (i=0; (i+shift)<=order; i++) |
| 262 |
gamma[shift] += 2*CONJ(ar_coeff[i+shift])*ar_coeff[i]; |
| 263 |
// printf("gamma[%d] = %f%+fi\n",shift,creal(gamma[shift]),cimag(gamma[shift])); |
| 264 |
} |
| 265 |
// for (i=0; i<=order; i++) |
| 266 |
// printf("ar_coeff[%d] = %f%+fi.\n",i,creal(ar_coeff[i]),cimag(ar_coeff[i])); |
| 267 |
|
| 268 |
|
| 269 |
// Find gaps. |
| 270 |
gap_idx = (int*)calloc(n,sizeof(int)); |
| 271 |
gap_idx_inv = (int*)calloc(n,sizeof(int)); |
| 272 |
block = (int*)calloc(n,sizeof(int)); |
| 273 |
n_data = 0; |
| 274 |
n_gap = 0; |
| 275 |
n_blocks = 0; |
| 276 |
maxblock = 0; |
| 277 |
for (i=0; i<n; i++) |
| 278 |
{ |
| 279 |
if (isgood[i]) |
| 280 |
{ |
| 281 |
n_data++; |
| 282 |
if (i>0 && !isgood[i-1]) |
| 283 |
n_blocks++; |
| 284 |
} |
| 285 |
else |
| 286 |
{ |
| 287 |
// printf("gap at %4d\n",i); |
| 288 |
gap_idx[i] = n_gap; |
| 289 |
gap_idx_inv[n_gap] = i; |
| 290 |
n_gap++; |
| 291 |
block[n_blocks]++; |
| 292 |
if (block[n_blocks] > maxblock) |
| 293 |
maxblock = block[n_blocks]; |
| 294 |
} |
| 295 |
} |
| 296 |
|
| 297 |
/* Set up right-hand side for block Toeplitz system. */ |
| 298 |
for (shift=0; (shift+order)<n; shift++) |
| 299 |
{ |
| 300 |
rhs[shift] = 0; |
| 301 |
rhs[shift+(n-order)] = 0; |
| 302 |
for (i=0; i<=order; i++) |
| 303 |
{ |
| 304 |
if (isgood[i+shift]) |
| 305 |
{ |
| 306 |
rhs[shift] -= ar_coeff[order-i]*data[i+shift]; |
| 307 |
rhs[shift+(n-order)] -= CONJ(ar_coeff[i])*data[i+shift]; |
| 308 |
} |
| 309 |
} |
| 310 |
} |
| 311 |
|
| 312 |
rhs2 = (CTYPE *)calloc(n_gap,sizeof(CTYPE)); |
| 313 |
for (shift=0; (shift+order)<n; shift++) |
| 314 |
{ |
| 315 |
for (i=0; i<=order; i++) |
| 316 |
{ |
| 317 |
if ( !isgood[i+shift] ) |
| 318 |
{ |
| 319 |
idx = gap_idx[i+shift]; |
| 320 |
// assert(idx>=0 && idx < n_gap); |
| 321 |
rhs2[idx] += CONJ(ar_coeff[order-i])*rhs[shift]; |
| 322 |
rhs2[idx] += ar_coeff[i]*rhs[shift+(n-order)]; |
| 323 |
|
| 324 |
} |
| 325 |
} |
| 326 |
} |
| 327 |
|
| 328 |
// printf("n_gap = %d\n",n_gap); |
| 329 |
//for (i=0; i<n_gap; i++) |
| 330 |
// printf("rhs2[%d] = %f%+fi.\n",i,creal(rhs2[i]),cimag(rhs2[i])); |
| 331 |
|
| 332 |
|
| 333 |
data_gap = (CTYPE*)calloc(n_gap,sizeof(CTYPE)); |
| 334 |
if (n_gap == 1) |
| 335 |
data_gap[0] = rhs2[0] / gamma[0]; |
| 336 |
else |
| 337 |
{ |
| 338 |
// Prepare PCG inputs defining the matrix. |
| 339 |
aarg[0] = (void *) ℴ |
| 340 |
aarg[1] = (void *) gap_idx_inv; |
| 341 |
aarg[2] = (void *) gamma; |
| 342 |
|
| 343 |
// Prepare PCG inputs defining the preconditioner. |
| 344 |
gamma2 = (CTYPE *) calloc(maxblock,sizeof(CTYPE)); |
| 345 |
scratch = (CTYPE *) calloc(maxblock,sizeof(CTYPE)); |
| 346 |
for (i=0; i<=min(order,maxblock-1); i++) |
| 347 |
gamma2[i] = gamma[i]; |
| 348 |
marg[0] = (void *) block; |
| 349 |
marg[1] = (void *) gamma2; |
| 350 |
marg[2] = (void *) scratch; |
| 351 |
|
| 352 |
// Solve normal equations using PCG with a block Toeplitz preconditioner. |
| 353 |
memset(data_gap, 0, n_gap*sizeof(CTYPE)); |
| 354 |
tol = sqrt(order)*PCG_TOL; |
| 355 |
iter = PCG(n_gap, PCG_MAXIT+2, tol, TOEPLITZ_MULT, |
| 356 |
PRECONDITION, rhs2, data_gap, &rnorm, aarg, marg); |
| 357 |
if (rnorm>tol) |
| 358 |
fprintf(stderr, "Warning: PCG did not converge. " |
| 359 |
"After %d iterations (maxit=%d) the relative" |
| 360 |
" residual was %e (tol=%e)\n", |
| 361 |
iter, PCG_MAXIT+2, rnorm,tol); |
| 362 |
|
| 363 |
/* printf("n_gap=%3d, maxblock=%3d, iter=%2d, rnorm=%e\n", |
| 364 |
n_gap, maxblock, iter,rnorm); */ |
| 365 |
// for (i=0; i<n_gap; i++) |
| 366 |
//printf("x[%d] = %f%+fi\n",i,creal(data_gap[i]),cimag(data_gap[i])); |
| 367 |
free(scratch); free(gamma2); |
| 368 |
} |
| 369 |
// Fill gap the original time series with extrapolated data. |
| 370 |
for (i=0; i<n_gap; i++) |
| 371 |
data[gap_idx_inv[i]] = data_gap[i]; |
| 372 |
|
| 373 |
free(block); free(gap_idx); free(gap_idx_inv); |
| 374 |
free(data_gap); free(gamma); free(rhs); free(rhs2); |
| 375 |
} |
| 376 |
|
| 377 |
|
| 378 |
|
| 379 |
static int GAPSTRUCTURE(int n, CTYPE *data, int *isgood, gapped_timeseries *ts, |
| 380 |
int **data_length, CTYPE ***good_data, int minpercentage) |
| 381 |
{ |
| 382 |
int i,j,first; |
| 383 |
|
| 384 |
// Count number of good and bad points and |
| 385 |
// number of good and bad intervals. |
| 386 |
ts->n_data = 0; ts->m_data = 0; |
| 387 |
ts->n_gap = 0; ts->m_gap = 0; |
| 388 |
for (i=0; i<n-1; i++) |
| 389 |
{ |
| 390 |
// assert(isgood[i]==1 || isgood[i]==0); |
| 391 |
if (isgood[i]==1) |
| 392 |
{ |
| 393 |
ts->n_data++; |
| 394 |
if (i==0) |
| 395 |
{ |
| 396 |
ts->first_is_good = 1; |
| 397 |
ts->m_data++; |
| 398 |
} |
| 399 |
if (isgood[i+1]==0) |
| 400 |
ts->m_gap++; |
| 401 |
} |
| 402 |
else if (isgood[i]==0) |
| 403 |
{ |
| 404 |
ts->n_gap++; |
| 405 |
if (i==0) |
| 406 |
{ |
| 407 |
ts->first_is_good = 0; |
| 408 |
ts->m_gap++; |
| 409 |
} |
| 410 |
if (isgood[i+1]==1) |
| 411 |
ts->m_data++; |
| 412 |
} |
| 413 |
else |
| 414 |
{ |
| 415 |
fprintf(stderr, "ERROR: problem in routine gapstructure: value of window function must be either 0 or 1\n"); |
| 416 |
exit(1); |
| 417 |
} |
| 418 |
} |
| 419 |
if (isgood[n-1]==1) |
| 420 |
ts->n_data++; |
| 421 |
else |
| 422 |
ts->n_gap++; |
| 423 |
|
| 424 |
|
| 425 |
#ifdef DEBUGOUT |
| 426 |
printf("n_data = %d, m_data = %d\nn_gap = %d, m_gap = %d\n", |
| 427 |
ts->n_data,ts->m_data,ts->n_gap,ts->m_gap); |
| 428 |
#endif |
| 429 |
/* |
| 430 |
assert((isgood[0]==0 && isgood[n-1]==0 && ts->m_gap == ts->m_data+1) || |
| 431 |
(isgood[0]==1 && isgood[n-1]==1 && ts->m_gap == ts->m_data-1) || |
| 432 |
(isgood[0]!=isgood[n-1] && ts->m_gap == ts->m_data)); |
| 433 |
*/ |
| 434 |
|
| 435 |
if (!((isgood[0]==0 && isgood[n-1]==0 && ts->m_gap == ts->m_data+1) || |
| 436 |
(isgood[0]==1 && isgood[n-1]==1 && ts->m_gap == ts->m_data-1) || |
| 437 |
(isgood[0]!=isgood[n-1] && ts->m_gap == ts->m_data))) |
| 438 |
{ |
| 439 |
fprintf(stderr, "ERROR: problem in routine gapstructure\n"); |
| 440 |
exit(1); |
| 441 |
} |
| 442 |
|
| 443 |
|
| 444 |
// Allocate space for data pointers and AR coefficients. |
| 445 |
*good_data = (CTYPE **) malloc(ts->m_data*sizeof(CTYPE *)); |
| 446 |
*data_length = (int *) malloc(ts->m_data*sizeof(int)); |
| 447 |
ts->gap_int = (interval *) malloc(ts->m_gap*sizeof(interval)); |
| 448 |
ts->data_int = (interval *) malloc(ts->m_data*sizeof(interval)); |
| 449 |
|
| 450 |
|
| 451 |
// Scan again and set up data structures describing gap structure. |
| 452 |
j = 0; |
| 453 |
first=0; |
| 454 |
if (isgood[0]==1) |
| 455 |
{ |
| 456 |
first=0; |
| 457 |
ts->data_int[0].first = 0; |
| 458 |
(*good_data)[0] = &data[0]; |
| 459 |
j++; |
| 460 |
} |
| 461 |
for (i=0; i<n-1; i++) |
| 462 |
{ |
| 463 |
// first of an interval |
| 464 |
if (isgood[i]==0 && isgood[i+1]==1) |
| 465 |
{ |
| 466 |
first = i+1; |
| 467 |
ts->data_int[j].first = first; |
| 468 |
(*good_data)[j] = &(data[i+1]); |
| 469 |
// printf("first[%d] = %d\n",j,i+1); |
| 470 |
j++; |
| 471 |
} |
| 472 |
// end of an interval |
| 473 |
if (isgood[i]==1 && isgood[i+1]==0) |
| 474 |
{ |
| 475 |
(*data_length)[j-1] = i-first+1; |
| 476 |
ts->data_int[j-1].last = i; |
| 477 |
} |
| 478 |
} |
| 479 |
if (isgood[n-1]==1) |
| 480 |
{ |
| 481 |
(*data_length)[j-1] = n-1-first+1; |
| 482 |
ts->data_int[j-1].last = n-1; |
| 483 |
} |
| 484 |
if (ts->first_is_good) |
| 485 |
{ |
| 486 |
for(j=0; j<ts->m_data;j++) |
| 487 |
{ |
| 488 |
if (j<ts->m_gap) |
| 489 |
ts->gap_int[j].first = ts->data_int[j].last+1; |
| 490 |
if (j>0) |
| 491 |
ts->gap_int[j-1].last = ts->data_int[j].first-1; |
| 492 |
} |
| 493 |
if (ts->m_data==ts->m_gap) |
| 494 |
ts->gap_int[ts->m_gap-1].last = n-1; |
| 495 |
} |
| 496 |
else |
| 497 |
{ |
| 498 |
ts->gap_int[0].first = 0; |
| 499 |
for(j=0; j<ts->m_data;j++) |
| 500 |
{ |
| 501 |
if (j<ts->m_gap) |
| 502 |
ts->gap_int[j].last = ts->data_int[j].first-1; |
| 503 |
if (j<ts->m_gap-1) |
| 504 |
ts->gap_int[j+1].first = ts->data_int[j].last+1; |
| 505 |
} |
| 506 |
if (ts->m_gap==ts->m_data+1) |
| 507 |
ts->gap_int[ts->m_gap-1].last = n-1; |
| 508 |
} |
| 509 |
return maxorder(ts, *data_length, minpercentage); |
| 510 |
} |
| 511 |
|
| 512 |
|
| 513 |
|
| 514 |
// Matrix-vector multiply with symmetric Toeplitz matrix. |
| 515 |
// y = T(r) x |
| 516 |
|
| 517 |
static void TOEPLITZ_MULT(int n, CTYPE x[n], CTYPE y[n], void **data) |
| 518 |
{ |
| 519 |
int i,j,order,idx,ii; |
| 520 |
int *gap_idx_inv; |
| 521 |
CTYPE *r; |
| 522 |
|
| 523 |
order = *((int *) data[0]); |
| 524 |
gap_idx_inv = (int *) data[1]; |
| 525 |
r = (CTYPE *) data[2]; |
| 526 |
|
| 527 |
for (i=0; i<n; i++) |
| 528 |
y[i] = r[0]*x[i]; |
| 529 |
|
| 530 |
for (i=0; i<n; i++) |
| 531 |
{ |
| 532 |
ii = gap_idx_inv[i]; |
| 533 |
for (j=i+1; j<n; j++) |
| 534 |
{ |
| 535 |
idx = ii-gap_idx_inv[j]; |
| 536 |
if (idx <= order && idx >= -order) |
| 537 |
{ |
| 538 |
if (idx>=0) |
| 539 |
{ |
| 540 |
y[i] += CONJ(r[idx])*x[j]; |
| 541 |
y[j] += r[idx]*x[i]; |
| 542 |
} |
| 543 |
else |
| 544 |
{ |
| 545 |
y[i] += r[-idx]*x[j]; |
| 546 |
y[j] += CONJ(r[-idx])*x[i]; |
| 547 |
} |
| 548 |
} |
| 549 |
} |
| 550 |
} |
| 551 |
} |
| 552 |
|
| 553 |
// Solve M y = x; |
| 554 |
static void PRECONDITION(int n, CTYPE *b, CTYPE *x, void **data) |
| 555 |
{ |
| 556 |
int i,first,*block; |
| 557 |
CTYPE *gamma, *scratch; |
| 558 |
block = (int *) data[0]; |
| 559 |
gamma = (CTYPE *) data[1]; |
| 560 |
scratch = (CTYPE *) data[2]; |
| 561 |
|
| 562 |
|
| 563 |
for(i=0, first = 0; first<n; first+=block[i], i++) |
| 564 |
{ |
| 565 |
memcpy(scratch,&b[first],block[i]*sizeof(CTYPE)); |
| 566 |
LEVINSON(block[i],gamma,scratch,&x[first]); |
| 567 |
} |
| 568 |
} |
| 569 |
|
| 570 |
|
| 571 |
#undef GAPSTRUCTURE |
| 572 |
#undef FAHLMAN_ULRYCH |
| 573 |
#undef BURG |
| 574 |
#undef FILL_GAPS |
| 575 |
#undef TOEPLITZ_MULT |
| 576 |
#undef PRECONDITION |
| 577 |
#undef PCG |
| 578 |
#undef LEVINSON |
| 579 |
#undef EPSILON |
| 580 |
#include "ctypes_undef.h" |
