subroutine bmp6(n,nz,value,row_ind,col_ptr,
     1	    row_perm, col_perm, nblocks, block_ptr, iwork,
     2	    dwork, new_row_ind,iflag )
c
c---------------------------------------------------------------------------
c	This is an implementation of P6 algorithm.
c---------------------------------------------------------------------------

c***********************************************************************
c	Arguments:
c
c	n - integer: the order of the matrix
c	nz - integer: the number of non-zero elements in the matrix
c	max_row - integer: the dimension of the row_ind array
c	max_col - integer: the dimension of the col_ind array
c	value - double precision(nz): the values for each element
c	    of the array
c	row_ind - integer(max_row): this array contains the row
c	    indices for each of the elements in the array
c	col_ind - integer(max_col): this array contains the column
c	    indices for each of the elements in the array
c	col_ptr - integer(n+1): contains a pointer to the start of
c	    each column
c	row_perm - integer(n): row permutation
c	col_perm - integer(n): column permutation
c	nblocks - integer: the number of blocks found
c	block_ptr - integer(n): pointer to the first row in each block
c	iwork - integer(n,13): this is also work space
c	dwork - double precision(n): this is double precision work space
c	iflag - integer: this is the error return flag
c***********************************************************************
	integer n, nz, max_row, max_col, max_block_size, bsize
	double precision value(nz)
	integer row_ind(nz), new_row_ind(nz)
	integer row_perm(n), col_perm(n), col_ptr(n+1)
	integer nblocks, block_ptr(n), border, offset, save
	integer iwork(n,10)
	double precision dwork(nz)
	integer iflag, block_start, min_blk_size, max_blk_size
	integer num_new, point, col, temp, row, new_border,i1
	logical get
c
c	For the information, here is the list of the usage of the
c	working interger space 
c
c	equivalence (stack(1), iwork(1,7)), 
c     1       (stack_inv(1), iwork(1,8)),
c     1       (row_len(1),iwork(1,8)),
c     1       (block_end(1), iwork(1,1)), 
c     1       (stack_pos(1), iwork(1,2)),
c     1       (prev_node(1),iwork(1,3)), 
c     1       (lowest(1), iwork(1,4)),
c     1       (jused(1), iwork(1,5)), 
c     1       (visited(1), iwork(1,6)),
c     1       (row_count(1), iwork(1,1)), 
c     1       (col_count(1), iwork(1,2)),
c     1       (min_col(1), iwork(1,3)), 
c     1       (candidate(1),iwork(1,4)),
c     1       (list_candidate(1), iwork(1,5)), 
c     1       (new_value(1),dwork(1)),
c     1       (new_col_ptr(1), iwork(1,1))
c     1       (border_ptr(1) , iwork(1,6))
c***********************************************************************
	
	do i=1, n
	    iwork(i,8)= col_ptr(i+1)-col_ptr(i)
	enddo
	nblocks = 0
	call tarjan( n, nz, nblocks, row_ind, col_ptr, iwork(1,8),
     1	    iwork(1,2), block_ptr, iwork(1,1), iwork(1,7), iwork(1,3),
     2	    iwork(1,4), iwork(1,5), iwork(1,6) )
c	print *,'tarjan finish'
c	print *,'number of blocks:',nblocks
c	print *,block_ptr(1:nblocks+1)
c	print *,iwork(1:nblocks+1,1)

c***********************************************************************
c	get stack_inv

	do i = 1, n
	    iwork(iwork(i,7),8) = i
	enddo

c***********************************************************************
c      cvd$ Cncall
	do i = 1, nblocks
	   call p5(n,nz, col_ptr, row_ind, col_perm, row_perm,
     1                i, block_ptr, iwork(1,7), iwork(1,8), 
     1                iwork(1,1), iwork(1,2), iwork(1,3), 
     1                iwork(1,5), iwork(1,4), iwork(1,6),iflag)
	   if (iflag .ne. 0 ) then
		write(0,*) 'there is errro in',i,'block',iflag
		iflag = 0
	   endif
	enddo

c***********************************************************************
c	permute matrix into bordered block triangular by moving
c	all the border rows from every blocks to the globe border.
c***********************************************************************

c	------------------------------------------------------------
c	permute the rows between border_ptr(i) and block_ptr(i+1)-1
c	the border sort is not needed. 
c	------------------------------------------------------------
	save = n
	border = n
	offset = 0
	do i = 1, nblocks
	    do j = block_ptr(i), block_ptr(i+1)-1
		if (iwork(i,6).le.row_perm(j) .and.
     1              row_perm(j).lt.block_ptr(i+1)) then
		    row_perm(j) = border
		    border = border-1
		else
		    row_perm(j) = row_perm(j)-offset
		endif
		if (iwork(i,6).le.j .and. 
     1              j.lt.block_ptr(i+1)) then
		    iwork(save,5) = col_perm(j)
		    save = save - 1
	 	else
		    col_perm(j-offset) = col_perm(j)
		endif
		if (iwork(j,7).gt.0) then
		    iwork(j,7) = iwork(j,7) - offset
		endif
	    enddo

	    if ( save .ne. border) then
		error = 1
		goto 1000
	    endif

	    block_ptr(i) = block_ptr(i) - offset
	    offset = n - border
	enddo

	do i = border+1, n
	    col_perm(i) = iwork(i,5)
	enddo

	j = 1
	do i = 1,n
	    if ( iwork(i,7).gt.0) then
		block_ptr(j) = iwork(i,7)
		j = j + 1
	    endif
	enddo
	block_ptr(j) = border+1
	nblocks = j - 1
	block_ptr(j+1) = n+1
c
c-------------------------------------------------------------------
c	adjust the block size to the min_blk_size.
c-------------------------------------------------------------------
c	
	min_blk_size = max(20,n/15)
	max_blk_size = max(1,n/10)
	point = 1
	sum = 0
	do i = 2, nblocks+1
	    sum = sum + (block_ptr(i)-block_ptr(i-1))
	    if (sum .ge. min_blk_size ) then
		block_ptr (point + 1) = block_ptr(point) + sum
		point = point + 1
		sum = 0
	    endif
	enddo
	if (sum .ne. 0) then
	block_ptr(point+1) = block_ptr(point) + sum
	nblocks = point
	else
	nblocks = point-1
	endif
c
c-------------------------------------------------------------------
c	sort the border
c-------------------------------------------------------------------
c	sort the columns.

	num_iter = 0
	do i = 1, n
	   iwork(i,1) = 0
	   iwork(i,2) = 0
	end do
	border = block_ptr(nblocks+1)
	do i = border, n
	    iwork(i,1) = 0
	    ii = col_perm(i)
	    do jj = col_ptr(ii),col_ptr(ii+1)-1
		j = abs(row_perm(iwork(row_ind(jj),8)))
		if (j.gt.iwork(i,1).and.j.lt.border) iwork(i,1) = j
	    enddo
	    k = i-1
	    do while (k.ge.border .and. iwork(k,1).lt.iwork(k+1,1))
	 	temp = col_perm(k+1)
		col_perm(k+1) = col_perm(k)
		col_perm(k) = temp
		temp = iwork(k+1,1)
		iwork(k+1,1) = iwork(k,1)
		iwork(k,1) = temp
		k = k - 1
	    enddo
	enddo

	jj = 1
	border = block_ptr(nblocks+1)
	do while (jj.le.n .and. border.le.n)
	    j = col_perm(jj)
	    do ii = col_ptr(j), col_ptr(j+1)-1
		i = row_perm(iwork(row_ind(ii),8))
		if (i.ge.block_ptr(nblocks+1) .and. i.le.n) then
		    row_perm(iwork(row_ind(ii),8)) = -border
		    iwork(border, 2) = jj
		    border = border + 1
		endif
	    enddo
	    jj = jj + 1
	enddo
	border = block_ptr(nblocks+1)
c
c	write (0,*) 'border sort is finished'
c	extand the diagnal blocks.
c
	num_new = 0
	border = block_ptr(nblocks+1)
	point = 1
	block_start = point

	    i1 = 1
	    col = n
	    do while (.not. get .and. iwork(col,1).lt.i1
     1                .and. col .ge.border .and.
     1                point-block_start .lt. max_blk_size)
		if (iwork(col,1).lt.0) goto 240
		j = col_perm(col)
		ii = col_ptr(j)
		do while(ii .lt. col_ptr(j+1))
		    i = abs(row_perm(iwork(row_ind(ii),8)))
		    if ( i.ge.border.and. i.le.n .and.
     1                   iwork(i,2).ge.0.and.
     1                   iwork(i,2).ne.0) then
			get = .true.
			row = i
			ii = col_ptr(j+1)
		    endif
		    ii = ii + 1
		enddo
	        if(get) then
		iwork(row,2) = 0
		iwork(col,1) = -1
		iwork(point,3) = col_perm(col)
		iwork(row,4) = point
		point = point + 1
		num_new = num_new+1
		j0 = col_perm(col)
		do j1 = col_ptr(j0),col_ptr(j0+1)-1
		    j2 = abs(row_perm(iwork(row_ind(j1),8)))
		    if (i1.le.iwork(j2,2).and.j2.ge.border
     1                         .and. j2.ne.row) 
     1                  iwork(j2,2) = -1
		enddo
		get = .false.
		endif
240		col = col -1
	    enddo

	do i0 = 1, nblocks
	    block_start = point
	    do i1 = block_ptr(i0),block_ptr(i0+1)-1
		iwork(point,3) = col_perm(i1)
		iwork(i1,4) = point
		iwork(i1,2) = 0
		point = point + 1
	    enddo
	    i1 = block_ptr(i0+1)
	    col = n
	    get = .false.
	    do while (.not. get .and. iwork(col,1).lt.i1
     1                .and. col .ge.border .and.
     1                point-block_start .lt. max_blk_size)
		if (iwork(col,1).lt.0) goto 40
		j = col_perm(col)
		ii = col_ptr(j)
		do while(ii .lt. col_ptr(j+1))
		    i = abs(row_perm(iwork(row_ind(ii),8)))
		    if ( i.ge.border.and. i.le.n .and.
     1                   iwork(i,2).ge.block_ptr(i0).and.
     1                   iwork(i,2).ne.0) then
			get = .true.
			row = i
			ii = col_ptr(j+1)
		    endif
		    ii = ii + 1
		enddo
	        if(get) then
		iwork(row,2) = 0
		iwork(col,1) = -1
		iwork(point,3) = col_perm(col)
		iwork(row,4) = point
		point = point + 1
		num_new = num_new+1
		j0 = col_perm(col)
		do j1 = col_ptr(j0),col_ptr(j0+1)-1
		    j2 = abs(row_perm(iwork(row_ind(j1),8)))
		    if (i1.le.iwork(j2,2).and.j2.ge.border
     1                         .and. j2.ne.row) 
     1                  iwork(j2,2) = max(1,block_ptr(i0))
		enddo
		get = .false.
		endif
40		col = col -1
	    enddo
	enddo
	new_border = point 
	nn = 0
	do i = 1,n
	    if (abs(row_perm(i)).ge.border .and. abs(row_perm(i))
     1         .le.n .and. iwork(abs(row_perm(i)),2).eq.0) then
		nn = nn+1
	    endif
	enddo
	point = 1
	if (iwork(point,3).eq.col_perm(1)) then
	    nb = -1
	else
	    nb = 0
	endif

	j = block_ptr(1)
	do i = 2,nblocks+1
	    do while (iwork(point,3).ne.col_perm(j))
	    	point= point + 1
	    enddo
	    j = block_ptr(i)
	    block_ptr(i+nb) = point
	enddo
	nblocks = nblocks+1+nb
	block_ptr(nblocks+1) = new_border
c
c	reset the col_perm.
c
	point = new_border
	do i = border,n
	    if(iwork(i,1) .ge. 0) then
		iwork(point,3) = col_perm(i)
		point = point + 1
	    endif
	enddo
	if (point .ne. n+1) then
		iflag = 10
		goto 1000
	endif
	do i = 1, point-1
	    col_perm(i) = iwork(i,3)
	enddo
c
c	reset row_perm.
c
	point = new_border 
	nn = 0
	do i = 1,n
	    if (abs(row_perm(i)).ge.border .and. abs(row_perm(i))
     1         .le.n .and. iwork(abs(row_perm(i)),2).ne.0) then
		iwork(abs(row_perm(i)),4) = point
		point = point + 1
	    elseif (abs(row_perm(i)).ge.border .and. abs(row_perm(i))
     1         .le.n .and. iwork(abs(row_perm(i)),2).eq.0) then
		nn = nn+1
	    endif
	enddo
	if (point .ne. n+1) then
		iflag = 11
		goto 1000
	endif
c	print *,'check the permutation'
c	do i=1,n
c	iwork(i,1)= 0
cc	iwork(i,2)= 0
c	enddo
c	do i=1,n
c	iwork(iwork(i,4),2)=iwork(iwork(i,4),2)+1
c	iwork(col_perm(i),1)=iwork(col_perm(i),1)+1
c	enddo
c	do i=1,n
c	if (iwork(i,1).ne.1) then
c	print *,'there is problem in ', i,'th row',iwork(i,1)
c	iflag = 12
c	goto 1000
c	endif
c	enddo

	do i = 1,n
	    row_perm(i) = iwork(abs(row_perm(i)),4)
	enddo


c	print *,'check the permutation'
c	do i=1,n
c	iwork(i,1)= 0
c	iwork(i,2)= 0
c	enddo
c	do i=1,n
c	iwork(row_perm(i),2)=iwork(row_perm(i),2)+1
c	iwork(iwork(i,3),1)=iwork(iwork(i,3),1)+1
c	enddo
c	do i=1,n
c	if (iwork(i,1).ne.1) then
c	print *,'there is problem in ', i,'th row',iwork(i,1)
c	iflag = 13
c	goto 1000
c	endif
c	enddo

c
c	write (0,*) 'sort the border again'
c
	jj = 1
	border = block_ptr(nblocks+1)
	do while (jj.le.n .and. border.le.n)
	    j = col_perm(jj)
	    do ii = col_ptr(j), col_ptr(j+1)-1
		i = row_perm(iwork(row_ind(ii),8))
		if (i.ge.block_ptr(nblocks+1) .and. i.le.n) then
		    row_perm(iwork(row_ind(ii),8)) = -border
		    iwork(border, 2) = jj
		    border = border + 1
		endif
	    enddo
	    jj = jj + 1
	enddo
c
c	permute the rows and columns to the position.
c
	point = 1
	do col = 1,n
	    iwork(col,1) = point
	    j = col_perm(col)
	    do i = col_ptr(j), col_ptr(j+1)-1
		dwork(point) = value(i)
		new_row_ind(point) = abs(row_perm(iwork(row_ind(i),8)))
		point = point + 1
	    enddo
	enddo
	if (point-1 .eq.nz) then
	    do i = 1, n
	    	col_ptr(i) = iwork(i,1)
	    enddo
	    do i = 1, nz
		row_ind(i) = new_row_ind(i)
		value(i) = dwork(i)
	    enddo
		col_ptr(n+1) = nz+1
	else
	    iflag = 1
	endif

1000	continue
	return
	end