Initial commit from Polytechnique Montreal

1 files changed, 222 insertions, 0 deletions

diff --git a/Trivac/src/TRICHP.f b/Trivac/src/TRICHP.f
new file mode 100755
index 0000000..b518acf
--- /dev/null
+++ b/Trivac/src/TRICHP.f
@@ -0,0 +1,222 @@
+*DECK TRICHP
+      SUBROUTINE TRICHP(IEL,LX,LY,LZ,L4,MAT,KN,MUX,MUY,MUZ,IPY,IPZ,
+     1 IMPX)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Primal finite element unknown numbering for ADI solution in a 3D
+* domain. Compute the storage info for ADI matrices in compressed
+* diagona storage mode. Compute the ADI permutation vectors.
+*
+*Copyright:
+* Copyright (C) 2002 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): A. Hebert
+*
+*Parameters: input
+* IMPX    print parameter.
+* LX      number of elements along the X axis.
+* LY      number of elements along the Y axis.
+* LZ      number of elements along the Z axis.
+* IEL     degree of the Lagrangian finite elements. =1 (linear);
+*         =2 (parabolic); =3 (cubic); =4 (quartic).
+* L4      total number of unknown (variational coefficients) per
+*         energy group (order of system matrices).
+* MAT     mixture index assigned to each element.
+* KN      element-ordered unknown list.
+*
+*Parameters: output
+* MUX     X-directed compressed diagonal storage mode indices.
+* MUY     Y-directed compressed diagonal storage mode indices.
+* MUZ     Z-directed compressed diagonal storage mode indices.
+* IPY     Y-directed permutation vectors.
+* IPZ     Z-directed permutation vectors.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER IEL,LX,LY,LZ,L4,MAT(LX*LY*LZ),KN(LX*LY*LZ*(IEL+1)**3),
+     1 MUX(L4),MUY(L4),MUZ(L4),IPY(L4),IPZ(L4),IMPX
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER IJ1(125),IJ2(125),IJ3(125)
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: IWRK
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(IWRK(LX*IEL+1,LY*IEL+1,LZ*IEL+1))
+*
+      LC=IEL+1
+      LL=LC*LC*LC
+      DO 5 L=1,LL
+      L1=1+MOD(L-1,LC)
+      L2=1+(L-L1)/LC
+      L3=1+MOD(L2-1,LC)
+      IJ1(L)=L1
+      IJ2(L)=L3
+      IJ3(L)=1+(L2-L3)/LC
+    5 CONTINUE
+*----
+*  JUXTAPOSITION OF A CHECKERBOARD OVER A PLANE IN THE REACTOR
+*----
+      L2M=0
+      LZTOT=LZ*(LC-1)+1
+      LYTOT=LY*(LC-1)+1
+      LXTOT=LX*(LC-1)+1
+      DO 12 K=1,LZTOT
+      DO 11 J=1,LYTOT
+      DO 10 I=1,LXTOT
+      IWRK(I,J,K)=0
+   10 CONTINUE
+   11 CONTINUE
+   12 CONTINUE
+      NUM1=0
+      KEL=0
+      DO 32 K0=1,LZ
+      LK0=(K0-1)*(LC-1)
+      DO 31 K1=1,LY
+      LK1=(K1-1)*(LC-1)
+      DO 30 K2=1,LX
+      KEL=KEL+1
+      IF(MAT(KEL).EQ.0) GO TO 30
+      L2M=L2M+1
+      LK2=(K2-1)*(LC-1)
+      L=0
+      DO 22 IK0=LK0+1,LK0+LC
+      DO 21 IK1=LK1+1,LK1+LC
+      DO 20 IK2=LK2+1,LK2+LC
+      L=L+1
+      IND1=KN(NUM1+L)
+      IF(IND1.EQ.0) GO TO 20
+      IF(IWRK(IK2,IK1,IK0).EQ.0) THEN
+         IWRK(IK2,IK1,IK0)=IND1
+      ELSE IF(IWRK(IK2,IK1,IK0).NE.IND1) THEN
+         CALL XABORT('TRICHP: FAILURE OF THE RENUMBERING ALGORITHM(1).')
+      ENDIF
+   20 CONTINUE
+   21 CONTINUE
+   22 CONTINUE
+      NUM1=NUM1+LL
+   30 CONTINUE
+   31 CONTINUE
+   32 CONTINUE
+*----
+*  CALCULATION OF PERMUTATION VECTORS IPY AND IPZ
+*----
+      DO 40 I=1,L4
+      IPY(I)=0
+      IPZ(I)=0
+   40 CONTINUE
+      INEW=0
+      DO 52 K0=1,LZTOT
+      DO 51 K2=1,LXTOT
+      IF(IWRK(K2,1,K0).EQ.IWRK(K2,LC,K0)) THEN
+         K1MIN=1+LC/2
+      ELSE
+         K1MIN=1
+      ENDIF
+      DO 50 K1=K1MIN,LYTOT
+      I=IWRK(K2,K1,K0)
+      IF(I.EQ.0) GO TO 50
+      IF(IPY(I).EQ.0) THEN
+         INEW=INEW+1
+         IPY(I)=INEW
+      ENDIF
+   50 CONTINUE
+   51 CONTINUE
+   52 CONTINUE
+      IF(INEW.NE.L4) THEN
+         CALL XABORT('TRICHP: FAILURE OF THE RENUMBERING ALGORITHM(2).')
+      ENDIF
+      INEW=0
+      DO 72 K1=1,LYTOT
+      DO 71 K2=1,LXTOT
+      IF(IWRK(K2,K1,1).EQ.IWRK(K2,K1,LC)) THEN
+         K0MIN=1+LC/2
+      ELSE
+         K0MIN=1
+      ENDIF
+      DO 70 K0=K0MIN,LZTOT
+      I=IWRK(K2,K1,K0)
+      IF(I.EQ.0) GO TO 70
+      IF(IPZ(I).EQ.0) THEN
+         INEW=INEW+1
+         IPZ(I)=INEW
+      ENDIF
+   70 CONTINUE
+   71 CONTINUE
+   72 CONTINUE
+      IF(INEW.NE.L4) THEN
+         CALL XABORT('TRICHP: FAILURE OF THE RENUMBERING ALGORITHM(3).')
+      ENDIF
+*----
+*  CALCULATION OF VECTORS MUX, MUY AND MUZ
+*----
+      DO 100 I=1,L4
+      MUX(I)=1
+      MUY(I)=1
+      MUZ(I)=1
+  100 CONTINUE
+      NUM1=0
+      DO 130 K=1,L2M
+      DO 120 I=1,LL
+      INX1=KN(NUM1+I)
+      IF(INX1.EQ.0) GO TO 120
+      INY1=IPY(INX1)
+      INZ1=IPZ(INX1)
+      DO 110 J=1,LL
+      INX2=KN(NUM1+J)
+      IF(INX2.EQ.0) GO TO 110
+      INY2=IPY(INX2)
+      INZ2=IPZ(INX2)
+      IF((IJ2(I).EQ.IJ2(J)).AND.(IJ3(I).EQ.IJ3(J)).AND.(INX2.LT.INX1))
+     1 THEN
+         MUX(INX1)=MAX0(MUX(INX1),INX1-INX2+1)
+      ELSE IF((IJ1(I).EQ.IJ1(J)).AND.(IJ3(I).EQ.IJ3(J)).AND.
+     1 (INY2.LT.INY1)) THEN
+         MUY(INY1)=MAX0(MUY(INY1),INY1-INY2+1)
+      ELSE IF((IJ1(I).EQ.IJ1(J)).AND.(IJ2(I).EQ.IJ2(J)).AND.
+     1 (INZ2.LT.INZ1)) THEN
+         MUZ(INZ1)=MAX0(MUZ(INZ1),INZ1-INZ2+1)
+      ENDIF
+  110 CONTINUE
+  120 CONTINUE
+      NUM1=NUM1+LL
+  130 CONTINUE
+*
+      MUXMAX=0
+      MUYMAX=0
+      MUZMAX=0
+      IIMAXX=0
+      IIMAXY=0
+      IIMAXZ=0
+      DO 140 I=1,L4
+      MUXMAX=MAX(MUXMAX,MUX(I))
+      MUYMAX=MAX(MUYMAX,MUY(I))
+      MUZMAX=MAX(MUZMAX,MUZ(I))
+      IIMAXX=IIMAXX+MUX(I)
+      MUX(I)=IIMAXX
+      IIMAXY=IIMAXY+MUY(I)
+      MUY(I)=IIMAXY
+      IIMAXZ=IIMAXZ+MUZ(I)
+      MUZ(I)=IIMAXZ
+  140 CONTINUE
+      IF(IMPX.GT.0) WRITE (6,500) MUXMAX,MUYMAX,MUZMAX
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(IWRK)
+      RETURN
+*
+  500 FORMAT(/52H TRICHP: MAXIMUM BANDWIDTH FOR X-ORIENTED MATRICES =,
+     1 I4/27X,25HFOR Y-ORIENTED MATRICES =,I4/27X,16HFOR Z-ORIENTED M,
+     2 9HATRICES =,I4)
+      END