1!----------------------------------------------------------------------- 2 subroutine eval(dimn,x,rslt) 3 integer dimn 4 real*8 x(3) 5 real*8 rslt 6 7 rslt=tanh(x(1)+0.1) 8 if (dimn>1) then 9 rslt=rslt+atan(x(2)+0.2) 10 endif 11 if (dimn>2) then 12 rslt=rslt+exp(-(x(3)+0.3)*(x(3)+0.3)) 13 endif 14 15 end 16!----------------------------------------------------------------------- 17 program test 18 implicit none 19 include 'ceedf.h' 20 21 integer ceed,err 22 integer x,xq,u,uq,ones,gtposeones 23 integer bxl,bug 24 integer dimn,d 25 integer i 26 integer p 27 integer q 28 parameter(p=8) 29 parameter(q=10) 30 integer maxdim 31 parameter(maxdim=3) 32 integer qdimnmax 33 parameter(qdimnmax=q**maxdim) 34 integer pdimnmax 35 parameter(pdimnmax=p**maxdim) 36 integer xdimmax 37 parameter(xdimmax=2**maxdim) 38 integer pdimn,qdimn,xdim 39 40 real*8 xx(xdimmax*maxdim) 41 real*8 xxx(maxdim) 42 real*8 xxq(pdimnmax*maxdim) 43 real*8 uuq(qdimnmax*maxdim) 44 real*8 uu(pdimnmax) 45 real*8 ggtposeones(pdimnmax) 46 real*8 sum1 47 real*8 sum2 48 integer dimxqdimn 49 integer*8 uoffset,xoffset,offset1,offset2,offset3 50 51 character arg*32 52 53 call getarg(1,arg) 54 call ceedinit(trim(arg)//char(0),ceed,err) 55 56 do dimn=1,maxdim 57 qdimn=q**dimn 58 pdimn=p**dimn 59 xdim=2**dimn 60 dimxqdimn=dimn*qdimn 61 sum1=0 62 sum2=0 63 64 do d=0,dimn-1 65 do i=1,xdim 66 if ((mod(i-1,2**(dimn-d))/(2**(dimn-d-1))).ne.0) then 67 xx(d*xdim+i)=1 68 else 69 xx(d*xdim+i)=-1 70 endif 71 enddo 72 enddo 73 74 call ceedvectorcreate(ceed,xdim*dimn,x,err) 75 xoffset=0 76 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,xx,xoffset,err) 77 call ceedvectorcreate(ceed,pdimn*dimn,xq,err) 78 call ceedvectorsetvalue(xq,0.d0,err) 79 call ceedvectorcreate(ceed,pdimn,u,err) 80 call ceedvectorcreate(ceed,qdimn*dimn,uq,err) 81 call ceedvectorsetvalue(uq,0.d0,err) 82 call ceedvectorcreate(ceed,qdimn*dimn,ones,err) 83 call ceedvectorsetvalue(ones,1.d0,err) 84 call ceedvectorcreate(ceed,pdimn,gtposeones,err) 85 call ceedvectorsetvalue(gtposeones,0.d0,err) 86 87 call ceedbasiscreatetensorh1lagrange(ceed,dimn,dimn,2,p,& 88 & ceed_gauss_lobatto,bxl,err) 89 call ceedbasisapply(bxl,1,ceed_notranspose,ceed_eval_interp,x,xq,err) 90 91 call ceedvectorgetarrayread(xq,ceed_mem_host,xxq,offset1,err) 92 do i=1,pdimn 93 do d=0,dimn-1 94 xxx(d+1)=xxq(d*pdimn+i+offset1) 95 enddo 96 call eval(dimn,xxx,uu(i)) 97 enddo 98 call ceedvectorrestorearrayread(xq,xxq,offset1,err) 99 uoffset=0 100 call ceedvectorsetarray(u,ceed_mem_host,ceed_use_pointer,uu,uoffset,err) 101 102 call ceedbasiscreatetensorh1lagrange(ceed,dimn,1,p,q,ceed_gauss,bug,err) 103 104 call ceedbasisapply(bug,1,ceed_notranspose,ceed_eval_grad,u,uq,err) 105 call ceedbasisapply(bug,1,ceed_transpose,ceed_eval_grad,ones,& 106 & gtposeones,err) 107 108 call ceedvectorgetarrayread(gtposeones,ceed_mem_host,ggtposeones,& 109 & offset1,err) 110 call ceedvectorgetarrayread(u,ceed_mem_host,uu,offset2,err) 111 call ceedvectorgetarrayread(uq,ceed_mem_host,uuq,offset3,err) 112 do i=1,pdimn 113 sum1=sum1+ggtposeones(i+offset1)*uu(i+offset2) 114 enddo 115 do i=1,dimxqdimn 116 sum2=sum2+uuq(i+offset3) 117 enddo 118 call ceedvectorrestorearrayread(gtposeones,ggtposeones,offset1,err) 119 call ceedvectorrestorearrayread(u,uu,offset2,err) 120 call ceedvectorrestorearrayread(uq,uuq,offset3,err) 121 if(abs(sum1-sum2) > 1.0D-10) then 122! LCOV_EXCL_START 123 write(*,'(A,I1,A,F12.6,A,F12.6)')'[',dimn,'] Error: ',sum1, ' != ',& 124 & sum2 125! LCOV_EXCL_STOP 126 endif 127 128 call ceedvectordestroy(x,err) 129 call ceedvectordestroy(xq,err) 130 call ceedvectordestroy(u,err) 131 call ceedvectordestroy(uq,err) 132 call ceedvectordestroy(ones,err) 133 call ceedvectordestroy(gtposeones,err) 134 call ceedbasisdestroy(bxl,err) 135 call ceedbasisdestroy(bug,err) 136 enddo 137 138 call ceeddestroy(ceed,err) 139 end 140!----------------------------------------------------------------------- 141