Initial commit from Polytechnique Montreal

1 files changed, 121 insertions, 0 deletions

diff --git a/Dragon/src/PIJCMP.f b/Dragon/src/PIJCMP.f
new file mode 100644
index 0000000..f62fb2e
--- /dev/null
+++ b/Dragon/src/PIJCMP.f
@@ -0,0 +1,121 @@
+*DECK PIJCMP
+      SUBROUTINE PIJCMP(NREG,NSOUT,NCOR,DPR,VOLSUR,LPIJK,PROB)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compression of PIJ matrices in symmetric format.
+*
+*Copyright:
+* Copyright (C) 1994 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): R. Roy
+*
+*Parameters: input
+* NREG    number of zones for geometry.
+* NSOUT   number of surfaces for geometry. 
+* NCOR    maximum number of corners.
+* DPR     collision probabilities.
+* VOLSUR  volumes.
+* LPIJK   pijk flag.
+*
+*Parameters: output
+* PROB    compress probability matrix.
+*
+*Comments:
+*  Format of compress probability matrix
+*    NPLEN=(NREG+NSOUT+2)*(NREG+NSOUT+1)/2
+*    IND(I,J)=MAX(I+NSOUT+1,J+NSOUT+1)
+*            *(MAX(I+NSOUT+1,J+NSOUT+1)-1)/2
+*            +MIN(I+NSOUT+1,J+NSOUT+1)
+*    IS=-NSOUT,-1; JS=-NSOUT,IS; I=IND(IS,JS)
+*      PROB(I)=VOLSUR(IS)*PSS(IS,JS)
+*    IV=1,NREG; JS=-NSOUT,-1;    I=IND(IV,JS)
+*      SIGT(IV).GT.0.0
+*        PROB(I)=SIGT(IV)*VOLSUR(IV)*PVS(IV,JS)
+*      SIGT(IV).EQ.0.0
+*        PROB(I)=VOLSUR(IV)*PVS(IV,JS)
+*     IV=1,NREG; JV=1,IV;       I=IND(IV,JV)
+*       SIGT(IV).GT.0.0 AND SIGT(JV).GT.0.0
+*         PROB(I)=SIGT(IV)*SIGT(JV)*VOLSUR(IV)*PVV(IV,JV)
+*       SIGT(IV).GT.0.0 AND SIGT(JV).EQ.0.0
+*         PROB(I)=SIGT(IV)*VOLSUR(IV)*PVV(IV,JV)
+*       SIGT(IV).EQ.0.0 AND SIGT(JV).GT.0.0
+*         PROB(I)=SIGT(JV)*VOLSUR(IV)*PVV(IV,JV)
+*       SIGT(IV).EQ.0.0 AND SIGT(JV).EQ.0.0
+*         PROB(I)=VOLSUR(IV)*PVV(IV,JV)
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*----
+* VARIABLES
+*----
+      INTEGER          NREG,NSOUT,NCOR,IPR,IL,JL,IVOL,IUN
+      DOUBLE PRECISION DPR(-NSOUT:NREG,-NSOUT:NREG),PROB(*),ZCOR,ZCOR1
+      REAL             VOLSUR(-NSOUT:NREG),COEF
+      LOGICAL          LPIJK
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: WORK
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(WORK((NREG+NSOUT+2)*(NREG+NSOUT+1)/2))
+*----
+*  SYMMETRIZE AND STORE IN PROB
+*----
+      IPR= 0
+      DO  150 IL    = -NSOUT, NREG
+         DO  160 JL = -NSOUT, IL
+            IPR= IPR+1
+            WORK(IPR)= DPR(IL,JL) + DPR(JL,IL)
+  160    CONTINUE
+  150 CONTINUE
+      IF( NCOR.EQ.1 )THEN
+         IPR= 0
+         DO  250 IL    = -NSOUT, NREG
+            DO  260 JL = -NSOUT, IL
+               IPR= IPR+1
+               PROB(IPR)= WORK(IPR)
+  260       CONTINUE
+  250    CONTINUE
+      ELSE
+         IPR= 0
+         ZCOR1= 1.0D0/DBLE(NCOR)
+         ZCOR=  1.0D0/DBLE(NCOR*NCOR)
+         DO  251 IL    = -NSOUT, NREG
+            IF( IL.GT.0 ) ZCOR= ZCOR1
+            DO  261 JL = -NSOUT, IL
+               IPR= IPR+1
+               IF( JL.GT.0 ) ZCOR= 1.0D0
+               PROB(IPR)= ZCOR * WORK(IPR)
+  261       CONTINUE
+  251    CONTINUE
+      ENDIF
+*----
+*  CHARGE VOLUMES IN THE PROB MATRIX
+*----
+      COEF=1.0
+      IVOL= NSOUT*(NSOUT+1)/2
+      DO 300 IUN= -NSOUT, NREG
+         IF( IUN.LE.0 )THEN
+            IVOL= IVOL+1
+            IF(LPIJK) COEF= 3./4.
+         ELSE
+            IVOL= IVOL+NSOUT+IUN
+            IF(LPIJK) COEF= 2./3.
+         ENDIF
+         IF( PROB(IVOL).NE.0.0 )THEN
+            CALL XABORT( 'PIJCMP: UNEXPECTED VALUE IN PROB MATRIX' )
+         ENDIF
+         PROB(IVOL) = VOLSUR(IUN)*COEF
+  300 CONTINUE
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(WORK)
+      RETURN
+      END