FORTRAN求解实数方阵的逆矩阵

PROGRAMMatrix_inv!定义动态数组REAL,ALLOCATABLE::X(:,:),L(:,:),U(:,:),y(:),z(:),inv_X(:,:)!定义矩阵的维数INTEGERN,I,JPRINT*,'请输入方阵的维数N'READ*,NALLOCATE(X(N,N))ALLOCATE(U(N,N))ALLOCATE(L(N,N))ALLOCATE(y(N))ALLOCATE(z(N))ALLOCATE(inv_X(N,N))!输入矩阵XPRINT*,'请确保矩阵可逆!'PRINT*,'请输入矩阵X',N,'行',N,'列'DOI=1,NPRINT*,'请输入第',I,'行'READ*,X(I,:)ENDDO!将矩阵X进行LU分解CALLMatrix_LU(X,N,L,U)!求解逆矩阵DOI=1,NCALLcal_y(L,y,N,I)CALLcal_z(U,y,z,N)DOJ=1,Ninv_X(J,I)=z(J)ENDDOENDDO!打印结果DOI=1,NPRINT*,inv_X(I,:)ENDDO!======================================================================!!各子程序!!======================================================================!!子程序1=======LU分解CONTAINSSUBROUTINEMatrix_LU(X,N,L,U)REALL(N,N),U(N,N),X(N,N)INTEGERI,N,J,K,MDOI=1,NU(1,I)=X(1,I)L(I,1)=X(I,1)/U(1,1)ENDDO!求出U的第K行和U的第K列(K=2,3...N)DOK=2,NDOJ=K,NU(K,J)=X(K,J)DOM=1,K-1U(K,J)=U(K,J)-L(K,M)*U(M,J)ENDDOENDDOIF(N==2)THENEXITENDIFDOI=K+1,NL(I,K)=X(I,K)DOM=1,K-1L(I,K)=L(I,K)-L(I,M)*U(M,K)ENDDOL(I,K)=L(I,K)/U(K,K)ENDDOENDDO!L矩阵的主对角线全为1DOI=1,NL(I,I)=1ENDDO!L为下三角,U为上三角DOI=1,NDOJ=1,NIF(IJ)THENL(I,J)=0ENDIFENDDOENDDODOI=1,NDOJ=1,NIF(IJ)THENU(I,J)=0ENDIFENDDOENDDOENDSUBROUTINEMatrix_LU!子程序2=======解Ly=b方程SUBROUTINEcal_y(L,y,N,K)REALy(N),L(N,N),B(N)INTEGERI,J,N,KDOI=1,NB(I)=0ENDDOB(K)=1y(1)=B(1)DOI=2,Ny(I)=B(I)DOJ=1,I-1y(I)=y(I)-L(I,J)*y(J)ENDDOENDDOENDSUBROUTINEcal_y!子程序3=======解Uz=y方程SUBROUTINEcal_z(U,y,z,N)REALy(N),z(N),U(N,N)INTEGERI,J,Nz(N)=y(N)/U(N,N)DOI=N-1,1z(I)=y(I)DOJ=I+1,Nz(I)=z(I)-U(I,J)*z(J)ENDDOz(I)=z(I)/U(I,I)ENDDOENDSUBROUTINEcal_zEND