!*************************************************************************** !* Inversion of a complex square matrix by LU decomposition * !* with dynamic allocations * !* * !* F90 version by J-P Moreau, Paris * !* (www.jpmoreau.fr) * !* ----------------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* Input file (inv_clu.dat): * !* * !* 4 * !* 8 0 2 0 3 0 12 1 * !* 2 0 4 0 7 0 0.25 2 * !* 3 0 7 0 3 0 5 3 * !* 12 0 0.25 0 5 0 2 4 * !* * !* Output file (inv_clu.lst): * !* * !* INVERSION OF A COMPLEX SQUARE MATRIX: * !* * !* N = 4 * !* * !* 8.00000 0.00000 2.00000 0.00000 3.00000 0.00000 12.00000 1.00000 * !* 2.00000 0.00000 4.00000 0.00000 7.00000 0.00000 0.25000 2.00000 * !* 3.00000 0.00000 7.00000 0.00000 3.00000 0.00000 5.00000 3.00000 * !* 12.00000 0.00000 0.25000 0.00000 5.00000 0.00000 2.00000 4.00000 * !* * !* Inverted matrix Y: * !* * !* -0.03022 -0.04162 -0.08426 -0.00641 0.05306 0.01466 0.10426 0.02515 * !* -0.07421 -0.04847 -0.06024 -0.00746 0.19813 0.01707 0.00998 0.02929 * !* 0.05443 0.00890 0.20187 0.00137 -0.12957 -0.00313 -0.03754 -0.00538 * !* 0.10431 0.02491 0.01606 0.00384 -0.03673 -0.00877 -0.06304 -0.01505 * !* * !* Verification A*Y = I: * !* * !* 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 * !* 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 * !* 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 * !* 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 * !* * !* Uses: module CLU. * !*************************************************************************** Program INVERSE_LU USE CLU COMPLEX, pointer :: A(:,:) !real matrix (n x n) COMPLEX, pointer :: A1(:,:) !copy of matrix A COMPLEX, pointer :: Y(:,:) !real matrix (n x n) COMPLEX, pointer :: temp(:) !real temporary vector (n) integer d, rc character*12 input, output character*8 s ! input data file name from screen print *,' ' write(*,50,advance='no'); read *, s J=0 DO I=1,LEN(s) IF(s(I:I)<>' ') J=J+1 !J=real length of string s ENDDO input=s(1:J)//'.dat' output =s(1:J)//'.lst' ! open input and output files open(unit=1,file=input,status='old') open(unit=2,file=output,status='unknown') read(1,*) n !size of given linear system !dynamic allocations allocate(A(n,n),stat=ialloc) allocate(A1(n,n),stat=ialloc) allocate(Y(n,n),stat=ialloc) allocate(temp(n),stat=ialloc) !Read data in file #1 and copy to file #2 write(2,*) ' ' write(2,*) ' INVERSION OF A COMPLEX SQUARE MATRIX:' write(2,60) n do i=1, n !read one line call ReadCVec(1,n,temp) do j=1, n A(i,j) = temp(j) A1(i,j) = temp(j) Y(i,j) = 0.d0 end do Y(i,i) = 1.d0 !print one line of matrix A call WriteCVec(2,n,temp) end do close(1) !call LU inversion routine call CINVERT_LU(A,n,Y,rc) !the inverse matrix is now in matrix Y !the original matrix A is destroyed !print results or error message if (rc.eq.1) then write(2,*) ' The matrix is singular, no solution !' else write(2,*) ' ' write(2,*) ' Inverted matrix Y:' write(2,*) call WriteCMat(2,n, n, Y) end if !verify A1 x Y = I (result put in A) call CMATMUL(A1,Y,A,n,n) !A should now contain identity matrix write(2,*) ' ' write(2,*) ' Verification A*Y = I:' write(2,*) call WriteCMat(2,n, n, A) ! section fin Close(2) print *,' ' print *,'Results in file ',output,'.' print *,' ' stop 50 format(' Input data file name (without .dat): ') 60 format(/' N = ',I2/) END !************************************************* ! Read a complex vector from unit u from index = 1 !************************************************* Subroutine ReadCVec(u, n, x) INTEGER u COMPLEX x(1:n) REAL tmp(1:2*n) read(u,*) (tmp(i),i=1,2*n) j=1 do i=1, n x(i) = CMPLX(tmp(j),tmp(j+1)) j=j+2 end do return End Subroutine WriteCMat(u,n, m, a) !*======================================================================* !* * !* Put out a complex m x n matrix in output file ( 10 items per line ) * !* (index begins at one) * !* * !* -------------------------------------------------------------------- * !* * !* Input parameters: * !* ================ * !* u integer u; # of unit to read * !* n int m; ( m > 0 ) * !* row number of matrix * !* m int n; ( n > 0 ) * !* column number of matrix * !* a COMPLEX a(n,n) * !* matrix to write * !* * !*======================================================================* INTEGER u COMPLEX a(n,m) do i=1, n write(u,10) (a(i,j),j=1,m) enddo return 10 format(8F9.5) End Subroutine WriteCVec(u, n, x) INTEGER u, n COMPLEX x(1:n) !*====================================================================* !* * !* Put out complex vector x of length n to output file (10 items per* !* line). Index starts from ONE. * !* * !*====================================================================* !* * !* Input parameters: * !* ================ * !* * !* n integer n; * !* lenght of vector * !* x COMPLEX x(n) * !* vector * !* * !*====================================================================* write(u,10) (x(i),i=1,n) return 10 format(8F9.5) End SUBROUTINE CMATMUL(A,B,C,N,M) !******************************************* !* MULTIPLY TWO COMPLEX MATRICES * !* --------------------------------------- * !* INPUTS: A MATRIX N*N * !* B MATRIX N*M * !* N INTEGER * !* M INTEGER * !* --------------------------------------- * !* OUTPUT: C MATRIX N*M, PRODUCT A*B * !* * !******************************************* COMPLEX A(N,N),B(N,M),C(N,M),SUM DO I=1,N DO J=1,M SUM=0. DO K=1,N SUM=SUM+A(I,K)*B(K,J) ENDDO C(I,J)=SUM ENDDO ENDDO RETURN END ! end of file tinv_clu.f90