编辑: 阿拉蕾 | 2017-09-16 |
D0)) real (rkind), parameter :: Zero=0.D0,One=1.D0,Two=2.D0,Three=3.D0, &
&
Four=4.D0,Five=5.D0,Six=6.D0,Seven=7.D0,Eight=8.D0,Nine=9.D0, &
&
Ten=10.D0 contains function matinv(A) result (B) real(rkind) ,intent (in)::A(:,:) !real(rkind) , allocatable::B(:,:) real(rkind) , pointer::B(:,:) integer(ikind):: N,I,J,K real(rkind)::D,T real(rkind), allocatable::IS(:),JS(:) N=size(A,dim=2) allocate(B(N,N)) allocate(IS(N));
allocate(JS(N)) B=A do K=1,N D=0.0D0 do I=K,N do J=K,N if(abs(B(I,J))>
D) then D=abs(B(I,J)) IS(K)=I JS(K)=J end if end do end do do J=1,N T=B(K,J) B(K,J)=B(IS(K),J) B(IS(K),J)=T end do do I=1,N T=B(I,K) B(I,K)=B(I,JS(K)) B(I,JS(K))=T end do B(K,K)=1/B(K,K) do J=1,N if(J.NE.K) then B(K,J)=B(K,J)*B(K,K) end if end do do I=1,N if(I.NE.K) then do J=1,N if(J.NE.K) then B(I,J)=B(I,J)-B(I,K)*B(K,J) end if end do end if end do do I=1,N if(I.NE.K) then B(I,K)=-B(I,K)*B(K,K) end if end do end do do K=N,1,-1 do J=1,N T=B(K,J) B(K,J)=B(JS(K),J) B(JS(K),J)=T end do do I=1,N T=B(I,K) B(I,K)=B(I,IS(K)) B(I,IS(K))=T end do end do return end function matinv subroutine IntSwap(a,b) integer(ikind),intent(in out)::a,b integer(ikind)::t t=a;
a=b;
b=t end subroutine IntSwap subroutine RealSwap(a,b) real(rkind),intent(in out)::a,b real(rkind)::t t=a;
a=b;
b=t end subroutine RealSwap subroutine matprint(A,n) real(rkind),intent(in)::A(:,:) integer(ikind)::n integer(ikind)::n1,n2 integer(ikind)::i,j character(10)::C n1=size(A,dim=1) n2=size(A,dim=2) C='
('
//trim(itoc(n2))//'
E'
//trim(itoc(n))//&
'
.'
//trim(itoc(n-7))//'
)'
do I=1,n1 write(*,C)(A(I,J),J=1,n2) end do end subroutine matprint function matdet(B) result(det) real(rkind),intent(in)::B(:,:) real(rkind)::det integer(ikind)::n,i,j,k,is,js real(rkind),pointer::A(:,:) real(rkind)::f,d,q n=size(B,dim=1) allocate (A(n,n)) A=B f=1.0D0;
det=1.0D0 do k=1,n-1 q=0.0D0 do i=k,n do j=k,n if(abs(a(i,j)).gt.q) then q=abs(a(i,j)) is=i js=j end if end do end do if(q+1.0D0.eq.1.0D0) then det=0.0d0 return end if if(is.ne.k) then f=-f do j=k,n d=a(k,j) a(k,j)=a(is,j) a(is,j)=d end do end if if(js.ne.k) then f=-f do i=k,n d=a(i,js) a(i,js)=a(i,k) a(i,k)=d end do end if det=det*a(k,k) do i=k+1,n d=a(i,k)/a(k,k) do j=k+1,n a(i,j)=a(i,j)-d*a(k,j) end do end do end do det=f*det*a(n,n) deallocate (a) return end function matdet function itoc(i1) result (c) integer(ikind),intent(in)::i1 character(len=2)::c real(rkind)::x integer(ikind) ::n,b,i,j i=i1 x=i c(1:2)='
'
x=log10(x) n=int(x)+2 do j=n-2,0,-1 b=mod(i,10**j) b=(i-b)/(10**j) i=i-b*(10**j) c(n-j-1:n-j-1)=achar(iachar('
0'
)+b) end do end function itoc subroutine Gauss(GStif,GLoad,GDisp) real (rkind),intent (in) :: GStif(:,:),GLoad(:) real (rkind),intent (out) :: GDisp(:) integer (ikind) :: i,j,k integer (ikind) :: N real (rkind) :: P,I1,X,Y real (rkind),allocatable :: A(:,:) N=size(GDisp,dim=1) allocate (A(N,N+1)) A(1:N,1:N)=GStif(1:N,1:N) A(1:N,N+1)=GLoad(1:N) DO j=1,N P=0.0D0 DO k=j,N IF(ABS(A(k,j)).LE.P) cycle P=ABS(A(k,j)) I1=k end do IF(P.GE.1E-15)GO TO
230 WRITE(22,'
(A)'
) '
NO UNIQUE SOLUTION'
RETURN
230 IF(I1.EQ.j)GO TO
280 DO
270 K=J,N+1 X=A(J,K) A(J,K)=A(I1,K)
270 A(I1,K)=X
280 Y=1.D0/A(J,J) DO
310 K=J,N+1
310 A(J,K)=Y*A(J,K) DO
380 I=1,N IF(I.EQ.J)GO TO
380 Y=-A(I,J) DO
370 K=J,N+1
370 A(I,K)=A(I,K)+Y*A(J,K)
380 CONTINUE