integral_dke_jd

C/C++ source

0001 #include <stdio.h>
0002 #include <stdlib.h>
0003 #include "Fokker_Planck_DKE_loop.h"
0004 #include <math.h>
0005 
0006 
0007 void integral_dke_jd(struct Tensor B, struct Tensor A,int d,struct Tensor *C){
0008 int i,j=0,k;
0009 /*A Tensor,B vector to contract with, C output*/
0010 
0011 /*char* file="/home/sccp/gsem/FT090307/Project_DKE/fred/Csource/A";
0012 char* file2="/home/sccp/gsem/FT090307/Project_DKE/fred/Csource/B";*/
0013 d-=1; /*for Consistency with Matlab indices starting at 1*/
0014 
0015 /* printf("start integral\n");*/
0016 if(A.sizedim[d]!=B.sizedim[1]){
0017 printf("dimension mismatch\n");
0018 printf("sizeA=%i\n",A.sizedim[d]);
0019 printf("sizeB=%i\n",B.sizedim[1]);
0020 return;}
0021 else{
0022 /* printf("sizeA=%i\n",A.sizedim[d]);*/
0023 }
0024 /*printf("sizeB=%i\n",B.sizedim[1]);*/
0025 int *depth=malloc((A.ndims+1)*sizeof(int));
0026 int ndims;
0027 int *sizedim;
0028 for(k=0;k<=A.ndims;k++){
0029 depth[k]=1;
0030    for(i=0;i/*printf("compd depth\n");*/
0034 
0035 if(A.ndims !=2){
0036    ndims=A.ndims-1;
0037    sizedim=malloc(ndims*sizeof(int));
0038    for(i=0;i<(*C).ndims;i++){
0039       if(i==d){j++;}
0040          *(sizedim+i)=*(A.sizedim+j);
0041    j++;
0042    }
0043 }
0044 
0045 /*Vectors are coded as (n,1) matrices*/
0046 else {
0047    ndims=2;
0048    sizedim=malloc(ndims*sizeof(int));
0049    sizedim[0]=1;
0050    sizedim[1]=A.sizedim[abs(d-1)];
0051 }
0052 Tensor_Create(ndims,sizedim,C);
0053 /*printf("computing C\n");*/
0054 for(i=0;ifor(j=0;jfor(k=0;k/*fprintf(stdout,"C.ndims=%i\n",(*C).ndims);
0064 for(i=0;i<(*C).ndims;i++){fprintf(stdout,"%i\n",(*C).sizedim[i]);}*/
0065 /* printf("integral done\n");*/
0066 }