! inverseC.f90  compute complex AI = A^-1  modified simeq.f90

subroutine inverseC(n, sz, A, AI)
  implicit none
  integer, intent(in) :: n  ! number of equations
  integer, intent(in) :: sz ! dimension of arrays
  double complex, dimension(sz,sz), intent(in) :: A
  double complex, dimension(sz,sz), intent(inout) :: AI

!      PURPOSE : COMPUTE INVERSEC WITH REAL COEFFICIENTS  |AI| = |A|^-1
!                                                                   
!      INPUT  : THE NUMBER OF ROWS  n
!               THE DIMENSION OF A, sz
!               THE REAL MATRIX  A
!      OUTPUT : THE REAL MATRIX  AI                                  

  integer, dimension(n) :: ROW             ! ROW INTERCHANGE INDICIES 
  integer, dimension(n) :: COL             ! COL INTERCHANGE INDICIES 
  double complex, dimension(n) :: TEMP   ! INTERCHANGE VECTOR
  integer :: HOLD , I_PIVOT, J_PIVOT       ! PIVOT INDICIES 
  double complex :: PIVOT                ! PIVOT ELEMENT VALUE 
  double precision :: ABS_PIVOT, NORM1
  integer :: i, j, k

  NORM1 = 0.0D0;
  ! BUILD WORKING DATA STRUCTURE 
  do i=1,n
    do j=1,n
      AI(i,j) = A(i,j)
      if( abs(AI(i,j)) > NORM1 ) then
        NORM1 = abs(AI(i,j))
      end if
    end do ! j
  end do ! i
  ! SET UP ROW AND COL  INTERCHANGE VECTORS 
  do k=1,n
    ROW(k) = k
    COL(k) = k
  end do ! k

  ! BEGIN MAIN REDUCTION LOOP 
  do k=1,n
    ! FIND LARGEST ELEMENT FOR PIVOT 
    PIVOT = AI(ROW(k), COL(k))
    I_PIVOT = k
    J_PIVOT = k
    do i=k,n
      do j=k,n
        ABS_PIVOT = abs(PIVOT)
        if( abs(AI(ROW(i), COL(j))) > ABS_PIVOT ) then
          I_PIVOT = i
          J_PIVOT = j
          PIVOT = AI(ROW(i), COL(j))
        end if
      end do ! j
    end do ! i
    ABS_PIVOT = abs(PIVOT)

    ! HAVE PIVOT, INTERCHANGE ROW, COL POINTERS 
    HOLD = ROW(k)
    ROW(k) = ROW(I_PIVOT)
    ROW(I_PIVOT) = HOLD
    HOLD = COL(k)
    COL(k) = COL(J_PIVOT)
    COL(J_PIVOT) = HOLD

    ! CHECK FOR NEAR SINGULAR 
    if( ABS_PIVOT < 1.0D-52*NORM1 ) then
      do j=1,n
        AI(ROW(k),j) = 0.0D0
      end do ! j
      do i=1,n
        AI(i,COL(k)) = 0.0D0
      end do ! i
      print *, 'redundant row (singular) ', ROW(k)
    else
      !                            REDUCE ABOUT PIVOT
      AI(ROW(k), COL(k)) = 1.0 / PIVOT
      do j=1,n
        if( j .ne. k ) then
          AI(ROW(k), COL(j)) = AI(ROW(k), COL(j)) * AI(ROW(k), COL(k))
        end if
      end do ! j
      !                            INNER REDUCTION LOOP
      do i=1,n
        if( k .ne. i ) then
          do j=1,n
            if( k .ne. j ) then
              AI(ROW(i), COL(j)) = AI(ROW(i), COL(j)) - &
                                   AI(ROW(i), COL(k)) * AI(ROW(k), COL(j))
            end if
          end do ! j
          AI(ROW(i), COL(k)) = - AI(ROW(i), COL(k)) * AI(ROW(k), COL(k))
        end if
      end do ! i
    end if
    ! FINISHED INNER REDUCTION 
  end do ! k
  ! END OF MAIN REDUCTION LOOP 

  !                              UNSCRAMBLE ROWS
  do j=1,n
    do i=1,n
      TEMP(COL(i)) = AI(ROW(i), j)
    end do ! i
    do i=1,n
      AI(i,j)= TEMP(i)
    end do !i
  end do ! j
  !                              UNSCRAMBLE COLUMNS
  do i=1,n
    do j=1,n
      TEMP(ROW(j)) = AI(i,COL(j))
    end do ! j
    do j=1,n
      AI(i,j)= TEMP(j)
    end do ! j
  end do ! i
end subroutine inverseC