SUBROUTINE cmpscscbm (maxoffnz, a, SchTol, SC)

USE m_dump
USE m_build

INTEGER							, INTENT(IN)    :: maxoffnz
DOUBLE PRECISION					, INTENT(IN)    :: SchTol
TYPE (scdematrix)					, POINTER       :: sc
TYPE (scbmmatrix)					, POINTER       :: a

!     Compute Schur-complement from SCBM-type matrix.
!     Only called from the subroutine  schurcmpl!

!     Arguments:
!     ==========
!     A%n    	i   Number of rows/columns in matrix A.
!     A%g     	i   Size of the upper partition.
!     MaxOffNz 	i   Maximum number of off-diagonal non-zeros that can
!                   be stored in new Schur-complement 'SC'.
!     A%a11d    i   The block-diagonal submatrix  inv(A_11).C
!     A%a12%beg i   A%a12%beg(i): index in 'A%a12%jco' and 'A%a12%co' of the first
!                   element in row i of submatrix A_12 of A.
!     A%a12%jco i   A%a12%jco(nz): column number of element A%a12%co(nz).
!     A%a12%co  i   A%a12%co(nz): value of nonzero element of A_12.
!     A%a21%beg i   A%a21%beg(i): index in 'A%a21%jco' and 'A%a21%co' of the first
!                   element in row i of submatrix A_21 of A.
!     A%a21%jco i   A%a21%jco(nz): column number of element A%a21%co(nz).
!     A%a21%co  i   A%a21%co(nz): value of nonzero element of A_21.
!     A%a22%beg i   A%a22%beg(i): index in 'A%a22%jco' and 'A%a22%co' of the first
!                   element in row i of submatrix A_22 of A.
!     A%a22%jco i   A%a22%jco(nz): column number of element A%a22%co(nz).
!     A%a22%co  i   A%a22%co(nz): value of nonzero element of A_22.
!     SchTol   	i   All non-zero off-diagonal elements of the Schur-
!                   complement, S, to be stored into 'SC' should be
!                   greater than 'SchTol'.
!     sc%offd%beg    o   sc%offd%beg(i): index in 'sc%offd%jco' and 'sc%offd%co' of the first
!                  off-diagonal element in row i of matrix SC.
!     sc%offd%jco    o   sc%offd%jco(nz): column number of off-diagonal element
!                  sc%offd%co(nz).
!     sc%offd%co     o   sc%offd%co(nz): value of nonzero off-diagonal element of SC.
!     sc%dia%com    o   The values of the elements in the diagonal blocks of
!                  matrix SC, stored in column major order.

!     Array arguments used as local variables:
!     ========================================
!     ColNr    ColNr(i), 1 <= i <= NnzRow <= A%Nschur, column number
!              of i-th element that should be stored in the new
!              Schur-complement.
!     StorCol  StorCol(i), 1 <= i <= A%Nschur,
!              = .TRUE.  Element in column 'i' of actual row should be
!                     stored in new Schur-complement.
!              = .FALSE.   Element in column 'i' of actual row should not be
!                     stored in new Schur-complement.
!     ValCol   ValCol(i), 1 <= i <= A%Nschur, value to be stored in
!              column 'i' of actual row of the new Schur-complement.

!     Local Parameters:
!     =================
CHARACTER (LEN=*), PARAMETER :: rounam = 'cmpscscbm'

!     Local Variables:
!     ================
!     Row   1 <= Row <= A%n-A%g, the number of the actual row.
!     NnzRow   Number of nonzeros in actual row 'Row'.
!     Sc%offd%nnz    Number of non-zeros in off-diagonal part of the new
!              Schur-complement 'SC'.

INTEGER 					:: scnnz, ier
INTEGER 					:: basecol, BaseRow, col
INTEGER 					:: Row, nnzrow
INTEGER 					:: nza12, nza21, nza22, nz
INTEGER 					:: cola21, rowa12
DOUBLE PRECISION 				:: fact
LOGICAL, ALLOCATABLE, DIMENSION(:)        	:: StorCol
INTEGER, ALLOCATABLE, DIMENSION(:)             	:: ColNr
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)    	:: ValCol

#ifdef DEBUG

!     TRACE INFORMATION
PRINT '(A, X, A)' , 'Entry:', rounam
#endif

ALLOCATE( StorCol(1:A%Nschur), STAT=ier )
IF (ier /= 0) CALL dump(__FILE__,__LINE__,'Allocation error')
ALLOCATE( ColNr(1:A%Nschur), STAT=ier )
IF (ier /= 0) CALL dump(__FILE__,__LINE__,'Allocation error')
ALLOCATE( ValCol(1:A%Nschur), STAT=ier )
IF (ier /= 0) CALL dump(__FILE__,__LINE__,'Allocation error')

!     Clear column search pointer:
StorCol = .false.

!     Initialise number of nonzeros and pointer to first non-zero in
!     'SC':
scnnz    = 0
sc%offd%beg(1) = scnnz + 1

!     Construct Schur-Complement, of size A%NSCHUR, in  'SC' and 'Sc%dia':

IF (sc%dia%BlkSiz > 1) THEN
  DO Row = 1, A%Nschur
    BaseRow = ((Row - 1) / sc%dia%BlkSiz) * sc%dia%BlkSiz
    
    nnzrow = 0
    
!        Copy elements in row 'Row' of A_22 into new row of SC:
    
    DO nza22 = a%a22%beg(Row), a%a22%beg(Row+1)-1
      nnzrow        = nnzrow + 1
      col          = a%a22%jco(nza22)
      ColNr(nnzrow) = col
      StorCol(col) = .true.
      ValCol(col)  = a%a22%co(nza22)
    END DO
    
!        Add zero block-diagonal elements of A_22 if necessary:
    DO col = BaseRow+1, BaseRow+sc%dia%BlkSiz
      IF (.not. StorCol(col)) THEN
        nnzrow        = nnzrow + 1
        ColNr(nnzrow) = col
        StorCol(col) = .true.
        ValCol(col)  = 0.0D0
      END IF
    END DO
    
    
!        Subtract product  A_21*(1/A_11)*A_12  from  A_22
    
!        For each non-zero in row of A_21:
    DO nza21 = a%a21%beg(Row), a%a21%beg(Row+1)-1
      cola21 = a%a21%jco(nza21)
      
      basecol = ((cola21 - 1) / sc%dia%BlkSiz) * sc%dia%BlkSiz
      DO rowa12 = basecol+1, basecol+sc%dia%BlkSiz
        
!              Common factor for all non-zeros in row of A_12
        fact = a%a21%co(nza21)*a%a11d%com(cola21 - basecol,rowa12)
        
!              For each non-zero in row of A_12:
        DO nza12 = a%a12%beg(rowa12), a%a12%beg(rowa12+1)-1
          
          col = a%a12%jco(nza12)
          
          IF (StorCol(col)) THEN
!                    An existing element:
            
            ValCol(col) = ValCol(col) - fact*a%a12%co(nza12)
          ELSE
!                    A new element, insert into new row:
            
            nnzrow        = nnzrow + 1
            ColNr(nnzrow) = col
            StorCol(col) = .true.
            ValCol(col)  = - (fact * a%a12%co(nza12))
          END IF
        END DO
      END DO
    END DO
    
    
    IF (scnnz + (nnzrow - sc%dia%BlkSiz) > maxoffnz) THEN
!           Not enough space left to store off-diagonals of 'Row'.
      PRINT '(A, 2X, A, A, /, 3X, A, I8, /, 3X, A, I8)' ,  &
          'Error in', rounam, '.  Not enough space to store Schur-complement!',  &
          'Detected when storing row:         ', Row,  &
          'Number of rows in Schur-complement:', a%Nschur
      CALL dump(__FILE__,__LINE__, 'Not enough space left to store off-diagonals!')
    END IF
    
!        Store elements:
    
    
!        Off-diagonal elements:
    
!        Store off-diagonal elements, with absolute value > 'SchTol', into
!        the approximate Schur-complement 'SC'.  Off-diagonal elements
!        with absolute value <= SchTol in row 'Row' are lumped onto the
!        diagonal element in this row.
    
!        For each non-zero in new row:
    DO nz = 1, nnzrow
      col = ColNr(nz)
      IF (col <= BaseRow  .Or. col > BaseRow+sc%dia%BlkSiz) THEN
        StorCol(col) = .false.
        
        IF (DABS(ValCol(col)) > SchTol) THEN
!                 Off-diagonal element cannot be lumped, store:
          scnnz        = scnnz + 1
          sc%offd%co(scnnz)  = ValCol(col)
          sc%offd%jco(scnnz) = col
        ELSE
!                 Lump off-diagonal onto main diagonal:
           ValCol(Row) = ValCol(Row) + ValCol(col)
        END IF
      END IF
    END DO
    
!        Pointer to first non-zero Off-diagonal in next row:
    sc%offd%beg(Row+1) = scnnz + 1
    
    
!        Block-diagonal elements:
    
!        Store the block part of the actual row into block diagonal part
!        of Schur complement 'Sc%dia':
    
    StorCol(BaseRow+1:BaseRow+sc%dia%BlkSiz) = .false.
    sc%dia%com(Row - BaseRow,BaseRow+1:BaseRow+sc%dia%BlkSiz) = ValCol(BaseRow+1:BaseRow+sc%dia%BlkSiz)
    
  END DO
ELSE
!     {  Sc%dia%BlkSiz = 1  }
  DO Row = 1, a%Nschur
    
    nnzrow = 0
    
!        Copy elements in row 'Row' of A_22 into new row SC:
    
    DO nza22 = a%a22%beg(Row), a%a22%beg(Row+1)-1
      nnzrow        = nnzrow + 1
      col          = a%a22%jco(nza22)
      ColNr(nnzrow) = col
      StorCol(col) = .true.
      ValCol(col)  = a%a22%co(nza22)
    END DO
    
!        Add zero block-diagonal element of A_22 if necessary:
    IF (.not. StorCol(Row)) THEN
      nnzrow          = nnzrow + 1
      ColNr(nnzrow)   = Row
      StorCol(Row) = .true.
      ValCol(Row)  = 0.0D0
    END IF
    
!        Subtract product  A_21*(1/A_11)*A_12  from  A_22
    
!        For each non-zero in row of A_21:
    DO nza21 = a%a21%beg(Row), a%a21%beg(Row+1)-1
      cola21 = a%a21%jco(nza21)
      
!           Common factor for all non-zeros in row of A_12
      fact = a%a21%co(nza21)*a%a11d%com(1,cola21)
      
!           For each non-zero in row of A_12:
      DO nza12 = a%a12%beg(cola21), a%a12%beg(cola21+1)-1
        
        col = a%a12%jco(nza12)
        
        IF (StorCol(col)) THEN
!                 An existing element:
          
          ValCol(col) = ValCol(col) - fact*a%a12%co(nza12)
        ELSE
!                 A new element, insert into new row:
          
          nnzrow        = nnzrow + 1
          ColNr(nnzrow) = col
          StorCol(col) = .true.
          ValCol(col)  = - (fact * a%a12%co(nza12))
        END IF
      END DO
    END DO
    
    IF (scnnz + (nnzrow - sc%dia%BlkSiz) > maxoffnz) THEN
!           Not enough space left to store off-diagonals of 'Row'.
      PRINT '(A, 2X, A, A, /, 3X, A, I8, /, 3X, A, I8)' ,  &
          'Error in', rounam, '.  Not enough space to store Schur-complement!',  &
          'Detected when storing row:         ', Row,  &
          'Number of rows in Schur-complement:', a%Nschur
      CALL dump(__FILE__,__LINE__, 'Not enough space left to store off-diagonals!')
    END IF
    
    
!        Off-diagonal elements:
    
!        Store off-diagonal elements, with absolute value > 'SchTol', into
!        the approximate Schur-complement 'SC'.  Off-diagonal elements
!        with absolute value <= SchTol in row 'Row' are lumped onto the
!        diagonal element in this row.
    
!        For each off-diagonal non-zero in new row:
    DO nz = 1, nnzrow
      col = ColNr(nz)
      IF (col /= Row) THEN
        StorCol(col) = .false.
        
        IF (DABS(ValCol(col)) > SchTol) THEN
!                 Element cannot be lumped, store off-diagonal:
          
          scnnz        = scnnz + 1
          sc%offd%co(scnnz)  = ValCol(col)
          sc%offd%jco(scnnz) = col
        ELSE
!                 Lump off-diagonal onto main diagonal:
          ValCol(Row) = ValCol(Row) + ValCol(col)
        END IF
      END IF
    END DO
    
!        Pointer to first non-zero Off-diagonal in next row:
    sc%offd%beg(Row+1) = scnnz + 1
    
    
!        Diagonal element:
    
!        Store the Diagonal element of the actual row into the diagonal
!        part of the Schur complement 'Sc%dia':
    
    StorCol(Row) = .false.
    sc%dia%com(1,Row)   = ValCol(Row)
  END DO
END IF
sc%offd%nnz = scnnz

DEALLOCATE( StorCol, STAT=ier )
IF (ier /= 0) CALL dump(__FILE__,__LINE__,'Deallocation error')
DEALLOCATE( ColNr, STAT=ier )
IF (ier /= 0) CALL dump(__FILE__,__LINE__,'Deallocation error')
DEALLOCATE( ValCol, STAT=ier )
IF (ier /= 0) CALL dump(__FILE__,__LINE__,'Deallocation error')

END SUBROUTINE cmpscscbm
