subroutine p6(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 integer num_new, point, col, temp, row, new_border,i1 logical get c c Usage of iwork 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 c c print *,'start tarjan' c 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*********************************************************************** c get stack_inv do i = 1, n iwork(iwork(i,7),8) = i enddo c*********************************************************************** 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 sort the border c------------------------------------------------------------------- c sort the columns. num_iter = 0 2000 continue do i= 1,n iwork(i,1) = 0 iwork(i,2) = 0 enddo 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 c c write (0,*) 'start to sort the border' 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 border = block_ptr(nblocks+1) if (num_iter .ge. 1) goto 3000 num_new = 0 border = block_ptr(nblocks+1) point = 1 do i1 = 1, border-1 10 col = n get = .false. do while (.not. get .and. iwork(col,1).lt.i1 1 .and. col .ge.border) 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.i1-1.and.iwork(i,2).ne.0) then get = .true. row = i ii = col_ptr(j+1) endif ii = ii + 1 enddo if (get) goto 100 40 col = col -1 enddo 100 continue 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,i1-1) enddo goto 10 else iwork(point,3) = col_perm(i1) iwork(i1,4) = point iwork(i1,2) = 0 point = point + 1 endif enddo new_border = point c c reset the block_ptr. 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).eq.0) then nn = nn+1 endif enddo point = 1 do i = 1,nblocks j = block_ptr(i) do while (iwork(point,3).ne.col_perm(j)) point= point + 1 enddo block_ptr(i) = point enddo 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 reset row_perm. 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 c 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 if (num_new .ne. 0 ) then num_iter = num_iter +1 c write(0,*) 'the num_border:', block_ptr(nblocks+1) goto 2000 endif c c permute the rows and columns to the position. c 3000 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