Civil Engineering Reference
In-Depth Information
!dir$ ivdep
elements_2a: DO iel=1,nels
km=storkm(:,:,iel); g_utemp=MATMUL(km,g_pmul)
u(g_g(:,iel))=u(g_g(:,iel))+g_utemp(:,iel)
END DO elements_2a
IF(fixed_freedoms/=0)u(no)=p(no)*store; up=DOT_PRODUCT(loads,d)
alpha=up/DOT_PRODUCT(p,u); xnew=x+p*alpha; loads=loads-u*alpha
d=diag_precon*loads; beta=DOT_PRODUCT(loads,d)/up; p=d+p*beta
big=zero; cg_converged=.TRUE.
DO i=1,neq; IF(ABS(xnew(i))>big)big=ABS(xnew(i)); END DO
DO i=1,neq; IF(ABS(xnew(i)-x(i))/big>cg_tol)cg_converged=.FALSE.; END DO
x=xnew; IF(cg_converged.OR.cg_iters==cg_limit)EXIT
END DO pcg
WRITE(11,'(A,I5)')" Number of cg iterations to convergence was",cg_iters
WRITE(11,'(/A)')" Node x-disp y-disp z-disp"; loads=xnew
DO k=1,nn; WRITE(11,'(I5,3E12.4)')k,loads(nf(:,k)); END DO
!-----------------------recover stresses at nip integrating point---------
nip=1; DEALLOCATE(points,weights); ALLOCATE(points(nip,ndim),weights(nip))
CALL sample(element,points,weights); loads(0)=zero
WRITE(11,'(/A,I2,A)')" The integration point (nip=",nip,") stresses are:"
WRITE(11,'(A,/,A)')"
Element
x-coord
y-coord
z-coord",
&
"
sig_x
sig_y
sig_z
tau_xy
tau_yz
tau_zx"
elements_4: DO iel=1,nels
CALL deemat(dee,prop(1,etype(iel)),prop(2,etype(iel))); num=g_num(:,iel)
coord=TRANSPOSE(g_coord(:,num)); g=g_g(:,iel); eld=loads(g)
gauss_pts_2: DO i=1,nip
CALL shape_der(der,points,i); CALL shape_fun(fun,points,i)
gc=MATMUL(fun,coord); jac=MATMUL(der,coord); CALL invert(jac)
deriv=MATMUL(jac,der); CALL beemat(bee,deriv)
sigma=MATMUL(dee,MATMUL(bee,eld)); WRITE(11,'(I8,4X,3E12.4)')iel,gc
WRITE(11,'(6E12.4)')sigma
END DO gauss_pts_2
END DO elements_4
STOP
END PROGRAM p56
New dynamic real arrays:
g pmul
“gathers” the
p
vectors for all elements
g utemp
holds all the products of
km
and
p
This program is used to solve exactly the same problem as was detailed for Programs 5.3
and 5.5, using the mesh-free strategy of Program 5.5. However, it is used as an example of
some of the issues which arise when programming for vector computers (see Section 1.4).
All practical vector computers enable code to be analysed to see where most time is
being used and where there are features of the program inhibiting most effective use of
vectorising compilers.
When Program 5.5 was run through such an analysis program, three main points arose:
a) There is potential “dependency” in the scatter operation
u(g)=u(g)+MATMUL(km,p(g))
and the compiler does not know whether there may be repeated entries in
g
and so
does not vectorise this statement.
