source: trunk/ADOL-C/src/drivers/taylor.c @ 48

Last change on this file since 48 was 48, checked in by awalther, 10 years ago

bug in taylor.c for computation of derivatives of inverse functions fixed

  • Property svn:keywords set to Author Date Id Revision
File size: 26.9 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     drivers/taylor.c
4 Revision: $Id: taylor.c 48 2009-08-31 14:37:29Z awalther $
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 <drivers/taylor.h>
16#include <interfaces.h>
17#include <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
121    if (k > n)
122        return 0;
123
124    if (k > n/2)
125        k = n-k;
126
127    long double accum = 1;
128    unsigned int i;
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 (%le)\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, 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        for (j=0; j<d; j++)
567        {
568            Xhelp[i][j] = 0;
569            X[i][j+1] = 0;
570            W[i][j] = 0;
571        }
572
573    while (--ii > 0) {
574        di = dd[ii-1]-1;
575        Di = dd[ii-1]-dd[ii]-1;
576        MINDEC(rc,hos_forward(tag,n,n,di,Di+1,xold,Xhelp,w,W));
577        MINDEC(rc,hov_reverse(tag,n,n,Di,n,I,A,nonzero));
578        da = dd[ii];
579        for (l=da; l<dd[ii-1]; l++) {
580            for (i=0; i<n; i++) {
581                if (l == 0)
582                    bi = w[i]-Y[i][0];
583                else
584                    bi = W[i][l-1]-Y[i][l];
585                for (j=0; j<n; j++)
586                    if (nonzero[i][j]>1) {
587                        Aij = A[i][j];
588                        int indj = j;
589                        indexA = l-1;
590                        Xj = X[j]+l;
591                        int indX = l;
592                        indexX = 1;
593                        if (da == l-1)
594                          {
595                            bi += (*(++Aij))*(*(--Xj));
596                          }
597                        else
598                          {
599                            for (q=da; q<l; q++)
600                              {
601                                bi += (*(++Aij))*(*(--Xj));
602                                bi += A[i][j][indexA]*X[j][indexX];
603                                indexA--;
604                                indexX++;
605                              }
606                          }
607                    }
608                b[i] = -bi;
609            }
610            MINDEC(rc,jac_solv(tag,n,xold,b,2));
611           if (rc == -3)
612                return -3;
613            for (i=0; i<n; i++) {
614              X[i][l] += b[i];
615              Xhelp[i][l-1] += b[i];
616            }
617        }
618    }
619    return rc;
620}
621
622
623/****************************************************************************/
624int inverse_tensor_eval( short tag, int n, int d, int p,
625                         double *x, double **tensor, double** S ) {
626    static int dim;
627    static int dold,pold;
628    static struct item *coeff_list;
629    int i,j,dimten;
630    int *it = (int*) malloc(d*sizeof(int));
631    double** X;
632    double** Y;
633    int *jm;
634    double *y = (double*) malloc(n*sizeof(double));
635    struct item *ptr;
636    int rc = 3;
637
638    dimten=binomi(p+d,d);
639    for(i=0;i<n;i++)
640        for(j=0;j<dimten;j++)
641            tensor[i][j] = 0;
642    MINDEC(rc,zos_forward(1,n,n,0,x,y));
643    if (d > 0) {
644        if ((d != dold) || (p != pold)) {
645            if (pold) { /* olvo 980728 */
646                dim = binomi(pold+dold-1,dold);
647                freecoefflist(dim,coeff_list);
648                free((char*) coeff_list);
649            }
650            dim = binomi(p+d-1,d);
651            coeff_list = (struct item *) malloc(sizeof(struct item)*dim);
652            coeff(p,d, coeff_list);
653            dold = d;
654            pold = p;
655        }
656        jm = (int *)malloc(sizeof(int)*p);
657        X = myalloc2(n,d+1);
658        Y = myalloc2(n,d+1);
659        for (i=0; i<n; i++) {
660            X[i][0] = x[i];
661            Y[i][0] = y[i];
662            for (j=1; j<d; j++)
663              {
664                X[i][j] = 0;
665                Y[i][j] = 0;
666              }
667        }
668        if (d == 1) {
669            it[0] = 0;
670            for (i=0; i<dim; i++)  /* sum over all multiindices jm with |jm| = d */
671            { it[0] = it[0]+1;
672                convert(p,d,it,jm);
673                ptr = &coeff_list[i];
674                multma2vec1(n,p,d,Y,S,jm);
675                MINDEC(rc,inverse_Taylor_prop(tag,n,d,Y,X));
676                if (rc == -3)
677                    return -3;
678                do {
679                    for(j=0;j<n;j++)
680                        tensor[j][ptr->a] += X[j][ptr->b]*ptr->c;
681                    ptr = ptr->next;
682                } while (ptr != NULL);
683            }
684        } else {
685            for (i=0; i<d-1; i++)
686                it[i] = 1;
687            it[d-1] = 0;
688            for (i=0; i<dim; i++)  /* sum over all multiindices jm with |jm| = d */
689            { it[d-1] = it[d-1]+1;
690                for (j=d-2; j>=0; j--)
691                    it[j] = it[j] + it[j+1]/(p+1);
692                for (j=1; j<d; j++)
693                    if (it[j] > p) it[j] = it[j-1];
694                convert(p,d,it,jm);
695                multma2vec1(n,p,d,Y,S,jm); /* Store S*jm in Y */
696                MINDEC(rc,inverse_Taylor_prop(tag,n,d,Y,X));
697                if (rc == -3)
698                    return -3;
699                ptr = &coeff_list[i];
700                do {
701                    for(j=0;j<n;j++)
702                      {
703                        tensor[j][ptr->a] += X[j][ptr->b]*ptr->c;
704                      }
705                    ptr = ptr->next;
706                } while (ptr != NULL);
707            }
708        }
709        free((char*) jm);
710        free((char*) *X);
711        free((char*) X);
712        free((char*) *Y);
713        free((char*) Y);
714    }
715    for(i=0;i<n;i++)
716        tensor[i][0] = x[i];
717    free((char*) y);
718    free((char*) it);
719    return rc;
720}
721
722
723/****************************************************************************/
724int tensor_eval( short tag, int m, int n, int d, int p,
725                 double* x, double **tensor, double **S ) {
726    static int bd,dim;
727    static int dold,pold;
728    static struct item *coeff_list;
729    int i,j,k,dimten,ctr;
730    int **jm, jmbd=0;
731    int *it = (int*) malloc(d*sizeof(int));
732    double *y = (double*) malloc(m*sizeof(double));
733    double*** X;
734    double*** Y;
735    struct item *ptr[10];
736    int rc = 3;
737
738    dimten=binomi(p+d,d);
739    for (i=0; i<m; i++)
740        for (j=0; j<dimten; j++)
741            tensor[i][j] = 0;
742    if (d == 0) {
743        MINDEC(rc,zos_forward(1,m,n,0,x,y));
744    } else {
745        if ((d != dold) || (p != pold)) {
746            if (pold) {
747                dim = binomi(pold+dold-1,dold);
748                freecoefflist(dim,coeff_list);
749                free((char*) coeff_list);
750            }
751            dim = binomi(p+d-1,d);
752            if (dim < 10)
753                bd = dim;
754            else
755                bd = 10;
756            coeff_list = (struct item *) malloc(sizeof(struct item)*dim);
757            coeff(p,d, coeff_list);
758            dold = d;
759            pold = p;
760        }
761        jmbd = bd;
762        jm = (int **) malloc(jmbd*sizeof(int*));
763        for (i=0; i<jmbd; i++)
764            jm[i] = (int *) malloc(p*sizeof(int));
765        if (d == 1) {
766            X = myalloc3(1,n,bd);
767            Y = myalloc3(1,m,bd);
768            ctr   = 0;
769            it[0] = 0;
770            for (i=0; i<dim; i++) /* sum over all multiindices jm with |jm| = d */
771            { it[0] = it[0]+1;
772                convert(p,d,it,jm[ctr]);
773                ptr[ctr] = &coeff_list[i];
774                if (ctr < bd-1)
775                    ctr += 1;
776                else {
777                    multma2vec2(n,p,bd,X[0],S,jm);
778                    MINDEC(rc,fov_forward(tag,m,n,bd,x,X[0],y,Y[0]));
779                    for (k=0; k<bd; k++)
780                        do {
781                            for (j=0; j<m; j++)
782                                tensor[j][ptr[k]->a] += Y[0][j][k]*ptr[k]->c;
783                            ptr[k] = ptr[k]->next;
784                        } while (ptr[k] != NULL);
785                    if (dim-i < bd)
786                        bd = dim-i-1;
787                    ctr = 0;
788                }
789            }
790        } else {
791            X = myalloc3(n,bd,d);
792            Y = myalloc3(m,bd,d);
793            ctr = 0;
794            for (i=0; i<d-1; i++)
795                it[i] = 1;
796            it[d-1] = 0;
797            for (i=0; i<dim; i++) /* sum over all multiindices jm with |jm| = d */
798            { it[d-1] = it[d-1]+1;
799                for (j=d-2; j>=0; j--)
800                    it[j] = it[j] + it[j+1]/(p+1);
801                for (j=1; j<d; j++)
802                    if (it[j] > p)
803                        it[j] = it[j-1];
804                convert(p,d,it,jm[ctr]);
805                ptr[ctr] = &coeff_list[i];
806                if (ctr < bd-1)
807                    ctr += 1;
808                else {
809                    multma3vec2(n,p,d,bd,X,S,jm);
810                    MINDEC(rc,hov_forward(tag,m,n,d,bd,x,X,y,Y));
811                    for (k=0; k<bd; k++)
812                        do {
813                            for (j=0; j<m; j++)
814                                tensor[j][ptr[k]->a] += Y[j][k][ptr[k]->b-1]*ptr[k]->c;
815                            ptr[k] = ptr[k]->next;
816                        } while (ptr[k] != NULL);
817                    if (dim-i < bd)
818                        bd = dim-i-1;
819                    ctr = 0;
820                }
821            }
822        }
823        for (i=0; i<jmbd; i++)
824            free((char*) *(jm+i));
825        free((char*) jm);
826        free((char*) **X);
827        free((char*) *X);
828        free((char*) X);
829        free((char*) **Y);
830        free((char*) *Y);
831        free((char*) Y);
832    }
833    for(i=0;i<m;i++)
834        tensor[i][0] = y[i];
835    bd = jmbd;
836    free((char*) y);
837    free((char*) it);
838    return rc;
839}
840
841
842/****************************************************************************/
843void tensor_value( int d, int m, double *y, double **tensor, int *multi ) {
844    int i, j, max, ind, add;
845    int *im = (int*) malloc(d*sizeof(int));
846
847    max = 0;
848    ind = d-1;
849    for (i=0; i<d; i++) {
850        if (multi[i] > max)
851            max = multi[i];
852        im[i] = 0;
853    }
854    for (i=0; i<d; i++) {
855        if (multi[i] == max)  /* olvo 980728 == instead of = */
856        { im[ind] = multi[i];
857            multi[i] = 0;
858            max = 0;
859            ind -= 1;
860            for (j=0; j<d; j++)
861                if (multi[j] > max)
862                    max = multi[j];
863        }
864    }
865    add = address(d,im);
866    for (i=0; i<m; i++)
867        y[i] = tensor[i][add];
868    free((char*) im);
869}
870
871END_C_DECLS
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
Note: See TracBrowser for help on using the repository browser.