1*752c3701SJeremy L Thompson !----------------------------------------------------------------------- 2*752c3701SJeremy L Thompson subroutine setup_mass(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,& 3*752c3701SJeremy L Thompson & u14,u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,& 4*752c3701SJeremy L Thompson & v16,ierr) 5*752c3701SJeremy L Thompson real*8 ctx 6*752c3701SJeremy L Thompson real*8 u1(1) 7*752c3701SJeremy L Thompson real*8 u2(1) 8*752c3701SJeremy L Thompson real*8 v1(1) 9*752c3701SJeremy L Thompson integer q,ierr 10*752c3701SJeremy L Thompson 11*752c3701SJeremy L Thompson do i=1,q 12*752c3701SJeremy L Thompson v1(i)=u2(i)*(u1(i+q*0)*u1(i+q*3)-u1(i+q*1)*u1(i+q*2)) 13*752c3701SJeremy L Thompson enddo 14*752c3701SJeremy L Thompson 15*752c3701SJeremy L Thompson ierr=0 16*752c3701SJeremy L Thompson end 17*752c3701SJeremy L Thompson !----------------------------------------------------------------------- 18*752c3701SJeremy L Thompson subroutine setup_diff(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,& 19*752c3701SJeremy L Thompson & u14,u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,& 20*752c3701SJeremy L Thompson & v16,ierr) 21*752c3701SJeremy L Thompson real*8 ctx 22*752c3701SJeremy L Thompson real*8 u1(1) 23*752c3701SJeremy L Thompson real*8 u2(1) 24*752c3701SJeremy L Thompson real*8 v1(1) 25*752c3701SJeremy L Thompson real*8 w 26*752c3701SJeremy L Thompson integer q,ierr 27*752c3701SJeremy L Thompson 28*752c3701SJeremy L Thompson do i=1,q 29*752c3701SJeremy L Thompson w=u2(i)/(u1(i+q*0)*u1(i+q*3)-u1(i+q*1)*u1(i+q*2)) 30*752c3701SJeremy L Thompson v1(i+q*0)=w*(u1(i+q*2)*u1(i+q*2)+u1(i+q*3)*u1(i+q*3)) 31*752c3701SJeremy L Thompson v1(i+q*1)=w*(u1(i+q*0)*u1(i+q*0)+u1(i+q*1)*u1(i+q*1)) 32*752c3701SJeremy L Thompson v1(i+q*2)=-w*(u1(i+q*0)*u1(i+q*2)+u1(i+q*2)*u1(i+q*3)) 33*752c3701SJeremy L Thompson enddo 34*752c3701SJeremy L Thompson 35*752c3701SJeremy L Thompson ierr=0 36*752c3701SJeremy L Thompson end 37*752c3701SJeremy L Thompson !----------------------------------------------------------------------- 38*752c3701SJeremy L Thompson subroutine apply(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 39*752c3701SJeremy L Thompson & u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 40*752c3701SJeremy L Thompson real*8 ctx 41*752c3701SJeremy L Thompson real*8 u1(1) 42*752c3701SJeremy L Thompson real*8 u2(1) 43*752c3701SJeremy L Thompson real*8 u3(1) 44*752c3701SJeremy L Thompson real*8 u4(1) 45*752c3701SJeremy L Thompson real*8 v1(1) 46*752c3701SJeremy L Thompson real*8 v2(1) 47*752c3701SJeremy L Thompson real*8 du0,du1 48*752c3701SJeremy L Thompson integer q,ierr 49*752c3701SJeremy L Thompson 50*752c3701SJeremy L Thompson do i=1,q 51*752c3701SJeremy L Thompson ! mass 52*752c3701SJeremy L Thompson v1(i) = u2(i)*u4(i) 53*752c3701SJeremy L Thompson ! diff 54*752c3701SJeremy L Thompson du0=u1(i+q*0) 55*752c3701SJeremy L Thompson du1=u1(i+q*1) 56*752c3701SJeremy L Thompson v2(i+q*0)=u3(i+q*0)*du0+u3(i+q*2)*du1 57*752c3701SJeremy L Thompson v2(i+q*1)=u3(i+q*2)*du0+u3(i+q*1)*du1 58*752c3701SJeremy L Thompson enddo 59*752c3701SJeremy L Thompson 60*752c3701SJeremy L Thompson ierr=0 61*752c3701SJeremy L Thompson end 62*752c3701SJeremy L Thompson !----------------------------------------------------------------------- 63*752c3701SJeremy L Thompson subroutine apply_lin(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,& 64*752c3701SJeremy L Thompson & u14,u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,& 65*752c3701SJeremy L Thompson & v16,ierr) 66*752c3701SJeremy L Thompson real*8 ctx 67*752c3701SJeremy L Thompson real*8 u1(1) 68*752c3701SJeremy L Thompson real*8 u2(1) 69*752c3701SJeremy L Thompson real*8 u3(1) 70*752c3701SJeremy L Thompson real*8 v1(1) 71*752c3701SJeremy L Thompson real*8 v2(1) 72*752c3701SJeremy L Thompson real*8 du0,du1 73*752c3701SJeremy L Thompson integer q,ierr 74*752c3701SJeremy L Thompson 75*752c3701SJeremy L Thompson do i=1,q 76*752c3701SJeremy L Thompson du0=u1(i+q*0) 77*752c3701SJeremy L Thompson du1=u1(i+q*1) 78*752c3701SJeremy L Thompson v1(i+q*0)=u2(i+q*0)*du0+u2(i+q*3)*du1+u2(i+q*6)*u3(i) 79*752c3701SJeremy L Thompson v2(i+q*0)=u2(i+q*1)*du0+u2(i+q*4)*du1+u2(i+q*7)*u3(i) 80*752c3701SJeremy L Thompson v2(i+q*1)=u2(i+q*2)*du0+u2(i+q*5)*du1+u2(i+q*8)*u3(i) 81*752c3701SJeremy L Thompson enddo 82*752c3701SJeremy L Thompson 83*752c3701SJeremy L Thompson ierr=0 84*752c3701SJeremy L Thompson end 85*752c3701SJeremy L Thompson !----------------------------------------------------------------------- 86