sd_ssf.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE sd_ssf
!                   *****************
!
     &(n,p,ip,ia,ja,iju,ju,iu,jumax,q,mark,jl,flag)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    SYMBOLIC UT-D-U FACTORISATION OF SPARSE SYMMETRIC MATRIX.
!+                COMPRESSED STORAGE OF SPARSE MATRICES.
!code
!+  IMPORTANT NOTE: INSPIRED FROM PACKAGE CMLIB3 - YALE UNIVERSITE-YSMP
!+
!+
!+    THE STRICT UPPER TRIANGULAR PORTION OF THE MATRIX U IS STORED IN
!+    (IA,JA,A) FORMAT USING THE ARRAYS IU, JU, AND U, EXCEPT THAT AN
!+    ADDITIONAL ARRAY IJU IS USED TO REDUCE THE STORAGE REQUIRED FOR JU
!+    BY ALLOWING SOME SEQUENCES OF COLUMN INDICES TO CORRESPOND TO MORE
!+    THAN ONE ROW.  FOR I
!+    ENTRY FOR THE I-TH ROW;  IJU(N) IS THE NUMBER OF ENTRIES IN JU.
!+    THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS  IU(I+1) - IU(I),
!+    THE NONZERO ENTRIES OF THE I-TH ROW ARE STORED CONSECUTIVELY IN
!+
!+        U(IU(I)),   U(IU(I)+1),   ..., U(IU(I+1)-1),
!+
!+    AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN
!+
!+        JU(IJU(I)), JU(IJU(I)+1), ..., JU(IJU(I)+IU(I+1)-IU(I)-1).
!+
!+    COMPRESSION IN JU OCCURS IN TWO WAYS.  FIRST, IF A ROW I WAS MERGED
!+    INTO ROW K, AND THE NUMBER OF ELEMENTS MERGED IN FROM (THE TAIL
!+    PORTION OF) ROW I IS THE SAME AS THE FINAL LENGTH OF ROW K, THEN
!+    THE KTH ROW AND THE TAIL OF ROW I ARE IDENTICAL AND IJU(K) POINTS
!+    TO THE START OF THE TAIL.  SECOND, IF SOME TAIL PORTION OF THE
!+    (K-1)ST ROW IS IDENTICAL TO THE HEAD OF THE KTH ROW, THEN IJU(K)
!+    POINTS TO THE START OF THAT TAIL PORTION.  FOR EXAMPLE, THE NONZERO
!+    STRUCTURE OF THE STRICT UPPER TRIANGULAR PART OF THE MATRIX
!+
!+             ( D 0 0 0 X X X )
!+             ( 0 D 0 X X 0 0 )
!+             ( 0 0 D 0 X X 0 )
!+         U = ( 0 0 0 D X X 0 )
!+             ( 0 0 0 0 D X X )
!+             ( 0 0 0 0 0 D X )
!+             ( 0 0 0 0 0 0 D )
!+
!+    WOULD BE STORED AS
!+
!+             \ 1  2  3  4  5  6  7  8
!+         ----+------------------------
!+          IU \ 1  4  6  8 10 12 13 13
!+          JU \ 5  6  7  4  5  6
!+         IJU \ 1  4  5  5  2  3  6           .
!+
!+    THE DIAGONAL ENTRIES OF U ARE EQUAL TO ONE AND ARE NOT STORED.
!
!history  E. RAZAFINDRAKOTO (LNH)
!+        20/11/06
!+        V5P7
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FLAG           |<--| LOGICAL VARIABLE;  IF DFLAG = .TRUE., THEN
!|                |   | STORE NONZERO DIAGONAL ELEMENTS AT THE
!|                |   | BEGINNING OF THE ROW
!| IA             |<--| INTEGER ONE-DIMENSIONAL ARRAY CONTAINING
!|                |   | POINTERS TO DELIMIT ROWS IN JA AND A;
!|                |   | DIMENSION = N+1
!| IJU            |---| INTEGER ONE-DIMENSIONAL ARRAY CONTAINING
!|                |   | POINTERS TO THE START OF EACH ROW IN JU;  DIMENSION = N
!| IP             |<--| INTEGER ONE-DIMENSIONAL ARRAY USED TO RETURN
!|                |   | THE INVERSE OF THE PERMUTATION RETURNED IN P;
!|                |   | DIMENSION = N
!| IU             |---| INTEGER ONE-DIMENSIONAL ARRAY CONTAINING
!|                |   | POINTERS TO DELIMIT ROWS IN U;  DIMENSION = N+1
!| JA             |<--| INTEGER ONE-DIMENSIONAL ARRAY CONTAINING THE
!|                |   | COLUMN INDICES CORRESPONDING TO THE ELEMENTS
!|                |   | OF A;  DIMENSION = NUMBER OF NONZERO ENTRIES
!|                |   | IN (THE UPPER TRIANGLE OF) M
!| JL             |-->| INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION N
!|                |   | DIMENSION = NUMBER OF NONZERO ENTRIES IN THE
!|                |   | UPPER TRIANGLE OF M. JL CONTAINS LISTS OF ROWS
!|                |   | TO BE MERGED INTO UNELIMINATED ROWS --
!|                |   | I GE K => JL(I) IS THE FIRST ROW TO BE
!|                |   | MERGED INTO ROW I
!|                |   | I LT K => JL(I) IS THE ROW FOLLOWING ROW I IN
!|                |   | SOME LIST OF ROWS
!|                |   | IN EITHER CASE, JL(I) = 0 INDICATES THE
!|                |   | END OF A LIST
!| JU             |---| INTEGER ONE-DIMENSIONAL ARRAY CONTAINING THE
!|                |   | COLUMN INDICES CORRESPONDING TO THE ELEMENTS
!|                |   | OF U;  DIMENSION = JUMAX
!| JUMAX          |---| DECLARED DIMENSION OF THE ONE-DIMENSIONAL
!|                |   | ARRAY JU;  JUMAX MUST BE AT LEAST THE SIZE
!|                |   | OF U MINUS COMPRESSION (IJU(N) AFTER THE
!|                |   | CALL TO SSF)
!| MARK           |---| INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION N
!|                |   | THE LAST ROW STORED IN JU FOR WHICH U(MARK(I),I)
!|                |   | NE 0
!| N              |-->| ORDER OF THE MATRIX
!| P              |<--| INTEGER ONE-DIMENSIONAL ARRAY USED TO RETURN
!|                |   | THE PERMUTATION OF THE ROWS AND COLUMNS OF M
!|                |   | CORRESPONDING TO THE MINIMUM DEGREE ORDERING;
!|                |   | DIMENSION = N
!| Q              |-->| INTEGER ONE-DIMENSIONAL WORK ARRAY, DIMENSION N
!|                |   | Q CONTAINS AN ORDERED LINKED LIST
!|                |   | REPRESENTATION OF THE NONZERO
!|                |   | STRUCTURE OF THE K-TH ROW OF U --
!|                |   | Q(K) IS THE FIRST COLUMN WITH A NONZERO ENTRY
!|                |   | Q(I) IS THE NEXT COLUMN WITH A NONZERO ENTRY
!|                |   | AFTER COLUMN I
!|                |   | IN EITHER CASE, Q(I) = N+1 INDICATES THE
!|                |   | END OF THE LIST
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_sd_ssf => sd_ssf
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: N,JUMAX
      INTEGER, INTENT(INOUT) :: P(n),IP(n),IA(n+1),JA(*),IJU(n),FLAG
      INTEGER, INTENT(INOUT) :: JU(jumax),IU(n+1),Q(n),MARK(n),JL(n)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,J,M,TAG,VJ,QM,JUMIN,JUPTR,K,LUK,LUI,JMIN,JMAX,LMAX
      LOGICAL CLIQUE
!
!-----------------------------------------------------------------------
!
!----INITIALISES
!
!    JUMIN AND JUPTR ARE THE INDICES IN JU OF THE FIRST AND LAST
!    ELEMENTS IN THE LAST ROW SAVED IN JU
!
!    LUK IS THE NUMBER OF NONZERO ENTRIES IN THE K-TH ROW
!
      jumin = 1
      juptr = 0
      iu(1) = 1
      DO k=1,n
        mark(k) = 0
        jl(k) = 0
      ENDDO
!
!----FOR EACH ROW K
!
      DO 19 k=1,n
        luk = 0
        q(k) = n+1
!
        tag = mark(k)
        clique = .false.
        IF(jl(k).NE.0)  clique = jl(jl(k)).EQ.0
!
!------INITIALISES NONZERO STRUCTURE OF K-TH ROW TO ROW P(K) OF M
!
        jmin = ia(p(k))
        jmax = ia(p(k)+1) - 1
        IF (jmin.GT.jmax)  GO TO 4
        DO 3 j=jmin,jmax
          vj = ip(ja(j))
          IF(vj.LE.k)  GO TO 3
!
          qm = k
2         m = qm
          qm = q(m)
          IF (qm.LT.vj)  GO TO 2
          IF (qm.EQ.vj)  GO TO 102
          luk = luk+1
          q(m) = vj
          q(vj) = qm
          IF (mark(vj).NE.tag)  clique = .false.
!
3       CONTINUE
!
!------IF EXACTLY ONE ROW IS TO BE MERGED INTO THE K-TH ROW AND THERE IS
!------A NONZERO ENTRY IN EVERY COLUMN IN THAT ROW IN WHICH THERE IS A
!------NONZERO ENTRY IN ROW P(K) OF M, THEN DOES NOT COMPUTE FILL-IN, JUST
!------USES THE COLUMN INDICES FOR THE ROW WHICH WAS TO HAVE BEEN MERGED
!
4       IF(.NOT.clique)  GO TO 5
        iju(k) = iju(jl(k)) + 1
        luk = iu(jl(k)+1) - (iu(jl(k))+1)
        GO TO 17
!
!------MODIFIES NONZERO STRUCTURE OF K-TH ROW BY COMPUTING FILL-IN
!------FOR EACH ROW I TO BE MERGED IN
!
5       lmax = 0
        iju(k) = juptr
!
        i = k
6       i = jl(i)
        IF (i.EQ.0)  GO TO 10
!
!--------MERGES ROW I INTO K-TH ROW
!
        lui = iu(i+1) - (iu(i)+1)
        jmin = iju(i) +  1
        jmax = iju(i) + lui
        qm = k
!
        DO 8 j=jmin,jmax
          vj = ju(j)
7         m = qm
          qm = q(m)
          IF (qm.LT.vj)  GO TO 7
          IF (qm.EQ.vj)  GO TO 8
          luk = luk+1
          q(m) = vj
          q(vj) = qm
          qm = vj
8       CONTINUE
!
!--------REMEMBERS LENGTH AND POSITION IN JU OF LONGEST ROW MERGED
!
        IF(lui.LE.lmax)  GO TO 6
        lmax = lui
        iju(k) = jmin
!
        GO TO 6
!
!------IF THE K-TH ROW IS THE SAME LENGTH AS THE LONGEST ROW MERGED,
!------THEN USES THE COLUMN INDICES FOR THAT ROW
!
10      IF (luk.EQ.lmax)  GO TO 17
!
!------IF THE TAIL OF THE LAST ROW SAVED IN JU IS THE SAME AS THE HEAD
!------OF THE K-TH ROW, THEN OVERLAPS THE TWO SETS OF COLUMN INDICES --
!--------SEARCHES LAST ROW SAVED FOR FIRST NONZERO ENTRY IN K-TH ROW ...
!
        i = q(k)
        IF (jumin.GT.juptr)  GO TO 12
        DO 11 jmin=jumin,juptr
          IF (ju(jmin) < i) GOTO 11
          IF (ju(jmin) == i) GOTO 13
          IF (ju(jmin) > i) GOTO 12
11      CONTINUE
12      GO TO 15
!
!--------... AND THEN TESTS WHETHER TAIL MATCHES HEAD OF K-TH ROW
!
13      iju(k) = jmin
        DO j=jmin,juptr
          IF (ju(j).NE.i)  GO TO 15
          i = q(i)
          IF (i.GT.n)  GO TO 17
        ENDDO
        juptr = jmin - 1
!
!------SAVES NONZERO STRUCTURE OF K-TH ROW IN JU
!
15      i = k
        jumin = juptr +  1
        juptr = juptr + luk
        IF (juptr.GT.jumax)  GO TO 106
        DO j=jumin,juptr
          i = q(i)
          ju(j) = i
          mark(i) = k
        ENDDO
        iju(k) = jumin
!
!------ADDS K TO ROW LIST FOR FIRST NONZERO ELEMENT IN K-TH ROW
!
17      IF(luk.LE.1)  GO TO 18
        i = ju(iju(k))
        jl(k) = jl(i)
        jl(i) = k
!
18    iu(k+1) = iu(k) + luk
!
19    CONTINUE
!
      flag = 0
      RETURN
!
! ** ERROR -- DUPLICATE ENTRY IN A
!
102   flag = 2*n + p(k)
      RETURN
!
! ** ERROR -- INSUFFICIENT STORAGE FOR JU
!
106   flag = 6*n + k
      RETURN
      END