Civil Engineering Reference
In-Depth Information
CALL deemat(dee,prop(1,etype(iel)),prop(2,etype(iel)))
num=g_num(:,iel); coord=TRANSPOSE(g_coord(:,num))
g=g_g(:,iel); km=zero; mm=zero
integrating_pts_1: DO i=1,nip
CALL shape_fun(fun,points,i); CALL shape_der(der,points,i)
jac=MATMUL(der,coord); det=determinant(jac); CALL invert(jac)
deriv=MATMUL(jac,der); call beemat(bee,deriv)
km=km+MATMUL(MATMUL(TRANSPOSE(bee),dee),bee)*det*weights(i)
CALL ecmat(ecm,fun,ndof,nodof)
mm=mm+ecm*det*weights(i)*prop(3,etype(iel))
END DO integrating_pts_1
DO i=1,ndof; diag(g(i))=diag(g(i))+SUM(mm(i,:)); END DO
END DO elements_2
!------------------------------find eigenvalues---------------------------
diag=one/SQRT(diag); diag(0)=zero ! diag holds l**(-1/2)
DO iters=1,lalfa
CALL lancz1(neq,el,er,acc,leig,lx,lalfa,lp,iflag,ua,va,eig,jeig,neig,x, &
del,nu,alfa,beta,v_store)
IF(iflag==0)EXIT
IF(iflag>1)THEN
WRITE(11,'(A,I5)')" Lancz1 is signalling failure, with iflag = ",iflag
STOP
END IF
!----- iflag = 1 therefore form u + a * v ( done element by element )----
vdiag=va; vdiag=vdiag*diag ! vdiag is l**(-1/2).va
udiag=zero; vdiag(0)=zero
elements_3: DO iel=1,nels
g=g_g(:,iel); udiag(g)=udiag(g)+MATMUL(km,vdiag(g))
END DO elements_3
udiag=udiag*diag; ua=ua+udiag
END DO
!-------------- iflag = 0 therefore write out the spectrum ---------------
WRITE(11,'(A,I4,A/)')" It took ",iters," iterations"
WRITE(11,'(3(A,E12.4))')" Eigenvalues in the range",el," to",er," are:"
WRITE(11,'(6E12.4)')eig(1:neig)
! calculate the eigenvectors
IF(neig>10)neig=10
CALL lancz2(neq,lalfa,lp,eig,jeig,neig,alfa,beta,lz,jflag,y,w1,z,v_store)
!------------------if jflag is zero calculate the eigenvectors ----------
IF(jflag==0)THEN
DO i=1,nmodes
udiag(:)=y(:,i)
udiag=udiag*diag
WRITE(11,'(" Eigenvector number",I4," is:")')i
WRITE(11,'(6E12.4)')udiag(1:)/MAXVAL(udiag(1:))
END DO
ELSE
! lancz2 fails
WRITE(11,'(A,I5)')" Lancz2 is signalling failure with jflag = ",jflag
END IF
STOP
END PROGRAM p104
New dynamic real arrays:
vdiag
used in element-by-element products
