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