!**************************************************************** !* calculate inverse of a complex square matrix by Gauss Method * !* with full pivoting. * !* ------------------------------------------------------------ * !* SAMPLE RUN: * !* Calculate inverse of complex square matrix: * !* * !* | (0, 1) ( 2,0) (-1,-1) | * !* A = | (1, 0) ( 1,1) ( 0,-3) | * !* | (0,-2) (-1,1) ( 3, 0) | * !* * !* Inverse matrix: * !* ( 0.0000, 0.0000) ( 0.3333, 0.0000) ( 0.0000, 0.3333) * !* ( 0.7500, 0.0000) (-0.1667,-0.0833) ( 0.3333, 0.0833) * !* ( 0.2500,-0.2500) (-0.0833, 0.2500) ( 0.2500,-0.0833) * !* * !* Product AM1*A: * !* ( 1.0000, 0.0000) ( 0.0000, 0.0000) ( 0.0000, 0.0000) * !* ( 0.0000, 0.0000) ( 1.0000, 0.0000) ( 0.0000, 0.0000) * !* ( 0.0000, 0.0000) ( 0.0000, 0.0000) ( 1.0000, 0.0000) * !* * !* ------------------------------------------------------------ * !* Ref.: "Algèbre, Algorithmes et programmes en Pascal * !* By Jean-Louis Jardrin, DUNOD Paris, 1988". * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !**************************************************************** Program InvMatC; integer, parameter :: NMAX = 20 integer i,it,j,m, N real eps Complex A(NMAX,NMAX), AM1(NMAX,NMAX), MI(NMAX,NMAX) N=3 !size of complex matrix A A(1,1) = CMPLX( 0.0, 1.0) A(1,2) = CMPLX( 2.0, 0.0) A(1,3) = CMPLX(-1.0,-1.0) A(2,1) = CMPLX(1.0, 0.0) A(2,2) = CMPLX(1.0, 1.0) A(2,3) = CMPLX(0.0,-3.0) A(3,1) = CMPLX( 0.0,-2.0) A(3,2) = CMPLX(-1.0, 1.0) A(3,3) = CMPLX( 3.0, 0.0) eps=1E-10 call ICGT(eps,N,A,it,AM1) print *,' ' print *,' Inverse matrix:' do i=1, N write(*,10) (AM1(i,j),j=1,N) end do ! Check AM1 * A = I call MATMUL(AM1, A, MI, N) print *,' ' print *,' Product AM1*A:' do i=1, N write(*,10) (MI(i,j),j=1,N) end do print *,' ' stop 10 format(3(' (',F7.4,',',F7.4,')')) END !***************************************************************** !* TSCGT procedure implements the triangularization algorithm of * !* Gauss with full pivoting at each step for a complex matrix, A * !* and saves the made transformations in KP and LP. * !* ------------------------------------------------------------- * !* INPUTS: * !* N: size of complex matrix A * !* A: complex matrix of size N x N * !* OUTPUTS; * !* it: =0 if A is singular, else =1. * !* C: contains the upper triangular matrix and the * !* multipliers used during the process. * !* KP: contains the column exchanges. * !* LP: contains the line exchanges. * !***************************************************************** Subroutine TSCGT(eps, N, A, it, C, KP, LP) integer, parameter :: NMAX = 20 Complex A(NMAX,NMAX), C(NMAX,NMAX) integer KP(NMAX), LP(NMAX) Complex C0,C1,P0 C = A it=1; K=1; do while (it==1.and.k ABS(P0)) then P0=C(i,j) l0=i; k0=j end if end do end do LP(k)=l0; KP(k)=k0 if (ABS(P0) < eps) then it=0 else if (l0<>k) then do j=k, N !Swap(C(k,j),C(l0,j)) C0=C(l0,j); C(l0,j)=C(k,j); C(k,j)=C0 end do end if if (k0<>k) then do i=1, N !Swap(C(i,k),C(i,k0)) C0=C(i,k0); C(i,k0)=C(i,k); C(i,k)=C0 end do end if do i=k+1, N C0=C(i,k) C(i,k) = C0/P0 do j=k+1, N C0=C(i,j) C1 = C(i,k) * C(k,j) C(i,j) = C0 - C1 end do end do k = k + 1 end if end do; if (it==1.and.ABS(C(N,N)) < eps) it=0 return End !TSCGT !************************************************************** !* Function BSCGT calculates the solution of upper triangular * !* system [S(n-1)] by the back substitution method and deduces* !* from it the solution of system [S]. The call must be made * !* only after a call to TSCGT and the matrix of [S] must be * !* regular. * !* ---------------------------------------------------------- * !* INPUTS: * !* N : size of matrix C * !* C: contains the upper triangular matrix and the * !* multipliers used during the triangularization * !* process (in output of TSCGT). . * !* W : contains the right-side coefficients * !* KP: contains the column exchanges. * !* LP: contains the line exchanges. * !* OUTPUT: * !* Z : system solution complex vector. * !************************************************************** Subroutine BSCGT(N, C, W, KP, LP, Z) integer, parameter :: NMAX = 20 Complex C(NMAX,NMAX), W(NMAX), Z(NMAX) integer KP(NMAX), LP(NMAX) Complex C0,C1,S,Z0 Z0=CMPLX(0.0,0.0) do k=1, N-1 l0=LP(k) if (l0.ne.k) then !Swap(W(k),W(l0)) C0=W(l0); W(l0)=W(k); W(k)=C0 end if do i=k+1, N C0=W(i) C1 = C(i,k) * W(k) W(i) = C0 - c1 end do end do Z(N) = W(N)/C(N,N) do i=N-1, 1, -1 S=Z0 do j=i+1, N C0=S C1 = C(i,j) * Z(j) S = C0 + C1 end do C0 = W(i) - S Z(i) = C0/C(i,i) end do do k=N-1, 1, -1 k0=KP(k) if (k0.ne.k) then !Swap(Z(k),Z(k0)) C0=Z(k0); Z(k0)=Z(k); Z(k)=C0 end if end do return End !BSCGT !**************************************************************** !* calculate inverse of a complex square matrix by Gauss Method * !* with full pivoting. * !* ------------------------------------------------------------ * !* INPUTS: * !* eps : required precision * !* N : size of complex matrix A * !* A : input complex matrix (NxN) * !* OUTPUTS: * !* it : =0, matrix is singular; =1, matris is regular * !* AM1 : inverse complex matrix (NxN). * !**************************************************************** Subroutine ICGT(eps, N, A, it, AM1) integer, parameter :: NMAX = 20 Complex A(NMAX,NMAX), AM1(NMAX,NMAX) Complex Z0, Z1 Complex C(NMAX,NMAX), W(NMAX), Z(NMAX) integer KP(NMAX), LP(NMAX) Z0=CMPLX(0.0, 0.0) Z1=CMPLX(1.0, 0.0) call TSCGT(eps,N,A,it,C,KP,LP) if (it==1) then do l=1, N do i=1, N if (i==l) then W(I)=Z1 else W(I)=Z0 end if end do call BSCGT(N,C,W,KP,LP,Z) do i=1, N AM1(i,l)=Z(i) end do end do end if return End !ICGT Subroutine MATMUL(A, B, C, N) !******************************************* !* MULTIPLICATION OF TWO SQUARE COMPLEX * !* MATRICES * !* --------------------------------------- * !* INPUTS: A MATRIX N*N * !* B MATRIX N*N * !* N INTEGER * !* --------------------------------------- * !* OUTPUTS: C MATRIX N*N PRODUCT A*B * !* * !******************************************* integer, parameter :: NMAX = 20 Complex A(NMAX,NMAX), B(NMAX,NMAX), C(NMAX,NMAX) Complex SUM, PROD do I=1, N do J=1, N SUM=CMPLX(0.0,0.0) do K=1, N SUM=SUM+A(I,K)*B(K,J) end do C(I,J)=SUM end do end do return END ! end of file Invmatc.f90