source: stable/2.4/ADOL-C/src/drivers/taylor.c @ 397

Last change on this file since 397 was 333, checked in by awalther, 7 years ago

fix error in taylor.c: number of directions left should be tested with <= instead of <

  • Property svn:keywords set to Author Date Id Revision
File size: 27.0 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     drivers/taylor.c
4 Revision: $Id: taylor.c 333 2012-07-04 14:00:42Z kulshres $
5 Contents: Easy to use drivers for the evaluation of higher order derivative
6           tensors and inverse/impicit function differentiation
7 
8 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz, Olaf Vogel
9 
10 This file is part of ADOL-C. This software is provided as open source.
11 Any use, reproduction, or distribution of the software constitutes
12 recipient's acceptance of the terms of the accompanying license file.
13 
14----------------------------------------------------------------------------*/
15#include <adolc/drivers/taylor.h>
16#include <adolc/interfaces.h>
17#include <adolc/adalloc.h>
18#include "taping_p.h"
19
20#include <math.h>
21
22BEGIN_C_DECLS
23
24
25/****************************************************************************/
26/*                                                              STRUCT ITEM */
27struct item {
28    int a;                 /* address in array of derivatives */
29    int b;                 /* absolute value of the correspondig multiindex i */
30    double c;              /* value of the coefficient c_{i,j} */
31    struct item *next;     /* next item */
32};
33
34/****************************************************************************/
35/*                                                     DEALLOCATE COEFFLIST */
36void freecoefflist( int dim, struct item *coeff_list ) {
37    int i;
38    struct item *ptr1;
39    struct item *ptr2;
40
41    for (i=0; i<dim; i++)  /* sum over all multiindices jm with |jm| = d */
42    { ptr1 = &coeff_list[i];
43        ptr1 = ptr1->next;
44        while (ptr1 != NULL) {
45            ptr2 = ptr1->next;
46            free((char *) ptr1);
47            ptr1 = ptr2;
48        }
49    }
50}
51
52
53/****************************************************************************/
54/*                                            ALLOCATE/DEALLOCATE COEFFLIST */
55double* tensoriglob;
56
57/*--------------------------------------------------------------------------*/
58/* Allcoate space for symmetric derivative tensors
59   of up to order d in n variables, derivatives are  */
60void* tensorpoint( int n, int d ) {
61    int i;
62    void* t;
63
64    if (d == 1) {
65        t = (void*) tensoriglob;
66        tensoriglob += n+1;
67    } else {
68        t = (void*) malloc((n+1)*sizeof(void*));
69        for (i=0; i<=n; i++)
70            ((void**)t)[i] = (void*) tensorpoint(i,d-1);
71    }
72    return t;
73}
74
75/*--------------------------------------------------------------------------*/
76void** tensorsetup( int m, int n, int d, double** tensorig ) {
77    int i;
78    void** t = (void**) malloc(m*sizeof(void*));
79
80    for (i=0; i<m; i++) {
81        tensoriglob = tensorig[i];
82        t[i] = (void*) tensorpoint(n,d);
83    }
84    return t;
85}
86
87/*--------------------------------------------------------------------------*/
88/* Deallocate space for symmetric derivative tensors
89   of up to order d in n variables  */
90void freetensorpoint( int n, int d, double** tensor ) {
91    int i;
92    double* t;
93
94    if (d > 2)
95        for(i=0;i<=n;i++) {
96            t = tensor[i];
97            freetensorpoint(i,d-1,(double **) t);
98            free((char *) t);
99        }
100}
101
102/*--------------------------------------------------------------------------*/
103void freetensor( int m, int n, int d, double** tensor ) {
104    int i;
105    double* t;
106
107    for (i=0; i<m; i++) {
108        t = tensor[i];
109        freetensorpoint(n,d,(double **)t);
110        free((char *) t);
111    }
112}
113
114
115/****************************************************************************/
116/*                                                           SOME UTILITIES */
117
118
119long binomi(int n, int k) {
120    long double accum = 1;
121    unsigned int i;
122
123    if (k > n)
124        return 0;
125
126    if (k > n/2)
127        k = n-k;
128
129    for (i = 1; i <= k; i++)
130         accum = accum * (n-k+i) / i;
131
132    return (long) accum + 0.5;
133}
134
135/*--------------------------------------------------------------------------*/
136double dbinomi( double a, int b ) {
137    int i;
138    double result = 1.0;
139
140    for (i=1; i<=b; i++)
141        result = result*(a-i+1)/i;
142    return result;
143}
144
145/*--------------------------------------------------------------------------*/
146double summand(int p, int d, int* jm, int* km, int order_im, int order_km, long binomiZ) {   /* calculates summation value for fixed j, i, k with terms used in the article.*/
147    int i;
148    double result, order_k_by_d;
149
150    order_k_by_d = order_km/(double)d;
151    result = 1.0;
152    for (i=0; i<order_im; i++) result *= order_k_by_d;     /* (|k|/d)^i */
153    if ((order_im+order_km)%2==1) result *= -1.0;             /* (-1)^(|i-k|) */
154    result *= binomiZ;
155    for (i=0; i<p; i++) result *= dbinomi(d*km[i]/(double)order_km, jm[i]);
156    return result;
157}
158
159/****************************************************************************/
160/*                                                    EVALUATE COEFFICIENTS */
161
162void coeff(int p, int d, struct item* coeff_list) {
163    int i, j, u, n, index_coeff_list, order_im, order_km, address;
164    int* jm = (int*) malloc(p*sizeof(int));     /* Multiindex j */
165    int* im = (int*) malloc(p*sizeof(int));     /* Multiindex i */
166    int* km = (int*) malloc(p*sizeof(int));     /* Multiindex k */
167    struct item* ptr;
168    double sum;
169    long binomiZ;           /* whole number binomial coefficient */
170
171    jm[0] = d;
172    for (i=1; i<p; i++) jm[i] = im[i] = 0;
173    for (i=0; i<p; i++) km[i] = 0;
174    order_km = 0;
175
176    for (index_coeff_list = 0; 1; index_coeff_list++) {  /* travers coeff_list, i.e. create all j with |j| = d. */
177        ptr = NULL;
178        for (order_im=1; order_im<=d; order_im++) {  /* travers all orders from 1 to d for i */
179            im[p-1]=0;
180            im[0] = order_im;
181            while (1) {   /* create all i with order order_im. */
182                sum = 0;
183                binomiZ = 1;
184                for (i=0; i<p; i++) /* check, whether i valid. */
185                    if ((jm[i]>0)&&(im[i]==0)) break;
186                if (i==p)
187                    while (1) {   /* create all k where 0<k<=i */
188                        for (i=p-1; i>=0; i--)
189                            if (km[i]<im[i]) {
190                                km[i]++;
191                                order_km++;
192                                binomiZ *= im[i]-km[i]+1;  /* for (i over k)/(i over k-1) = (i-k+1)/k */
193                                binomiZ /= km[i];
194                                break;
195                            } else {
196                                /* binomiZ doesn't change, for (i over k) = 1 if k=0 and also if k=i */
197                                order_km -= km[i];
198                                km[i] = 0;
199                            };
200                        if (i==-1) break;
201
202                        sum += summand(p,d,jm,km,order_im,order_km,binomiZ);
203                    };
204
205                if (fabs(sum) > 0) { /* Store coefficient */
206                    if (ptr==NULL)
207                        ptr = &coeff_list[index_coeff_list];
208                    else {
209                        ptr->next = (struct item*) malloc(sizeof(struct item));
210                        ptr = ptr->next;
211                    };
212
213                    address = 0; /* calculate address for ptr->a */
214                    j = d-order_im+1;
215                    for (u=0; u<p; u++)  /* It is sum(binomial(i+k,j+k),k=0..n-1) = */
216                        if (im[u]!=0)      /* = ((j+n)*binomial(i+n,j+n)-j*binomial(i,j))/(1+i-j) */
217                        {
218                            i = u+j;
219                            n = im[u];
220                            address += ((j+n)*binomi(i+n,j+n)-binomi(i,j)*j)/(1+i-j);
221                            j += n;
222                        };
223                    ptr->a = address;
224                    ptr->b = order_im;
225                    ptr->c = sum;
226                };
227
228                if ((im[p-1]==order_im)||(p==1)) break;
229                for (i=p-2; im[i]==0; i--); /* find first nonvanishing entry on the right. */
230                im[i]--;
231                im[i+1] = im[p-1]+1;
232                if (i!=p-2) im[p-1] = 0;
233            };
234        };
235
236        ptr->next = NULL; /* mark end of queue. */
237
238        if ((jm[p-1]==d)||(p==1)) break;
239        for (i=p-2; jm[i]==0; i--); /* find first nonvanishing entry on the right. */
240        jm[i]--;
241        jm[i+1] = jm[p-1]+1;
242        if (i!=p-2) jm[p-1] = 0;
243    };
244
245    free((char*) jm);
246    free((char*) im);
247    free((char*) km);
248}
249
250
251/*--------------------------------------------------------------------------*/
252void convert( int p, int d, int *im, int *multi ) {
253    int i;
254
255    for (i=0; i<p; i++)
256        multi[i] = 0;
257    for (i=0; i<d; i++)
258        if (im[i]) /* olvo/walther 980804 new tl */
259            multi[im[i]-1] += 1;
260}
261
262/*--------------------------------------------------------------------------*/
263int address( int d, int* im ) {
264    int i,
265    add = 0;
266
267    for (i=0; i<d; i++)
268        add += binomi(im[i]+i,i+1);
269    return add;
270}
271
272
273
274/****************************************************************************/
275/*                                                           MORE UTILITIES */
276
277/*--------------------------------------------------------------------------*/
278void multma3vec2( int n, int p, int d, int bd,
279                  double ***X, double **S, int **jm ) {
280    int i,j,k;
281    double sum;
282
283    for (i=0; i<n; i++)
284        for (k=0; k<bd; k++) {
285            sum = 0;
286            for (j=0; j<p; j++)
287                sum += S[i][j]*jm[k][j];
288            X[i][k][0] = sum;
289            for (j=1; j<d; j++)
290                X[i][k][j] = 0;
291        }
292}
293
294/*--------------------------------------------------------------------------*/
295void multma2vec2( int n, int p, int bd, double **X, double **S, int **jm ) {
296    int i,j,k;
297    double sum;
298
299    for (i=0; i<n; i++)
300        for (k=0; k<bd; k++) {
301            sum = 0;
302            for (j=0; j<p; j++)
303                sum += S[i][j]*jm[k][j];
304            X[i][k] = sum;
305        }
306}
307
308/*--------------------------------------------------------------------------*/
309void multma2vec1( int n, int p, int d, double **X, double **S, int *jm ) {
310    int i,j;
311    double sum;
312
313    for (i=0; i<n; i++) {
314        sum = 0;
315        for (j=0; j<p; j++)
316            sum += S[i][j]*jm[j];
317        X[i][1] = sum;
318        for (j=2; j<d; j++)
319            X[i][j] = 0;
320    }
321}
322
323
324/****************************************************************************/
325
326/* test if zero */
327#define ZERO 1.0E-15
328
329/*--------------------------------------------------------------------------*/
330int LUFactorization( double** J, int n, int* RI, int* CI ) {
331    int i, j, k, cIdx, rIdx, h;
332    double v;
333
334    for (i=0; i<n; i++)
335        RI[i]=i;
336    for (j=0; j<n; j++)
337        CI[j]=j;
338    /* n Gausz-steps with full Pivoting */
339    for (k=0; k<n; k++) {
340        v=0.0;
341        cIdx=rIdx=0;
342        /* Pivotsearch */
343        for (i=k; i<n; i++)
344            for (j=k; j<n; j++)
345                if (fabs(J[RI[i]][CI[j]])>v) {
346                    v=fabs(J[RI[i]][CI[j]]);
347                    rIdx=i;
348                    cIdx=j;
349                }
350        if (ZERO > v) {
351            fprintf(DIAG_OUT,
352                    "Error:LUFactorisation(..): no Pivot in step %d (%E)\n",k+1,v);
353            return -(k+1);
354        }
355        /* row and column change resp. */
356        if (rIdx > k) {
357            h=RI[k];
358            RI[k]=RI[rIdx];
359            RI[rIdx]=h;
360        }
361        if (cIdx > k) {
362            h=CI[k];
363            CI[k]=CI[cIdx];
364            CI[cIdx]=h;
365        }
366        /* Factorisation step */
367        for (i=k+1; i<n; i++) { /* L-part */
368            J[RI[i]][CI[k]]/=J[RI[k]][CI[k]];
369            /* R-part */
370            for (j=k+1; j<n; j++)
371                J[RI[i]][CI[j]]-=J[RI[i]][CI[k]]*J[RI[k]][CI[j]];
372        }
373    }
374    return k;
375}
376
377/*--------------------------------------------------------------------------*/
378void GauszSolve( double** J, int n, int* RI, int* CI, double* b ) {
379    double* tmpZ;
380    int i,j;
381
382    tmpZ = myalloc1(n);
383    for (i=0; i<n; i++) {
384        tmpZ[i]=b[RI[i]];
385        for (j=0; j<i; j++)
386            tmpZ[i]-=J[RI[i]][CI[j]]*tmpZ[j];
387    }
388    for (i=n-1; i>=0; i--) {
389        b[CI[i]]=tmpZ[i];
390        for (j=i+1; j<n; j++)
391            b[CI[i]]-=J[RI[i]][CI[j]]*b[CI[j]];
392        b[CI[i]]/=J[RI[i]][CI[i]];
393    }
394    free(tmpZ);
395}
396
397
398/****************************************************************************/
399int jac_solv( unsigned short tag, int n, const double* x, double* b, unsigned short mode ) {
400    double *y;
401    int i, newX = 0;
402    int rc = 3;
403    ADOLC_OPENMP_THREAD_NUMBER;
404
405    ADOLC_OPENMP_GET_THREAD_NUMBER;
406    y = myalloc1(n);
407    if (n != ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_nax) {
408        if (ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_nax) {
409            free(ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_ci);
410            free(ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_ri);
411            myfree1(ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_xold);
412            myfreeI2(ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_nax,
413                    ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_I);
414            myfree2(ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_J);
415        }
416        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_J = myalloc2(n,n);
417        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_I = myallocI2(n);
418        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_xold = myalloc1(n);
419        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_ri =
420            (int*)malloc(n*sizeof(int));
421        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_ci =
422            (int*)malloc(n*sizeof(int));
423
424        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_modeold = 0;
425        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_nax = n;
426    }
427    for (i = 0; i < n; ++i)
428        if (x[i] != ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_xold[i]) {
429            ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_xold[i] = x[i];
430            newX = 1;
431        }
432    switch(mode) {
433        case 0:
434            MINDEC(rc,zos_forward(tag, n, n, 1, x, y));
435            MINDEC(rc,fov_reverse(tag, n, n, n,
436                        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_I,
437                        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_J));
438            break;
439        case 1:
440            MINDEC(rc,zos_forward(tag, n, n, 1, x, y));
441            MINDEC(rc,fov_reverse(tag, n, n, n,
442                            ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_I,
443                            ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_J));
444            if (LUFactorization(
445                        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_J, n,
446                        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_ri,
447                        ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_ci) < 0)
448            {
449                rc = -3;
450                break;
451            }
452            ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_modeold = 1;
453            break;
454        case 2:
455            if ((ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_modeold < 1) ||
456                    (newX == 1))
457            {
458                MINDEC(rc,zos_forward(tag, n, n, 1, x, y));
459                MINDEC(rc,fov_reverse(tag, n, n, n,
460                            ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_I,
461                            ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_J));
462                if (LUFactorization(
463                            ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_J, n,
464                            ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_ri,
465                            ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_ci) < 0)
466                {
467                    rc = -3;
468                    break;
469                }
470            }
471            GauszSolve(ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_J, n,
472                    ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_ri,
473                    ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_ci, b);
474            ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.jacSolv_modeold = 2;
475            break;
476    }
477    myfree1(y);
478    return rc;
479}
480
481
482/****************************************************************************/
483int inverse_Taylor_prop( short tag, int n, int d,
484                         double** Y, double** X ) {
485    int i,j,l,q;
486    static double **I;
487    register double bi;
488    static double** Xhelp;
489    static double** W;
490    static double* xold;
491    static double ***A;
492    static double *w;
493    static int *dd;
494    static double *b;
495    static int nax,dax,bd,cgd;
496    static short **nonzero;
497    short* nz;
498    double* Aij;
499    double* Xj;
500    int ii, di, da, Di, indexA, indexX;
501    int rc = 3;
502
503    /* Re/Allocation Stuff */
504    if ((n != nax) || (d != dax))
505       {
506        if (nax)
507         {
508            free(**A);
509            free(*A);
510            free(A);
511            free(*I);
512            free(I);
513            free(*W);
514            free(W);
515            free(*Xhelp);
516            free(Xhelp);
517            free(w);
518            free(xold);
519            free(*nonzero);
520            free(nonzero);
521            free(dd);
522            free(b);
523        }
524        A = myalloc3(n,n,d+1);
525        I = myalloc2(n,n);
526        W = myalloc2(n,d);
527        Xhelp = myalloc2(n,d);
528        w = myalloc1(n);
529        dd = (int*)malloc((d+1)*sizeof(int));
530        b  = (double*)malloc(n*sizeof(double));
531        xold = (double*)malloc(n*sizeof(double));
532        nonzero = (short**)malloc(n*sizeof(short*));
533        nz = (short*)malloc(n*n*sizeof(short));
534        for (i=0; i<n; i++) {
535            nonzero[i] = nz;
536            nz = nz + n;
537            xold[i] = 0;
538            for (j=0; j<n; j++)
539                I[i][j]=(i==j)?1.0:0.0;
540        }
541        cgd = 1;
542        nax=n;
543        dax=d;
544        dd[0] = d+1;
545        i = -1;
546        while(dd[++i] > 1)
547            dd[i+1] = (int)ceil(dd[i]*0.5);
548        bd = i+1;
549    }
550    if (cgd == 0)
551        for (i=0; i<n; i++)
552            if (X[i][0] != xold[i])
553                cgd = 1;
554    for(i=0;i<n;i++)
555      b[i] = 0;
556    if (cgd == 1) {
557        cgd = 0;
558        for (i=0; i<n; i++)
559            xold[i] = X[i][0];
560        MINDEC(rc,jac_solv(tag,n,xold,b,1));
561        if (rc == -3)
562            return -3;
563    }
564    ii = bd;
565    for (i=0; i<n; i++)
566      {
567        for (j=0; j<d; j++)
568        {
569            Xhelp[i][j] = 0;
570            X[i][j+1] = 0;
571            W[i][j] = 0;
572        }
573        for (j=0; j<n; j++)
574          for (l=0; l<=d; l++)
575            A[i][j][l] = 0;
576      }
577
578    while (--ii > 0) {
579        di = dd[ii-1]-1;
580        Di = dd[ii-1]-dd[ii]-1;
581        MINDEC(rc,hos_forward(tag,n,n,di,Di+1,xold,Xhelp,w,W));
582        MINDEC(rc,hov_reverse(tag,n,n,Di,n,I,A,nonzero));
583        da = dd[ii];
584        for (l=da; l<dd[ii-1]; l++) {
585            for (i=0; i<n; i++) {
586                if (l == 0)
587                    bi = w[i]-Y[i][0];
588                else
589                    bi = W[i][l-1]-Y[i][l];
590                for (j=0; j<n; j++)
591                    if (nonzero[i][j]>1) {
592                        Aij = A[i][j];
593                        indexA = l-1;
594                        Xj = X[j]+l;
595                        indexX = 1;
596                        if (da == l-1)
597                          {
598                            bi += (*(++Aij))*(*(--Xj));
599                          }
600                        else
601                          {
602                            for (q=da; q<l; q++)
603                              {
604                                bi += (*(++Aij))*(*(--Xj));
605                                bi += A[i][j][indexA]*X[j][indexX];
606                                indexA--;
607                                indexX++;
608                              }
609                          }
610                    }
611                b[i] = -bi;
612            }
613            MINDEC(rc,jac_solv(tag,n,xold,b,2));
614           if (rc == -3)
615                return -3;
616            for (i=0; i<n; i++) {
617              X[i][l] += b[i];
618              Xhelp[i][l-1] += b[i];
619            }
620        }
621    }
622    return rc;
623}
624
625
626/****************************************************************************/
627int inverse_tensor_eval( short tag, int n, int d, int p,
628                         double *x, double **tensor, double** S ) {
629    static int dim;
630    static int dold,pold;
631    static struct item *coeff_list;
632    int i,j,dimten;
633    int *it = (int*) malloc(d*sizeof(int));
634    double** X;
635    double** Y;
636    int *jm;
637    double *y = (double*) malloc(n*sizeof(double));
638    struct item *ptr;
639    int rc = 3;
640
641    dimten=binomi(p+d,d);
642    for(i=0;i<n;i++)
643        for(j=0;j<dimten;j++)
644            tensor[i][j] = 0;
645    MINDEC(rc,zos_forward(tag,n,n,0,x,y));
646    if (d > 0) {
647        if ((d != dold) || (p != pold)) {
648            if (pold) { /* olvo 980728 */
649                dim = binomi(pold+dold-1,dold);
650                freecoefflist(dim,coeff_list);
651                free((char*) coeff_list);
652            }
653            dim = binomi(p+d-1,d);
654            coeff_list = (struct item *) malloc(sizeof(struct item)*dim);
655            coeff(p,d, coeff_list);
656            dold = d;
657            pold = p;
658        }
659        jm = (int *)malloc(sizeof(int)*p);
660        X = myalloc2(n,d+1);
661        Y = myalloc2(n,d+1);
662        for (i=0; i<n; i++) {
663            X[i][0] = x[i];
664            Y[i][0] = y[i];
665            for (j=1; j<d; j++)
666              {
667                X[i][j] = 0;
668                Y[i][j] = 0;
669              }
670        }
671        if (d == 1) {
672            it[0] = 0;
673            for (i=0; i<dim; i++)  /* sum over all multiindices jm with |jm| = d */
674            { it[0] = it[0]+1;
675                convert(p,d,it,jm);
676                ptr = &coeff_list[i];
677                multma2vec1(n,p,d,Y,S,jm);
678                MINDEC(rc,inverse_Taylor_prop(tag,n,d,Y,X));
679                if (rc == -3)
680                    return -3;
681                do {
682                    for(j=0;j<n;j++)
683                        tensor[j][ptr->a] += X[j][ptr->b]*ptr->c;
684                    ptr = ptr->next;
685                } while (ptr != NULL);
686            }
687        } else {
688            for (i=0; i<d-1; i++)
689                it[i] = 1;
690            it[d-1] = 0;
691            for (i=0; i<dim; i++)  /* sum over all multiindices jm with |jm| = d */
692            { it[d-1] = it[d-1]+1;
693                for (j=d-2; j>=0; j--)
694                    it[j] = it[j] + it[j+1]/(p+1);
695                for (j=1; j<d; j++)
696                    if (it[j] > p) it[j] = it[j-1];
697                convert(p,d,it,jm);
698                multma2vec1(n,p,d,Y,S,jm); /* Store S*jm in Y */
699                MINDEC(rc,inverse_Taylor_prop(tag,n,d,Y,X));
700                if (rc == -3)
701                    return -3;
702                ptr = &coeff_list[i];
703                do {
704                    for(j=0;j<n;j++)
705                      {
706                        tensor[j][ptr->a] += X[j][ptr->b]*ptr->c;
707                      }
708                    ptr = ptr->next;
709                } while (ptr != NULL);
710            }
711        }
712        free((char*) jm);
713        free((char*) *X);
714        free((char*) X);
715        free((char*) *Y);
716        free((char*) Y);
717    }
718    for(i=0;i<n;i++)
719        tensor[i][0] = x[i];
720    free((char*) y);
721    free((char*) it);
722    return rc;
723}
724
725
726/****************************************************************************/
727int tensor_eval( short tag, int m, int n, int d, int p,
728                 double* x, double **tensor, double **S ) {
729    static int bd,dim;
730    static int dold,pold;
731    static struct item *coeff_list;
732    int i,j,k,dimten,ctr;
733    int **jm, jmbd=0;
734    int *it = (int*) malloc(d*sizeof(int));
735    double *y = (double*) malloc(m*sizeof(double));
736    double*** X;
737    double*** Y;
738    struct item *ptr[10];
739    int rc = 3;
740
741    dimten=binomi(p+d,d);
742    for (i=0; i<m; i++)
743        for (j=0; j<dimten; j++)
744            tensor[i][j] = 0;
745    if (d == 0) {
746        MINDEC(rc,zos_forward(tag,m,n,0,x,y));
747    } else {
748        if ((d != dold) || (p != pold)) {
749            if (pold) {
750                dim = binomi(pold+dold-1,dold);
751                freecoefflist(dim,coeff_list);
752                free((char*) coeff_list);
753            }
754            dim = binomi(p+d-1,d);
755            if (dim < 10)
756                bd = dim;
757            else
758                bd = 10;
759            coeff_list = (struct item *) malloc(sizeof(struct item)*dim);
760            coeff(p,d, coeff_list);
761            dold = d;
762            pold = p;
763        }
764        jmbd = bd;
765        jm = (int **) malloc(jmbd*sizeof(int*));
766        for (i=0; i<jmbd; i++)
767            jm[i] = (int *) malloc(p*sizeof(int));
768        if (d == 1) {
769            X = myalloc3(1,n,bd);
770            Y = myalloc3(1,m,bd);
771            ctr   = 0;
772            it[0] = 0;
773            for (i=0; i<dim; i++) /* sum over all multiindices jm with |jm| = d */
774            { it[0] = it[0]+1;
775                convert(p,d,it,jm[ctr]);
776                ptr[ctr] = &coeff_list[i];
777                if (ctr < bd-1)
778                    ctr += 1;
779                else {
780                    multma2vec2(n,p,bd,X[0],S,jm);
781                    MINDEC(rc,fov_forward(tag,m,n,bd,x,X[0],y,Y[0]));
782                    for (k=0; k<bd; k++)
783                        do {
784                            for (j=0; j<m; j++)
785                              {
786                                tensor[j][ptr[k]->a] += Y[0][j][k]*ptr[k]->c;
787                              }
788                           ptr[k] = ptr[k]->next;
789                        } while (ptr[k] != NULL);
790                    if (dim-i <= bd)
791                        bd = dim-i-1;
792                    ctr = 0;
793                }
794            }
795        } else {
796            X = myalloc3(n,bd,d);
797            Y = myalloc3(m,bd,d);
798            ctr = 0;
799            for (i=0; i<d-1; i++)
800                it[i] = 1;
801            it[d-1] = 0;
802            for (i=0; i<dim; i++) /* sum over all multiindices jm with |jm| = d */
803            { it[d-1] = it[d-1]+1;
804                for (j=d-2; j>=0; j--)
805                    it[j] = it[j] + it[j+1]/(p+1);
806                for (j=1; j<d; j++)
807                    if (it[j] > p)
808                        it[j] = it[j-1];
809                convert(p,d,it,jm[ctr]);
810                ptr[ctr] = &coeff_list[i];
811                if (ctr < bd-1)
812                    ctr += 1;
813                else {
814                    multma3vec2(n,p,d,bd,X,S,jm);
815                    MINDEC(rc,hov_forward(tag,m,n,d,bd,x,X,y,Y));
816                    for (k=0; k<bd; k++)
817                        do {
818                            for (j=0; j<m; j++)
819                                tensor[j][ptr[k]->a] += Y[j][k][ptr[k]->b-1]*ptr[k]->c;
820                            ptr[k] = ptr[k]->next;
821                        } while (ptr[k] != NULL);
822                    if (dim-i <= bd)
823                        bd = dim-i-1;
824                    ctr = 0;
825                }
826            }
827        }
828        for (i=0; i<jmbd; i++)
829            free((char*) *(jm+i));
830        free((char*) jm);
831        free((char*) **X);
832        free((char*) *X);
833        free((char*) X);
834        free((char*) **Y);
835        free((char*) *Y);
836        free((char*) Y);
837    }
838    for(i=0;i<m;i++)
839        tensor[i][0] = y[i];
840    bd = jmbd;
841    free((char*) y);
842    free((char*) it);
843    return rc;
844}
845
846
847/****************************************************************************/
848void tensor_value( int d, int m, double *y, double **tensor, int *multi ) {
849    int i, j, max, ind, add;
850    int *im = (int*) malloc(d*sizeof(int));
851
852    max = 0;
853    ind = d-1;
854    for (i=0; i<d; i++) {
855        if (multi[i] > max)
856            max = multi[i];
857        im[i] = 0;
858    }
859    for (i=0; i<d; i++) {
860        if (multi[i] == max)  /* olvo 980728 == instead of = */
861        { im[ind] = multi[i];
862            multi[i] = 0;
863            max = 0;
864            ind -= 1;
865            for (j=0; j<d; j++)
866                if (multi[j] > max)
867                    max = multi[j];
868        }
869    }
870    add = address(d,im);
871    for (i=0; i<m; i++)
872        y[i] = tensor[i][add];
873    free((char*) im);
874}
875
876END_C_DECLS
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
Note: See TracBrowser for help on using the repository browser.