|
c1 = 1.0/bta/dt/dt
c2 = 1.0/bta/dt
c3 = gma/bta/dt
c4 = 1.0 - gma/bta
c5 = 1.0 - 0.5*gma/bta
c6 = 0.5/bta - 1.0
c7 = c1*M_m + c3*M_c
c8 = c2*M_m - c4*M_c
c9 = c6*M_m - dt*c5*M_c
do i=1,length-1,1
ph = -ME*ag(i+1) + matmul(c7,M_u(i,:)) + matmul(c8,M_v(i,:)) + matmul(c9,M_a(i,:))
temp=M_kh
temp1=ph
call dgetrf(ndof,ndof,temp,ndof,ipiv,info)
call dgetrs("N",ndof,ndof,temp,ndof,ipiv,ph,ndof,info)
M_u(i+1,:)=ph
M_v(i+1,:) = c3*(M_u(i+1,:)-M_u(i,:)) + c4*M_v(i,:) + dt*c5*M_a(i,:)
M_a(i+1,:) = c1*(M_u(i+1,:)-M_u(i,:)) - c2*M_v(i,:) - c6*M_a(i,:)
end do
在编写newmark方法中,执行MKL的dgetrs后,发现c9的数值发生了变化。如将dgetrs函数中的ph改为temp1,则一切正常。请大家帮忙分析下是何原因,感谢!!
|
|