Initial commit from Polytechnique Montreal

1 files changed, 475 insertions, 0 deletions

diff --git a/Trivac/src/FLDSMB.f b/Trivac/src/FLDSMB.f
new file mode 100755
index 0000000..f479a9a
--- /dev/null
+++ b/Trivac/src/FLDSMB.f
@@ -0,0 +1,475 @@
+*DECK FLDSMB
+      SUBROUTINE FLDSMB (IPTRK,IPSYS,IPFLUX,LL4,ITY,NUN,NGRP,ICL1,ICL2,
+     1 IMPX,IMPH,TITR,EPS2,MAXOUT,MAXINR,EPSINR,EVECT,FKEFF)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solution of a multigroup eigenvalue system for the calculation of the
+* direct neutron flux in BIVAC. Use the preconditionned power method
+* with a two-parameter SVAT acceleration technique.
+*
+*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
+* IPTRK   L_TRACK pointer to the BIVAC tracking information.
+* IPSYS   L_SYSTEM pointer to system matrices.
+* IPFLUX  L_FLUX pointer to the solution.
+* LL4     order of the system matrices.
+* ITY     type of algorithm: 1: Diffusion theory; 11: Simplified PN
+*         approximation.
+* NUN     number of unknowns in each energy group.
+* NGRP    number of energy groups.
+* ICL1    number of free iterations in one cycle of the inverse power
+*         method.
+* ICL2    number of accelerated iterations in one cycle.
+* IMPX    print parameter: =0: no print; =1: minimum printing;
+*         =2: iteration history is printed.
+* IMPH    type of histogram processing:
+*         =0: no action is taken;
+*         =1: the flux is compared to a reference flux stored on LCM;
+*         =2: the convergence histogram is printed;
+*         =3: the convergence histogram is printed with axis and
+*            titles. The plotting file is completed;
+*         =4: the convergence histogram is printed with axis, acce-
+*            leration factors and titles. The plotting file is
+*            completed.
+* TITR    title.
+* EPS2    convergence criteria for the flux.
+* MAXOUT  maximum number of outer iterations.
+* MAXINR  maximum number of thermal iterations.
+* EPSINR  thermal iteration epsilon.
+* EVECT   initial estimate of the unknown vector.
+*
+*Parameters: output
+* EVECT   converged unknown vector.
+* FKEFF   effective multiplication factor.
+*
+*Reference:
+* A. H\'ebert, 'Preconditioning the power method for reactor
+* calculations', Nucl. Sci. Eng., 94, 1 (1986).
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK,IPSYS,IPFLUX
+      CHARACTER*72 TITR
+      INTEGER LL4,ITY,NUN,NGRP,ICL1,ICL2,IMPX,IMPH,MAXOUT,MAXINR
+      REAL FKEFF,EPS2,EPSINR,EVECT(NUN,NGRP)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (MMAXX=250,EPS1=1.0E-5)
+      CHARACTER*12 TEXT12
+      LOGICAL LOGTES
+      DOUBLE PRECISION AEAE,AEAG,AEAH,AGAG,AGAH,AHAH,BEBE,BEBG,BEBH,
+     1 BGBG,BGBH,BHBH,AEBE,AEBG,AEBH,AGBE,AGBG,AGBH,AHBE,AHBG,AHBH,
+     2 X,DXDA,DXDB,Y,DYDA,DYDB,Z,DZDA,DZDB,F,D2F(2,3),EVAL,ALP,BET,
+     3 FMIN
+      INTEGER ITITR(18)
+      REAL ERR(MMAXX),ALPH(MMAXX),BETA(MMAXX)
+      DOUBLE PRECISION, PARAMETER :: ALP_TAB(24) = (/ 0.2, 0.4, 0.6,
+     1  0.8, 1.0, 1.2, 1.5, 2.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0,
+     2  40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0 /)
+      DOUBLE PRECISION, PARAMETER :: BET_TAB(11) = (/ -1.0, -0.8, -0.6,
+     1 -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 /)
+      REAL, DIMENSION(:), ALLOCATABLE :: WORK
+      REAL, DIMENSION(:,:), ALLOCATABLE :: GRAD1,GRAD2,GAR1,GAR2,GAR3,
+     1 GAF1,GAF2,GAF3
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(GRAD1(NUN,NGRP),GRAD2(NUN,NGRP),GAR1(NUN,NGRP),
+     1 GAR2(NUN,NGRP),GAR3(NUN,NGRP),GAF1(NUN,NGRP),GAF2(NUN,NGRP),
+     2 GAF3(NUN,NGRP),WORK(NUN))
+*
+*     TKT : CPU TIME FOR THE SOLUTION OF LINEAR SYSTEMS.
+*     TKB : CPU TIME FOR BILINEAR PRODUCT EVALUATIONS.
+      TKT=0.0
+      TKB=0.0
+      CALL KDRCPU(TK1)
+*----
+*  PRECONDITIONED POWER METHOD
+*----
+      EVAL=1.0
+      VVV=0.0
+      ISTART=1
+      TEST=0.0
+      IF(IMPX.GE.1) WRITE (6,600)
+      IF(IMPX.GE.2) WRITE (6,610)
+      DO 25 IGR=1,NGRP
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EVECT(1,IGR),GAR1(1,IGR))
+      DO 20 JGR=1,NGRP
+      IF(JGR.EQ.IGR) GO TO 20
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+      IF(ILONG.EQ.0) GO TO 20
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EVECT(1,JGR),WORK(1))
+      DO 10 I=1,LL4
+      GAR1(I,IGR)=GAR1(I,IGR)-WORK(I)
+   10 CONTINUE
+   20 CONTINUE
+   25 CONTINUE
+      CALL KDRCPU(TK2)
+      TKB=TKB+(TK2-TK1)
+*
+      M=0
+   30 M=M+1
+*----
+*  EIGENVALUE EVALUATION
+*----
+      CALL KDRCPU(TK1)
+      AEBE=0.0D0
+      BEBE=0.0D0
+      DO 75 IGR=1,NGRP
+      DO 40 I=1,LL4
+      GAF1(I,IGR)=0.0
+   40 CONTINUE
+      DO 60 JGR=1,NGRP
+      WRITE(TEXT12,'(1HB,2I3.3)') IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+      IF(ILONG.EQ.0) GO TO 60
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EVECT(1,JGR),WORK(1))
+      DO 50 I=1,LL4
+      GAF1(I,IGR)=GAF1(I,IGR)+WORK(I)
+   50 CONTINUE
+   60 CONTINUE
+      DO 70 I=1,LL4
+      AEBE=AEBE+GAR1(I,IGR)*GAF1(I,IGR)
+      BEBE=BEBE+GAF1(I,IGR)**2
+   70 CONTINUE
+   75 CONTINUE
+      EVAL=AEBE/BEBE
+      CALL KDRCPU(TK2)
+      TKB=TKB+(TK2-TK1)
+*----
+*  DIRECTION EVALUATION
+*----
+      DO 110 IGR=1,NGRP
+      CALL KDRCPU(TK1)
+      DO 80 I=1,LL4
+      GRAD1(I,IGR)=REAL(EVAL)*GAF1(I,IGR)-GAR1(I,IGR)
+   80 CONTINUE
+      DO 100 JGR=1,IGR-1
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+      IF(ILONG.EQ.0) GO TO 100
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,JGR),WORK(1))
+      DO 90 I=1,LL4
+      GRAD1(I,IGR)=GRAD1(I,IGR)+WORK(I)
+   90 CONTINUE
+  100 CONTINUE
+      CALL KDRCPU(TK2)
+      TKB=TKB+(TK2-TK1)
+*
+      CALL KDRCPU(TK1)
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+      CALL MTLDLS(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,IGR))
+      CALL KDRCPU(TK2)
+      TKT=TKT+(TK2-TK1)
+  110 CONTINUE
+*----
+*  PERFORM THERMAL (UP-SCATTERING) ITERATIONS
+*----
+      KTER=0
+      NADI=5 ! used with SPN approximations
+      IF(MAXINR.GT.1) THEN
+         CALL FLDBHR(IPTRK,IPSYS,.FALSE.,LL4,ITY,NUN,NGRP,ICL1,ICL2,
+     1   IMPX,NADI,MAXINR,EPSINR,KTER,TKT,TKB,GRAD1)
+      ENDIF
+*----
+*  DISPLACEMENT EVALUATION
+*----
+      F=0.0D0
+      DELS=ABS(REAL((EVAL-VVV)/EVAL))
+      VVV=REAL(EVAL)
+      CALL KDRCPU(TK1)
+*----
+*  EVALUATION OF THE TWO ACCELERATION PARAMETERS ALP AND BET
+*----
+      ALP=1.0D0
+      BET=0.0D0
+      N=0
+      AEAE=0.0D0
+      AEAG=0.0D0
+      AEAH=0.0D0
+      AGAG=0.0D0
+      AGAH=0.0D0
+      AHAH=0.0D0
+      BEBG=0.0D0
+      BEBH=0.0D0
+      BGBG=0.0D0
+      BGBH=0.0D0
+      BHBH=0.0D0
+      AEBG=0.0D0
+      AEBH=0.0D0
+      AGBE=0.0D0
+      AGBG=0.0D0
+      AGBH=0.0D0
+      AHBE=0.0D0
+      AHBG=0.0D0
+      AHBH=0.0D0
+      DO 165 IGR=1,NGRP
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,IGR),GAR2(1,IGR))
+      DO 120 I=1,LL4
+      GAF2(I,IGR)=0.0
+  120 CONTINUE
+      DO 160 JGR=1,NGRP
+      IF(JGR.EQ.IGR) GO TO 140
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+      IF(ILONG.EQ.0) GO TO 140
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,JGR),WORK(1))
+      DO 130 I=1,LL4
+      GAR2(I,IGR)=GAR2(I,IGR)-WORK(I)
+  130 CONTINUE
+  140 WRITE(TEXT12,'(1HB,2I3.3)') IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+      IF(ILONG.EQ.0) GO TO 160
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,JGR),WORK(1))
+      DO 150 I=1,LL4
+      GAF2(I,IGR)=GAF2(I,IGR)+WORK(I)
+  150 CONTINUE
+  160 CONTINUE
+  165 CONTINUE
+      IF(1+MOD(M-ISTART,ICL1+ICL2).GT.ICL1) THEN
+         DO 175 IGR=1,NGRP
+         DO 170 I=1,LL4
+*        COMPUTE (A ,A )
+         AEAE=AEAE+GAR1(I,IGR)**2
+         AEAG=AEAG+GAR1(I,IGR)*GAR2(I,IGR)
+         AEAH=AEAH+GAR1(I,IGR)*GAR3(I,IGR)
+         AGAG=AGAG+GAR2(I,IGR)**2
+         AGAH=AGAH+GAR2(I,IGR)*GAR3(I,IGR)
+         AHAH=AHAH+GAR3(I,IGR)**2
+*        COMPUTE (B ,B )
+         BEBG=BEBG+GAF1(I,IGR)*GAF2(I,IGR)
+         BEBH=BEBH+GAF1(I,IGR)*GAF3(I,IGR)
+         BGBG=BGBG+GAF2(I,IGR)**2
+         BGBH=BGBH+GAF2(I,IGR)*GAF3(I,IGR)
+         BHBH=BHBH+GAF3(I,IGR)**2
+*        COMPUTE (A ,B )
+         AEBG=AEBG+GAR1(I,IGR)*GAF2(I,IGR)
+         AEBH=AEBH+GAR1(I,IGR)*GAF3(I,IGR)
+         AGBE=AGBE+GAR2(I,IGR)*GAF1(I,IGR)
+         AGBG=AGBG+GAR2(I,IGR)*GAF2(I,IGR)
+         AGBH=AGBH+GAR2(I,IGR)*GAF3(I,IGR)
+         AHBE=AHBE+GAR3(I,IGR)*GAF1(I,IGR)
+         AHBG=AHBG+GAR3(I,IGR)*GAF2(I,IGR)
+         AHBH=AHBH+GAR3(I,IGR)*GAF3(I,IGR)
+  170    CONTINUE
+  175    CONTINUE
+*
+  180    N=N+1
+         IF(N.GT.10) GO TO 185
+*        COMPUTE X(M+1)
+         X=BEBE+ALP*ALP*BGBG+BET*BET*BHBH+2.0D0*(ALP*BEBG+BET*BEBH
+     1   +ALP*BET*BGBH)
+         DXDA=2.0D0*(BEBG+ALP*BGBG+BET*BGBH)
+         DXDB=2.0D0*(BEBH+ALP*BGBH+BET*BHBH)
+*        COMPUTE Y(M+1)
+         Y=AEAE+ALP*ALP*AGAG+BET*BET*AHAH+2.0D0*(ALP*AEAG+BET*AEAH
+     1   +ALP*BET*AGAH)
+         DYDA=2.0D0*(AEAG+ALP*AGAG+BET*AGAH)
+         DYDB=2.0D0*(AEAH+ALP*AGAH+BET*AHAH)
+*        COMPUTE Z(M+1)
+         Z=AEBE+ALP*ALP*AGBG+BET*BET*AHBH+ALP*(AEBG+AGBE)
+     1   +BET*(AEBH+AHBE)+ALP*BET*(AGBH+AHBG)
+         DZDA=AEBG+AGBE+2.0D0*ALP*AGBG+BET*(AGBH+AHBG)
+         DZDB=AEBH+AHBE+ALP*(AGBH+AHBG)+2.0D0*BET*AHBH
+*        COMPUTE F(M+1)
+         F=X*Y-Z*Z
+         D2F(1,1)=2.0D0*(BGBG*Y+DXDA*DYDA+X*AGAG-DZDA**2-2.0D0*Z*AGBG)
+         D2F(1,2)=2.0D0*BGBH*Y+DXDA*DYDB+DXDB*DYDA+2.0D0*X*AGAH
+     1   -2.0D0*DZDA*DZDB-2.0D0*Z*(AGBH+AHBG)
+         D2F(2,2)=2.0D0*(BHBH*Y+DXDB*DYDB+X*AHAH-DZDB**2-2.0D0*Z*AHBH)
+         D2F(2,1)=D2F(1,2)
+         D2F(1,3)=DXDA*Y+X*DYDA-2.0D0*Z*DZDA
+         D2F(2,3)=DXDB*Y+X*DYDB-2.0D0*Z*DZDB
+*        SOLUTION OF A LINEAR SYSTEM.
+         CALL ALSBD(2,1,D2F,IER,2)
+         IF(IER.NE.0) GO TO 185
+         ALP=ALP-D2F(1,3)
+         BET=BET-D2F(2,3)
+         IF(ALP.GT.100.0) GO TO 185
+         IF((ABS(D2F(1,3)).LE.1.0D-4).AND.(ABS(D2F(2,3)).LE.1.0D-4))
+     1   GO TO 190
+         GO TO 180
+*
+*        alternative algorithm in case of Newton-Raphton failure
+  185    IF(IMPX.GT.0) WRITE(6,'(/30H FLDSMB: FAILURE OF THE NEWTON,
+     1   55H-RAPHTON ALGORIHTHM FOR COMPUTING THE OVERRELAXATION PA,
+     2   9HRAMETERS.)')
+         IAMIN=999
+         IBMIN=999
+         FMIN=HUGE(FMIN)
+         DO IA=1,SIZE(ALP_TAB)
+           ALP=ALP_TAB(IA)
+           DO IB=1,SIZE(BET_TAB)
+             BET=BET_TAB(IB)
+*            COMPUTE X
+             X=BEBE+ALP*ALP*BGBG+BET*BET*BHBH+2.0D0*(ALP*BEBG+BET*BEBH
+     1       +ALP*BET*BGBH)
+*            COMPUTE Y
+             Y=AEAE+ALP*ALP*AGAG+BET*BET*AHAH+2.0D0*(ALP*AEAG+BET*AEAH
+     1       +ALP*BET*AGAH)
+*            COMPUTE Z
+             Z=AEBE+ALP*ALP*AGBG+BET*BET*AHBH+ALP*(AEBG+AGBE)
+     1       +BET*(AEBH+AHBE)+ALP*BET*(AGBH+AHBG)
+*            COMPUTE F
+             F=X*Y-Z*Z
+             IF(F.LT.FMIN) THEN
+               IAMIN=IA
+               IBMIN=IB
+               FMIN=F
+             ENDIF
+           ENDDO
+         ENDDO
+         ALP=ALP_TAB(IAMIN)
+         BET=BET_TAB(IBMIN)
+  190    BET=BET/ALP
+         IF((ALP.LT.1.0D0).AND.(ALP.GT.0.0D0)) THEN
+            ALP=1.0D0
+            BET=0.0D0
+         ELSE IF(ALP.LE.0.0D0) THEN
+            ISTART=M+1
+            ALP=1.0D0
+            BET=0.0D0
+         ENDIF
+         DO 205 IGR=1,NGRP
+         DO 200 I=1,LL4
+         GRAD1(I,IGR)=REAL(ALP)*(GRAD1(I,IGR)+REAL(BET)*GRAD2(I,IGR))
+         GAR2(I,IGR)=REAL(ALP)*(GAR2(I,IGR)+REAL(BET)*GAR3(I,IGR))
+         GAF2(I,IGR)=REAL(ALP)*(GAF2(I,IGR)+REAL(BET)*GAF3(I,IGR))
+  200    CONTINUE
+  205    CONTINUE
+      ENDIF
+      CALL KDRCPU(TK2)
+      TKB=TKB+(TK2-TK1)
+*
+      LOGTES=(M.LT.ICL1).OR.(MOD(M-ISTART,ICL1+ICL2).EQ.ICL1-1)
+      IF(LOGTES.AND.(DELS.LE.EPS1)) THEN
+         DELT=0.0
+         DO 220 IGR=1,NGRP
+         DELN=0.0
+         DELD=0.0
+         DO 210 I=1,LL4
+         EVECT(I,IGR)=EVECT(I,IGR)+GRAD1(I,IGR)
+         GAR1(I,IGR)=GAR1(I,IGR)+GAR2(I,IGR)
+         GAF1(I,IGR)=GAF1(I,IGR)+GAF2(I,IGR)
+         GRAD2(I,IGR)=GRAD1(I,IGR)
+         GAR3(I,IGR)=GAR2(I,IGR)
+         GAF3(I,IGR)=GAF2(I,IGR)
+         DELN=MAX(DELN,ABS(GAF2(I,IGR)))
+         DELD=MAX(DELD,ABS(GAF1(I,IGR)))
+  210    CONTINUE
+         IF(DELD.NE.0.0) DELT=MAX(DELT,DELN/DELD)
+  220    CONTINUE
+         IF(IMPX.GE.2) WRITE (6,615) M,AEAE,AEAG,AEAH,AGAG,AGAH,AHAH,
+     1   BEBE,ALP,BET,EVAL,F,DELS,DELT,N,BEBG,BEBH,BGBG,BGBH,BHBH,
+     2   AEBE,AEBG,AEBH,AGBE,AGBG,AGBH,AHBE,AHBG,AHBH
+*        COMPUTE THE CONVERGENCE HISTOGRAM.
+         IF((IMPH.GE.1).AND.(M.LE.MMAXX)) THEN
+            CALL FLDXCO(IPFLUX,LL4,NUN,EVECT(1,NGRP),.TRUE.,ERR(M))
+            ALPH(M)=REAL(ALP)
+            BETA(M)=REAL(BET)
+         ENDIF
+         IF(DELT.LE.EPS2) GO TO 240
+      ELSE
+         DO 235 IGR=1,NGRP
+         DO 230 I=1,LL4
+         EVECT(I,IGR)=EVECT(I,IGR)+GRAD1(I,IGR)
+         GAR1(I,IGR)=GAR1(I,IGR)+GAR2(I,IGR)
+         GAF1(I,IGR)=GAF1(I,IGR)+GAF2(I,IGR)
+         GRAD2(I,IGR)=GRAD1(I,IGR)
+         GAR3(I,IGR)=GAR2(I,IGR)
+         GAF3(I,IGR)=GAF2(I,IGR)
+  230    CONTINUE
+  235    CONTINUE
+         IF(IMPX.GE.2) WRITE (6,620) M,AEAE,AEAG,AEAH,AGAG,AGAH,AHAH,
+     1   BEBE,ALP,BET,EVAL,F,DELS,N,BEBG,BEBH,BGBG,BGBH,BHBH,AEBE,
+     2   AEBG,AEBH,AGBE,AGBG,AGBH,AHBE,AHBG,AHBH
+*        COMPUTE THE CONVERGENCE HISTOGRAM.
+         IF((IMPH.GE.1).AND.(M.LE.MMAXX)) THEN
+            CALL FLDXCO(IPFLUX,LL4,NUN,EVECT(1,NGRP),.TRUE.,ERR(M))
+            ALPH(M)=REAL(ALP)
+            BETA(M)=REAL(BET)
+         ENDIF
+      ENDIF
+      IF(M.EQ.1) TEST=DELS
+      IF((M.GT.5).AND.(DELS.GT.TEST)) CALL XABORT('FLDSMB: CONVERGENCE'
+     1 //' FAILURE.')
+      IF(M.GE.MAXOUT) CALL XABORT('FLDSMB: MAXIMUM NUMBER OF ITERATION'
+     1 //'S REACHED.')
+      GO TO 30
+*----
+*  SOLUTION EDITION
+*----
+  240 FKEFF=REAL(1.0D0/EVAL)
+      IF(IMPX.EQ.1) WRITE (6,640) M
+      IF(IMPX.GE.1) THEN
+         WRITE (6,650) TKT,TKB,TKT+TKB
+         WRITE (6,670) FKEFF
+      ENDIF
+      IF(IMPX.EQ.3) THEN
+         DO 250 IGR=1,NGRP
+         WRITE (6,680) IGR,(EVECT(I,IGR),I=1,LL4)
+  250    CONTINUE
+      ENDIF
+      IF(IMPH.EQ.1) THEN
+         CALL LCMLEN(IPFLUX,'REF',ILONG,ITYLCM)
+         IF(ILONG.EQ.0) THEN
+            WRITE(6,'(40H FLDSMB: STORE A REFERENCE THERMAL FLUX.)')
+            CALL LCMPUT(IPFLUX,'REF',NUN,2,EVECT(1,NGRP))
+         ENDIF
+      ELSE IF(IMPH.GE.2) THEN
+         IGRAPH=0
+  260    IGRAPH=IGRAPH+1
+         WRITE (TEXT12,'(5HHISTO,I3)') IGRAPH
+         CALL LCMLEN (IPFLUX,TEXT12,ILENG,ITYLCM)
+         IF(ILENG.EQ.0) THEN
+            MM=MIN(M,MMAXX)
+            READ (TITR,'(18A4)') ITITR
+            CALL LCMSIX (IPFLUX,TEXT12,1)
+            CALL LCMPUT (IPFLUX,'HTITLE',18,3,ITITR)
+            CALL LCMPUT (IPFLUX,'ALPHA',MM,2,ALPH)
+            CALL LCMPUT (IPFLUX,'BETA',MM,2,BETA)
+            CALL LCMPUT (IPFLUX,'ERROR',MM,2,ERR)
+            CALL LCMPUT (IPFLUX,'IMPH',1,1,IMPH)
+            CALL LCMSIX (IPFLUX,' ',2)
+         ELSE
+            GO TO 260
+         ENDIF
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(GRAD1,GRAD2,GAR1,GAR2,GAR3,GAF1,GAF2,GAF3,WORK)
+      RETURN
+*
+  600 FORMAT(1H1/50H FLDSMB: ITERATIVE PROCEDURE BASED ON PRECONDITION,
+     1 16HED POWER METHOD./9X,16HDIRECT EQUATION.)
+  610 FORMAT(//5X,17HBILINEAR PRODUCTS,48X,5HALPHA,3X,4HBETA,3X,
+     1 12HEIGENVALUE..,12X,8HACCURACY,11(1H.),2X,1HN)
+  615 FORMAT(1X,I3,1P,7E9.1,0P,2F8.3,E14.6,3E10.2,I4/(4X,1P,7E9.1))
+  620 FORMAT(1X,I3,1P,7E9.1,0P,2F8.3,E14.6,2E10.2,10X,I4/(4X,1P,7E9.1))
+  640 FORMAT(/23H FLDSMB: CONVERGENCE IN,I4,12H ITERATIONS.)
+  650 FORMAT(/53H FLDSMB: CPU TIME USED TO SOLVE THE TRIANGULAR LINEAR,
+     1 10H SYSTEMS =,F10.3/23X,34HTO COMPUTE THE BILINEAR PRODUCTS =,
+     2 F10.3,20X,16HTOTAL CPU TIME =,F10.3)
+  670 FORMAT(//42H FLDSMB: EFFECTIVE MULTIPLICATION FACTOR =,1P,E17.10/)
+  680 FORMAT(//47H FLDSMB: EIGENVECTOR CORRESPONDING TO THE GROUP,I4
+     1 //(5X,1P,8E14.5))
+      END