!********************************************************************
!* 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<N)
      P0=C(k,k)
      l0=k; k0=k
      do i=k, N
        do j=k, N
          if (ABS(C(i,j)) > 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