sd_nnfc.f

Go to the documentation of this file.

!                   ******************
                    SUBROUTINE sd_nnfc
!                   ******************
!
     &(n,r,c,ic,ia,ja,a,z,b,lmax,il,jl,ijl,l,d,umax,iu,ju,iju,u,
     & row,tmp,irl,jrl,flag)
!
!***********************************************************************
! BIEF   V6P3                                   30/06/2013
!***********************************************************************
!
!brief NUMERICAL LDU-FACTORIZATION OF SPARSE NONSYMMETRIC MATRIX AND
!+      SOLUTION OF SYSTEM OF LINEAR EQUATIONS (COMPRESSED POINTER
!+      STORAGE)
!+
!+
!+       INPUT VARIABLES..  N, R, C, IC, IA, JA, A, B,
!+                          IL, JL, IJL, LMAX, IU, JU, IJU, UMAX
!+       OUTPUT VARIABLES.. Z, L, D, U, FLAG
!+
!+ FIA   \ ROW   - HOLDS INTERMEDIATE VALUES IN CALCULATION OF  U AND L.
!+       \           SIZE = N.
!+ FIA   \ TMP   - HOLDS NEW RIGHT-HAND SIDE  B*  FOR SOLUTION OF THE
!+       \           EQUATION UX = B*.
!+       \           SIZE = N.
!+
!+
!note     IMPORTANT : INSPIRED FROM PACKAGE CMLIB3 - YALE UNIVERSITE-YSMP
!
!         DON'T HESITATE TO CHANGE IN/OUTPUT VARIABLES COMMENTS
!         FOR CLARITY
!
!history  C. PEYRARD (LNHE)
!+        30/06/13
!+        V6P3
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| NONZERO ENTRIES OF THE COEFFICIENT MATRIX M,
!|                |   | STORED BY ROWS
!| B              |-->| RIGHT-HAND SIDE B ;
!| C              |-->| ORDERING OF THE COLUMNS OF MATRIX
!| D              |<--| DIAGONAL FACTORIZED OF MATRIX
!| FLAG           |<--| INDICATOR ERROR :
!|                |   |= 4*N + 1:INSUFFICIENT STORAGE FOR L
!|                |   |= 7*N + 1:INSUFFICIENT STORAGE FOR U
!| IA, JA         |-->| STRUCTURE OF A NONSYMMETRICAL MATRIX
!| IA, JA         |-->| STRUCTURE OF A NONSYMMETRICAL MATRIX
!| IC             |-->| INVERSE OF THE ORDERING OF THE COLUMNS OF MATRIX
!| IL, JL         |-->| STRUCTURE OF LOWER FACTORISED TRIANGULAR MATRIX
!| IU, JU         |-->| STRUCTURE OF UPPER FACTORISED TRIANGULAR MATRIX
!| IJU,IJL        |-->| USED TO COMPRESS STORAGE OF JU and JL
!| IRL,JRL        |-->| VECTORS USED TO FINDTHE ROWS OF L, AT THE kth STEP
!| L              |<--| LOWER FACTORIZED TRIANGULAR MATRIX
!| LMAX           |-->| PREVISIONAL MAXIMUM DIMENSION OF JL
!| N              |-->| RANK OF MATRIX
!| R              |-->| ORDERING OF THE ROWS OF MATRIX
!| ROW            |---| REAL ONE-DIMENSIONAL WORK ARRAY
!| TMP            |---| REAL ONE-DIMENSIONAL WORK ARRAY
!| U              |<--| UPPER FACTORIZED TRIANGULAR MATRIX
!| UMAX           |-->| PREVISIONAL MAXIMUM DIMENSION OF JU
!| Z              |<--| SOLUTION X
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER RK,UMAX
      INTEGER  R(*), C(*), IC(*), IA(*), JA(*), IL(*), JL(*), IJL(*)
      INTEGER  IU(*), JU(*), IJU(*), IRL(*), JRL(*), FLAG
      INTEGER  N,K,I1,I,I2,JMIN,JMAX,J,MU,IJLB,LMAX
      DOUBLE PRECISION  A(*), L(*), D(*), U(*), Z(*), B(*), ROW(*)
      DOUBLE PRECISION  TMP(*), LKI, SOMME, DK
!
!  ******  INITIALIZE POINTERS AND TEST STORAGE  ***********************
!
      IF(il(n+1)-1 .GT. lmax) GO TO 104
      IF(iu(n+1)-1 .GT. umax) GO TO 107
      DO 1 k=1,n
        irl(k) = il(k)
        jrl(k) = 0
   1    CONTINUE
!
!  ******  FOR EACH ROW  ***********************************************
      DO 19 k=1,n
!  ******  REVERSE JRL AND ZERO ROW WHERE KTH ROW OF L WILL FILL IN  ***
        row(k) = 0
        i1 = 0
        IF (jrl(k) .EQ. 0) GO TO 3
        i = jrl(k)
   2    i2 = jrl(i)
        jrl(i) = i1
        i1 = i
        row(i) = 0
        i = i2
        IF (i .NE. 0) GO TO 2
!  ******  SET ROW TO ZERO WHERE U WILL FILL IN  ***********************
   3    jmin = iju(k)
        jmax = jmin + iu(k+1) - iu(k) - 1
        IF (jmin .GT. jmax) GO TO 5
        DO 4 j=jmin,jmax
   4      row(ju(j)) = 0
!  ******  PLACE KTH ROW OF A IN ROW  **********************************
   5    rk = r(k)
        jmin = ia(rk)
        jmax = ia(rk+1) - 1
        DO 6 j=jmin,jmax
          row(ic(ja(j))) = a(j)
   6      CONTINUE
!  ******  INITIALIZE SOMME, AND LINK THROUGH JRL  *********************
        somme = b(rk)
        i = i1
        IF (i .EQ. 0) GO TO 10
!  ******  ASSIGN THE KTH ROW OF L AND ADJUST ROW, SOMME  **************
   7      lki = -row(i)
!  ******  IF L IS NOT REQUIRED, THEN COMMENT OUT THE FOLLOWING LINE  **
          l(irl(i)) = -lki
          somme = somme + lki * tmp(i)
          jmin = iu(i)
          jmax = iu(i+1) - 1
          IF (jmin .GT. jmax) GO TO 9
          mu = iju(i) - jmin
          DO 8 j=jmin,jmax
   8        row(ju(mu+j)) = row(ju(mu+j)) + lki * u(j)
   9      i = jrl(i)
          IF (i .NE. 0) GO TO 7
!
!  ******  ASSIGN KTH ROW OF U AND DIAGONAL D, SET TMP(K)  *************
  10    IF (row(k) .EQ. 0.0d0) GO TO 108
        dk = 1.0d0 / row(k)
        d(k) = dk
        tmp(k) = somme * dk
        IF (k .EQ. n) GO TO 19
        jmin = iu(k)
        jmax = iu(k+1) - 1
        IF (jmin .GT. jmax)  GO TO 12
        mu = iju(k) - jmin
        DO 11 j=jmin,jmax
  11      u(j) = row(ju(mu+j)) * dk
  12    CONTINUE
!
!  ******  UPDATE IRL AND JRL, KEEPING JRL IN DECREASING ORDER  ********
        i = i1
        IF (i .EQ. 0) GO TO 18
  14    irl(i) = irl(i) + 1
        i1 = jrl(i)
        IF (irl(i) .GE. il(i+1)) GO TO 17
        ijlb = irl(i) - il(i) + ijl(i)
        j = jl(ijlb)
  15    IF (i .GT. jrl(j)) GO TO 16
          j = jrl(j)
          GO TO 15
  16    jrl(i) = jrl(j)
        jrl(j) = i
  17    i = i1
        IF (i .NE. 0) GO TO 14
  18    IF (irl(k) .GE. il(k+1)) GO TO 19
        j = jl(ijl(k))
        jrl(k) = jrl(j)
        jrl(j) = k
  19    CONTINUE
!
!  ******  SOLVE  UX = TMP  BY BACK SUBSTITUTION  **********************
      k = n
      DO 22 i=1,n
        somme =  tmp(k)
        jmin = iu(k)
        jmax = iu(k+1) - 1
        IF (jmin .GT. jmax)  GO TO 21
        mu = iju(k) - jmin
        DO 20 j=jmin,jmax
  20      somme = somme - u(j) * tmp(ju(mu+j))
  21    tmp(k) =  somme
        z(c(k)) =  somme
  22    k = k-1
      flag = 0
      RETURN
!
! ** ERROR.. INSUFFICIENT STORAGE FOR L
 104  flag = 4*n + 1
      RETURN
! ** ERROR.. INSUFFICIENT STORAGE FOR U
 107  flag = 7*n + 1
      RETURN
! ** ERROR.. ZERO PIVOT
 108  flag = 8*n + k
      RETURN
      END