C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
C     File:        ZGAUSSJDRV.F (Fortran 77 source for the ZGAUSSDRV prog)
C     Author:      Fredrik Jonsson <fj@phys.soton.ac.uk>
C     Date:        December 29, 2005
C     Last change: December 29, 2006
C     Description: The ZGAUSSJDRV program is simply a driver for the
C                  ZGAUSSJ routine for solving complex-valued systems by
C                  Gauss-Jordan elimination.  See the documentation of
C                  the ZGAUSSJ routine, enclosed in file ZGAUSSJ.F, for
C                  further information.
C
C     Copyright (C) 2006, Fredrik Jonsson
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
C     PROGRAM MAIN
      IMPLICIT LOGICAL (A-Z)
      EXTERNAL zgaussj
      INTEGER J,K,L
      COMPLEX*8 A(3,3),AINV(3,3),T(3,3),B(3),X(3),D(3)

C     *****************************************************************
C     Define the complex-valued test matrix A of the system A*x=b.
C     *****************************************************************
      A(1,1)=CMPLX(1.4d0,3.3d0)
      A(1,2)=CMPLX(-2.3d0,1.7d0)
      A(1,3)=CMPLX(2.1d0,2.9d0)
      A(2,1)=CMPLX(-3.4d0,2.3d0)
      A(2,2)=CMPLX(1.3d0,-1.7d0)
      A(2,3)=CMPLX(-3.2d0,-1.4d0)
      A(3,1)=CMPLX(1.3d0,-2.3d0)
      A(3,2)=CMPLX(2.8d0,1.1d0)
      A(3,3)=CMPLX(-1.7d0,-1.2d0)

C     *****************************************************************
C     Define the complex-valued right-hand side b of the system A*x=b.
C     *****************************************************************
      B(1)=1.3
      B(2)=-2.1
      B(3)=0.7

C     *****************************************************************
C     Calculate the matrix inverse inv(A) using the zgaussj routine.
C     *****************************************************************
      DO 100,J=1,3
         X(J)=B(J)
         DO 110,K=1,3
            AINV(J,K)=A(J,K)
 110     CONTINUE
 100  CONTINUE
      CALL zgaussj(AINV,3,3,X,1,1)

C     *****************************************************************
C     Display the complex-valued system matrix A.
C     *****************************************************************
      WRITE (*,*) 'A:'
      DO 130,J=1,3
         WRITE(6,30) A(J,1),A(J,2),A(J,3)
 130  CONTINUE

C     *****************************************************************
C     Display the inverted matrix inv(A).
C     *****************************************************************
      WRITE (*,*) ' '
      WRITE (*,*) 'inv(A):'
      DO 140,J=1,3
         WRITE(6,30) AINV(J,1),AINV(J,2),AINV(J,3)
 140  CONTINUE

C     *****************************************************************
C     As a check on the validity of the output of the zgaussj routine,
C     verify that the matrix A times its inverse inv(A) yields the
C     identity matrix.
C     *****************************************************************
      DO 150,J=1,3
         DO 160,K=1,3
            T(J,K)=CMPLX(0.0,0.0)
            DO 170,L=1,3
               T(J,K)=T(J,K)+A(J,L)*AINV(L,K)
 170        CONTINUE
 160     CONTINUE
 150  CONTINUE
      WRITE (*,*) ' '
      WRITE (*,*) 'A*inv(A):'
      DO 180,J=1,3
         WRITE(6,30) T(J,1),T(J,2),T(J,3)
 180  CONTINUE

C     *****************************************************************
C     Also perform the analogous check of inv(a)*A.
C     *****************************************************************
      DO 190,J=1,3
         DO 200,K=1,3
            T(J,K)=CMPLX(0.0,0.0)
            DO 210,L=1,3
               T(J,K)=T(J,K)+AINV(J,L)*A(L,K)
 210        CONTINUE
 200     CONTINUE
 190  CONTINUE
      WRITE (*,*) ' '
      WRITE (*,*) 'inv(A)*A:'
      DO 220,J=1,3
         WRITE(6,30) T(J,1),T(J,2),T(J,3)
 220  CONTINUE

C     *****************************************************************
C     Check the obtained solution x of A*x=b by computing the
C     difference d=inv(A)*b-x.
C     *****************************************************************
      DO 230,J=1,3
         D(J)=-X(J)
         DO 240,K=1,3
            D(J)=D(J)+AINV(J,K)*B(K)
 240     CONTINUE
 230  CONTINUE
      WRITE (*,*) ' '
      WRITE (*,*) 'inv(A)*b-x:'
      WRITE(6,30) D(1),D(2),D(3)
      WRITE (*,*) ' '

 30   FORMAT(2X,'(',F11.8,',',F11.8,') (',F11.8,',',F11.8,') (',
     +    F11.8,',',F11.8,')')

      STOP
      END