1*52bfb9bbSJeremy L Thompson!----------------------------------------------------------------------- 2*52bfb9bbSJeremy L Thompson! 3*52bfb9bbSJeremy L Thompson! Header with common subroutine 4*52bfb9bbSJeremy L Thompson! 5*52bfb9bbSJeremy L Thompson include 't320-basis-f.h' 6*52bfb9bbSJeremy L Thompson!----------------------------------------------------------------------- 7*52bfb9bbSJeremy L Thompson subroutine feval(x1,x2,val) 8*52bfb9bbSJeremy L Thompson real*8 x1,x2,val 9*52bfb9bbSJeremy L Thompson 10*52bfb9bbSJeremy L Thompson val=x1*x1+x2*x2+x1*x2+1 11*52bfb9bbSJeremy L Thompson 12*52bfb9bbSJeremy L Thompson end 13*52bfb9bbSJeremy L Thompson!----------------------------------------------------------------------- 14*52bfb9bbSJeremy L Thompson subroutine dfeval(x1,x2,val) 15*52bfb9bbSJeremy L Thompson real*8 x1,x2,val 16*52bfb9bbSJeremy L Thompson 17*52bfb9bbSJeremy L Thompson val=2*x1+x2 18*52bfb9bbSJeremy L Thompson 19*52bfb9bbSJeremy L Thompson end 20*52bfb9bbSJeremy L Thompson!----------------------------------------------------------------------- 21*52bfb9bbSJeremy L Thompson program test 22*52bfb9bbSJeremy L Thompson 23*52bfb9bbSJeremy L Thompson include 'ceedf.h' 24*52bfb9bbSJeremy L Thompson 25*52bfb9bbSJeremy L Thompson integer ceed,err 26*52bfb9bbSJeremy L Thompson integer input,output 27*52bfb9bbSJeremy L Thompson integer p,q,d 28*52bfb9bbSJeremy L Thompson parameter(p=6) 29*52bfb9bbSJeremy L Thompson parameter(q=4) 30*52bfb9bbSJeremy L Thompson parameter(d=2) 31*52bfb9bbSJeremy L Thompson 32*52bfb9bbSJeremy L Thompson real*8 qref(d*q) 33*52bfb9bbSJeremy L Thompson real*8 qweight(q) 34*52bfb9bbSJeremy L Thompson real*8 interp(p*q) 35*52bfb9bbSJeremy L Thompson real*8 grad(d*p*q) 36*52bfb9bbSJeremy L Thompson real*8 xq(d*q) 37*52bfb9bbSJeremy L Thompson real*8 xr(d*p) 38*52bfb9bbSJeremy L Thompson real*8 iinput(p) 39*52bfb9bbSJeremy L Thompson real*8 ooutput(d*q) 40*52bfb9bbSJeremy L Thompson real*8 val,diff 41*52bfb9bbSJeremy L Thompson real*8 x1,x2 42*52bfb9bbSJeremy L Thompson integer*8 ioffset,ooffset 43*52bfb9bbSJeremy L Thompson 44*52bfb9bbSJeremy L Thompson integer b 45*52bfb9bbSJeremy L Thompson 46*52bfb9bbSJeremy L Thompson character arg*32 47*52bfb9bbSJeremy L Thompson 48*52bfb9bbSJeremy L Thompson xq=(/2.d-1,6.d-1,1.d0/3.d0,2.d-1,2.d-1,2.d-1, 1.d0/3.d0,6.d-1/) 49*52bfb9bbSJeremy L Thompson xr=(/0.d0,5.d-1,1.d0,0.d0,5.d-1,0.d0,0.d0,0.d0, 0.d0,5.d-1,5.d-1,1.d0/) 50*52bfb9bbSJeremy L Thompson 51*52bfb9bbSJeremy L Thompson call getarg(1,arg) 52*52bfb9bbSJeremy L Thompson 53*52bfb9bbSJeremy L Thompson call buildmats(qref,qweight,interp,grad) 54*52bfb9bbSJeremy L Thompson 55*52bfb9bbSJeremy L Thompson call ceedinit(trim(arg)//char(0),ceed,err) 56*52bfb9bbSJeremy L Thompson 57*52bfb9bbSJeremy L Thompson call ceedbasiscreateh1(ceed,ceed_triangle,1,p,q,interp,grad,qref,qweight,& 58*52bfb9bbSJeremy L Thompson & b,err) 59*52bfb9bbSJeremy L Thompson 60*52bfb9bbSJeremy L Thompson do i=1,p 61*52bfb9bbSJeremy L Thompson x1=xr(0*p+i) 62*52bfb9bbSJeremy L Thompson x2=xr(1*p+i) 63*52bfb9bbSJeremy L Thompson call feval(x1,x2,val) 64*52bfb9bbSJeremy L Thompson iinput(i)=val 65*52bfb9bbSJeremy L Thompson enddo 66*52bfb9bbSJeremy L Thompson 67*52bfb9bbSJeremy L Thompson call ceedvectorcreate(ceed,p,input,err) 68*52bfb9bbSJeremy L Thompson ioffset=0 69*52bfb9bbSJeremy L Thompson call ceedvectorsetarray(input,ceed_mem_host,ceed_use_pointer,iinput,& 70*52bfb9bbSJeremy L Thompson & ioffset,err) 71*52bfb9bbSJeremy L Thompson call ceedvectorcreate(ceed,q*d,output,err) 72*52bfb9bbSJeremy L Thompson call ceedvectorsetvalue(output,0.d0,err) 73*52bfb9bbSJeremy L Thompson 74*52bfb9bbSJeremy L Thompson call ceedbasisapply(b,1,ceed_notranspose,ceed_eval_grad,input,output,err) 75*52bfb9bbSJeremy L Thompson 76*52bfb9bbSJeremy L Thompson call ceedvectorgetarrayread(output,ceed_mem_host,ooutput,ooffset,err) 77*52bfb9bbSJeremy L Thompson do i=1,q 78*52bfb9bbSJeremy L Thompson x1=xq(0*q+i) 79*52bfb9bbSJeremy L Thompson x2=xq(1*q+i) 80*52bfb9bbSJeremy L Thompson call dfeval(x1,x2,val) 81*52bfb9bbSJeremy L Thompson diff=val-ooutput(0*q+i+ooffset) 82*52bfb9bbSJeremy L Thompson if (abs(diff)>1.0d-10) then 83*52bfb9bbSJeremy L Thompson! LCOV_EXCL_START 84*52bfb9bbSJeremy L Thompson write(*,'(A,I1,A,F12.8,A,F12.8)') '[',i,'] ',ooutput(i+ooffset),& 85*52bfb9bbSJeremy L Thompson & ' != ',val 86*52bfb9bbSJeremy L Thompson! LCOV_EXCL_STOP 87*52bfb9bbSJeremy L Thompson endif 88*52bfb9bbSJeremy L Thompson call dfeval(x2,x1,val) 89*52bfb9bbSJeremy L Thompson diff=val-ooutput(1*q+i+ooffset) 90*52bfb9bbSJeremy L Thompson if (abs(diff)>1.0d-10) then 91*52bfb9bbSJeremy L Thompson! LCOV_EXCL_START 92*52bfb9bbSJeremy L Thompson write(*,'(A,I1,A,F12.8,A,F12.8)') '[',i,'] ',ooutput(i+ooffset),& 93*52bfb9bbSJeremy L Thompson & ' != ',val 94*52bfb9bbSJeremy L Thompson! LCOV_EXCL_STOP 95*52bfb9bbSJeremy L Thompson endif 96*52bfb9bbSJeremy L Thompson enddo 97*52bfb9bbSJeremy L Thompson call ceedvectorrestorearrayread(output,ooutput,ooffset,err) 98*52bfb9bbSJeremy L Thompson 99*52bfb9bbSJeremy L Thompson call ceedvectordestroy(input,err) 100*52bfb9bbSJeremy L Thompson call ceedvectordestroy(output,err) 101*52bfb9bbSJeremy L Thompson call ceedbasisdestroy(b,err) 102*52bfb9bbSJeremy L Thompson call ceeddestroy(ceed,err) 103*52bfb9bbSJeremy L Thompson 104*52bfb9bbSJeremy L Thompson end 105*52bfb9bbSJeremy L Thompson!----------------------------------------------------------------------- 106