!******************************************************************** !* Solve a Complex Linear System By Gauss Method with full pivoting * !* and correction process. * !* ---------------------------------------------------------------- * !* SAMPLE RUN: * !* Solve Complex Linear System AX = B with: * !* * !* ( 47,-15) ( 62,5) ( 0,-72) (61, 20) ( 629,988) * !* A = ( 6, 14) (-17,3) (-102,91) ( 7,-12) and B = (-180,825) * !* ( 13, 55) ( 32,8) ( 41, 7) (25, 1) ( 877,441) * !* (111,25) ( 40,0) ( 12,-82) (58,-30) ( 734,-88) * !* * !* System solutions: * !* ( -1.0000, 3.0000) * !* ( 2.0000, 10.0000) * !* ( 7.0000, -5.0000) * !* ( 17.0000, 6.0000) * !* * !* ---------------------------------------------------------------- * !* Reference: "Algebre, Algorithmes et programmes en Pascal * !* By Jean-Louis Jardrin, DUNOD Paris, 1988". * !* * !* F90 Release By J-P MOreau, Paris. * !* (www.jpmoreau.fr) * !******************************************************************** Program CSysLin integer, parameter :: NMAX = 20 Complex A(NMAX,NMAX) Complex B(NMAX), X(NMAX) N = 4 !Size of system A(1,1) = CMPLX(47.0,-15.0) A(1,2) = CMPLX(62.0, 5.0) A(1,3) = CMPLX( 0.0,-72.0) A(1,4) = CMPLX(61.0, 20.0) A(2,1) = CMPLX( 6.0, 14.0) A(2,2) = CMPLX( -17.0, 3.0) A(2,3) = CMPLX(-102.0, 91.0) A(2,4) = CMPLX( 7.0,-12.0) A(3,1) = CMPLX(13.0,-55.0) A(3,2) = CMPLX(32.0, 8.0) A(3,3) = CMPLX(41.0, 7.0) A(3,4) = CMPLX(25.0, 1.0) A(4,1) = CMPLX(111.0, 25.0) A(4,2) = CMPLX( 40.0, 0.0) A(4,3) = CMPLX( 12.0,-82.0) A(4,4) = CMPLX( 58.0,-30.0) B(1) = CMPLX( 629.0,988.0) B(2) = CMPLX(-180.0,825.0) B(3) = CMPLX( 877.0,441.0) B(4) = CMPLX( 734.0,-88.0) eps=1.E-10; dta=1.E-10; M=1 call RSLCGTC(eps,dta,M,N,A,B,it,X) print *,' ' print *,' System solutions:' do I=1, N write(*,10) X(I) end do print *,' ' stop 10 format(' (',F8.4,',',F8.4,')') END Subroutine CSwap(Z1, Z2) COMPLEX Z1, Z2, C C=Z1; Z1=Z2; Z2=C return 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 * !* (modified in the process). * !* 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 !************************************************************* !* Solve a Complex Linear System AX = B By Gauss Method with * !* full pivoting and a correction process. * !* --------------------------------------------------------- * !* INPUTS: * !* eps, dta : absolute and relative precisions * !* M : maximum number of iterations * !* N : size of linear system * !* A : complex matrix * !* B : right-hand side (complex vector) * !* OUTPUTS: * !* it: flag, =-1 if no convergence, =0 if matrix A * !* is singular, =1 convergence ok. * !* X : contains system solution (X1,X2,...,Xn) * !* * !************************************************************* Subroutine RSLCGTC(eps, dta, M, N, A, B, it, X) integer, parameter :: NMAX = 20 Complex A(NMAX,NMAX), B(NMAX), X(NMAX) COMPLEX C0,C1,S,Z0; Complex C(NMAX,NMAX), B1(NMAX), W(NMAX), Z(NMAX) integer KP(NMAX), LP(NMAX) Z0 = CMPLX(0.0,0.0) call TSCGT(eps,N,A,it,C,KP,LP) if (it==1) then ! Save B in B1 before BSCGT do J=1, N B1(J) = B(J) end do call BSCGT(N,C,B,KP,LP,X) ! Restore B after BSCGT do J=1, N B(J) = B1(J) end do it=-1; L=1 do while (it==-1.and.L<=M) phi1=0.0 do I=1, N if (ABS(X(I)) > phi1) phi1=ABS(X(I)) end do do I=1, N S=Z0 do J=1, N C0=S C1 = A(I,J) * X(J) S = C0 + C1 end do W(I) = B(I) - S end do call BSCGT(N,C,W,KP,LP,Z) do I=1, N C0=X(I) X(I) = C0 + Z(I) end do phi2=0.0 do I=1, N if (ABS(Z(I)) > phi2) phi2=ABS(Z(I)) end do if (phi2/phi1 < dta) then it=1 else L = L + 1 end if ; end do end if return End ! RSLCGTC ! end of file trslcgtc.f90