C--------------------- Program: BT.FOR --------------------------------
C
C
C                    Dr. Neculai Andrei
C           Research Institute for Informatics         
C      Center for Advanced Modeling and Optimization
C
C
C
C Date created: November 2, 1992
C------------------------------------------------------------
C
C
C
	CHARACTER*30 NAMF
	DIMENSION A(10000),INDC(10000),NZL(1000)
        DIMENSION IP(1000),IQ(1000)
	DIMENSION NZND(1000),IZ(1000,4),IZL(1000,3),LOC(2)
	COMMON/RP02F/NO,NBLOC,MBLOC,NRNZ,STOPAT
C
C-------------------------------------------------------------
C---Main program for Block Triangularization of a large-scale
C   sparse matrix.
C
C   Remark: This is a sample program to illustrate the using
C           and the running of the subroutines for BT (Block
C           Triangular) permutation of a matrix. 
C
C	    The user may create itself its main program in 
C	    different circumstances.
C	    The dimensions allocated in this variant of the
C	    main program are for matrices up to 1000 rows (or
C	    columns) and 10000 non-zero elements.  Increasing  
C	    the above dimensions of the arrays, larger matrices
C	    can be considered.
C-------------------------------------------------------------
C   The following subroutines are invokated:
C   -RP02A
C      -RP02B
C         -RP02D
C      -RP02C
C         -RP02E
C--------------------------------------------------------------
C
C   In this example to use the subroutines for BT computation
C   it is necessary to create an input data file with the
C   structure and the non-zero elements of the matrix.
C   Such a file will be presented below.
C--------------------------------------------------------------
C
C
9	WRITE(*,10)
10	FORMAT(' Input data file name ? ',$)
	READ(*,11)NAMF
11	FORMAT(A30)
	OPEN(5,FILE=NAMF,STATUS='OLD',ERR=12)
	GO TO 14
12	WRITE(*,13)
13	FORMAT(' *** ERR: Open failure in input data file !')
	GO TO 9
14	CONTINUE  

C
C--- Output file
C                                       

      OPEN(UNIT=6,FILE='BT3.OUT',STATUS='UNKNOWN')   

C
C---Read the input data: 
C-----------------------
C   N   Dimension of the matrix.
C   IA  The length of the A and INDC arrays.
C
	READ(5,1) N,IA
1	FORMAT(2I3)
C
C---Read the structure and the non-zeros of the matrix:
C------------------------------------------------------
C   NZL   Array with the number of non-zero elements
C	  on rows.
C   INDC  Array with the column indices of non-zeros in rows.
C   A     Array with the non-zeros of the matrix.
C	  These two arrays INDC and A are of length IA.
C
C-----------------------------------------------------------------
	READ(5,2) (NZL(I),I=1,N)
2	FORMAT(40I2)
C
	NZ=0
	DO 3 I=1,N
	NZ=NZ+NZL(I)
3	CONTINUE
C
	READ(5,2) (INDC(I),I=1,NZ)
	READ(5,4) (A(I),I=1,NZ)
4	FORMAT(10F8.1)
C
C-------------------------------------------------------------
C---In order to illustrate the structure and the contents of
C   the input data file we shall consider a small example:
C
C---Example:
C-----------
C   Let us consider the following matrix for which we will
C   compute the BT form:
C
C      20  14   0   0   0   7   0  28
C      23   0   0  12   0  34   0  30
C       6   0   0  11   0  31  26  25
C      17  35  15  27   3   1  22   8
C      18   0   0   0   0   5   0   0
C      21  33   9  19   4  13   0  29
C       0   0   0  32   0  16   0  24
C       2   0   0   0   0  10   0   0
C
C   In this example we have:
C-  N = 8
C-  IA =100 (This is the length considered for this particular
C	     example. For larger applications the user may
C	     change the dimensions of the arrays in this
C	     program, and changing accordingly IA value.
C	     Please remember: this main program is taylored up
C	     to 1000 rows or columns and 10000 non-zeros into
C	     the matrix, but only IA locations from A and INDC
C	     arrays will be used.)
C
C-  The NZL array is as follows:
C
C   NZL = (4  4  5  8  2  7  3  2)
C
C   NZL array contains the number of non-zero elements in
C       each row of the matrix. As we can see the first row
C       has 4 non-zero elements, the second one 4, and so on.
C
C-  The INDC array is as follows:
C  
C   INDC = (1  2  6  8  1  4  6  8  1  4  6  7  8  1  2  3  
C           4  5  6  7  8  1  6  1  2  3  4  5  6  8  4  6
C           8  1  6)
C
C   INDC array contains the column indices of non-zero
C        elements from each row, given in rows order.
C        So, the first row has non-zeros in columns: 1,2,6,8;
C        the second row has non-zeros in columns: 1,4,6,8;
C        and so on.
C
C-  The A array is as follows:
C
C   A = (20.,14.,7.,28.,23.,12.,34.,30.,6.,11.,31.,26.,25.,
C        17.,35.,15.,27.,3.,1.,22.,8.,18.,5.,21.,33.,9.,19.,
C        4.,13.,29.,32.,16.,24.,2.,10.,)
C
C   A array contains the value of non-zero elements from the
C     matrix in accordance with the INDC array.
C
C   The user must create these arrays and introduce them in
C   an user file in above specified formats.
C   (Please see the above READ instructions.)
C
C-----------------------------------------------------------------
C
C---Permutation of the matrix to Block Triangular form.
C   Calling of the main subroutine RP02A.
C------------------------------------------------------
C
	CALL RP02A(N,IA,INDC,A,NZL,IP,IQ,NZND,LOC,IZ,IZL,
     +  MESAJ)

C
      WRITE(6,41) NBLOC, MBLOC, NRNZ
41    FORMAT(3I4)      

C
C-----------------------------------------------------------------
C---Print out the permutation IP(Rows) and IQ(Columns).
C------------------------------------------------------
	WRITE(6,5) (IP(I),IQ(I),I=1,N)
5	FORMAT(6X,I4,5X,I4)
C
C---Print out the non-zeros of the permuted matrix.
C--------------------------------------------------
	WRITE(6,6) (I,A(I),INDC(I),I=1,IA)
6	FORMAT(6X,I4,3X,E10.3,3X,I4)
C
C---Print out some other information of the permuted matrix.
C-----------------------------------------------------------
	WRITE(6,7) NBLOC,MBLOC,NRNZ
7	FORMAT(/8X,'THE NUMBER OF DIAGONAL BLOCKS: ',I4,/
     +      8X,'THE DIMENSION OF THE LARGEST BLOCK: ',I4,/
     +      8X,'THE MAXIMAL NUMBER OF NON-ZEROS ON DIAGONAL: ',I4)
C
	CLOSE(1)
C
C----------------------------------------------------------------
C
C   For the above example the IP and IQ arrays, containing
C   the rows permutation and the columns permutation respectively,
C   are as follows:
C
C   IP array with the rows permutation:
C   IP = (5  -8  2  -7  -1  -3  4  6)
C
C   IQ array with the columns permutation:
C   IQ = (6  1  4  8  2  7  5  3)
C
C   The matrix permuted accordind with the above permutations
C   is as follows:
C
C
C		 5  18   0   0   0   0   0   0
C		10   2   0   0   0   0   0   0
C		34  23  12  30   0   0   0   0
C		16   0  32  24   0   0   0   0
C		 7  20   0  28  14   0   0   0
C		31   6  11  25   0  26   0   0
C		 1  17  27   8  35  22   3  15
C		13  21  19  29  33   0   4   9
C
C
C   As we can see there are 5 bloc diagonal submatrices. The
C   largest one has the dimension 2. The main diagonal is zero
C   free.
C   Some more comments of the BT process as well as of the
C   arrays A, INDC and NZL are given in the subroutines for
C   BT computation. (Please see the comments in the head of
C   RP02*.FOR subroutines.)
C
C
C   Remark: BT program is accompanied by the RP02A, RP02B,
C           RP02C, RP02D and RP02E subroutines.
C--------------------------------------------------------------
C
	STOP
	END


C   Dr. Neculai Andrei
C   Research Institute for Informatics - Bucharest, Romania
C   8-10 Averescu Avenue, Sector 1
C   Center for Advanced Modeling and Optimization
C---------------------------------------------------- Last line.









C----------Subroutine RP02A.FOR
C------------------------------
C November 2, 1992
C
C
	SUBROUTINE RP02A(N,IA,INDC,A,NZL,IP,IQ,NZND,LOC,
     +  IZ,IZL,MESAJ)
C
C  ************************************************************
C
C0. PURPOSE.
C   *******
C   This subroutine, for a given large-scale sparse matrix A,
C   row packed, if it is possible, permute its rows and columns
C   to the lower-block-triangular form. 
C   The matrix is permuted to the form PAQ, according to P and
C   Q permutation matrices, so that the non-zero elements in the
C   off-diagonal blocks preceed those in the diagonal blocks,
C   which are in order.
C   If the user introduce his matrix by columns, then RP02A
C   subroutine will produce the upper-block-triangular form
C   and the corresponding permutation matrices.
C
C   For more details please see:
C     1.  F.GUSTAVSON., Finding the block lower triangular
C         form of a sparse matrix.
C         in J.R.BUNCH & D.J.ROSE (Ed.), Sparse matrix
C         computation. Academic Press, (1976), pp.275-289.
C     2.  I.S.DUFF., MA28 - A set of FORTRAN subroutines for
C         sparse unsymmetric linear equations.
C         Report AERE - R.8730, Harwell (1977).
C     3.  N.ANDREI, & C.RASTURNOIU., Sparse matrices and their
C         applications. Technical Press, Bucharest 1983.
C
C
C1. ARGUMENT LIST.
C   *************
C   CALL RP02A(N,IA,INDC,A,NZL,IP,IQ,NZND,LOC,IZ,IZL,MESAJ)
C
C   INPUTS:
C   ------
C   N      Integer must be set by the user to the order of the
C          matrix. It is not altered by the RP02A.
C   IA     Integer must be set by the user to the length of the
C          arrays A and INDC. It is not altered by RP02A.
C   INDC   is an integer array of length IA. On entry to RP02A,
C          it must be set by the user to contain the column indices
C          of the non-zero elements of the rows of the matrix.
C          The elements belonging to a single row must be
C          contiguous, but unordered on columns.
C          The matrix is row packed, i.e. the non-zeros of row I
C          preceed those of rows I+1, I=1,2,...,N-1, and no wasted
C          space is allowed between rows.
C          On exit from RP02A, the permuted column indices of the
C          rows in the diagonal part of the matrix are held in INDC
C          array from LOC(2) to IA locations, with the rows in
C          permuted order, without any space wasted between rows.
C          If there is more than one diagonal block, the original 
C          column indices of the parts of the rows in off-diagonal
C          blocks are in original row orded in INDC array from
C          1 to LOC(1)-1 locations without any space wasted between 
C          the rows.
C          The space in INDC array from LOC(1) to LOC(2)-1
C          locations is wasted.
C   A      is a real array of length IA. It is used only in RP02A
C          subroutine. On entry to RP02A it must be set by the user
C          to the values of the non-zeros of the matrix in the
C          corresponding positions from the array INDC.
C          On exit from RP02A, it contains the value of non-zeros
C          of the matrix reordered similarly to INDC.
C   NZL    is an integer array of length N. On entry to RP02A,
C          NZL(I) must be set by the user to the number of non-zeros
C          in row I, I=1,2,...,N.    On exit from RP02A,
C          NZL(I), I=1,2,...N contains the number of non-zeros
C          in the part of row I (of the permuted matrix) in its
C          diagonal block.
C
C
C   OUTPUTS:
C   -------
C   IP     is an integer array of length N. On exit from RP02A
C          modulus of IP(I), I=1,2,...N is the original row
C          index of the row which is the I-th in the permuted form.
C          The last row in the permuted form of each diagonal block
C          (except the last) is indicated by a negative value for IP.
C   IQ     is an integer array of length N.   On exit from RP02A
C          IQ(I), I=1,2,...N contains the original column index
C          of the column which is the I-th in the permuted form.
C   NZND   is an integer array of length N. On exit from RP02A
C          NZND(I), I=1,2,...,N contains the number of non-zeros
C          in the part of row I in the off-diagonal blocks.
C          If the matrix has only one diagonal block then
C          NZND(1)=-1, i.e. the matrix is irreductible.
C   LOC    is an integer array of length 2. On exit from RP02A,
C          LOC(1) gives the position in arrays A and INDC of the
C          first location after the off-diagonal blocks, and LOC(2)
C          gives the position in the same arrays of the first non-zero
C          in the diagonal blocks.
C  IZ      Is an integer work array of dimension (N,4).
C  IZL     Is an integer work array of dimension (N,2).
C
C
C   PARAMETERS:
C   ----------
C   MESAJ  integer used to transmit the error diagnostics.
C          A non-negative value on exit indicates success of RP02A.
C          After an unsuccessful entry a message is output on the line
C          printer and MESAJ is set negative to indicate one of the
C          following conditions:
C          -1  the matrix is structurally singular,
C          -2  IA is too small.
C
C
C2. USE OF COMMON.
C   *************
C   RP02A contains the following COMMON BLOCK:
C
C   COMMON/RP02F/NO,NBLOC,MBLOC,NRNZ,STOPAT
C
C   The following variables must be set by DATA in COMMON:
C   NO     is the index of the line printer.
C   STOPAT is a logical variable.  IF STOPAT=.TRUE. the 
C          process of ordering is suppressed in conditions
C          in which the matrix is structurally singular.
C
C   The following variables are returned by RP02A in COMMON
C   BLOCK to indicate:
C   NBLOC  The number of the diagonal blocks in the permuted
C          matrix.
C   MBLOC  The order of the largest diagonal block in the permuted 
C          matrix.
C   NRNZ   The number of non-zeros of the matrix which can be 
C          permuted onto the diagonal of the matrix.
C          If NRNZ < N, the matrix is structurally singular.
C
C
C3. OTHER ROUTINES.
C   **************
C   The following subroutines are called by RP02A:
C   RP02B  to find a row permutation to make diagonal of the
C          matrix zero-free.
C   RP02C  to find a symmetric permutation to the lower block 
C          triangular form.
C
C
C   Dr. Neculai Andrei
C   Research Institute for Informatics - Bucharest, Romania
C   8-10 Averescu Avenue, Sector 1
C   Center for Advanced Modeling and Optimization
C  ****************************************************************
C 
	REAL A(IA)
	INTEGER INDC(IA),NZL(N),IP(N),IQ(N),NZND(N)
	INTEGER IZ(N,4),IZL(N,2),LOC(2)
	LOGICAL STOPAT
	COMMON/RP02F/NO,NBLOC,MBLOC,NRNZ,STOPAT
C
	MESAJ=0
	IZL(1,1)=1
	NZND(1)=NZL(1)
	IF(N.EQ.1) GO TO 2
	DO 1 I=2,N
	NZND(I)=NZL(I)
1	IZL(I,1)=IZL(I-1,1)+NZL(I-1)
2	LOC(1)=IZL(N,1)+NZL(N)
C---FIND THE ROW PERMUTATION IP TO MAKE DIAGONAL
C   ZERO-FREE OF THE SPARSE  MATRIX.
	CALL RP02B(N,IA,INDC,IZL,NZL,IP,NRNZ,IZ)
	IF(NRNZ.NE.N.AND.STOPAT) GO TO 17
	DO 3 I=1,N
	II=IP(I)
	IZL(I,2)=IZL(II,1)
3	NZL(I)=NZND(II)
C---FIND THE SYMMETRIC PERMUTATION IQ TO BLOCK LOWER
C   TRIANGULAR FORM OF THE SPARSE MATRIX.
	CALL RP02C(N,IA,INDC,IZL(1,2),NZL,IQ,IZ(1,4),NBLOC,IZ)
	IF(NBLOC.NE.1) GO TO 6
C---THE MATRIX IS IRREDUCTIBLE.
C   THE PATTERN OF THE MATRIX IS MOVED TO THE 
C   END OF THE STORAGE.
	DO 4 I=1,N
	NZL(I)=NZND(I)
	IP(I)=I
4	IQ(I)=I
	NZND(1)=-1
	NZ=LOC(1)-1
	LOC(1)=1
	LOC(2)=IA-NZ+1
	MBLOC=N
	IF(NZ.EQ.IA) RETURN
	DO 5 K=1,NZ
	I=NZ-K+1
	II=IA-K+1
	A(II)=A(I)
5	INDC(II)=INDC(I)
	RETURN
C---REORDERING OF THE DATA STRUCTURE.
6	DO 7 I=1,N
	II=IQ(I)
7	IZ(I,1)=IP(II)
	DO 8 I=1,N
8	IP(I)=IZ(I,1)
	KP=IA+1
	KU=IA+1
	MBLOC=0
	DO 15 K=1,NBLOC
	IB=NBLOC-K+1
C---I1 IS THE FIRST ROW (IN THE PERMUTED FORM) 
C   OF THE IB BLOCK.
C---I2 IS THE LAST ROW (IN THE PERMUTED FORM)
C   OF THE IB BLOCK.
	I1=IZ(IB,4)
	I2=N
	IF(K.NE.1) I2=IZ(IB+1,4)-1
	MBLOC=MAX0(MBLOC,I2-I1+1)
C---PROCESS OF THE IB BLOCK IN REVERSE ORDER.
	DO 14 I=I1,I2
C---IN IS THE NEW ROW IN PERMUTED FORM.
C   IV IS THE OLD ROW IN ORIGINAL MATRIX.
	IN=I2+I1-I
	IV=IP(IN)
	IF(KU-LOC(1).GE.NZND(IV)) GO TO 11
	NPOZ=KP
	NEL=LOC(1)-1
	IF(NEL.LT.KP) GO TO 19
C---COMPRESS THE INDC AND A FILE.
C   MOVE THE OFF-DIAGONAL ELEMENTS AND UNTREATED ROWS
C   TO FRONT OF STORAGE.
	DO 9 J=KP,NEL
	IF(INDC(J).EQ.0) GO TO 9
	INDC(NPOZ)=INDC(J)
	A(NPOZ)=A(J)
	NPOZ=NPOZ+1
9	CONTINUE
	LOC(1)=NPOZ
	IF(KU-NPOZ.LT.NZND(IV)) GO TO 19
	KP=IA+1
C---RESET POINTERS TO THE BEGINNING OF THE ROWS.
	DO 10 II=2,N
10	IZL(II,1)=IZL(II-1,1)+NZND(II-1)
11	LB=IZL(IV,1)
	NZI=0
	LE=LB+NZND(IV)-1
	IF(LE.LT.LB) GO TO 13
	DO 12 J=LB,LE
	II=LE-J+LB
	JV=INDC(II)
	JN=IZ(JV,2)
	IF(JN.LT.I1) GO TO 12
	KU=KU-1
	A(KU)=A(II)
	INDC(KU)=JN
	KP=MIN0(KP,II)
	INDC(II)=0
	NZI=NZI+1
12	CONTINUE
	NZND(IV)=NZND(IV)-NZI
13	NZL(IN)=NZI
14	CONTINUE
	IP(I2)=-IP(I2)
15	CONTINUE
C---RESET IP(N) TO POSITIVE VALUE.
	IP(N)=-IP(N)
	LOC(2)=KU
	IF(KP.GT.IA) RETURN
	NPOZ=KP
	NEL=LOC(1)-1
	DO 16 J=KP,NEL
	IF(INDC(J).EQ.0) GO TO 16
	INDC(NPOZ)=INDC(J)
	A(NPOZ)=A(J)
	NPOZ=NPOZ+1
16	CONTINUE
	LOC(1)=NPOZ
	RETURN
C---THE FOLLOWING INSTRUCTIONS IMPLEMENT THE
C   FAILURE EXITS.
17	IF(NO.GT.0) WRITE(NO,18) NRNZ
18	FORMAT(3X,'THE MATRIX IS STRUCTURALLY SINGULAR',
     +  'RANK= ',I4)
	MESAJ=-1
	GO TO 21
19	IF(NO.GT.0) WRITE(NO,20)
20	FORMAT(3X,'IA IS TOO SMALL')
	MESAJ=-2
21	IF(NO.GT.0) WRITE(NO,22)
22	FORMAT(3X,'ERROR IN RP02A',//)
	RETURN
	END
C




	BLOCK DATA
C  **************************************************** 
	LOGICAL STOPAT
	COMMON/RP02F/NO,NBLOC,MBLOC,NRNZ,STOPAT
	DATA NO/6/,STOPAT/.FALSE./
C
C---------------------------------------------- End of RP02A.FOR
C
	END
C
C   Dr. Neculai Andrei
C   Research Institute for Informatics - Bucharest, Romania
C   8-10 Averescu Avenue, Sector 1
C   Center for Advanced Modeling and Optimization
C----------------------------------------------------------------




	


C---------- Subroutine RP02B.FOR
C-------------------------------
C November 2, 1992
C
C
	SUBROUTINE RP02B(N,IA,INDC,IPEL,NZL,IPER,NRNZ,IZ)
C
C  **********************************************************
C
C0. PURPOSE.
C   *******
C   This subroutine, for a given pattern of non-zeros of a 
C   sparse matrix, finds a row permutation that makes the
C   matrix have a maximum number of non-zero elements on its
C   diagonal, using a depth first search with look ahead
C   technique.
C
C   For more details please see:
C     1.  I.S.DUFF,. On algorithms for obtaining a maximum
C         transversal.
C         ACM-TOMS, Vol.7, Nr.3, (1981), pp.315-330.
C     2.  I.S.DUFF,. Algorithm 575. Permutations for a zero
C         free diagonal (F1).
C         ACM-TOMS, Vol.7, Nr.3, (1981), pp.387-390.
C     3.  N.ANDREI, & C.RASTURNOIU,. Sparse matrices and their
C         applications. Technical Press, Bucharest 1983.
C
C1. ARGUMENT LIST.
C   *************
C   CALL RP02B(N,IA,INDC,IPEL,NZL,IPER,NRNZ,IZ)
C
C   INPUTS:
C   ------
C   N      integer must be set by the user to the order of the
C          matrix. It is not altered by RP02B.
C   IA     integer must be set by the user to the length of array
C          INDC. It is not altered by RP02B.
C   INDC   is an integer array of length IA. On entry to and exit 
C          from RP02B, it must be set by the user to contain the
C          column indices of the non-zero elements.
C	   The elements belonging to a single row must be contigous, 
C          but unordered on columns within each row.
C          The matrix is row packed, in any order, wasted space 
C          between rows being permitted.
C   IPEL   is an integer array of length N. On entry to and exit
C          from RP02B, IPEL(I) contains the point to the start of
C          row I, I=1,2,...,N.
C   NZL    is an integer array of length N. On entry to and exit 
C          from RP02B, NZL(I) contains the number of non-zeros in
C          row I, I=1,2,...,N.
C
C
C   OUTPUTS:
C   -------
C   IPER   is an integer array of length N. On exit from RP02B
C          IPER(I) contains the position in the original matrix
C	   of row I in the permuted matrix, I=1,2,...,N.
C   NRNZ   integer which gives the number of non-zeros on the
C          diagonal of the permuted matrix. If NRNZ < N, then the
C	   matrix is structurally singular, and in this case
C          N-NRNZ of the elements (IPER(I),I) will be zero.
C   IZ     is an integer work array of dimension (N,4).
C
C
C2. OTHER ROUTINES.
C   **************
C   The following subroutine is called by RP02B:
C   RP02D   to process the row permutation for a zero-free
C           diagonal of a sparse matrix. This is the nucleus
C           for the RP02B subroutine.
C
C
C   Dr. Neculai Andrei
C   Research Institute for Informatics - Bucharest, Romania
C   8-10 Averescu Avenue, Sector 1
C   Center for Advanced Modeling and Optimization
C  ***************************************************************
C
	INTEGER INDC(IA),IPEL(N),NZL(N),IPER(N),IZ(N,4)
	CALL RP02D(N,IA,INDC,IPEL,NZL,IPER,NRNZ,IZ(1,1),
     +  IZ(1,2),IZ(1,3),IZ(1,4))
	RETURN
C
C------------------------------------------------ End of RP02B.FOR
C
	END
C
C   Dr. Neculai Andrei
C   Research Institute for Informatics - Bucharest, Romania
C   8-10 Averescu Avenue, Sector 1
C   Center for Advanced Modeling and Optimization
C--------------------------------------------------------------------






C---------- Subroutine RP02C.FOR
C-------------------------------
C November 2, 1992
C
C
	SUBROUTINE RP02C(N,IA,INDC,IPEL,NZL,IPER,LIB,NBLOC,IZL)
C
C  **************************************************************
C
C0. PURPOSE.
C   *******
C   This subroutine, for a given pattern of non-zeros of a sparse 
C   matrix, finds a symmetric permutation that makes the matrix
C   lower block triangular, using the Tarjan's depth first search
C   algorithm.
C   To obtain the best results, the user must first permute the
C   structure of the matrix so that it has a zero-free diagonal.
C   This can be done using RP02B subroutine.
C
C   For more details please see:
C     1.  I.S.DUFF, J.K.REID., An implementation of Tarjan's
C         algorithm for the block triangularization of a 
C         matrix. 
C         ACM-TOMS., Vol.4, Nr.2, (1978), pp.137-147.
C     2.  I.S.DUFF, J.K.REID., Algorithm 529, permutation to
C         block triangular form (F1).
C         ACM-TOMS., Vol.4, Nr.2, (1978), pp.189-192.
C     3.  N.ANDREI, & C.RASTURNOIU,. Sparse matrices and their
C         applications. Technical Press, Bucharest, 1983.
C 
C
C1. ARGUMENT LIST.
C   *************
C   CALL RP02C(N,IA,INDC,IPEL,NZL,IPER,LIB,NBLOC,IZL)
C
C   INPUTS:
C   ------
C   N      integer must be set by the user to the order of the
C	   matrix. It is not altered by RP02C.
C   IA     integer must be set by the user to the length of
C	   array INDC. It is not altered by RP02C.
C   INDC   is an integer array of length IA. On entry to and
C	   exit from RP02C, it must be set by the user to contain
C          the column indices of the non-zero elements.
C          The elements belonging to a single row must be contigous,
C	   but unordered on columns within each row.
C	   The matrix is row packed, in any order, wasted space 
C	   between rows being permitted.
C   IPEL   is an integer array of length N. On entry to and exit
C	   from RP02C, IPEL(I) contains the point to the start of
C          row I, I=1,2,...,N.
C   NZL    is an integer array of length N. On entry to and exit
C          from RP02C, NZL(I) contains the number of non-zeros in
C          row I, I=1,2,...,N.
C
C   OUTPUTS:
C   -------
C   IPER   is an integer array of length N. On exit from RP02C,
C          IPER(I) contains the position in the original matrix of 
C          row or column which is in position I in the permuted
C          matrix, I=1,2,...,N.
C   LIB    is an integer array of length N. On exit from RP02C,
C          LIB(I) contains the row indices in the permuted matrix
C          of the beginning of the I-th block, I=1,2,...,NBLOC.
C   NBLOC  integer which records the number of blocks found in the
C          permuted matrix.
C   IZL    is an integer work array of dimension (N,3).
C
C
C2. OTHER ROUTINES.
C   **************
C   The following subroutine is called by RP02C:
C   RP02E   To process the rows and columns permutation (symmetric
C           permutation) of a sparse matrix to lower block
C           triangular form. This is the nucleus for RP02C 
C	    subroutine.
C
C  *****************************************************************
C
	INTEGER INDC(IA),IPEL(N),NZL(N),IPER(N),LIB(N),IZL(N,3)
	CALL RP02E(N,IA,INDC,IPEL,NZL,IPER,LIB,NBLOC,IZL(1,1),
     +  IZL(1,2),IZL(1,3))
	RETURN
C
C-------------------------------------------------- End of RP02C.FOR
C
	END


C   Dr. Neculai Andrei
C   Research Institute for Informatics - Bucharest, Romania
C   8-10 Averescu Avenue, Sector 1
C   Center for Advanced Modeling and Optimization
C------------------------------------------------------------------







C---------- Subroutine RP02D.FOR
C-------------------------------
C November 2, 1992
C
C
	SUBROUTINE RP02D(N,IA,INDC,IPEL,NZL,IPER,NRNZ,LIA,
     +  AI,COLV,NEX)
C--------------------------------------------------------------
	INTEGER INDC(IA),IPEL(N),NZL(N),IPER(N)
	INTEGER LIA(N),COLV(N),AI(N),NEX(N)
C---INITIALIZATION.
	NRNZ=0
	DO 1 I=1,N
	AI(I)=NZL(I)-1
	COLV(I)=0
1	IPER(I)=0
C---MAIN LOOP FOR A NEW ASSIGNAMENT OR FOR A ROW
C   WITH NO ASSIGNMENT.
	DO 10 I=1,N
	J=I
	LIA(J)=-1
C---LOOK FOR A CHEAP ASSIGNMENT.
	DO 7 K=1,I
	KP=AI(J)
	IF(KP.LT.0) GO TO 3
	KU=IPEL(J)+NZL(J)-1
	KP=KU-KP
	DO 2 KK=KP,KU
	IC=INDC(KK)
	IF(IPER(IC).EQ.0) GO TO 8
2	CONTINUE
C---NO CHEAP ASSIGNMENT IN ROW.
	AI(J)=-1
C---ASSIGNMENT CHAIN STARTING WITH ROW J.
3	NEX(J)=NZL(J)-1
C---EXTENDS CHAIN BY ONE OR BACKTRACKS.
	DO 6 II=1,I
	KP=NEX(J)
	IF(KP.LT.0) GO TO 5
	KU=IPEL(J)+NZL(J)-1
	KP=KU-KP
C---FORWARD SCAN FOR ASSIGNMENT OF ROW J IN
C   COLUMN IC
	DO 4 KK=KP,KU
	IC=INDC(KK)
	IF(COLV(IC).EQ.I) GO TO 4
C---COLUMN IC HAS NOT YET BEEN ACCESSED DURING
C   THIS PASS.
	JJ=J
	J=IPER(IC)
	COLV(IC)=I
	LIA(J)=JJ
	NEX(JJ)=KU-KK-1
C C---LOOK FOR A CHEAP ASSIGNMENT IN ROW IPER(IC)=J.
	GO TO 7
4	CONTINUE
C---BACKTRACKING STEP.
5	J=LIA(J)
	IF(J.EQ.-1) GO TO 10
6	CONTINUE
7	CONTINUE
C---A NEW ASSIGNMENT IS MADE.
8	IPER(IC)=J
	AI(J)=KU-KK-1
	NRNZ=NRNZ+1
C---REASSIGNMENT CHAIN.
	DO 9 K=1,I
	J=LIA(J)
	IF(J.EQ.-1) GO TO 10
	KK=IPEL(J)+NZL(J)-NEX(J)-2
	IC=INDC(KK)
C---COLUMN IC WILL BE COLUMN J IN THE PERMUTED MATRIX.
	IPER(IC)=J
9	CONTINUE
10	CONTINUE
	IF(NRNZ.EQ.N) RETURN
C---THE MATRIX IS STRUCTURALLY SINGULAR.
C   COMPLETE THE PERMUTATION IPER.	
	DO 11 I=1,N
11	AI(I)=0
	K=0
	DO 13 I=1,N
	IF(IPER(I).NE.0) GO TO 12
	K=K+1
	NEX(K)=I
	GO TO 13
12	J=IPER(I)
	AI(J)=I
13	CONTINUE
	K=0
	DO 14 I=1,N
	IF(AI(I).NE.0) GO TO 14
	K=K+1
	II=NEX(K)
	IPER(II)=I
14	CONTINUE
	RETURN
C----------------------------------------------- End of RP02D.FOR
	END


C   Dr. Neculai Andrei
C   Research Institute for Informatics - Bucharest, Romania
C   8-10 Averescu Avenue, Sector 1
C   Center for Advanced Modeling and Optimization
C----------------------------------------------------------------







C---------- Subroutine RP02E.FOR
C-------------------------------
C November 2, 1992
C
C
	SUBROUTINE RP02E(N,IA,INDC,IPEL,NZL,NAR,LIB,NBLOC,
     +  LEGINF,NUMB,ANTE)
C--------------------------------------------------------------
	INTEGER INDC(IA),IPEL(N),NZL(N),LIB(N)
	INTEGER LEGINF(N),NUMB(N),ANTE(N),NAR(N)
C---INITIALIZATION.
	NN=0
	NBLOC=0
	N2=N+N-1
	DO 1 J=1,N
	NUMB(J)=0
1	NAR(J)=NZL(J)-1
C---MAIN LOOP FOR ORDERING OF NODES.
	DO 9 I=1,N
C---LOOK FOR A STARTING NODE.
	IF(NUMB(I).NE.0) GO TO 9
	IV=I
C---NS IS THE NUMBER OF NODE ON THE STACK.
	NS=1
C---SET THE NODE IV AT BEGINNING OF THE STACK.
	LEGINF(IV)=1
	NUMB(IV)=1
	LIB(N)=IV
C---INTERNAL LOOP TO PUT A NEW NODE ON THE STACK
C   OR BACKTRACKING.
	DO 8 J=1,N2
	KP=NAR(IV)
	IF(KP.LT.0) GO TO 3
	KU=IPEL(IV)+NZL(IV)-1
	KP=KU-KP
C---LOOK AT EDGES LEAVING NODE IV UNTIL ONE ENTERS
C   A NEW NODE OR ALL EDGES ARE EXHAUSTED.
	DO 2 KK=KP,KU
	IW=INDC(KK)
	IF(NUMB(IW).EQ.0) GO TO 7
	IF(LEGINF(IW).LT.LEGINF(IV)) LEGINF(IV)=LEGINF(IW)
2	CONTINUE
	NAR(IV)=-1
C---IF NODE IV IS NOT THE ROOT OF A BLOCK GO TO 6.
3	IF(LEGINF(IV).LT.NUMB(IV)) GO TO 6
C---OBTAIN STRONG COMPONENT FROM TOP OF STACK.
	NBLOC=NBLOC+1
	NS1=N+1-NS
	NN1=NN+1
C---REMOVE NODES FROM PATH AND STACK.
	DO 4 K=NS1,N
	IW=LIB(K)
	LEGINF(IW)=N+1
	NN=NN+1
	NUMB(IW)=NN
	IF(IW.EQ.IV) GO TO 5
4	CONTINUE
5	NS=N-K
	LIB(NBLOC)=NN1
C---IS PATH NULL ?
	IF(NS.NE.0) GO TO 6
C---HAVE ALL THE NODES BEEN ORDERED ?
	IF(NN.LT.N) GO TO 9
	GO TO 10
C---BACKTRACK TO PREVIOUS NODE ON PATH.
6	IW=IV
	IV=ANTE(IV)
C---ADJUST VALUE OF LEGINF(IV) IF NECESSARY.
	IF(LEGINF(IW).LT.LEGINF(IV)) LEGINF(IV)=LEGINF(IW)
	GO TO 8
C---PUT NEW NODE ON THE STACK AND PATH.
7	NAR(IV)=KU-KK-1
	ANTE(IW)=IV
	IV=IW
	NS=NS+1
	LEGINF(IV)=NS
	NUMB(IV)=NS
	II=N+1-NS
	LIB(II)=IV
8	CONTINUE
9	CONTINUE
C---SET THE PERMUTATION IN REQUIRED FORM WHICH PUTS
C   THE MATRIX IN BLOCK LOWER TRIANGULAR FORM.
10	DO 11 I=1,N
	II=NUMB(I)
11	NAR(II)=I
	RETURN
C-------------------------------------------------- End of RP02E.FOR
	END
	

c############################################################  
* Date assembled: February 8, 2010.
* 
* Dr. Neculai Andrei
* Research Institute for Informatics         
* Center for Advanced Modeling and Optimization
* 8-10, Averescu Avenue, Bucharest 1, Romania  
* E-mail: nandrei@ici.ro          
* Voice 666.58.74   0744-777748
c############################################################  
*
* Line Final