1*8980d4a7Sjeremylt!----------------------------------------------------------------------- 2*8980d4a7Sjeremylt subroutine polyeval(x,n,p,uq) 3*8980d4a7Sjeremylt real*8 x,y 4*8980d4a7Sjeremylt integer n,i 5*8980d4a7Sjeremylt real*8 p(1) 6*8980d4a7Sjeremylt real*8 uq 7*8980d4a7Sjeremylt 8*8980d4a7Sjeremylt y=p(n) 9*8980d4a7Sjeremylt 10*8980d4a7Sjeremylt do i=n-1,1,-1 11*8980d4a7Sjeremylt y=y*x+p(i) 12*8980d4a7Sjeremylt enddo 13*8980d4a7Sjeremylt 14*8980d4a7Sjeremylt uq=y 15*8980d4a7Sjeremylt 16*8980d4a7Sjeremylt end 17*8980d4a7Sjeremylt!----------------------------------------------------------------------- 18*8980d4a7Sjeremylt program test 19*8980d4a7Sjeremylt 20*8980d4a7Sjeremylt include 'ceedf.h' 21*8980d4a7Sjeremylt 22*8980d4a7Sjeremylt integer ceed,err 23*8980d4a7Sjeremylt integer x,xq,u,uq 24*8980d4a7Sjeremylt integer bxl,bul,bxg,bug 25*8980d4a7Sjeremylt integer i 26*8980d4a7Sjeremylt integer q 27*8980d4a7Sjeremylt parameter(q=6) 28*8980d4a7Sjeremylt 29*8980d4a7Sjeremylt real*8 p(6) 30*8980d4a7Sjeremylt real*8 xx(2) 31*8980d4a7Sjeremylt real*8 xxq(q) 32*8980d4a7Sjeremylt real*8 uuq(q) 33*8980d4a7Sjeremylt real*8 px 34*8980d4a7Sjeremylt integer*8 offset1,offset2 35*8980d4a7Sjeremylt 36*8980d4a7Sjeremylt character arg*32 37*8980d4a7Sjeremylt 38*8980d4a7Sjeremylt data p/1,2,3,4,5,6/ 39*8980d4a7Sjeremylt data xx/-1,1/ 40*8980d4a7Sjeremylt 41*8980d4a7Sjeremylt call getarg(1,arg) 42*8980d4a7Sjeremylt call ceedinit(trim(arg)//char(0),ceed,err) 43*8980d4a7Sjeremylt 44*8980d4a7Sjeremylt call ceedvectorcreate(ceed,2,x,err) 45*8980d4a7Sjeremylt call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,xx,err) 46*8980d4a7Sjeremylt call ceedvectorcreate(ceed,q,xq,err) 47*8980d4a7Sjeremylt call ceedvectorsetvalue(xq,0.d0,err) 48*8980d4a7Sjeremylt call ceedvectorcreate(ceed,q,u,err) 49*8980d4a7Sjeremylt call ceedvectorsetvalue(u,0.d0,err) 50*8980d4a7Sjeremylt call ceedvectorcreate(ceed,q,uq,err) 51*8980d4a7Sjeremylt call ceedvectorsetvalue(uq,0.d0,err) 52*8980d4a7Sjeremylt 53*8980d4a7Sjeremylt call ceedbasiscreatetensorh1lagrange(ceed,1,1,2,q,ceed_gauss_lobatto,& 54*8980d4a7Sjeremylt & bxl,err) 55*8980d4a7Sjeremylt call ceedbasiscreatetensorh1lagrange(ceed,1,1,q,q,ceed_gauss_lobatto,& 56*8980d4a7Sjeremylt & bul,err) 57*8980d4a7Sjeremylt 58*8980d4a7Sjeremylt call ceedbasisapply(bxl,1,ceed_notranspose,ceed_eval_interp,x,xq,err) 59*8980d4a7Sjeremylt 60*8980d4a7Sjeremylt call ceedvectorgetarrayread(xq,ceed_mem_host,xxq,offset1,err) 61*8980d4a7Sjeremylt do i=1,q 62*8980d4a7Sjeremylt call polyeval(xxq(i+offset1),6,p,uuq(i)) 63*8980d4a7Sjeremylt enddo 64*8980d4a7Sjeremylt call ceedvectorrestorearrayread(xq,xxq,offset1,err) 65*8980d4a7Sjeremylt call ceedvectorsetarray(uq,ceed_mem_host,ceed_use_pointer,uuq,err) 66*8980d4a7Sjeremylt 67*8980d4a7Sjeremylt call ceedbasisapply(bul,1,ceed_transpose,ceed_eval_interp,uq,u,err) 68*8980d4a7Sjeremylt 69*8980d4a7Sjeremylt call ceedbasiscreatetensorh1lagrange(ceed,1,1,2,q,ceed_gauss,bxg,err) 70*8980d4a7Sjeremylt call ceedbasiscreatetensorh1lagrange(ceed,1,1,q,q,ceed_gauss,bug,err) 71*8980d4a7Sjeremylt 72*8980d4a7Sjeremylt call ceedbasisapply(bxg,1,ceed_notranspose,ceed_eval_interp,x,xq,err) 73*8980d4a7Sjeremylt call ceedbasisapply(bug,1,ceed_notranspose,ceed_eval_interp,u,uq,err) 74*8980d4a7Sjeremylt 75*8980d4a7Sjeremylt call ceedvectorgetarrayread(xq,ceed_mem_host,xxq,offset1,err) 76*8980d4a7Sjeremylt call ceedvectorgetarrayread(uq,ceed_mem_host,uuq,offset2,err) 77*8980d4a7Sjeremylt do i=1,q 78*8980d4a7Sjeremylt call polyeval(xxq(i+offset1),6,p,px) 79*8980d4a7Sjeremylt if (abs(uuq(i+offset2)-px) > 1e-14) then 80*8980d4a7Sjeremylt write(*,*) uuq(i+offset2),' not eqaul to ',px,'=p(',xxq(i+offset1),')' 81*8980d4a7Sjeremylt endif 82*8980d4a7Sjeremylt enddo 83*8980d4a7Sjeremylt call ceedvectorrestorearrayread(xq,xxq,offset1,err) 84*8980d4a7Sjeremylt call ceedvectorrestorearrayread(uq,uuq,offest2,err) 85*8980d4a7Sjeremylt 86*8980d4a7Sjeremylt call ceedvectordestroy(x,err) 87*8980d4a7Sjeremylt call ceedvectordestroy(xq,err) 88*8980d4a7Sjeremylt call ceedvectordestroy(u,err) 89*8980d4a7Sjeremylt call ceedvectordestroy(uq,err) 90*8980d4a7Sjeremylt call ceedbasisdestroy(bxl,err) 91*8980d4a7Sjeremylt call ceedbasisdestroy(bul,err) 92*8980d4a7Sjeremylt call ceedbasisdestroy(bxg,err) 93*8980d4a7Sjeremylt call ceedbasisdestroy(bug,err) 94*8980d4a7Sjeremylt call ceeddestroy(ceed,err) 95*8980d4a7Sjeremylt end 96*8980d4a7Sjeremylt!----------------------------------------------------------------------- 97