SUBROUTINE ZSUR(A,INT,X,N,NC,EMACH)  
C     ------------------------------------------------------------------  
C     ZSUR IS A STANDARD BACK-SUBSTITUTION SUBROUTINE USING THE   
C     OUTPUT OF ZGER TO CALCULATE X TIMES A-INVERSE, RETURNED IN X.
C     IT HAS ONLY 21 PROGRAM LINES.
C     FOR THE MATRIX PROBLEM WHERE THE MATRIX 'A' MULTIPLIES A 
C     VECTOR X ON THE RIGHT:
C
C                         X*A=B
C
C     INT RECORDS PIVOTING DETAILS OF THE GAUSS ELIMINATION 
C         PERFOMED IN ZGER.
C     EMACH IS A CUTOFF ON THE MATRIX ELEMENTS.
C     ------------------------------------------------------------------  
      IMPLICIT NONE 
C  
C ..  SCALAR ARGUMENTS  ..  
C  
      INTEGER N,NC  
      REAL*8 EMACH  
C  
C ..  ARRAY ARGUMENTS  ..  
C  
      INTEGER    INT(NC)  
      COMPLEX*16 A(NC,NC),X(NC)  
C  
C ..  LOCAL SCALARS  ..  
C  
      INTEGER    I,II,IN,J,IJ  
      COMPLEX*16 DUM  
C  
C ..  INTRINSIC FUNCTIONS  ..  
C  
*      INTRINSIC ABS  
C     ------------------------------------------------------------------  
C 
C Rearanging the "right-hand side" using pivoting information from
C ZGER:
 
      DO 5 II=2,N  
      I=II-1               !from 1 to N-1
      IF(INT(I)-I)1,2,1    !If INT(I).NE.I exchange the I-th and INT(I)-th 
                           !elements of the vector X
   1  IN=INT(I)  
      DUM=X(IN)  
      X(IN)=X(I)  
      X(I)=DUM  
*
C Rearanging the "right-hand side" by mirroring the alegbaric 
C operations performed on columns of 'A' matrix in ZGER:
* If A(I,J) is sufficiently large, X(J)=X(J)-X(I)*A(I,J)
*
   2  DO 4 J=II,N                 !from I+1 to N
      IF(ABS(A(I,J))-EMACH)4,4,3  
   3  X(J)=X(J)-X(I)*A(I,J)
   4  CONTINUE  

   5  CONTINUE               !the I-th row of A multiplied by X(I) 
                             !subtracted from X                    
*
* BACKSUBSTITUTION:
*
      DO 10 II=1,N  
      I=N-II+1         !from N   to 1
      IJ=I+1           !from N+1 to 2
      IF(I-N)6,8,6  
   6  DO 7 J=IJ,N  
   7  X(I)=X(I)-X(J)*A(J,I)

* Imposing a cutoff on the possible underflow.
* You can change the factor 1.0D-7 below to some
* other number, depending on the precision of your
* computer

   8  IF(ABS(A(I,I))-EMACH*1.0D-7)9,10,10  
   9  A(I,I)=EMACH*1.0D-7*(1.D0,1.D0)  
  10  X(I)=X(I)/A(I,I)  
      RETURN  
      END  

C (C) Copr. 03/2003  Alexander Moroz