Replaced STOP by RETURN before error print statements in src/msspec/spec/fortran/eig/mi/do_main.f

In src/msspec/spec/fortran/eig/common/, modified eig_mat_ms.f to call subroutines in new
files diagonalize_matrix.f and renormalization.f to implement renormalization in the
eigenvalue "spectroscopy"

src/msspec/spec/fortran/eig/common/diagonalize_matrix.f

@@ -0,0 +1,53 @@
+c
+c=======================================================================
+c
+c   This version: Kevin Dunseath, 9 December 2019
+c
+      subroutine diag_mat (n, a, lda, w, info)
+c
+      use outunits_mod, only: iuo1
+c
+      implicit none
+c
+      integer,    intent(in)  :: n, lda
+      integer,    intent(out) :: info
+      complex*16, intent(in)  :: a(lda,*)
+      complex*16, intent(out) :: w(*)
+c
+c Local variables
+c
+      integer    :: lwork
+      complex*16 :: wquery
+      complex*16 :: vl(1,1), vr(1,1)
+c
+      real*8,     allocatable :: rwork(:)
+      complex*16, allocatable :: work(:)
+c
+c
+      info = 0
+c
+      allocate(rwork(2*n))
+c
+c Get optimal workspace
+c
+      lwork = -1
+      call zgeev('n','n',n,a,lda,w,vl,1,vr,1,wquery,lwork,rwork,info)
+c
+      if (info.ne.0) then
+         write(iuo1,*) '                        '
+         write(iuo1,*) '    --->  work(1),info =',wquery,info
+         write(iuo1,*) '                        '
+      end if
+c
+      lwork = int(wquery)
+      allocate(work(lwork))
+c
+      call zgeev('n','n',n,a,lda,w,vl,1,vr,1,work,lwork,rwork,info)
+c
+      deallocate(work,rwork)
+c
+      return
+      end subroutine diag_mat
+c
+c=======================================================================
+c

src/msspec/spec/fortran/eig/common/eig_mat_ms.f

@@ -17,18 +17,22 @@ C
       USE OUTFILES_MOD
       USE OUTUNITS_MOD
       USE TRANS_MOD
+CKMD
+      USE RENORM_MOD
+CKMD
 C
 C!      PARAMETER(NLTWO=2*NL_M) !Moved to DIM_MOD
 C
       CHARACTER*24 OUTFILE,PATH
 C
       COMPLEX*16 HL1(0:NLTWO),SM(LINMAX*NATCLU_M,LINMAX*NATCLU_M)
-      COMPLEX*16 SUM_L,IC,ZEROC,WORK(32*LINMAX*NATCLU_M)
+CKMD      COMPLEX*16 SUM_L,IC,ZEROC,WORK(32*LINMAX*NATCLU_M)
+      COMPLEX*16 SUM_L,IC,ZEROC
       COMPLEX*16 YLM(0:NLTWO,-NLTWO:NLTWO),TLK,EXPKJ
       COMPLEX*16 W(LINMAX*NATCLU_M)
-      COMPLEX*16 VL(1,1),VR(1,1)
+CKMD      COMPLEX*16 VL(1,1),VR(1,1)
 C
-      DOUBLE PRECISION RWORK(2*LINMAX*NATCLU_M)
+CKMD      DOUBLE PRECISION RWORK(2*LINMAX*NATCLU_M)
 C
       REAL*8 PI,ATTKJ,GNT(0:N_GAUNT),XKJ,YKJ,ZKJ,RKJ,ZDKJ,KRKJ
 C
@@ -64,6 +68,8 @@ C
 C  Construction of the multiple scattering kernel matrix G_o T.
 C       Elements are stored using a linear index LINJ
 C                    representing (J,LJ)
+CKMD
+      SM = CMPLX(0.0D0, 0.0D0)
 C
       JLIN=0
       DO JTYP=1,N_PROT
@@ -139,16 +145,33 @@ C
       ENDDO
 C
       N_DIM=LINMAX*NATCLU_M
+C
+      IF (I_REN.gt.0) THEN
+C
+CKMD Renormalize the matrix SM
+C
+          CALL RENORM_MATRIX(JLIN,SM,N_DIM)
+C
+CKMD SM now contains the renormalized matrix
+C
+      END IF
 C
 C   Eigenvalues of the kernel multiple scattering matrix SM
 C
-      CALL ZGEEV('N','N',JLIN,SM,N_DIM,W,VL,1,VR,1,WORK,34*N_DIM,RWORK,
+CKMDC      CALL ZGEEV('N','N',JLIN,SM,N_DIM,W,VL,1,VR,1,WORK,32*N_DIM,RWORK,
-     &INFO)
+CKMD          LWORK = 32*N_DIM
-      IF(INFO.NE.0) THEN
+CKMD          CALL ZGEEV('N','N',JLIN,SM,N_DIM,W,VL,1,VR,1,WORK,LWORK,
-          WRITE(IUO1,*) '                        '
+CKMD     &RWORK,INFO)
-          WRITE(IUO1,*) '    --->  WORK(1),INFO =',WORK(1),INFO
+CKMD
-          WRITE(IUO1,*) '                        '
+      CALL DIAG_MAT(JLIN,SM,N_DIM,W,INFO)
-      ENDIF
+CKMD
+CKMD SM has been overwritten here
+C
+CKMD      IF(INFO.NE.0) THEN
+CKMD          WRITE(IUO1,*) '                        '
+CKMD          WRITE(IUO1,*) '    --->  WORK(1),INFO =',WORK(1),INFO
+CKMD          WRITE(IUO1,*) '                        '
+CKMD      ENDIF
 C
       N_EIG=0
 C
@@ -182,6 +205,10 @@ C
       CALL ORDRE(JLIN,W1,NFIN,W2)
 C
 C
+CKMD
+      WRITE(IUO1,10)
+      WRITE(IUO1,12) JLIN
+CKMD
       WRITE(IUO1,10)
       WRITE(IUO1,10)
       WRITE(IUO1,15) W2(1)
@@ -213,8 +240,19 @@ C
       ENDDO
       WRITE(IUO1,10)
       WRITE(IUO1,10)
-      WRITE(IUO1,45) W2(1)
+CKMD
-      WRITE(IUO2,*) E_KIN,W2(1)
+      IF (I_REN.NE.0) THEN
+         WRITE(IUO1,46) REN_R, REN_I
+         WRITE(IUO1,47) W2(1)
+      ELSE
+         WRITE(IUO1,45) W2(1)
+      ENDIF
+CKMD
+      IF (I_REN.NE.0) THEN
+         WRITE(IUO2,*) E_KIN,W2(1),REN_R,REN_I
+      ELSE
+         WRITE(IUO2,*) E_KIN,W2(1)
+      ENDIF
       IF(N_EIG.EQ.0) THEN
         WRITE(IUO1,50)
       ELSE
@@ -226,9 +264,13 @@ C
 C
       RETURN
 C
-   5  FORMAT(/,11X,'-----------------  EIGENVALUE ANALYSIS ','---------
+CKMD
-     &--------')
+   5  FORMAT(/,11X,'-----------------  EIGENVALUE ANALYSIS ',
+     &'-----------------')
   10  FORMAT(11X,'-',54X,'-')
+CKMD
+  12  FORMAT(11X,'-',14X,'MATRIX DIMENSION : ',I8,13X,'-')
+CKMD
   15  FORMAT(11X,'-',14X,'MAXIMUM MODULUS : ',F9.6,13X,'-')
   20  FORMAT(11X,'-',14X,'MINIMUM MODULUS : ',F9.6,13X,'-')
   25  FORMAT(11X,'-',6X,'1 EIGENVALUE IS > 1 ON A TOTAL OF ',I8,6X,'-')
@@ -238,12 +280,18 @@ C
   40  FORMAT(11X,'-',6X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,5X,'-')
   45  FORMAT(11X,'-',5X,'SPECTRAL RADIUS OF THE KERNEL MATRIX :',F6.3,
      &5X,'-')
+CKMD
+  46  FORMAT(11X,'-',16X,'OMEGA = (',F6.3,',',F6.3,')',15X,'-')
+  47  FORMAT(11X,'-',2X,'SPECTRAL RADIUS OF THE RENORMALIZED MATRIX: ',
+     &F6.3,2X,'-')
+CKMD
   50  FORMAT(11X,'-',5X,'---> THE MULTIPLE SCATTERING SERIES ',
      &'CONVERGES',4X,'-')
   55  FORMAT(11X,'-',10X,'---> NO CONVERGENCE OF THE MULTIPLE',9X,'-',/
      &,11X,'-',18X,'SCATTERING SERIES',19X,'-')
-  60  FORMAT(11X,'----------------------------------------','----------
+CKMD
-     &------',/)
+  60  FORMAT(11X,'----------------------------------------',
+     &'----------------',/)
   65  FORMAT(11X,'-',5X,'    LABEL OF LARGEST EIGENVALUE   : ',I5,8X,'-
      &')
   70  FORMAT(11X,'-',5X,'    LARGEST EIGENVALUE  : ','(',F6.3,',',F6.3,

src/msspec/spec/fortran/eig/common/renormalization.f

@@ -0,0 +1,140 @@
+c
+c=======================================================================
+c
+      subroutine renorm_matrix (n, a, lda)
+c
+c  This subroutine computes the renormalized matrices square matrices
+c     B_n : G_n, Sigma_n, Z_n, Pi_1.
+c
+c  The general renormalization scheme is given by:
+c
+c           (I - A)^{-1} = (I - M)^{-1} N
+c
+c     where matrix N is given by:
+c
+c     I_REN = 1 : N = REN2 * I with REN2 = REN**N_REN
+c           = 2 : N = REN2 * I with REN2 = (ONEC-REN**(N_REN+1))/(DFLOAT(N_REN+1)*(ONEC-REN))
+c           = 3 : N = REN2 * I with REN2 = -(REN-ONEC)**(N_REN+1)
+c           = 4 : N = I + REN * A
+c
+c
+c  Input parameters:
+c
+c       * N        : size of A
+c       * A        : original matrix
+c       * A2       : A * A
+c
+c  Input parameters:
+c
+c       * A        : renormalized matrix
+c
+c
+c  COMMON /RENORM/:
+c
+c     I_REN = 1 : renormalization in terms of the B_n = G_n matrices (n : N_REN)
+c           = 2 : renormalization in terms of the B_n = Sigma_n matrices
+c           = 3 : renormalization in terms of the B_n = Z_n matrices
+c           = 4 : renormalization in terms of the B_n = Pi_1 matrix
+c
+c     N_REN = n
+c
+c     REN   = REN_R+IC*REN_I
+c
+c
+c  Using the renormalization coefficient REN = REN_R + i REN_I, they
+c    are defined as
+c
+c         I_REN = 1-3 : (I - A)^{-1} = REN2 * (I - B_n)^{-1}
+c
+c         I_REN = 4   : (I - A)^{-1} = (I - B_n)^{-1} * (I + REN * A)
+c
+c    which in turn implies
+c
+c         I_REN = 1-3 :  B_n = (1 - REN) * I + REN * A
+c
+c         I_REN = 4   :  B_n = I - (1 - REN) * A - REN * A2
+c
+c
+c   Author :  D. Sébilleau
+c
+c                                         Last modified : 23 Apr 2019
+c
+c   This version: Kevin Dunseath, 9 December 2019
+c
+      use renorm_mod, only: i_ren, n_ren, ren_r, ren_i
+c
+      implicit none
+c
+      integer,    intent(in)    :: n, lda
+      complex*16, intent(inout) :: a(lda,*)
+c
+c Local variables
+c
+      complex*16, parameter :: zero = (0.0d0,0.0d0)
+      complex*16, parameter :: onec = (1.0d0,0.0d0)
+      complex*16, parameter :: ic =   (0.0d0,1.0d0)
+c
+      integer    :: i, j
+      complex*16 :: ren, ren2
+c
+      complex*16, allocatable :: a2(:,:)
+c
+c
+      ren = dble(ren_r) + ic*dble(ren_i)
+c
+c  Computing the modified renormalization parameter REN2 (g_n,s_n,zeta_n)
+c
+      select case (i_ren)
+      case (1)
+c
+c.....(g_n,G_n) renormalization:  g_n = omega^n
+c
+         ren2 = ren**n_ren
+c
+      case (2)
+c
+c.....(s_{n},Sigma_n) renormalization:
+c
+         if (abs(onec-ren).lt.1.0d-10) then
+            ren2 = onec
+         else
+            ren2 = (onec-ren**(n_ren+1))/(dfloat(n_ren+1)*(onec-ren))
+         end if
+c
+      case (3)
+c
+c.....(zeta_{n},Z_n) renormalization
+c
+         ren2 = -(ren-onec)**(n_ren+1)
+c
+      case (4)
+c
+         ren2 = ren
+c
+      end select
+c
+c  Calculation of the renormalized matrix
+c
+      if (i_ren.le.3) then
+         do j=1,n
+            do i=1,n
+               a(i,j) = ren2*a(i,j)
+            end do
+            a(j,j) = a(j,j) + (onec-ren2)
+         end do
+      else if (i_ren.eq.4) then
+         allocate(a2(n,n))
+         call zgemm('n','n',n,n,n,onec,a,lda,a,lda,zero,a2,n)
+         do j=1,n
+            do i=1,n
+               a(i,j) = (onec-ren2)*a(i,j) + ren2*a2(i,j)
+            end do
+         end do
+         deallocate(a2)
+      end if
+c
+      return
+      end
+c
+c=======================================================================
+c

src/msspec/spec/fortran/eig/mi/do_main.f

@@ -1231,7 +1231,8 @@ c          ENDIF
 C
 C!      IF((ISOM.NE.0).OR.(NFICHLEC.EQ.1)) CLOSE(IUO1)
       IF(ISOM.NE.0) CLOSE(IUO2)
-      STOP
+CKMD      STOP
+      return
 C
    1  WRITE(IUO1,60)
       STOP