LCOV - code coverage report
Current view: top level - hybrid - wrapper.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 1292 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 84 0.0 %

          Line data    Source code
       1             : # define test
       2             : 
       3             :       module m_wrapper
       4             : 
       5             :       interface blockmat
       6             :       module procedure  blockmat_d,blockmat_z
       7             :       end interface
       8             : 
       9             :       interface packmat
      10             :       module procedure  packmat_d,packmat_z
      11             :       end interface
      12             : 
      13             :       interface packmatcoul
      14             :       module procedure  packmatcoul_d,packmatcoul_z
      15             :       end interface
      16             : 
      17             :       interface unpackmat
      18             :       module procedure  unpackmat_d,unpackmat_z
      19             :       end interface
      20             : 
      21             :       interface dotprod
      22             :       module procedure  dotprod_dd, dotprod_dz, dotprod_zd, dotprod_zz
      23             :       end interface
      24             : 
      25             :       interface matvec
      26             :       module procedure  matvec_dpd, matvec_dpz, matvec_zpd, matvec_zpz
      27             :       end interface
      28             : 
      29             :       interface matmat
      30             :       module procedure  
      31             :      &     matmat_dpdp, matmat_dpzp, matmat_zpdp, matmat_zpzp,
      32             :      &                  matmat_dmzp, matmat_dpzm, matmat_dmzm,
      33             :      &                  matmat_zmdp, matmat_zpdm, matmat_zmdm,
      34             :      &                  matmat_zmzp, matmat_zpzm, matmat_zmzm
      35             :       end interface
      36             : 
      37             :       interface matmatmatd
      38             :       module procedure  matmatmatd_ddd, matmatmatd_dzd
      39             :       end interface
      40             : 
      41             :       interface diagonalize
      42             : !       module procedure  diagonalize_de,  diagonalize_dv,  diagonalize_dpe,  diagonalize_dpv,
      43             : !      &                  diagonalize_ze,  diagonalize_zv,  diagonalize_zpe,  diagonalize_zpv,
      44             : !      &                  diagonalize_deo, diagonalize_dvo, diagonalize_dpeo, diagonalize_dpvo,
      45             : !      &                  diagonalize_zeo, diagonalize_zvo, diagonalize_zpeo, diagonalize_zpvo
      46             :       module procedure  
      47             :      +     diagonalize_de,    diagonalize_dv,    diagonalize_dpe,
      48             :      +     diagonalize_dpv,   diagonalize_ze,    diagonalize_zv,
      49             :      +     diagonalize_zpe,   diagonalize_zpv,   diagonalize_deo,
      50             :      +     diagonalize_dvo,   diagonalize_dpeo,  diagonalize_dpvo,
      51             :      +     diagonalize_zeo,   diagonalize_zvo,   diagonalize_zpeo, 
      52             :      +     diagonalize_zpvo,  diagonalize_dvs,   diagonalize_dvos,
      53             :      +     diagonalize_dpvs,  diagonalize_dpvos, diagonalize_zvs,
      54             :      +     diagonalize_zvos,  diagonalize_zpvs,  diagonalize_zpvos,
      55             :      +     diagonalize_dvx,   diagonalize_dvox,  diagonalize_dpvx,
      56             :      +     diagonalize_dpvox, diagonalize_zvx,   diagonalize_zvox,
      57             :      +     diagonalize_zpvx,  diagonalize_zpvox
      58             :       end interface
      59             : 
      60             :       interface geteigen
      61             :       module procedure  geteigen_zpvo
      62             :       end interface
      63             : 
      64             :       interface inverse
      65             :       module procedure  inverse_d,  inverse_dp,  inverse_z,  inverse_zp,
      66             :      &                  inverse_d1, inverse_dp1, inverse_z1, inverse_zp1
      67             :       end interface
      68             : 
      69             :       interface sqrtmat
      70             :       module procedure  sqrtmat_d,  sqrtmat_dp,  sqrtmat_z,  sqrtmat_zp,
      71             :      &                  sqrtmat_d1, sqrtmat_dp1, sqrtmat_z1, sqrtmat_zp1
      72             :       end interface
      73             :       
      74             :       contains
      75             : 
      76             : c     --------
      77             : 
      78           0 :       function identity(n)
      79             :       implicit none
      80             :       integer, intent(in)  :: n
      81             :       integer              :: identity(n,n)
      82             :       integer              :: i,j
      83           0 :       identity = 0
      84           0 :       do i = 1,n
      85           0 :         identity(i,i) = 1
      86             :       enddo
      87             :       end function identity
      88             : 
      89             : c     --------
      90             : 
      91           0 :       function blockmat_d(a,b)
      92             :       implicit none
      93             :       real, intent(in) :: a(:,:),b(:,:)
      94             :       real             :: blockmat_d(size(a,1)+size(b,1),
      95             :      +                                  size(a,1)+size(b,1))
      96             :       integer             :: na,nb
      97           0 :       na = size(a,1) ; nb = size(b,1)
      98           0 :       if(size(a,2).ne.na) 
      99           0 :      +  stop 'blockmat_d: dimensions of first array differ.'
     100           0 :       if(size(b,2).ne.nb)
     101           0 :      +  stop 'blockmat_d: dimensions of second array differ.'
     102           0 :       blockmat_d              = 0.0
     103           0 :       blockmat_d(  :na,  :na) = a
     104           0 :       blockmat_d(na+1:,na+1:) = b
     105             :       end function blockmat_d
     106             : 
     107           0 :       function blockmat_z(a,b)
     108             :       implicit none
     109             :       complex, intent(in) :: a(:,:),b(:,:)
     110             :       complex             :: blockmat_z(size(a,1)+size(b,1),
     111             :      +                                     size(a,1)+size(b,1))
     112             :       integer                :: na,nb
     113           0 :       na = size(a,1) ; nb = size(b,1)
     114           0 :       if(size(a,2).ne.na)
     115           0 :      +  stop 'blockmat_z: dimensions of first array differ.'
     116           0 :       if(size(b,2).ne.nb)
     117           0 :      +  stop 'blockmat_z: dimensions of second array differ.'
     118           0 :       blockmat_z              = 0.0
     119           0 :       blockmat_z(  :na,  :na) = a
     120           0 :       blockmat_z(na+1:,na+1:) = b
     121             :       end function blockmat_z
     122             : 
     123             : c     --------
     124             : 
     125           0 :       function packmat_d(mat)
     126             :       implicit none
     127             :       real, intent(in)  :: mat(:,:)
     128             :       real              :: packmat_d(size(mat,1)*(size(mat,1)+1)/2)
     129             :       integer              :: n,nn,i,j,k
     130           0 :       n = size(mat,1) ; nn = n*(n+1)/2
     131           0 :       if(size(mat,2).ne.n) stop 'packmat_d: array dimensions differ.'
     132           0 :       k = 0
     133           0 :       do j = 1,n
     134           0 :         do i = 1,j
     135           0 :           k            = k + 1
     136           0 :           packmat_d(k) = mat(i,j)
     137             : # ifdef test
     138           0 :           if(abs(mat(j,i)-mat(i,j)).gt.1e-8)
     139           0 :      +       STOP 'packmat_d: input matrix not symmetric'
     140             : # endif
     141             :         enddo
     142             :       enddo
     143             :       end function packmat_d
     144             :       
     145           0 :       function packmatcoul_d(mat)
     146             :       implicit none
     147             :       real, intent(in)  :: mat(:,:)
     148             :       real              :: packmatcoul_d(
     149             :      +                                    size(mat,1)*(size(mat,1)+1)/2)
     150             :       integer              :: n,nn,i,j,k
     151           0 :       n = size(mat,1) ; nn = n*(n+1)/2
     152           0 :       if(size(mat,2).ne.n) stop 'packmat_d: array dimensions differ.'
     153           0 :       k = 0
     154           0 :       do j = 1,n
     155           0 :         do i = 1,j
     156           0 :           k            = k + 1
     157             : 
     158           0 :           packmatcoul_d(k) = (mat(i,j)+mat(j,i))/2.
     159             : 
     160             : !           if(abs(mat(j,i)-mat(i,j)).gt.1e-6) then
     161             : !             write(*,*) 'packmatcoul_d: input matrix not symmetric; deviation .gt. 1E-06'         
     162             : !           endif
     163             :         enddo
     164             :       enddo
     165             :       end function packmatcoul_d
     166             : 
     167             :             
     168           0 :       function unpackmat_d(mat)
     169             :       implicit none
     170             :       real, intent(in)  :: mat(:)
     171             :       real              :: unpackmat_d(
     172             :      +                             nint(sqrt(0.25+2*size(mat))-0.5),
     173             :      &                             nint(sqrt(0.25+2*size(mat))-0.5))
     174             :       integer              :: n,nn,i,j,k
     175           0 :       nn = size(mat) ; n = nint(sqrt(0.25+2*nn)-0.5)
     176           0 :       k  = 0
     177           0 :       do j = 1,n
     178           0 :         do i = 1,j
     179           0 :           k                = k + 1
     180           0 :           unpackmat_d(i,j) = mat(k)
     181           0 :           unpackmat_d(j,i) = mat(k)
     182             :         enddo
     183             :       enddo
     184             :       end function unpackmat_d
     185             : 
     186             :       
     187             : 
     188           0 :       function packmat_z(mat)
     189             :       implicit none
     190             :       complex, intent(in) :: mat(:,:)
     191             :       complex             :: packmat_z(size(mat,1)*(size(mat,1)+1)/2)
     192             :       integer                :: n,nn,i,j,k
     193           0 :       n = size(mat,1) ; nn = n*(n+1)/2
     194           0 :       if(size(mat,2).ne.n) stop 'packmat_z: array dimensions differ.'
     195             :       k = 0
     196           0 :       do j = 1,n
     197           0 :         do i = 1,j
     198           0 :           k            = k + 1
     199           0 :           packmat_z(k) = mat(i,j)
     200             : # ifdef test
     201           0 :           if(abs(conjg(mat(j,i))-mat(i,j)).gt.1e-8) 
     202           0 :      +                    stop 'packmat_z: input matrix not Hermitian.'
     203             : # endif
     204             :         enddo
     205             :       enddo
     206             :       end function packmat_z
     207             : 
     208           0 :       function packmatcoul_z(mat)
     209             :       implicit none
     210             :       complex, intent(in)  :: mat(:,:)
     211             :       complex              :: packmatcoul_z(
     212             :      +                                    size(mat,1)*(size(mat,1)+1)/2)
     213             :       integer                 :: n,nn,i,j,k
     214           0 :       n = size(mat,1) ; nn = n*(n+1)/2
     215           0 :       if(size(mat,2).ne.n) stop 'packmat_z: array dimensions differ.'
     216             :       k = 0
     217           0 :       do j = 1,n
     218           0 :         do i = 1,j
     219           0 :           k            = k + 1
     220           0 :           packmatcoul_z(k) = (mat(i,j) + conjg(mat(j,i)))/2.
     221             : 
     222           0 :           if(abs(conjg(mat(j,i))-mat(i,j)).gt.1e-4) then
     223             :             stop 'packmatcoul_z: input matrix not Hermitian; &
     224           0 :      & deviation .gt. 1E-04.'
     225             :           endif
     226             :         enddo
     227             :       enddo
     228             :       end function packmatcoul_z
     229             :   
     230           0 :       function unpackmat_z(mat)
     231             :       implicit none
     232             :       complex, intent(in)  :: mat(:)
     233             :       complex              :: unpackmat_z(
     234             :      +                             nint(sqrt(0.25+2*size(mat))-0.5),
     235             :      &                             nint(sqrt(0.25+2*size(mat))-0.5))
     236             :       integer                 :: n,nn,i,j,k
     237           0 :       nn = size(mat) ; n = nint(sqrt(0.25+2*nn)-0.5)
     238           0 :       k  = 0
     239           0 :       do j = 1,n
     240           0 :         do i = 1,j
     241           0 :           k                = k + 1
     242           0 :           unpackmat_z(i,j) = mat(k)
     243           0 :           unpackmat_z(j,i) = conjg(mat(k))
     244             :         enddo
     245             :       enddo
     246             :       end function unpackmat_z
     247             : 
     248             : c     --------
     249             : 
     250           0 :       function dotprod_dd(vec1,vec2)
     251             :       implicit none
     252             :       real, intent(in) :: vec1(:),vec2(:)
     253             :       real             :: dotprod_dd
     254             :       integer             :: n
     255             :       real             :: ddot
     256           0 :       n = size(vec1) 
     257           0 :       if(size(vec2).ne.n) 
     258           0 :      +  stop 'dotprod_dd: sizes of input vectors differ.'
     259           0 :       dotprod_dd = ddot(n,vec1,1,vec2,1)
     260           0 :       end function dotprod_dd
     261             : 
     262           0 :       function dotprod_dz(vec1,vec2)
     263             :       implicit none
     264             :       real,    intent(in) :: vec1(:)
     265             :       complex, intent(in) :: vec2(:)
     266             :       complex             :: dotprod_dz
     267             :       integer                :: n
     268             :       real                :: ddot
     269           0 :       n = size(vec1) 
     270           0 :       if(size(vec2).ne.n) 
     271           0 :      +  stop 'dotprod_dz: sizes of input vectors differ.'
     272             :       dotprod_dz = ddot(n,vec1,1,real(vec2),1)
     273           0 :      +             +(0.0,1.0)*ddot(n,vec1,1,aimag(vec2),1)
     274           0 :       end function dotprod_dz
     275             : 
     276           0 :       function dotprod_zd(vec1,vec2)
     277             :       implicit none
     278             :       complex, intent(in) :: vec1(:)
     279             :       real,    intent(in) :: vec2(:)
     280             :       complex             :: dotprod_zd
     281             :       integer                :: n
     282             :       real                :: ddot
     283           0 :       n = size(vec1) 
     284           0 :       if(size(vec2).ne.n)
     285           0 :      +  stop 'dotprod_zd: sizes of input vectors differ.'
     286             :       dotprod_zd = ddot(n,real(vec1),1,vec2,1)
     287           0 :      +             -(0.0,1.0)*ddot(n,aimag(vec1),1,vec2,1)
     288           0 :       end function dotprod_zd
     289             : 
     290           0 :       function dotprod_zz(vec1,vec2)
     291             :       implicit none
     292             :       complex, intent(in) :: vec1(:),vec2(:)
     293             :       complex             :: dotprod_zz
     294             :       integer                :: n
     295             :       complex             :: zdotc
     296           0 :       n = size(vec1)
     297           0 :       if(size(vec2).ne.n)
     298           0 :      +  stop 'dotprod_zz: sizes of input vectors differ.'
     299           0 :       dotprod_zz = zdotc(n,vec1,1,vec2,1)
     300           0 :       end function dotprod_zz
     301             : 
     302             : c     --------
     303             : 
     304           0 :       function matvec_dpd(mat,vec)
     305             :       implicit none
     306             :       real, intent(in)  :: mat(:),vec(:)
     307             :       real              :: matvec_dpd(size(vec))
     308             :       integer              :: nn,n
     309           0 :       n  = size(vec)
     310           0 :       nn = n*(n+1)/2 
     311           0 :       if(size(mat).ne.nn) stop 'matvec_dpd: input array has wrong size.'
     312           0 :       call dspmv('U',n,1.0,mat,vec,1,0.0,matvec_dpd,1)
     313             :       end function matvec_dpd
     314             : 
     315           0 :       function matvec_dpz(mat,vec)
     316             :       implicit none
     317             :       real,    intent(in) :: mat(:)
     318             :       complex, intent(in) :: vec(:)
     319             :       complex             :: matvec_dpz(size(vec))
     320           0 :       real,   allocatable :: vecr(:),veci(:)
     321             :       integer                :: nn,n
     322           0 :       n  = size(vec) ; allocate ( vecr(n),veci(n) )
     323           0 :       nn = n*(n+1)/2 
     324           0 :       if(size(mat).ne.nn) stop 'matvec_dpz: input array has wrong size.'
     325           0 :       call dspmv('U',n,1.0,mat,real(vec),1,0.0,vecr,1)
     326           0 :       call dspmv('U',n,1.0,mat,aimag(vec),1,0.0,veci,1)
     327           0 :       matvec_dpz = vecr + (0.0,1.0) * veci
     328           0 :       deallocate ( vecr,veci )
     329             :       end function matvec_dpz
     330             : 
     331           0 :       function matvec_zpd(mat,vec)
     332             :       implicit none
     333             :       complex, intent(in) :: mat(:)
     334             :       real,    intent(in) :: vec(:)
     335             :       complex             :: matvec_zpd(size(vec))
     336           0 :       real,   allocatable :: vecr(:),veci(:)
     337             :       integer                :: nn,n
     338           0 :       n  = size(vec) ; allocate ( vecr(n),veci(n) )
     339           0 :       nn = n*(n+1)/2 
     340           0 :       if(size(mat).ne.nn) stop 'matvec_zpd: input array has wrong size.'
     341           0 :       call dspmv('U',n,1.0,real(mat),vec,1,0.0,vecr,1)
     342           0 :       call dspmv('U',n,1.0,aimag(mat),vec,1,0.0,veci,1)
     343           0 :       matvec_zpd = vecr + (0.0,1.0) * veci
     344           0 :       deallocate ( vecr,veci )
     345             :       end function matvec_zpd
     346             :       
     347           0 :       function matvec_zpz(mat,vec)
     348             :       implicit none
     349             :       complex, intent(in)  :: mat(:),vec(:)
     350             :       complex              :: matvec_zpz(size(vec))
     351             :       integer                 :: nn,n
     352           0 :       n  = size(vec)
     353           0 :       nn = n*(n+1)/2 
     354           0 :       if(size(mat).ne.nn) stop 'matvec_zpz: input array has wrong size.'
     355           0 :       call zhpmv('U',n,(1.0,0.0),mat,vec,1,(0.0,0.0),matvec_zpz,1)
     356             :       end function matvec_zpz
     357             : 
     358             : c     --------
     359             : 
     360           0 :       function matmat_dpdp(mat1,mat2)
     361             :       implicit none
     362             :       real, intent(in)  :: mat1(:),mat2(:)
     363             :       real              :: matmat_dpdp(
     364             :      +                            nint(sqrt(0.25+2*size(mat1))-0.5),
     365             :      &                            nint(sqrt(0.25+2*size(mat1))-0.5))
     366           0 :       real, allocatable :: vec(:),vec2(:)
     367             :       integer              :: nn,n,k1,i,j,k
     368           0 :       nn = size(mat1)
     369           0 :       n  = nint(sqrt(0.25+2*nn)-0.5) ; allocate ( vec(n),vec2(n) )
     370           0 :       if(size(mat2).ne.nn)
     371           0 :      +  stop 'matmat_dpdp: second input array has wrong size.'
     372             :       k  = 0
     373           0 :       do i = 1,n
     374           0 :         vec2(:i) = mat2(k+1:k+i)
     375           0 :         k1       = k+2*i
     376           0 :         do j = i+1,n
     377           0 :           vec2(j) = mat2(k1)
     378           0 :           k1      = k1 + j
     379             :         enddo
     380           0 :         call dspmv('U',n,1.0,mat1,vec2,1,0.0,vec,1)
     381           0 :         matmat_dpdp(:,i) = vec
     382           0 :         k = k + i
     383             :       enddo
     384           0 :       deallocate ( vec,vec2 )
     385             :       end function matmat_dpdp
     386             : 
     387           0 :       function matmat_dpzp(mat1,mat2)
     388             :       implicit none
     389             :       real,    intent(in)  :: mat1(:)
     390             :       complex, intent(in)  :: mat2(:)
     391             :       complex              :: matmat_dpzp(
     392             :      +                            nint(sqrt(0.25+2*size(mat1))-0.5),
     393             :      &                            nint(sqrt(0.25+2*size(mat1))-0.5))
     394           0 :       real,    allocatable :: vecr(:),veci(:)
     395           0 :       complex, allocatable :: vec2(:)
     396             :       integer                 :: nn,n,k1,i,j,k
     397           0 :       nn = size(mat1)
     398           0 :       n  = nint(sqrt(0.25+2*nn)-0.5) 
     399           0 :       allocate ( vecr(n),veci(n),vec2(n) )
     400           0 :       if(size(mat2).ne.nn)
     401           0 :      +  stop 'matmat_dpzp: second input array has wrong size.'
     402             :       k  = 0
     403           0 :       do i = 1,n
     404           0 :         vec2(:i) = mat2(k+1:k+i)
     405           0 :         k1       = k+2*i
     406           0 :         do j = i+1,n
     407           0 :           vec2(j) = conjg(mat2(k1))
     408           0 :           k1      = k1 + j
     409             :         enddo
     410           0 :         call dspmv('U',n,1.0,mat1,real(vec2),1,0.0,vecr,1)
     411           0 :         call dspmv('U',n,1.0,mat1,aimag(vec2),1,0.0,veci,1)
     412           0 :         matmat_dpzp(:,i) = vecr + (0.0,1.0) * veci
     413           0 :         k = k + i
     414             :       enddo
     415           0 :       deallocate ( vecr,veci,vec2 )
     416             :       end function matmat_dpzp
     417             : 
     418           0 :       function matmat_zpdp(mat1,mat2)
     419             :       implicit none
     420             :       complex, intent(in)  :: mat1(:)
     421             :       real,    intent(in)  :: mat2(:)
     422             :       complex              :: matmat_zpdp(
     423             :      +                            nint(sqrt(0.25+2*size(mat1))-0.5),
     424             :      &                            nint(sqrt(0.25+2*size(mat1))-0.5))
     425           0 :       real,    allocatable :: vecr(:),veci(:)
     426           0 :       complex, allocatable :: vec1(:)
     427             :       integer                 :: nn,n,k1,i,j,k
     428           0 :       nn = size(mat1)
     429           0 :       n  = nint(sqrt(0.25+2*nn)-0.5) 
     430           0 :       allocate ( vecr(n),veci(n),vec1(n) )
     431           0 :       if(size(mat2).ne.nn)
     432           0 :      +  stop 'matmat_zpdp: second input array has wrong size.'
     433             :       k  = 0
     434           0 :       do i = 1,n
     435           0 :         vec1(:i) = conjg(mat1(k+1:k+i))
     436           0 :         k1       = k+2*i
     437           0 :         do j = i+1,n
     438           0 :           vec1(j) = mat1(k1)
     439           0 :           k1      = k1 + j
     440             :         enddo
     441           0 :         call dspmv('U',n,1.0,mat2,real(vec1),1,0.0,vecr,1)
     442           0 :         call dspmv('U',n,1.0,mat2,aimag(vec1),1,0.0,veci,1)
     443           0 :         matmat_zpdp(i,:) = vecr + (0.0,1.0) * veci
     444           0 :         k = k + i
     445             :       enddo
     446           0 :       deallocate ( vecr,veci,vec1 )
     447             :       end function matmat_zpdp
     448             : 
     449           0 :       function matmat_zpzp(mat1,mat2)
     450             :       implicit none
     451             :       complex, intent(in)  :: mat1(:),mat2(:)
     452             :       complex              :: matmat_zpzp(
     453             :      +                            nint(sqrt(0.25+2*size(mat1))-0.5),
     454             :      &                            nint(sqrt(0.25+2*size(mat1))-0.5))
     455           0 :       complex, allocatable :: vec(:),vec2(:)
     456             :       integer                 :: nn,n,k1,i,j,k
     457           0 :       nn = size(mat1)
     458           0 :       n  = nint(sqrt(0.25+2*nn)-0.5) ; allocate ( vec(n),vec2(n) )
     459           0 :       if(size(mat2).ne.nn) 
     460           0 :      +  stop 'matmat_zpzp: second input array has wrong size.'
     461             :       k  = 0
     462           0 :       do i = 1,n
     463           0 :         vec2(:i) = mat2(k+1:k+i)
     464           0 :         k1       = k+2*i
     465           0 :         do j = i+1,n
     466           0 :           vec2(j) = conjg(mat2(k1))
     467           0 :           k1      = k1 + j
     468             :         enddo
     469           0 :         call zhpmv('U',n,(1.0,0.0),mat1,vec2,1,(0.0,0.0),vec,1)
     470           0 :         matmat_zpzp(:,i) = vec
     471           0 :         k = k + i
     472             :       enddo
     473           0 :       deallocate ( vec,vec2 )
     474             :       end function matmat_zpzp
     475             : 
     476           0 :       function matmat_dpdm(mat1,mat2)
     477             :       implicit none
     478             :       real, intent(in)  :: mat1(:),mat2(:,:)
     479             :       real              :: matmat_dpdm(size(mat2,1),size(mat2,1))
     480           0 :       real, allocatable :: vec(:),vec2(:)
     481             :       integer              :: nn,n,k1,i
     482           0 :       n = size(mat2,1) ; nn = n*(n+1)/2 ; allocate ( vec(n),vec2(n) )
     483           0 :       if(size(mat2,2).ne.n)  
     484           0 :      +  stop 'matmat_dpdm: dimensions of second input array differ.'
     485           0 :       if(size(mat1)  .ne.nn) 
     486           0 :      +  stop 'matmat_dpdm: first input array has wrong size.'
     487           0 :       do i = 1,n
     488           0 :         vec2 = mat2(:,i)
     489           0 :         call dspmv('U',n,1.0,mat1,vec2,1,0.0,vec,1)
     490           0 :         matmat_dpdm(:,i) = vec
     491             :       enddo
     492           0 :       deallocate ( vec,vec2 )
     493             :       end function matmat_dpdm
     494             : 
     495           0 :       function matmat_dmdp(mat1,mat2)
     496             :       implicit none
     497             :       real, intent(in)  :: mat1(:,:),mat2(:)
     498             :       real              :: matmat_dmdp(size(mat1,1),size(mat1,1))
     499           0 :       real, allocatable :: vec(:),vec2(:)
     500             :       integer              :: nn,n,k1,i
     501           0 :       n = size(mat1,1) ; nn = n*(n+1)/2 ; allocate ( vec(n),vec2(n) )
     502           0 :       if(size(mat1,2).ne.n)  
     503           0 :      +  stop 'matmat_dmdp: dimensions of first input array differ.'
     504           0 :       if(size(mat2)  .ne.nn) 
     505           0 :      +  stop 'matmat_dmdp: second input array has wrong size.'
     506           0 :       do i = 1,n
     507           0 :         vec2 = mat1(i,:)
     508           0 :         call dspmv('U',n,1.0,mat2,vec2,1,0.0,vec,1)
     509           0 :         matmat_dmdp(i,:) = vec
     510             :       enddo
     511           0 :       deallocate ( vec,vec2 )
     512             :       end function matmat_dmdp
     513             : 
     514           0 :       function matmat_dmdm(mat1,mat2)
     515             :       implicit none
     516             :       real, intent(in) :: mat1(:,:),mat2(:,:)
     517             :       real             :: matmat_dmdm(size(mat1,1),size(mat1,1))
     518             :       integer             :: n
     519           0 :       n = size(mat1,1)
     520           0 :       if(size(mat1,2).ne.n)
     521           0 :      +  stop 'matmat_dmdm: dimensions of first input array differ.'
     522           0 :       if(size(mat2,1).ne.n)
     523           0 :      +  stop 'matmat_dmdm: second input array has wrong dimensions.'
     524           0 :       if(size(mat2,2).ne.n)
     525           0 :      +  stop 'matmat_dmdm: dimensions of second input array differ.'
     526           0 :       call dgemm('N','N',n,n,n,1.0,mat1,n,mat2,n,0.0,matmat_dmdm,n)
     527             :       end function matmat_dmdm
     528             : 
     529           0 :       function matmat_dpzm(mat1,mat2)
     530             :       implicit none
     531             :       real,    intent(in)  :: mat1(:)
     532             :       complex, intent(in)  :: mat2(:,:)
     533             :       complex              :: matmat_dpzm(size(mat2,1),size(mat2,1))
     534           0 :       real,    allocatable :: vecr(:),veci(:)
     535           0 :       complex, allocatable :: vec2(:)
     536             :       integer                 :: nn,n,k1,i
     537           0 :       n = size(mat2,1) 
     538           0 :       nn = n*(n+1)/2 ; allocate ( vecr(n),veci(n),vec2(n) )
     539           0 :       if(size(mat2,2).ne.n)
     540           0 :      +  stop 'matmat_dpzm: dimensions of second input array differ.'
     541           0 :       if(size(mat1)  .ne.nn)
     542           0 :      +  stop 'matmat_dpzm: first input array has wrong size.'
     543           0 :       do i = 1,n
     544           0 :         vec2 = mat2(:,i)
     545           0 :         call dspmv('U',n,1.0,mat1,real(vec2),1,0.0,vecr,1)
     546           0 :         call dspmv('U',n,1.0,mat1,aimag(vec2),1,0.0,veci,1)
     547           0 :         matmat_dpzm(:,i) = vecr + (0.0,1.0) * veci
     548             :       enddo
     549           0 :       deallocate ( vecr,veci,vec2 )
     550             :       end function matmat_dpzm
     551             : 
     552           0 :       function matmat_dmzp(mat1,mat2)
     553             :       implicit none
     554             :       real,    intent(in)  :: mat1(:,:)
     555             :       complex, intent(in)  :: mat2(:)
     556             :       complex              :: matmat_dmzp(size(mat1,1),size(mat1,1))
     557           0 :       complex, allocatable :: vec1(:),vec(:)
     558             :       integer                 :: nn,n,k1,i
     559           0 :       n = size(mat1,1) ; nn = n*(n+1)/2 ; allocate ( vec(n),vec1(n) )
     560           0 :       if(size(mat1,2).ne.n)
     561           0 :      +  stop 'matmat_dmzp: dimensions of first input array differ.'
     562           0 :       if(size(mat2)  .ne.nn)
     563           0 :      +  stop 'matmat_dmzp: second input array has wrong size.'
     564           0 :       do i = 1,n
     565           0 :         vec1 = mat1(i,:)
     566           0 :         call zhpmv('U',n,(1.0,0.0),mat2,vec1,1,(0.0,0.0),vec,1)
     567           0 :         matmat_dmzp(i,:) = conjg(vec)
     568             :       enddo
     569           0 :       deallocate ( vec,vec1 )
     570             :       end function matmat_dmzp
     571             : 
     572           0 :       function matmat_dmzm(mat1,mat2)
     573             :       implicit none
     574             :       real,    intent(in) :: mat1(:,:)
     575             :       complex, intent(in) :: mat2(:,:)
     576             :       complex             :: matmat_dmzm(size(mat1,1),size(mat2,2))
     577           0 :       real                :: matr(size(mat1,1),size(mat2,2)),
     578           0 :      +                          mati(size(mat1,1),size(mat2,2))
     579             :       integer                :: n,n1,n2
     580           0 :       n1 = size(mat1,1)
     581           0 :       n  = size(mat1,2)
     582           0 :       n2 = size(mat2,2)
     583           0 :       if(size(mat2,1).ne.n)
     584           0 :      +  stop 'matmat_dmzm: dimensions of matrices are inconsistent.'
     585           0 :       call dgemm('N','N',n1,n2,n,1.0,mat1,n1,real(mat2),n,0.0,matr,n1)
     586           0 :       call dgemm('N','N',n1,n2,n,1.0,mat1,n1,aimag(mat2),n,0.0,mati,n1)
     587           0 :       matmat_dmzm = matr + (0.0,1.0) * mati
     588             :       end function matmat_dmzm
     589             : 
     590           0 :       function matmat_zpdm(mat1,mat2)
     591             :       implicit none
     592             :       complex, intent(in)  :: mat1(:)
     593             :       real,    intent(in)  :: mat2(:,:)
     594             :       complex              :: matmat_zpdm(size(mat2,1),size(mat2,1))
     595           0 :       complex, allocatable :: vec(:),vec2(:)
     596             :       integer                 :: nn,n,k1,i
     597           0 :       n = size(mat2,1) ; nn = n*(n+1)/2 ; allocate ( vec(n),vec2(n) )
     598           0 :       if(size(mat2,2).ne.n)
     599           0 :      +  stop 'matmat_zpdm: dimensions of second input array differ.'
     600           0 :       if(size(mat1)  .ne.nn)
     601           0 :      +  stop 'matmat_zpdm: first input array has wrong size.'
     602           0 :       do i = 1,n
     603           0 :         vec2 = mat2(:,i)
     604           0 :         call zhpmv('U',n,(1.0,0.0),mat1,vec2,1,(0.0,0.0),vec,1)
     605           0 :         matmat_zpdm(:,i) = vec
     606             :       enddo
     607           0 :       deallocate ( vec,vec2 )
     608             :       end function matmat_zpdm
     609             : 
     610           0 :       function matmat_zmdp(mat1,mat2)
     611             :       implicit none
     612             :       complex, intent(in)  :: mat1(:,:)
     613             :       real,    intent(in)  :: mat2(:)
     614             :       complex              :: matmat_zmdp(size(mat1,1),size(mat1,1))
     615           0 :       complex, allocatable :: vec1(:)
     616           0 :       real,    allocatable :: vecr(:),veci(:)
     617             :       integer                 :: nn,n,k1,i
     618           0 :       n = size(mat1,1) ; nn = n*(n+1)/2 
     619           0 :       allocate ( vecr(n),veci(n),vec1(n) )
     620           0 :       if(size(mat1,2).ne.n)
     621           0 :      +  stop 'matmat_zmdp: dimensions of first input array differ.'
     622           0 :       if(size(mat2)  .ne.nn)
     623           0 :      +  stop 'matmat_zmdp: second input array has wrong size.'
     624           0 :       do i = 1,n
     625           0 :         vec1 = conjg(mat1(i,:))
     626           0 :         call dspmv('U',n,1.0,mat2,real(vec1),1,0.0,vecr,1)
     627           0 :         call dspmv('U',n,1.0,mat2,aimag(vec1),1,0.0,veci,1)
     628           0 :         matmat_zmdp(i,:) = vecr - (0.0,1.0) * veci
     629             :       enddo
     630           0 :       deallocate ( vecr,veci,vec1 )
     631             :       end function matmat_zmdp
     632             : 
     633           0 :       function matmat_zmdm(mat1,mat2)
     634             :       implicit none
     635             :       complex, intent(in) :: mat1(:,:)
     636             :       real,    intent(in) :: mat2(:,:)
     637             :       complex             :: matmat_zmdm(size(mat1,1),size(mat2,2))
     638           0 :       real                :: matr(size(mat1,1),size(mat2,2)),
     639           0 :      +                          mati(size(mat1,1),size(mat2,2)) 
     640             :       integer                :: n,n1,n2
     641           0 :       n1 = size(mat1,1)
     642           0 :       n  = size(mat1,2)
     643           0 :       n2 = size(mat2,2)
     644           0 :       if(size(mat2,1).ne.n)
     645           0 :      +  stop 'matmat_zmdm: dimensions of matrices are inconsistent.'
     646           0 :       call dgemm('N','N',n1,n2,n,1.0,real(mat1),n1,mat2,n,0.0,matr,n1)
     647           0 :       call dgemm('N','N',n1,n2,n,1.0,aimag(mat1),n1,mat2,n,0.0,mati,n1)
     648           0 :       matmat_zmdm = matr + (0.0,1.0) * mati
     649             :       end function matmat_zmdm
     650             : 
     651           0 :       function matmat_zpzm(mat1,mat2)
     652             :       implicit none
     653             :       complex, intent(in)  :: mat1(:),mat2(:,:)
     654             :       complex              :: matmat_zpzm(size(mat2,1),size(mat2,2))
     655           0 :       complex, allocatable :: vec(:),vec2(:)
     656             :       integer                 :: nn,n,k1,i,n2
     657           0 :       n  = size(mat2,1) ; nn = n*(n+1)/2 ; allocate ( vec(n),vec2(n) )
     658           0 :       n2 = size(mat2,2)
     659           0 :       if(size(mat1).ne.nn)
     660           0 :      +  stop 'matmat_zpzm: first input array has wrong size.'
     661           0 :       do i = 1,n2
     662           0 :         vec2 = mat2(:,i)
     663           0 :         call zhpmv('U',n,(1.0,0.0),mat1,vec2,1,(0.0,0.0),vec,1)
     664           0 :         matmat_zpzm(:,i) = vec
     665             :       enddo
     666           0 :       deallocate ( vec,vec2 )
     667             :       end function matmat_zpzm
     668             : 
     669           0 :       function matmat_zmzp(mat1,mat2)
     670             :       implicit none
     671             :       complex, intent(in)  :: mat1(:,:),mat2(:)
     672             :       complex              :: matmat_zmzp(size(mat1,1),size(mat1,1))
     673           0 :       complex, allocatable :: vec(:),vec2(:)
     674             :       integer                 :: nn,n,k1,i
     675           0 :       n = size(mat1,1) ; nn = n*(n+1)/2 ; allocate ( vec(n),vec2(n) )
     676           0 :       if(size(mat1,2).ne.n)
     677           0 :      +  stop 'matmat_zmzp: dimensions of first input array differ.'
     678           0 :       if(size(mat2)  .ne.nn)
     679           0 :      +  stop 'matmat_zmzp: second input array has wrong size.'
     680           0 :       do i = 1,n
     681           0 :         vec2 = conjg(mat1(i,:))
     682           0 :         call zhpmv('U',n,(1.0,0.0),mat2,vec2,1,(0.0,0.0),vec,1)
     683           0 :         matmat_zmzp(i,:) = conjg(vec)
     684             :       enddo
     685           0 :       deallocate ( vec,vec2 )
     686             :       end function matmat_zmzp
     687             : 
     688           0 :       function matmat_zmzm(mat1,mat2)
     689             :       implicit none
     690             :       complex, intent(in) :: mat1(:,:),mat2(:,:)
     691             :       complex             :: matmat_zmzm(size(mat1,1),size(mat2,2))
     692             :       integer                :: n1,n,n2
     693             :       complex, parameter     :: one = (1,0), zero = 0
     694           0 :       n1 = size(mat1,1)
     695           0 :       n  = size(mat1,2)
     696           0 :       n2 = size(mat2,2)
     697           0 :       if(size(mat2,1).ne.n)
     698           0 :      +  stop 'matmat_zmzm: dimensions of matrices are inconsistent.'
     699           0 :       call zgemm('N','N',n1,n2,n,one,mat1,n1,mat2,n,zero,matmat_zmzm,n1)
     700             :       end function matmat_zmzm
     701             : 
     702             : c     --------
     703             : 
     704           0 :       function matmatmatd_ddd(diag1,mat,diag2)
     705             :       implicit none
     706             :       real, intent(in) :: diag1(:),mat(:,:),diag2(:)
     707             :       real             :: matmatmatd_ddd(size(diag1),size(diag1))
     708             :       integer             :: n,i
     709           0 :       n = size(diag1)
     710           0 :       if(size(mat,1).ne.n) 
     711           0 :      +  stop 'matmatmatd_ddd: input matrix has wrong size.'
     712           0 :       if(size(mat,2).ne.n)
     713           0 :      +  stop 'matmatmatd_ddd: dimensions of input matrix differ.'
     714           0 :       if(size(diag2).ne.n)
     715           0 :      +  stop 'matmatmatd_ddd: second diagonal matrix has wrong size.'
     716           0 :       do i = 1,n
     717           0 :         matmatmatd_ddd(:,i) = diag2(i) * mat(:,i)
     718             :       enddo
     719           0 :       do i = 1,n
     720           0 :         matmatmatd_ddd(i,:) = diag1(i) * matmatmatd_ddd(i,:)
     721             :       enddo
     722             :       end function matmatmatd_ddd
     723             : 
     724           0 :       function matmatmatd_dzd(diag1,mat,diag2)
     725             :       implicit none
     726             :       real,    intent(in) :: diag1(:),diag2(:)
     727             :       complex, intent(in) :: mat(:,:)
     728             :       complex             :: matmatmatd_dzd(size(diag1),size(diag1))
     729             :       integer                :: n,i
     730           0 :       n = size(diag1)
     731           0 :       if(size(mat,1).ne.n)
     732           0 :      +  stop 'matmatmatd_ddd: input matrix has wrong size.'
     733           0 :       if(size(mat,2).ne.n)
     734           0 :      +  stop 'matmatmatd_ddd: dimensions of input matrix differ.'
     735           0 :       if(size(diag2).ne.n)
     736           0 :      +  stop 'matmatmatd_ddd: second diagonal matrix has wrong size.'
     737           0 :       do i = 1,n
     738           0 :         matmatmatd_dzd(:,i) = diag2(i) * mat(:,i)
     739             :       enddo
     740           0 :       do i = 1,n
     741           0 :         matmatmatd_dzd(i,:) = diag1(i) * matmatmatd_dzd(i,:)
     742             :       enddo
     743             :       end function matmatmatd_dzd
     744             : 
     745             : c     --------
     746             :       
     747           0 :       subroutine diagonalize_de(eval,mat)
     748             :       implicit none
     749             :       real, intent(out) :: eval(:)
     750             :       real, intent(in)  :: mat(:,:)
     751           0 :       real, allocatable :: mat1(:,:),work(:)
     752             :       integer              :: n,info
     753           0 :       n = size(mat,1)
     754           0 :       if(n.eq.0)
     755           0 :      +  stop 'diagonalize_de: zero dimension in eigenvalue problem.'
     756           0 :       if(size(mat,2).ne.n)
     757           0 :      +  stop 'diagonalize_de: dimensions of input matrix differ.'
     758           0 :       if(size(eval) .ne.n)
     759           0 :      +  stop 'diagonalize_de: eigenvalue array has wrong size.'
     760           0 :       allocate ( mat1(n,n),work(3*n) ) ; mat1 = mat
     761           0 :       call dsyev('N','U',n,mat1,n,eval,work,3*n,info) 
     762           0 :       if(info.ne.0) stop 'diagonalize_de: dsyev failed.'
     763           0 :       deallocate ( mat1,work )
     764           0 :       end subroutine diagonalize_de
     765             : 
     766           0 :       subroutine diagonalize_dv(evec,eval,mat)
     767             :       implicit none
     768             :       real, intent(out) :: eval(:),evec(:,:)
     769             :       real, intent(in)  :: mat(:,:)
     770           0 :       real, allocatable :: work(:)
     771             :       integer              :: n,info
     772           0 :       n = size(mat,1)
     773           0 :       if(n.eq.0)
     774           0 :      +  stop 'diagonalize_dv: zero dimension in eigenvalue problem.'
     775           0 :       if(size(mat,2) .ne.n)
     776           0 :      +  stop 'diagonalize_dv: dimensions of input matrix differ.'
     777           0 :       if(size(eval)  .ne.n)
     778           0 :      +  stop 'diagonalize_dv: eigenvalue array has wrong size.'
     779           0 :       if(size(evec,1).ne.n)
     780           0 :      +  stop 'diagonalize_dv: eigenvector array has wrong dimensions.'
     781           0 :       if(size(evec,2).ne.n)
     782           0 :      +  stop 'diagonalize_dv: dimensions of eigenvector array differ.'
     783           0 :       allocate ( work(3*n) ) ; evec = mat
     784           0 :       call dsyev('V','U',n,evec,n,eval,work,3*n,info) 
     785           0 :       if(info.ne.0) stop 'diagonalize_dv: dsyev failed.'
     786           0 :       deallocate ( work )
     787           0 :       end subroutine diagonalize_dv
     788             : 
     789           0 :       subroutine diagonalize_dpe(eval,mat)
     790             :       implicit none
     791             :       real, intent(out) :: eval(:)
     792             :       real, intent(in)  :: mat(:)
     793           0 :       real, allocatable :: mat1(:),work(:)
     794             :       integer              :: n,nn,info
     795           0 :       n = size(eval,1) ; nn = n*(n+1)/2
     796           0 :       if(n.eq.0)
     797           0 :      +  stop 'diagonalize_dpe: zero dimension in eigenvalue problem.'
     798           0 :       if(size(mat).ne.nn)
     799           0 :      +  stop 'diagonalize_dpe: input matrix has wrong size.'
     800           0 :       allocate ( mat1(nn),work(3*n) ) ; mat1 = mat
     801           0 :       call dspev('N','U',n,mat1,eval,work,n,work,info) 
     802           0 :       if(info.ne.0) stop 'diagonalize_dpe: dspev failed.'
     803           0 :       deallocate ( mat1,work )
     804           0 :       end subroutine diagonalize_dpe
     805             : 
     806           0 :       subroutine diagonalize_dpv(evec,eval,mat)
     807             :       implicit none
     808             :       real, intent(out) :: eval(:),evec(:,:)
     809             :       real, intent(in)  :: mat(:)
     810           0 :       real, allocatable :: mat1(:),work(:)
     811             :       integer              :: n,nn,info
     812           0 :       n = size(eval,1) ; nn = n*(n+1)/2
     813           0 :       if(n.eq.0)
     814           0 :      +  stop 'diagonalize_dpv: zero dimension in eigenvalue problem.'
     815           0 :       if(size(mat)   .ne.nn)
     816           0 :      +  stop 'diagonalize_dpv: input matrix has wrong size.'
     817           0 :       if(size(evec,1).ne.n)
     818           0 :      +  stop 'diagonalize_dpv: eigenvector array has wrong dimensions.'
     819           0 :       if(size(evec,2).ne.n)
     820           0 :      +  stop 'diagonalize_dpv: dimensions of eigenvector array differ.'
     821           0 :       allocate ( mat1(nn),work(3*n) ) ; mat1 = mat
     822           0 :       call dspev('V','U',n,mat1,eval,evec,n,work,info) 
     823           0 :       if(info.ne.0) stop 'diagonalize_dpv: dspev failed.'
     824           0 :       deallocate ( mat1,work )
     825           0 :       end subroutine diagonalize_dpv
     826             : 
     827           0 :       subroutine diagonalize_ze(eval,mat)
     828             :       implicit none
     829             :       real,    intent(out) :: eval(:)
     830             :       complex, intent(in)  :: mat(:,:)
     831           0 :       complex, allocatable :: mat1(:,:),work(:)
     832           0 :       real,    allocatable :: rwork(:)
     833             :       integer                 :: n,info
     834           0 :       n = size(mat,1)
     835           0 :       if(n.eq.0)
     836           0 :      +  stop 'diagonalize_ze: zero dimension in eigenvalue problem.'
     837           0 :       if(size(mat,2).ne.n)
     838           0 :      +  stop 'diagonalize_ze: dimensions of input matrix differ.'
     839           0 :       if(size(eval) .ne.n)
     840           0 :      +  stop 'diagonalize_ze: eigenvalue array has wrong size.'
     841           0 :       allocate ( mat1(n,n),work(3*n),rwork(3*n) ) ; mat1 = mat
     842           0 :       call zheev('N','U',n,mat1,n,eval,work,3*n,rwork,info) 
     843           0 :       if(info.ne.0) stop 'diagonalize_ze: zheev failed.'
     844           0 :       deallocate ( mat1,work,rwork )
     845           0 :       end subroutine diagonalize_ze
     846             : 
     847           0 :       subroutine diagonalize_zv(evec,eval,mat)
     848             :       implicit none
     849             :       real,    intent(out) :: eval(:)
     850             :       complex, intent(out) :: evec(:,:)
     851             :       complex, intent(in)  :: mat(:,:)
     852           0 :       complex, allocatable :: work(:)
     853           0 :       real,    allocatable :: rwork(:)
     854             :       integer                 :: n,info
     855           0 :       n = size(mat,1)
     856           0 :       if(n.eq.0)
     857           0 :      +  stop 'diagonalize_zv: zero dimension in eigenvalue problem.'
     858           0 :       if(size(mat,2) .ne.n)
     859           0 :      +  stop 'diagonalize_zv: dimensions of input matrix differ.'
     860           0 :       if(size(eval)  .ne.n)
     861           0 :      +  stop 'diagonalize_zv: eigenvalue array has wrong size.'
     862           0 :       if(size(evec,1).ne.n)
     863           0 :      +  stop 'diagonalize_zv: eigenvector array has wrong dimensions.'
     864           0 :       if(size(evec,2).ne.n)
     865           0 :      +  stop 'diagonalize_zv: dimensions of eigenvector array differ.'
     866           0 :       allocate ( work(3*n),rwork(3*n) ) ; evec = mat
     867           0 :       call zheev('V','U',n,evec,n,eval,work,3*n,rwork,info) 
     868           0 :       if(info.ne.0) stop 'diagonalize_zv: zheev failed.'
     869           0 :       deallocate ( work,rwork )
     870           0 :       end subroutine diagonalize_zv
     871             : 
     872           0 :       subroutine diagonalize_zpe(eval,mat)
     873             :       implicit none
     874             :       real,    intent(out) :: eval(:)
     875             :       complex, intent(in)  :: mat(:)
     876           0 :       complex, allocatable :: mat1(:),work(:)
     877           0 :       real,    allocatable :: rwork(:)
     878             :       integer                 :: n,nn,info
     879           0 :       n = size(eval,1) ; nn = n*(n+1)/2
     880           0 :       if(n.eq.0)
     881           0 :      +  stop 'diagonalize_zpe: zero dimension in eigenvalue problem.'
     882           0 :       if(size(mat).ne.nn)
     883           0 :      +  stop 'diagonalize_zpe: input matrix has wrong size.'
     884           0 :       allocate ( mat1(nn),work(3*n),rwork(3*n) ) ; mat1 = mat
     885           0 :       call zhpev('N','U',n,mat1,eval,work,n,work,rwork,info) 
     886           0 :       if(info.ne.0) stop 'diagonalize_zpe: zhpev failed.'
     887           0 :       deallocate ( mat1,work,rwork )
     888           0 :       end subroutine diagonalize_zpe
     889             : 
     890           0 :       subroutine diagonalize_zpv(evec,eval,mat)
     891             :       implicit none
     892             :       real,    intent(out) :: eval(:)
     893             :       complex, intent(out) :: evec(:,:)
     894             :       complex, intent(in)  :: mat(:)
     895           0 :       complex, allocatable :: mat1(:),work(:)
     896           0 :       real,    allocatable :: rwork(:)
     897             :       integer                 :: n,nn,info
     898           0 :       n = size(eval,1) ; nn = n*(n+1)/2
     899           0 :       if(n.eq.0)
     900           0 :      +  stop 'diagonalize_zpv: zero dimension in eigenvalue problem.'
     901           0 :       if(size(mat)   .ne.nn)
     902           0 :      +  stop 'diagonalize_zpv: input matrix has wrong size.'
     903           0 :       if(size(evec,1).ne.n)
     904           0 :      +  stop 'diagonalize_zpv: eigenvector array has wrong dimensions.'
     905           0 :       if(size(evec,2).ne.n)
     906           0 :      +  stop 'diagonalize_zpv: dimensions of eigenvector array differ.'
     907           0 :       allocate ( mat1(nn),work(3*n),rwork(3*n) ) ; mat1 = mat 
     908           0 :       call zhpev('V','U',n,mat1,eval,evec,n,work,rwork,info) 
     909           0 :       if(info.ne.0) stop 'diagonalize_zpv: zhpev failed.'
     910           0 :       deallocate ( mat1,work,rwork )
     911           0 :       end subroutine diagonalize_zpv
     912             : 
     913           0 :       subroutine diagonalize_deo(eval,mat,olap)
     914             :       implicit none
     915             :       real, intent(out) :: eval(:)
     916             :       real, intent(in)  :: mat(:,:),olap(:,:)
     917           0 :       real, allocatable :: mat1(:,:),olap1(:,:),work(:)
     918             :       integer              :: n,info
     919           0 :       n = size(mat,1)
     920           0 :       if(n.eq.0)
     921           0 :      +  stop 'diagonalize_deo: zero dimension in eigenvalue problem.'
     922           0 :       if(size(mat,2) .ne.n)
     923           0 :      +  stop 'diagonalize_deo: dimensions of input matrix differ.'
     924           0 :       if(size(eval)  .ne.n)
     925           0 :      +  stop 'diagonalize_deo: eigenvalue array has wrong size.'
     926           0 :       if(size(olap,1).ne.n)
     927           0 :      +  stop 'diagonalize_deo: overlap matrix has wrong size.'
     928           0 :       if(size(olap,2).ne.n)
     929           0 :      +  stop 'diagonalize_deo: dimensions of overlap matrix differ.'
     930           0 :       allocate ( mat1(n,n),olap1(n,n),work(3*n) ) 
     931           0 :       mat1 = mat ; olap1 = olap
     932           0 :       call dsygv(1,'N','U',n,mat1,n,olap1,n,eval,work,3*n,info)
     933           0 :       if(info.ne.0) stop 'diagonalize_deo: dsygv failed.'
     934           0 :       deallocate ( mat1,olap1,work )
     935           0 :       end subroutine diagonalize_deo
     936             : 
     937           0 :       subroutine diagonalize_dvo(evec,eval,mat,olap)
     938             :       implicit none
     939             :       real, intent(out) :: eval(:),evec(:,:)
     940             :       real, intent(in)  :: mat(:,:),olap(:,:)
     941           0 :       real, allocatable :: olap1(:,:),work(:)
     942             :       integer              :: n,info
     943           0 :       n = size(mat,1)
     944           0 :       if(n.eq.0)
     945           0 :      +  stop 'diagonalize_dvo: zero dimension in eigenvalue problem.'
     946           0 :       if(size(mat,2) .ne.n)
     947           0 :      +  stop 'diagonalize_dvo: dimensions of input matrix differ.'
     948           0 :       if(size(eval)  .ne.n)
     949           0 :      +  stop 'diagonalize_dvo: eigenvalue array has wrong size.'
     950           0 :       if(size(evec,1).ne.n)
     951           0 :      +  stop 'diagonalize_dvo: eigenvector array has wrong dimensions.'
     952           0 :       if(size(evec,2).ne.n)
     953           0 :      +  stop 'diagonalize_dvo: dimensions of eigenvector array differ.'
     954           0 :       if(size(olap,1).ne.n)
     955           0 :      +  stop 'diagonalize_dvo: overlap matrix has wrong dimensions.'
     956           0 :       if(size(olap,2).ne.n)
     957           0 :      +  stop 'diagonalize_dvo: dimensions of overlap matrix differ.'
     958           0 :       allocate ( olap1(n,n),work(3*n) ) ; evec = mat ; olap1 = olap
     959           0 :       call dsygv(1,'V','U',n,evec,n,olap1,n,eval,work,3*n,info) 
     960           0 :       if(info.ne.0) stop 'diagonalize_dvo: dsygv failed.'
     961           0 :       deallocate ( olap1,work )
     962           0 :       end subroutine diagonalize_dvo
     963             : 
     964           0 :       subroutine diagonalize_dpeo(eval,mat,olap)
     965             :       implicit none
     966             :       real, intent(out) :: eval(:)
     967             :       real, intent(in)  :: mat(:),olap(:)
     968           0 :       real, allocatable :: mat1(:),olap1(:),work(:)
     969             :       integer              :: n,nn,info
     970           0 :       n = size(eval,1) ; nn = n*(n+1)/2
     971           0 :       if(n.eq.0)
     972           0 :      +  stop 'diagonalize_dpeo: zero dimension in eigenvalue problem.'
     973           0 :       if(size(mat) .ne.nn)
     974           0 :      +  stop 'diagonalize_dpeo: input matrix has wrong size.'
     975           0 :       if(size(olap).ne.nn)
     976           0 :      +  stop 'diagonalize_dpeo: overlap matrix has wrong size.'
     977           0 :       allocate ( mat1(nn),olap1(nn),work(3*n) ) 
     978           0 :       mat1 = mat ; olap1 = olap
     979           0 :       call dspgv(1,'N','U',n,mat1,olap1,eval,work,n,work,info) 
     980           0 :       if(info.ne.0) stop 'diagonalize_dpeo: dspgv failed.'
     981           0 :       deallocate ( mat1,olap1,work )
     982           0 :       end subroutine diagonalize_dpeo
     983             : 
     984           0 :       subroutine diagonalize_dpvo(evec,eval,mat,olap)
     985             :       implicit none
     986             :       real, intent(out) :: eval(:),evec(:,:)
     987             :       real, intent(in)  :: mat(:),olap(:)
     988           0 :       real, allocatable :: mat1(:),olap1(:),work(:)
     989             :       integer              :: n,nn,info
     990           0 :       n = size(eval,1) ; nn = n*(n+1)/2
     991           0 :       if(n.eq.0)
     992           0 :      +  stop 'diagonalize_dpvo: zero dimension in eigenvalue problem.'
     993           0 :       if(size(mat)   .ne.nn)
     994           0 :      +  stop 'diagonalize_dpvo: input matrix has wrong size.'
     995           0 :       if(size(olap)  .ne.nn)
     996           0 :      +  stop 'diagonalize_dpvo: overlap matrix has wrong size.'
     997           0 :       if(size(evec,1).ne.n)
     998           0 :      +  stop 'diagonalize_dpvo: eigenvector array has wrong dimensions.'
     999           0 :       if(size(evec,2).ne.n)
    1000           0 :      +  stop 'diagonalize_dpvo: dimensions of eigenvector array differ.'
    1001           0 :       allocate ( mat1(nn),olap1(nn),work(3*n) ) 
    1002           0 :       mat1 = mat ; olap1 = olap
    1003           0 :       call dspgv(1,'V','U',n,mat1,olap1,eval,evec,n,work,info) 
    1004           0 :       if(info.ne.0) stop 'diagonalize_dpvo: dspgv failed.'
    1005           0 :       deallocate ( mat1,olap1,work )
    1006           0 :       end subroutine diagonalize_dpvo
    1007             : 
    1008           0 :       subroutine diagonalize_zeo(eval,mat,olap)
    1009             :       implicit none
    1010             :       real,    intent(out) :: eval(:)
    1011             :       complex, intent(in)  :: mat(:,:),olap(:,:)
    1012           0 :       complex, allocatable :: mat1(:,:),olap1(:,:),work(:)
    1013           0 :       real,    allocatable :: rwork(:)
    1014             :       integer                 :: n,info
    1015           0 :       n = size(mat,1)
    1016           0 :       if(n.eq.0)
    1017           0 :      +  stop 'diagonalize_zeo: zero dimension in eigenvalue problem.'
    1018           0 :       if(size(mat,2) .ne.n)
    1019           0 :      +  stop 'diagonalize_zeo: dimensions of input matrix differ.'
    1020           0 :       if(size(eval)  .ne.n)
    1021           0 :      +  stop 'diagonalize_zeo: eigenvalue array has wrong size.'
    1022           0 :       if(size(olap,1).ne.n)
    1023           0 :      +  stop 'diagonalize_zeo: overlap matrix has wrong size.'
    1024           0 :       if(size(olap,2).ne.n)
    1025           0 :      +  stop 'diagonalize_zeo: dimensions of overlap matrix differ.'
    1026           0 :       allocate ( mat1(n,n),olap1(n,n),work(3*n),rwork(3*n) ) 
    1027           0 :       mat1 = mat ; olap1 = olap
    1028           0 :       call zhegv(1,'N','U',n,mat1,n,olap1,n,eval,work,3*n,rwork,info) 
    1029           0 :       if(info.ne.0) stop 'diagonalize_zeo: zhegv failed.'
    1030           0 :       deallocate ( mat1,olap1,work,rwork )
    1031           0 :       end subroutine diagonalize_zeo
    1032             : 
    1033           0 :       subroutine diagonalize_zvo(evec,eval,mat,olap)
    1034             :       implicit none
    1035             :       real,    intent(out) :: eval(:)
    1036             :       complex, intent(out) :: evec(:,:)
    1037             :       complex, intent(in)  :: mat(:,:),olap(:,:)
    1038           0 :       complex, allocatable :: olap1(:,:),work(:)
    1039           0 :       real,    allocatable :: rwork(:)
    1040             :       integer                 :: n,info
    1041           0 :       n = size(mat,1)
    1042           0 :       if(n.eq.0)
    1043           0 :      +  stop 'diagonalize_zvo: zero dimension in eigenvalue problem.'
    1044           0 :       if(size(mat,2) .ne.n)
    1045           0 :      +  stop 'diagonalize_zvo: dimensions of input matrix differ.'
    1046           0 :       if(size(eval)  .ne.n)
    1047           0 :      +  stop 'diagonalize_zvo: eigenvalue array has wrong size.'
    1048           0 :       if(size(evec,1).ne.n)
    1049           0 :      +  stop 'diagonalize_zvo: eigenvector array has wrong dimensions.'
    1050           0 :       if(size(evec,2).ne.n)
    1051           0 :      +  stop 'diagonalize_zvo: dimensions of eigenvector array differ.'
    1052           0 :       if(size(olap,1).ne.n)
    1053           0 :      +  stop 'diagonalize_zvo: overlap matrix has wrong dimensions.'
    1054           0 :       if(size(olap,2).ne.n)
    1055           0 :      +  stop 'diagonalize_zvo: dimensions of overlap matrix differ.'
    1056           0 :       allocate ( olap1(n,n),work(3*n),rwork(3*n) ) 
    1057           0 :       evec = mat ; olap1 = olap
    1058           0 :       call zhegv(1,'V','U',n,evec,n,olap1,n,eval,work,3*n,rwork,info) 
    1059           0 :       if(info.ne.0) stop 'diagonalize_zvo: zhegv failed.'
    1060           0 :       deallocate ( olap1,work,rwork )
    1061           0 :       end subroutine diagonalize_zvo
    1062             : 
    1063           0 :       subroutine diagonalize_zpeo(eval,mat,olap)
    1064             :       implicit none
    1065             :       real,    intent(out) :: eval(:)
    1066             :       complex, intent(in)  :: mat(:),olap(:)
    1067           0 :       complex, allocatable :: mat1(:),olap1(:),work(:)
    1068           0 :       real,    allocatable :: rwork(:)
    1069             :       integer                 :: n,nn,info
    1070           0 :       n = size(eval,1) ; nn = n*(n+1)/2
    1071           0 :       if(n.eq.0)
    1072           0 :      +  stop 'diagonalize_zpeo: zero dimension in eigenvalue problem.'
    1073           0 :       if(size(mat) .ne.nn)
    1074           0 :      +  stop 'diagonalize_zpeo: input matrix has wrong size.'
    1075           0 :       if(size(olap).ne.nn)
    1076           0 :      +  stop 'diagonalize_zpeo: overlap matrix has wrong size.'
    1077           0 :       allocate ( mat1(nn),olap1(nn),work(3*n),rwork(3*n) ) 
    1078           0 :       mat1 = mat ; olap1 = olap
    1079           0 :       call zhpgv(1,'N','U',n,mat1,olap1,eval,work,n,work,rwork,info) 
    1080           0 :       if(info.ne.0) stop 'diagonalize_zpeo: zhpev failed.'
    1081           0 :       deallocate ( mat1,olap1,work,rwork )
    1082           0 :       end subroutine diagonalize_zpeo
    1083             : 
    1084           0 :       subroutine diagonalize_zpvo(evec,eval,mat,olap)
    1085             :       implicit none
    1086             :       real,    intent(out) :: eval(:)
    1087             :       complex, intent(out) :: evec(:,:)
    1088             :       complex, intent(in)  :: mat(:),olap(:)
    1089           0 :       complex, allocatable :: mat1(:),olap1(:),work(:)
    1090           0 :       real,    allocatable :: rwork(:)
    1091             :       integer                 :: n,nn,info
    1092           0 :       n = size(eval,1) ; nn = n*(n+1)/2
    1093           0 :       if(n.eq.0)
    1094           0 :      +  stop 'diagonalize_zpvo: zero dimension in eigenvalue problem.'
    1095           0 :       if(size(mat)   .ne.nn)
    1096           0 :      +  stop 'diagonalize_zpvo: input matrix has wrong size.'
    1097           0 :       if(size(olap)  .ne.nn)
    1098           0 :      +  stop 'diagonalize_zpvo: overlap matrix has wrong size.'
    1099           0 :       if(size(evec,1).ne.n)
    1100           0 :      +  stop 'diagonalize_zpvo: eigenvector array has wrong dimensions.'
    1101           0 :       if(size(evec,2).ne.n)
    1102           0 :      +  stop 'diagonalize_zpvo: dimensions of eigenvector array differ.'
    1103           0 :       allocate ( mat1(nn),olap1(nn),work(3*n),rwork(3*n) ) 
    1104           0 :       mat1 = mat ; olap1 = olap
    1105           0 :       call zhpgv(1,'V','U',n,mat1,olap1,eval,evec,n,work,rwork,info) 
    1106           0 :       if(info.ne.0) stop 'diagonalize_zpvo: zhpgv failed.'
    1107           0 :       deallocate ( mat1,olap1,work,rwork )
    1108           0 :       end subroutine diagonalize_zpvo
    1109             : 
    1110           0 :       subroutine diagonalize_dvs(evec,eval,mat,m)
    1111             :       implicit none
    1112             :       real, intent(out) :: eval(:),evec(:,:)
    1113             :       real, intent(in)  :: mat(:,:)
    1114           0 :       real, allocatable :: work(:),mat1(:,:)
    1115             :       real              :: abstol,dlamch
    1116             :       integer, intent(in)  :: m
    1117           0 :       integer, allocatable :: iwork(:),ifail(:)
    1118             :       integer              :: n,ma,idum,info
    1119           0 :       ma = abs(m)
    1120           0 :       n  = size(mat,1)
    1121           0 :       if(n.eq.0)
    1122           0 :      +  stop 'diagonalize_dvs: zero dimension in eigenvalue problem.'
    1123           0 :       if(ma.gt.n) stop 'diagonalize_dvs: number of selected eigenvalues&
    1124           0 :      & exceeds maximal number.'
    1125           0 :       if(size(mat,2) .ne.n)
    1126           0 :      +  stop 'diagonalize_dvs: dimensions of input matrix differ.'
    1127           0 :       if(size(eval)  .ne.n)
    1128           0 :      +  stop 'diagonalize_dvs: eigenvalue array has wrong size.'
    1129           0 :       if(size(evec,1).ne.n)
    1130           0 :      +  stop 'diagonalize_dvs: first dimension of eigenvector is wrong.'
    1131           0 :       if(size(evec,2).ne.n)
    1132           0 :      + stop 'diagonalize_dvs: second dimension of eigenvector is wrong.'
    1133           0 :       allocate ( work(8*n),iwork(5*n),mat1(n,n),ifail(n) ) ; mat1 = mat
    1134           0 :       abstol = 2 * dlamch('S')
    1135           0 :       if(m.gt.0) then
    1136             :         call dsyevx('V','I','U',n,mat1,n,0.0,0.0,n-ma+1,n,abstol,idum,
    1137           0 :      +              eval,evec,n,work,8*n,iwork,ifail,info)
    1138             :       else
    1139             :         call dsyevx('V','I','U',n,mat1,n,0.0,0.0,    1,ma,abstol,idum,
    1140           0 :      +              eval,evec,n,work,8*n,iwork,ifail,info)
    1141             :       endif
    1142           0 :       if(info.ne.0) stop 'diagonalize_dvs: dsyevx failed.'
    1143           0 :       deallocate ( work,iwork,mat1,ifail )
    1144           0 :       end subroutine diagonalize_dvs
    1145             : 
    1146           0 :       subroutine diagonalize_dvos(evec,eval,mat,olap,m)
    1147             :       implicit none
    1148             :       real, intent(out) :: eval(:),evec(:,:)
    1149             :       real, intent(in)  :: mat(:,:),olap(:,:)
    1150           0 :       real, allocatable :: work(:),mat1(:,:),olap1(:,:)
    1151             :       real              :: abstol,dlamch
    1152             :       integer, intent(in)  :: m
    1153           0 :       integer, allocatable :: iwork(:),ifail(:)
    1154             :       integer              :: n,ma,idum,info
    1155           0 :       ma = abs(m)
    1156           0 :       n  = size(mat,1)
    1157           0 :       if(n.eq.0) 
    1158           0 :      +  stop 'diagonalize_dvos: zero dimension in eigenvalue problem.'
    1159           0 :       if(ma.gt.n) stop 'diagonalize_dvos: number of selected&
    1160           0 :      & eigenvalues exceeds maximal number.'
    1161           0 :       if(size(mat,2) .ne.n) stop 'diagonalize_dvos: dimensions of&
    1162           0 :      & input matrix differ.'
    1163           0 :       if(size(eval)  .ne.n)
    1164           0 :      +  stop 'diagonalize_dvos: eigenvalue array has wrong size.'
    1165           0 :       if(size(evec,1).ne.n)
    1166           0 :      + stop 'diagonalize_dvos: first dimension of eigenvector is wrong.'
    1167           0 :       if(size(evec,2).ne.n)
    1168           0 :      +stop 'diagonalize_dvos: second dimension of eigenvector is wrong.'
    1169           0 :       if(size(olap,1).ne.n) stop 'diagonalize_dvos: first dimension &
    1170           0 :      &of overlap matrix is wrong.'
    1171           0 :       if(size(olap,2).ne.n) stop 'diagonalize_dvos: second dimension of&
    1172           0 :      & overlap matrix is wrong.'
    1173           0 :       allocate ( work(8*n),iwork(5*n),mat1(n,n),olap1(n,n),ifail(n) ) 
    1174           0 :       mat1 = mat ; olap1 = olap
    1175           0 :       abstol = 2 * dlamch('S')
    1176           0 :       if(m.gt.0) then
    1177             :         call dsygvx(1,'V','I','U',n,mat1,n,olap1,n,0.0,0.0,n-ma+1,n,
    1178           0 :      +              abstol,idum,eval,evec,n,work,8*n,iwork,ifail,info)
    1179             :       else
    1180             :         call dsygvx(1,'V','I','U',n,mat1,n,olap1,n,0.0,0.0,    1,ma,
    1181           0 :      +              abstol,idum,eval,evec,n,work,8*n,iwork,ifail,info)
    1182             :       endif
    1183           0 :       if(info.ne.0) stop 'diagonalize_dvos: dsygvx failed.'
    1184           0 :       deallocate ( work,iwork,mat1,olap1,ifail )
    1185           0 :       end subroutine diagonalize_dvos
    1186             : 
    1187           0 :       subroutine diagonalize_dpvs(evec,eval,mat,m)
    1188             :       implicit none
    1189             :       real, intent(out) :: eval(:),evec(:,:)
    1190             :       real, intent(in)  :: mat(:)
    1191           0 :       real, allocatable :: work(:),mat1(:)
    1192             :       real              :: abstol,dlamch
    1193             :       integer, intent(in)  :: m
    1194           0 :       integer, allocatable :: iwork(:),ifail(:)
    1195             :       integer              :: n,nn,ma,idum,info
    1196           0 :       ma = abs(m)
    1197           0 :       n  = size(eval)
    1198           0 :       nn = n*(n+1)/2
    1199           0 :       if(n.eq.0)
    1200           0 :      +  stop 'diagonalize_dpvs: zero dimension in eigenvalue problem.'
    1201           0 :       if(ma.gt.n) stop 'diagonalize_dpvs: number of selected &
    1202           0 :      &eigenvalues exceeds maximal number.'
    1203           0 :       if(size(mat)   .ne.nn)
    1204           0 :      +  stop 'diagonalize_dpvs: input matrix has wrong size.'
    1205           0 :       if(size(evec,1).ne.n)  stop 'diagonalize_dpvs: first dimension &
    1206           0 :      &of eigenvector is wrong.'
    1207           0 :       if(size(evec,2).ne.n) 
    1208           0 :      +stop 'diagonalize_dpvs: second dimension of eigenvector is wrong.'
    1209           0 :       allocate ( work(8*n),iwork(5*n),mat1(nn),ifail(n) ) ; mat1 = mat
    1210           0 :       abstol = 2 * dlamch('S')
    1211           0 :       if(m.gt.0) then
    1212             :         call dspevx('V','I','U',n,mat1,0.0,0.0,n-ma+1,n,abstol,idum,
    1213           0 :      +              eval,evec,n,work,iwork,ifail,info)
    1214             :       else
    1215             :         call dspevx('V','I','U',n,mat1,0.0,0.0,    1,ma,abstol,idum,
    1216           0 :      +              eval,evec,n,work,iwork,ifail,info)
    1217             :       endif
    1218           0 :       if(info.ne.0) stop 'diagonalize_dpvs: dspevx failed.'
    1219           0 :       deallocate ( work,iwork,mat1,ifail )
    1220           0 :       end subroutine diagonalize_dpvs
    1221             : 
    1222           0 :       subroutine diagonalize_dpvos(evec,eval,mat,olap,m)
    1223             :       implicit none
    1224             :       real, intent(out) :: eval(:),evec(:,:)
    1225             :       real, intent(in)  :: mat(:),olap(:)
    1226           0 :       real, allocatable :: work(:),mat1(:),olap1(:)
    1227             :       real              :: abstol,dlamch
    1228             :       integer, intent(in)  :: m
    1229           0 :       integer, allocatable :: iwork(:),ifail(:)
    1230             :       integer              :: n,nn,ma,idum,info
    1231           0 :       ma = abs(m)
    1232           0 :       n  = size(eval)
    1233           0 :       nn = n*(n+1)/2
    1234           0 :       if(n.eq.0)
    1235           0 :      +  stop 'diagonalize_dpvos: zero dimension in eigenvalue problem.'
    1236           0 :       if(ma.gt.n) stop 'diagonalize_dpvos: number of selected &
    1237           0 :      &eigenvalues exceeds maximal number.'
    1238           0 :       if(size(mat)   .ne.nn)
    1239           0 :      + stop 'diagonalize_dpvos: input matrix has wrong size.'
    1240           0 :       if(size(olap)  .ne.nn)
    1241           0 :      +  stop 'diagonalize_dpvos: overlap matrix has wrong size.'
    1242           0 :       if(size(evec,1).ne.n)  stop 'diagonalize_dpvos: first dimension&
    1243           0 :      & of eigenvector is wrong.'
    1244           0 :       if(size(evec,2).ne.n)  stop 'diagonalize_dpvos: second dimension&
    1245           0 :      & of eigenvector is wrong.'
    1246           0 :       allocate ( work(8*n),iwork(5*n),mat1(nn),olap1(nn),ifail(n) ) 
    1247           0 :       mat1 = mat ; olap1 = olap
    1248           0 :       abstol = 2 * dlamch('S')
    1249           0 :       if(m.gt.0) then
    1250             :         call dspgvx(1,'V','I','U',n,mat1,olap1,0.0,0.0,n-ma+1,n,abstol,
    1251           0 :      +              idum,eval,evec,n,work,iwork,ifail,info)
    1252             :       else
    1253             :         call dspgvx(1,'V','I','U',n,mat1,olap1,0.0,0.0,    1,ma,abstol,
    1254           0 :      +              idum,eval,evec,n,work,iwork,ifail,info)
    1255             :       endif
    1256           0 :       if(info.ne.0) stop 'diagonalize_dpvos: dspgvx failed.'
    1257           0 :       deallocate ( work,iwork,mat1,olap1,ifail )
    1258           0 :       end subroutine diagonalize_dpvos
    1259             : 
    1260           0 :       subroutine diagonalize_zvs(evec,eval,mat,m)
    1261             :       implicit none
    1262             :       real,    intent(out) :: eval(:)
    1263             :       complex, intent(out) :: evec(:,:)
    1264             :       complex, intent(in)  :: mat(:,:)
    1265           0 :       complex, allocatable :: work(:),mat1(:,:)
    1266           0 :       real,    allocatable :: rwork(:)
    1267             :       real                 :: abstol,dlamch
    1268             :       integer,    intent(in)  :: m
    1269           0 :       integer,    allocatable :: iwork(:),ifail(:)
    1270             :       integer                 :: n,ma,idum,info
    1271           0 :       ma = abs(m)
    1272           0 :       n  = size(mat,1)
    1273           0 :       if(n.eq.0)
    1274           0 :      +  stop 'diagonalize_zvs: zero dimension in eigenvalue problem.'
    1275           0 :       if(ma.gt.n) stop 'diagonalize_zvs: number of selected eigenvalues&
    1276           0 :      & exceeds maximal number.'
    1277           0 :       if(size(mat,2) .ne.n)
    1278           0 :      +  stop 'diagonalize_zvs: dimensions of input matrix differ.'
    1279           0 :       if(size(eval)  .ne.n)
    1280           0 :      +  stop 'diagonalize_zvs: eigenvalue array has wrong size.'
    1281           0 :       if(size(evec,1).ne.n)
    1282           0 :      +  stop 'diagonalize_zvs: first dimension of eigenvector is wrong.'
    1283           0 :       if(size(evec,2).ne.n)
    1284           0 :      + stop 'diagonalize_zvs: second dimension of eigenvector is wrong.'
    1285           0 :       allocate ( work(2*n),rwork(7*n),iwork(5*n),mat1(n,n),ifail(n) ) 
    1286           0 :       mat1 = mat
    1287           0 :       abstol = 2 * dlamch('S')
    1288           0 :       if(m.gt.0) then
    1289             :         call zheevx('V','I','U',n,mat1,n,00.,0.0,n-ma+1,n,abstol,idum,
    1290           0 :      +              eval,evec,n,work,2*n,rwork,iwork,ifail,info)
    1291             :       else
    1292             :         call zheevx('V','I','U',n,mat1,n,0.0,0.0,    1,ma,abstol,idum,
    1293           0 :      +              eval,evec,n,work,2*n,rwork,iwork,ifail,info)
    1294             :       endif
    1295           0 :       if(info.ne.0) stop 'diagonalize_zvs: zheevx failed.'
    1296           0 :       deallocate ( work,rwork,iwork,mat1,ifail )
    1297           0 :       end subroutine diagonalize_zvs
    1298             : 
    1299           0 :       subroutine diagonalize_zvos(evec,eval,mat,olap,m)
    1300             :       implicit none
    1301             :       real,    intent(out) :: eval(:)
    1302             :       complex, intent(out) :: evec(:,:)
    1303             :       complex, intent(in)  :: mat(:,:),olap(:,:)
    1304           0 :       complex, allocatable :: work(:),mat1(:,:),olap1(:,:)
    1305           0 :       real,    allocatable :: rwork(:)
    1306             :       real                 :: abstol,dlamch
    1307             :       integer,    intent(in)  :: m
    1308           0 :       integer,    allocatable :: iwork(:),ifail(:)
    1309             :       integer                 :: n,ma,idum,info
    1310           0 :       ma = abs(m)
    1311           0 :       n  = size(mat,1)
    1312           0 :       if(n.eq.0)
    1313           0 :      +  stop 'diagonalize_zvos: zero dimension in eigenvalue problem.'
    1314           0 :       if(ma.gt.n) stop 'diagonalize_zvos: number of selected &
    1315           0 :      &eigenvalues exceeds maximal number.'
    1316           0 :       if(size(mat,2) .ne.n) stop 'diagonalize_zvos: dimensions of input&
    1317           0 :      & matrix differ.'
    1318           0 :       if(size(eval)  .ne.n) stop 'diagonalize_zvos: eigenvalue array&
    1319           0 :      & has wrong size.'
    1320           0 :       if(size(evec,1).ne.n)
    1321           0 :      + stop 'diagonalize_zvos: first dimension of eigenvector is wrong.'
    1322           0 :       if(size(evec,2).ne.n)
    1323           0 :      +stop 'diagonalize_zvos: second dimension of eigenvector is wrong.'
    1324           0 :       if(size(olap,1).ne.n) stop 'diagonalize_zvos: first dimension of&
    1325           0 :      & overlap matrix is wrong.'
    1326           0 :       if(size(olap,2).ne.n) stop 'diagonalize_zvos: second dimension of&
    1327           0 :      & overlap matrix is wrong.'
    1328             :       allocate ( work(2*n),rwork(7*n),iwork(5*n),mat1(n,n),olap1(n,n),
    1329           0 :      +           ifail(n) ) ; mat1 = mat ; olap1 = olap
    1330           0 :       abstol = 2 * dlamch('S')
    1331           0 :       if(m.gt.0) then
    1332             :         call zhegvx(1,'V','I','U',n,mat1,n,olap1,n,0.0,0.0,n-ma+1,n,
    1333           0 :      +          abstol,idum,eval,evec,n,work,2*n,rwork,iwork,ifail,info)
    1334             :       else
    1335             :         call zhegvx(1,'V','I','U',n,mat1,n,olap1,n,0.0,0.0,    1,ma,
    1336           0 :      +          abstol,idum,eval,evec,n,work,2*n,rwork,iwork,ifail,info)
    1337             :       endif
    1338           0 :       if(info.ne.0) stop 'diagonalize_zvos: zhegvx failed.'
    1339           0 :       deallocate ( work,rwork,iwork,mat1,olap1,ifail )
    1340           0 :       end subroutine diagonalize_zvos
    1341             : 
    1342           0 :       subroutine diagonalize_zpvs(evec,eval,mat,m)
    1343             :       implicit none
    1344             :       real,    intent(out) :: eval(:)
    1345             :       complex, intent(out) :: evec(:,:)
    1346             :       complex, intent(in)  :: mat(:)
    1347           0 :       complex, allocatable :: work(:),mat1(:)
    1348           0 :       real,    allocatable :: rwork(:)
    1349             :       real                 :: abstol,dlamch
    1350             :       integer,    intent(in)  :: m
    1351           0 :       integer,    allocatable :: iwork(:),ifail(:)
    1352             :       integer                 :: n,nn,ma,idum,info
    1353           0 :       ma = abs(m)
    1354           0 :       n  = size(eval)
    1355           0 :       nn = n*(n+1)/2
    1356           0 :       if(n.eq.0)
    1357           0 :      +  stop 'diagonalize_zpvs: zero dimension in eigenvalue problem.'
    1358           0 :       if(ma.gt.n)            stop 'diagonalize_zpvs: number of selected&
    1359           0 :      & eigenvalues exceeds maximal number.'
    1360           0 :       if(size(mat)   .ne.nn)
    1361           0 :      +  stop 'diagonalize_zpvs: input matrix has wrong size.'
    1362           0 :       if(size(evec,1).ne.n)
    1363           0 :      + stop 'diagonalize_zpvs: first dimension of eigenvector is wrong.'
    1364           0 :       if(size(evec,2).ne.n)
    1365           0 :      +stop 'diagonalize_zpvs: second dimension of eigenvector is wrong.'
    1366           0 :       allocate ( work(2*n),rwork(7*n),iwork(5*n),mat1(nn),ifail(n) ) 
    1367           0 :       mat1 = mat
    1368           0 :       abstol = 2 * dlamch('S')
    1369           0 :       if(m.gt.0) then
    1370             :         call zhpevx('V','I','U',n,mat1,0.0,0.0,n-ma+1,n,abstol,idum,
    1371           0 :      +              eval,evec,n,work,iwork,ifail,info)
    1372             :       else
    1373             :         call zhpevx('V','I','U',n,mat1,0.0,0.0,    1,ma,abstol,idum,
    1374           0 :      +              eval,evec,n,work,iwork,ifail,info)
    1375             :       endif
    1376           0 :       if(info.ne.0) stop 'diagonalize_zpvs: zhpevx failed.'
    1377           0 :       deallocate ( work,rwork,iwork,mat1,ifail )
    1378           0 :       end subroutine diagonalize_zpvs
    1379             : 
    1380           0 :       subroutine diagonalize_zpvos(evec,eval,mat,olap,m)
    1381             :       implicit none
    1382             :       real,    intent(out) :: eval(:)
    1383             :       complex, intent(out) :: evec(:,:)
    1384             :       complex, intent(in)  :: mat(:),olap(:)
    1385           0 :       complex, allocatable :: work(:),mat1(:),olap1(:)
    1386           0 :       real,    allocatable :: rwork(:)
    1387             :       real                 :: abstol,dlamch
    1388             :       integer,    intent(in)  :: m
    1389           0 :       integer,    allocatable :: iwork(:),ifail(:)
    1390             :       integer                 :: n,nn,ma,idum,info
    1391           0 :       ma = abs(m)
    1392           0 :       n  = size(eval)
    1393           0 :       nn = n*(n+1)/2
    1394           0 :       if(n.eq.0)
    1395           0 :      +  stop 'diagonalize_zpvos: zero dimension in eigenvalue problem.'
    1396           0 :       if(ma.gt.n)            stop 'diagonalize_zpvos: number of&
    1397           0 :      & selected eigenvalues exceeds maximal number.'
    1398           0 :       if(size(mat)   .ne.nn) 
    1399           0 :      + stop 'diagonalize_zpvos: input matrix has wrong size.'
    1400           0 :       if(size(olap)  .ne.nn) stop 'diagonalize_zpvos: overlap matrix&
    1401           0 :      + has wrong size.'
    1402           0 :       if(size(evec,1).ne.n)  stop 'diagonalize_zpvos: first dimension&
    1403           0 :      + of eigenvector is wrong.'
    1404           0 :       if(size(evec,2).ne.n)  stop 'diagonalize_zpvos: second dimension&
    1405           0 :      + of eigenvector is wrong.'
    1406             :       allocate ( work(2*n),rwork(7*n),iwork(5*n),mat1(nn),olap1(nn),
    1407           0 :      +           ifail(n) ) ; mat1 = mat ; olap1 = olap
    1408           0 :       abstol = 2 * dlamch('S')
    1409           0 :       if(m.gt.0) then
    1410             :         call zhpgvx(1,'V','I','U',n,mat1,olap1,0.0,0.0,n-ma+1,n,abstol,
    1411           0 :      +              idum,eval,evec,n,work,rwork,iwork,ifail,info)
    1412             :       else
    1413             :         call zhpgvx(1,'V','I','U',n,mat1,olap1,0.0,0.0,    1,ma,abstol,
    1414           0 :      +              idum,eval,evec,n,work,rwork,iwork,ifail,info)
    1415             :       endif
    1416           0 :       if(info.ne.0) stop 'diagonalize_zpvos: zhpgvx failed.'
    1417           0 :       deallocate ( work,rwork,iwork,mat1,olap1,ifail )
    1418           0 :       end subroutine diagonalize_zpvos
    1419             : 
    1420             :       ! routines for diagonalization: eigenvalue range [r1,r2) or index range [ir1,ir2].
    1421             :       ! the number of actually found eigenvectors is returned in ir2.
    1422             : 
    1423           0 :       subroutine diagonalize_dvx(evec,eval,mat,ir1,ir2,r1,r2)
    1424             :       implicit none
    1425             :       real, intent(out)   :: eval(:),evec(:,:)
    1426             :       real, intent(in)    :: mat(:,:)
    1427           0 :       real, allocatable   :: work(:),mat1(:,:)
    1428             :       real                :: abstol,dlamch
    1429             :       real, intent(in)    :: r1,r2
    1430             :       integer, intent(in)    :: ir1
    1431             :       integer, intent(inout) :: ir2
    1432           0 :       integer, allocatable   :: iwork(:),ifail(:)
    1433             :       integer                :: n,m,idum,info
    1434           0 :       n = size(mat,1)
    1435           0 :       m = ir2 - ir1 + 1
    1436           0 :       if(n.eq.0) 
    1437           0 :      +  stop 'diagonalize_dvx: zero dimension in eigenvalue problem.'
    1438           0 :       if(m.lt.0)
    1439           0 :      +  stop 'diagonalize_dvx: negative index range.'
    1440           0 :       if(m.gt.n)            stop 'diagonalize_dvx: number of selected&
    1441           0 :      & eigenvalues exceeds maximal number.'
    1442           0 :       if(size(mat,2) .ne.n)
    1443           0 :      +  stop 'diagonalize_dvx: dimensions of input matrix differ.'
    1444           0 :       if(size(eval)  .lt.m)
    1445           0 :      +  stop 'diagonalize_dvx: eigenvalue array too small.'
    1446           0 :       if(size(evec,1).ne.n) stop 'diagonalize_dvx: first dimension of&
    1447           0 :      & eigenvector is wrong.'
    1448           0 :       if(size(evec,2).lt.m) stop 'diagonalize_dvx: second dimension of&
    1449           0 :      & eigenvector too small.'
    1450           0 :       allocate ( work(8*n),iwork(5*n),mat1(n,n),ifail(n) ) ; mat1 = mat
    1451           0 :       abstol = 2 * dlamch('S')
    1452           0 :       if(r1.lt.r2) then
    1453             :         call dsyevx('V','V','U',n,mat1,n,r1,r2,0,0,      abstol,idum,
    1454           0 :      +              eval,evec,n,work,8*n,iwork,ifail,info)
    1455             :       else
    1456             :         call dsyevx('V','I','U',n,mat1,n,0.0,0.0,ir1,ir2,abstol,idum,
    1457           0 :      +              eval,evec,n,work,8*n,iwork,ifail,info)
    1458             :       endif
    1459           0 :       ir2 = idum
    1460           0 :       if(info.ne.0) stop 'diagonalize_dvx: dsyevx failed.'
    1461           0 :       deallocate ( work,iwork,mat1,ifail )
    1462           0 :       end subroutine diagonalize_dvx
    1463             : 
    1464           0 :       subroutine diagonalize_dvox(evec,eval,mat,olap,ir1,ir2,r1,r2)
    1465             :       implicit none
    1466             :       real, intent(out)   :: eval(:),evec(:,:)
    1467             :       real, intent(in)    :: mat(:,:),olap(:,:)
    1468           0 :       real, allocatable   :: work(:),mat1(:,:),olap1(:,:)
    1469             :       real                :: abstol,dlamch
    1470             :       real, intent(in)    :: r1,r2
    1471             :       integer, intent(in)    :: ir1
    1472             :       integer, intent(inout) :: ir2
    1473           0 :       integer, allocatable   :: iwork(:),ifail(:)
    1474             :       integer                :: n,m,ma,idum,info
    1475           0 :       n = size(mat,1)
    1476           0 :       m = ir2 - ir1 + 1 
    1477           0 :       if(n.eq.0) 
    1478           0 :      +  stop 'diagonalize_dvox: zero dimension in eigenvalue problem.'
    1479           0 :       if(m.lt.0)
    1480           0 :      +  stop 'diagonalize_dvox: negative index range.'
    1481           0 :       if(m.gt.n) stop 'diagonalize_dvox: number of selected eigenvalues&
    1482           0 :      & exceeds maximal number.'
    1483           0 :       if(size(mat,2) .ne.n) stop 'diagonalize_dvox: dimensions of input&
    1484           0 :      & matrix differ.'
    1485           0 :       if(size(eval)  .lt.m)
    1486           0 :      +  stop 'diagonalize_dvox: eigenvalue array too small.'
    1487           0 :       if(size(evec,1).ne.n)
    1488           0 :      + stop 'diagonalize_dvox: first dimension of eigenvector is wrong.'
    1489           0 :       if(size(evec,2).lt.m) stop 'diagonalize_dvox: second dimension of&
    1490           0 :      & eigenvector too small.'
    1491           0 :       if(size(olap,1).ne.n) stop 'diagonalize_dvox: first dimension of&
    1492           0 :      & overlap matrix is wrong.'
    1493           0 :       if(size(olap,2).ne.n) stop 'diagonalize_dvox: second dimension of&
    1494           0 :      & overlap matrix is wrong.'
    1495           0 :       allocate ( work(8*n),iwork(5*n),mat1(n,n),olap1(n,n),ifail(n) ) 
    1496           0 :       mat1 = mat ; olap1 = olap
    1497           0 :       abstol = 2 * dlamch('S')
    1498           0 :       if(r1.lt.r2) then
    1499             :         call dsygvx(1,'V','V','U',n,mat1,n,olap1,n,r1,r2,0,0,
    1500           0 :      +              abstol,idum,eval,evec,n,work,8*n,iwork,ifail,info)
    1501             :       else
    1502             :         call dsygvx(1,'V','I','U',n,mat1,n,olap1,n,0.0,0.0,ir1,ir2,
    1503           0 :      +              abstol,idum,eval,evec,n,work,8*n,iwork,ifail,info)
    1504             :       endif
    1505           0 :       ir2 = idum
    1506           0 :       if(info.ne.0) stop 'diagonalize_dvos: dsygvx failed.'
    1507           0 :       deallocate ( work,iwork,mat1,olap1,ifail )
    1508           0 :       end subroutine diagonalize_dvox
    1509             : 
    1510           0 :       subroutine diagonalize_dpvx(evec,eval,mat,ir1,ir2,r1,r2)
    1511             :       implicit none
    1512             :       real, intent(out)   :: eval(:),evec(:,:)
    1513             :       real, intent(in)    :: mat(:)
    1514           0 :       real, allocatable   :: work(:),mat1(:)
    1515             :       real                :: abstol,dlamch
    1516             :       real, intent(in)    :: r1,r2
    1517             :       integer, intent(in)    :: ir1
    1518             :       integer, intent(inout) :: ir2
    1519           0 :       integer, allocatable   :: iwork(:),ifail(:)
    1520             :       integer                :: n,m,nn,idum,info
    1521           0 :       n  = size(evec,1)
    1522           0 :       nn = n*(n+1)/2
    1523           0 :       m  = ir2 - ir1 + 1
    1524           0 :       if(n.eq.0)
    1525           0 :      +  stop 'diagonalize_dpvx: zero dimension in eigenvalue problem.'
    1526           0 :       if(m.lt.0)
    1527           0 :      +  stop 'diagonalize_dpvx: negative index range.'
    1528           0 :       if(m.gt.n)             stop 'diagonalize_dpvx: number of selected&
    1529           0 :      & eigenvalues exceeds maximal number.'
    1530           0 :       if(size(mat)   .ne.nn)
    1531           0 :      +  stop 'diagonalize_dpvx: input matrix has wrong size.'
    1532           0 :       if(size(eval)  .lt.m)
    1533           0 :      +  stop 'diagonalize_dpvx: eigenvalue array too small.'
    1534           0 :       if(size(evec,2).lt.m) stop 'diagonalize_dpvx: second dimension of&
    1535           0 :      & eigenvector too small.'
    1536           0 :       allocate ( work(8*n),iwork(5*n),mat1(nn),ifail(n) ) ; mat1 = mat
    1537           0 :       abstol = 2 * dlamch('S')
    1538           0 :       if(r1.lt.r2) then
    1539             :         call dspevx('V','V','U',n,mat1,r1,r2,0,0,      abstol,idum,eval,
    1540           0 :      +              evec,n,work,iwork,ifail,info)
    1541             :       else
    1542             :         call dspevx('V','I','U',n,mat1,0.0,0.0,ir1,ir2,abstol,idum,eval,
    1543           0 :      +              evec,n,work,iwork,ifail,info)
    1544             :       endif
    1545           0 :       ir2 = idum
    1546           0 :       if(info.ne.0) stop 'diagonalize_dpvx: dspevx failed.'
    1547           0 :       deallocate ( work,iwork,mat1,ifail )
    1548           0 :       end subroutine diagonalize_dpvx
    1549             : 
    1550           0 :       subroutine diagonalize_dpvox(evec,eval,mat,olap,ir1,ir2,r1,r2)
    1551             :       implicit none
    1552             :       real, intent(out)   :: eval(:),evec(:,:)
    1553             :       real, intent(in)    :: mat(:),olap(:)
    1554           0 :       real, allocatable   :: work(:),mat1(:),olap1(:)
    1555             :       real                :: abstol,dlamch
    1556             :       real, intent(in)    :: r1,r2
    1557             :       integer, intent(in)    :: ir1
    1558             :       integer, intent(inout) :: ir2
    1559           0 :       integer, allocatable   :: iwork(:),ifail(:)
    1560             :       integer                :: n,nn,m,idum,info
    1561           0 :       n  = size(evec,1)
    1562           0 :       nn = n*(n+1)/2
    1563           0 :       m  = ir2 - ir1 + 1
    1564           0 :       if(n.eq.0)
    1565           0 :      +  stop 'diagonalize_dpvox: zero dimension in eigenvalue problem.'
    1566           0 :       if(m.lt.0)
    1567           0 :      +  stop 'diagonalize_dpvox: negative index range.'
    1568           0 :       if(m.gt.n)            stop 'diagonalize_dpvox: number of selected&
    1569           0 :      & eigenvalues exceeds maximal number.'
    1570           0 :       if(size(mat)   .ne.nn) 
    1571           0 :      +  stop 'diagonalize_dpvox: input matrix has wrong size.'
    1572           0 :       if(size(olap)  .ne.nn) 
    1573           0 :      +  stop 'diagonalize_dpvox: overlap matrix has wrong size.'
    1574           0 :       if(size(eval)  .lt.m)
    1575           0 :      +  stop 'diagonalize_dpvox: eigenvalue array too small.'
    1576           0 :       if(size(evec,2).lt.m)  stop 'diagonalize_dpvox: second dimension&
    1577           0 :      & of eigenvector too small.'
    1578           0 :       allocate ( work(8*n),iwork(5*n),mat1(nn),olap1(nn),ifail(n) ) 
    1579           0 :       mat1 = mat ; olap1 = olap
    1580           0 :       abstol = 2 * dlamch('S')
    1581           0 :       if(r1.lt.r2) then
    1582             :         call dspgvx(1,'V','V','U',n,mat1,olap1,r1,r2,0,0,      abstol,
    1583           0 :      +              idum,eval,evec,n,work,iwork,ifail,info)
    1584             :       else
    1585             :         call dspgvx(1,'V','I','U',n,mat1,olap1,0.0,0.0,ir1,ir2,abstol,
    1586           0 :      +              idum,eval,evec,n,work,iwork,ifail,info)
    1587             :       endif
    1588           0 :       ir2 = idum
    1589           0 :       if(info.ne.0) stop 'diagonalize_dpvox: dspgvx failed.'
    1590           0 :       deallocate ( work,iwork,mat1,olap1,ifail )
    1591           0 :       end subroutine diagonalize_dpvox
    1592             : 
    1593           0 :       subroutine diagonalize_zvx(evec,eval,mat,ir1,ir2,r1,r2)
    1594             :       implicit none
    1595             :       real,    intent(out)   :: eval(:)
    1596             :       complex, intent(out)   :: evec(:,:)
    1597             :       complex, intent(in)    :: mat(:,:)
    1598           0 :       complex, allocatable   :: work(:),mat1(:,:)
    1599           0 :       real,    allocatable   :: rwork(:)
    1600             :       real                   :: abstol,dlamch
    1601             :       real,    intent(in)    :: r1,r2
    1602             :       integer,    intent(in)    :: ir1
    1603             :       integer,    intent(inout) :: ir2
    1604           0 :       integer,    allocatable   :: iwork(:),ifail(:)
    1605             :       integer                   :: n,m,idum,info
    1606           0 :       n = size(mat,1)
    1607           0 :       m = ir2 - ir1 + 1
    1608           0 :       if(n.eq.0)
    1609           0 :      +  stop 'diagonalize_zvx: zero dimension in eigenvalue problem.'
    1610           0 :       if(m.lt.0)           stop 'diagonalize_zvx: negative index range.'
    1611           0 :       if(m.gt.n)           stop 'diagonalize_zvx: number of selected&
    1612           0 :      + eigenvalues exceeds maximal number.'
    1613           0 :       if(size(mat,2) .ne.n)
    1614           0 :      +  stop 'diagonalize_zvx: dimensions of input matrix differ.'
    1615           0 :       if(size(eval)  .lt.m)
    1616           0 :      +  stop 'diagonalize_zvx: eigenvalue array too small.'
    1617           0 :       if(size(evec,1).ne.n)
    1618           0 :      +  stop 'diagonalize_zvx: first dimension of eigenvector is wrong.'
    1619           0 :       if(size(evec,2).lt.m) 
    1620           0 :      +stop 'diagonalize_zvx: second dimension of eigenvector too small.'
    1621           0 :       allocate ( work(2*n),rwork(7*n),iwork(5*n),mat1(n,n),ifail(n) ) 
    1622           0 :       mat1 = mat
    1623           0 :       abstol = 2 * dlamch('S')
    1624           0 :       if(r1.lt.r2) then
    1625             :         call zheevx('V','V','U',n,mat1,n,r1,r2,0,0,      abstol,idum,
    1626           0 :      +              eval,evec,n,work,2*n,rwork,iwork,ifail,info)
    1627             :       else
    1628             :         call zheevx('V','I','U',n,mat1,n,0.0,0.0,ir1,ir2,abstol,idum,
    1629           0 :      +              eval,evec,n,work,2*n,rwork,iwork,ifail,info)
    1630             :       endif
    1631           0 :       ir2 = idum
    1632           0 :       if(info.ne.0) stop 'diagonalize_zvx: zheevx failed.'
    1633           0 :       deallocate ( work,rwork,iwork,mat1,ifail )
    1634           0 :       end subroutine diagonalize_zvx
    1635             : 
    1636           0 :       subroutine diagonalize_zvox(evec,eval,mat,olap,ir1,ir2,r1,r2)
    1637             :       implicit none
    1638             :       real,    intent(out)   :: eval(:)
    1639             :       complex, intent(out)   :: evec(:,:)
    1640             :       complex, intent(in)    :: mat(:,:),olap(:,:)
    1641           0 :       complex, allocatable   :: work(:),mat1(:,:),olap1(:,:)
    1642           0 :       real,    allocatable   :: rwork(:)
    1643             :       real                   :: abstol,dlamch
    1644             :       real,    intent(in)    :: r1,r2
    1645             :       integer,    intent(in)    :: ir1
    1646             :       integer,    intent(inout) :: ir2
    1647           0 :       integer,    allocatable   :: iwork(:),ifail(:)
    1648             :       integer                   :: n,m,idum,info
    1649           0 :       n = size(mat,1)
    1650           0 :       m = ir2 - ir1 + 1
    1651           0 :       if(n.eq.0)
    1652           0 :      +  stop 'diagonalize_zvox: zero dimension in eigenvalue problem.'
    1653           0 :       if(m.lt.0)
    1654           0 :      +  stop 'diagonalize_zvox: negative index range.'
    1655           0 :       if(m.gt.n)            stop 'diagonalize_zvox: number of selected&
    1656           0 :      & eigenvalues exceeds maximal number.'
    1657           0 :       if(size(mat,2) .ne.n)
    1658           0 :      +  stop 'diagonalize_zvox: dimensions of input matrix differ.'
    1659           0 :       if(size(eval)  .lt.m) stop 'diagonalize_zvox: eigenvalue array&
    1660           0 :      & too small.'
    1661           0 :       if(size(evec,1).ne.n) stop 'diagonalize_zvox: first dimension of&
    1662           0 :      & eigenvector is wrong.'
    1663           0 :       if(size(evec,2).lt.m) stop 'diagonalize_zvox: second dimension of&
    1664           0 :      & eigenvector too small.'
    1665           0 :       if(size(olap,1).ne.n) stop 'diagonalize_zvox: first dimension of&
    1666           0 :      & overlap matrix is wrong.'
    1667           0 :       if(size(olap,2).ne.n) stop 'diagonalize_zvox: second dimension of&
    1668           0 :      & overlap matrix is wrong.'
    1669             :       allocate ( work(2*n),rwork(7*n),iwork(5*n),mat1(n,n),olap1(n,n),
    1670           0 :      +           ifail(n) ) ; mat1 = mat ; olap1 = olap
    1671           0 :       abstol = 2 * dlamch('S')
    1672           0 :       if(r1.lt.r2) then
    1673             :         call zhegvx(1,'V','V','U',n,mat1,n,olap1,n,r1,r2,0,0,
    1674           0 :      +          abstol,idum,eval,evec,n,work,2*n,rwork,iwork,ifail,info)
    1675             :       else
    1676             :         call zhegvx(1,'V','I','U',n,mat1,n,olap1,n,0.0,0.0,ir1,ir2,
    1677           0 :      +          abstol,idum,eval,evec,n,work,2*n,rwork,iwork,ifail,info)
    1678             :       endif
    1679           0 :       ir2 = idum
    1680           0 :       if(info.ne.0) stop 'diagonalize_zvox: zhegvx failed.'
    1681           0 :       deallocate ( work,rwork,iwork,mat1,olap1,ifail )
    1682           0 :       end subroutine diagonalize_zvox
    1683             : 
    1684           0 :       subroutine diagonalize_zpvx(evec,eval,mat,ir1,ir2,r1,r2)
    1685             :       implicit none
    1686             :       real,    intent(out)   :: eval(:)
    1687             :       complex, intent(out)   :: evec(:,:)
    1688             :       complex, intent(in)    :: mat(:)
    1689           0 :       complex, allocatable   :: work(:),mat1(:)
    1690           0 :       real,    allocatable   :: rwork(:)
    1691             :       real                   :: abstol,dlamch
    1692             :       real,    intent(in)    :: r1,r2
    1693             :       integer,    intent(in)    :: ir1
    1694             :       integer,    intent(inout) :: ir2
    1695           0 :       integer,    allocatable   :: iwork(:),ifail(:)
    1696             :       integer                   :: n,nn,m,idum,info
    1697           0 :       n  = size(evec,1)
    1698           0 :       nn = n*(n+1)/2
    1699           0 :       m  = ir2 - ir1 + 1
    1700           0 :       if(n.eq.0)
    1701           0 :      +  stop 'diagonalize_zpvx: zero dimension in eigenvalue problem.'
    1702           0 :       if(m.lt.0)
    1703           0 :      +  stop 'diagonalize_zpvx: negative index range.'
    1704           0 :       if(m.gt.n)             stop 'diagonalize_zpvx: number of selected&
    1705           0 :      & eigenvalues exceeds maximal number.'
    1706           0 :       if(size(mat)   .ne.nn)
    1707           0 :      +  stop 'diagonalize_zpvx: input matrix has wrong size.'
    1708           0 :       if(size(eval)  .lt.m)
    1709           0 :      +  stop 'diagonalize_zpvx: eigenvalue array too small.'
    1710           0 :       if(size(evec,2).lt.m)  stop 'diagonalize_zpvx: second dimension&
    1711           0 :      +  of eigenvector too small.'
    1712           0 :       allocate ( work(2*n),rwork(7*n),iwork(5*n),mat1(nn),ifail(n) ) 
    1713           0 :       mat1 = mat
    1714           0 :       abstol = 2 * dlamch('S')
    1715           0 :       if(r1.lt.r2) then
    1716             :         call zhpevx('V','V','U',n,mat1,r1,r2,0,0,      abstol,idum,
    1717           0 :      +              eval,evec,n,work,iwork,ifail,info)
    1718             :       else
    1719             :         call zhpevx('V','I','U',n,mat1,0.0,0.0,ir1,ir2,abstol,idum,
    1720           0 :      +              eval,evec,n,work,iwork,ifail,info)
    1721             :       endif
    1722           0 :       ir2 = idum
    1723           0 :       if(info.ne.0) stop 'diagonalize_zpvx: zhpevx failed.'
    1724           0 :       deallocate ( work,rwork,iwork,mat1,ifail )
    1725           0 :       end subroutine diagonalize_zpvx
    1726             : 
    1727           0 :       subroutine diagonalize_zpvox(evec,eval,mat,olap,ir1,ir2,r1,r2)
    1728             :       implicit none
    1729             :       real,    intent(out)   :: eval(:)
    1730             :       complex, intent(out)   :: evec(:,:)
    1731             :       complex, intent(in)    :: mat(:),olap(:)
    1732           0 :       complex, allocatable   :: work(:),mat1(:),olap1(:)
    1733           0 :       real,    allocatable   :: rwork(:)
    1734             :       real                   :: abstol,dlamch
    1735             :       real,    intent(in)    :: r1,r2
    1736             :       integer,    intent(in)    :: ir1
    1737             :       integer,    intent(inout) :: ir2
    1738           0 :       integer,    allocatable   :: iwork(:),ifail(:)
    1739             :       integer                   :: n,nn,m,idum,info
    1740           0 :       n  = size(evec,1)
    1741           0 :       nn = n*(n+1)/2
    1742           0 :       m  = ir2 - ir1 + 1
    1743           0 :       if(n.eq.0)
    1744           0 :      +  stop 'diagonalize_zpvox: zero dimension in eigenvalue problem.'
    1745           0 :       if(m.gt.n)            stop 'diagonalize_zpvox: number of selected&
    1746           0 :      & eigenvalues exceeds maximal number.'
    1747           0 :       if(m.lt.0)
    1748           0 :      +  stop 'diagonalize_zpvox: negative index range.'
    1749           0 :       if(size(mat)   .ne.nn)
    1750           0 :      +  stop 'diagonalize_zpvox: input matrix has wrong size.'
    1751           0 :       if(size(olap)  .ne.nn)
    1752           0 :      +  stop 'diagonalize_zpvox: overlap matrix has wrong size.'
    1753           0 :       if(size(eval)  .lt.m)  stop 'diagonalize_zpvox: first dimension&
    1754           0 :      & of eigenvector too small.'
    1755           0 :       if(size(evec,2).lt.m)  stop 'diagonalize_zpvox: second dimension&
    1756           0 :      & of eigenvector too small.'
    1757             :       allocate ( work(2*n),rwork(7*n),iwork(5*n),mat1(nn),olap1(nn),
    1758           0 :      +           ifail(n) ) ; mat1 = mat ; olap1 = olap
    1759           0 :       abstol = 2 * dlamch('S')
    1760           0 :       if(r1.lt.r2) then
    1761             :         call zhpgvx(1,'V','V','U',n,mat1,olap1,r1,r2,0,0,      abstol,
    1762           0 :      +              idum,eval,evec,n,work,rwork,iwork,ifail,info)
    1763             :       else
    1764             :         call zhpgvx(1,'V','I','U',n,mat1,olap1,0.0,0.0,ir1,ir2,abstol,
    1765           0 :      +              idum,eval,evec,n,work,rwork,iwork,ifail,info)
    1766             :       endif
    1767           0 :       ir2 = idum
    1768           0 :       if(info.ne.0) stop 'diagonalize_zpvox: zhpgvx failed.'
    1769           0 :       deallocate ( work,rwork,iwork,mat1,olap1,ifail )
    1770           0 :       end subroutine diagonalize_zpvox
    1771             : 
    1772             : c     --------
    1773             : 
    1774           0 :       subroutine geteigen_zpvo(evec,eval,mat,olap)
    1775             :       implicit none
    1776             :       real,    intent(out) :: eval
    1777             :       complex, intent(out) :: evec(:)
    1778             :       complex, intent(in)  :: mat(:),olap(:)
    1779           0 :       complex, allocatable :: mat1(:),olap1(:),work(:),evec1(:,:)
    1780           0 :       real,    allocatable :: eval1(:),rwork(:)
    1781           0 :       integer,    allocatable :: iwork(:),ifail(:)
    1782             :       integer                 :: n,nn,info
    1783             :       real                 :: dlamch
    1784           0 :       stop 'geteigen: disabled!'
    1785             :       n = size(evec) ; nn = n*(n+1)/2
    1786             :       if(size(mat) .ne.nn) 
    1787             :      +  stop 'diagonalize_zpvo: input matrix has wrong size.'
    1788             :       if(size(olap).ne.nn)
    1789             :      +  stop 'diagonalize_zpvo: overlap matrix has wrong size.'
    1790             :       allocate ( mat1(nn),olap1(nn),eval1(n),evec1(n,n),work(2*n),
    1791             :      +           rwork(7*n),iwork(5*n),ifail(n) ) 
    1792             :       mat1 = mat ; olap1 = olap
    1793             : c      call zhpgvx(1,'V','I','U',n,mat1,olap1,0.0,0.0,n,n,2*dlamch('S'),1,eval1,evec1,n,work,rwork,iwork,ifail,info)
    1794             :       if(info.ne.0) stop 'geteigen_zpvo: zhpgvx failed.'
    1795             :       evec = evec1(:,1)
    1796             :       eval = eval1(1)
    1797             :       deallocate ( mat1,olap1,eval1,evec1,work,rwork,iwork,ifail )
    1798             :       end subroutine geteigen_zpvo
    1799             : 
    1800             : c     --------
    1801             : 
    1802           0 :       subroutine inverse_d(mati,mat)
    1803             :       implicit none
    1804             :       real, intent(out) :: mati(:,:)
    1805             :       real, intent(in)  :: mat(:,:)
    1806             :       integer              :: n,info,i,j
    1807           0 :       n = size(mat,1)
    1808           0 :       if(size(mat,2) .ne.n)
    1809           0 :      +  stop 'inverse_d: dimensions of input array differ.'
    1810           0 :       if(size(mati,1).ne.n)
    1811           0 :      +  stop 'inverse_d: output array has wrong dimensions.'
    1812           0 :       if(size(mati,2).ne.n)
    1813           0 :      +  stop 'inverse_d: dimensions of output array differ.'
    1814           0 :       mati = mat
    1815           0 :       call dpotrf('U',n,mati,n,info) 
    1816           0 :       if(info.ne.0) stop 'inverse_d: dpotrf failed.'
    1817           0 :       call dpotri('U',n,mati,n,info) 
    1818           0 :       if(info.ne.0) stop 'inverse_d: dpotri failed.'
    1819           0 :       do i = 1,n ; do j = 1,i ; mati(i,j) = mati(j,i) ; enddo ; enddo
    1820           0 :       end subroutine inverse_d
    1821             : 
    1822           0 :       subroutine inverse_dp(mati,mat)
    1823             :       implicit none
    1824             :       real, intent(out) :: mati(:)
    1825             :       real, intent(in)  :: mat(:)
    1826             :       integer              :: n,nn,info
    1827           0 :       nn = size(mat,1) ; n = nint(sqrt(0.25+2*nn)-0.5)
    1828           0 :       if(size(mati).ne.nn)
    1829           0 :      +  stop 'inverse_dp: output array has wrong size.'
    1830           0 :       mati = mat
    1831           0 :       call dpptrf('U',n,mati,info) 
    1832           0 :       if(info.ne.0) stop 'inverse_dp: dpptrf failed.'
    1833           0 :       call dpptri('U',n,mati,info) 
    1834           0 :       if(info.ne.0) stop 'inverse_dp: dpptri failed.'
    1835           0 :       end subroutine inverse_dp
    1836             : 
    1837           0 :       subroutine inverse_z(mati,mat)
    1838             :       implicit none
    1839             :       complex, intent(out) :: mati(:,:)
    1840             :       complex, intent(in)  :: mat(:,:)
    1841             :       integer                 :: n,info,i,j
    1842           0 :       n = size(mat,1)
    1843           0 :       if(size(mat,2) .ne.n)
    1844           0 :      +  stop 'inverse_z: dimensions of input array differ.'
    1845           0 :       if(size(mati,1).ne.n)
    1846           0 :      +  stop 'inverse_z: output array has wrong dimensions.'
    1847           0 :       if(size(mati,2).ne.n)
    1848           0 :      +  stop 'inverse_z: dimensions of output array differ.'
    1849           0 :       mati = mat
    1850           0 :       call zpotrf('U',n,mati,n,info) 
    1851           0 :       if(info.ne.0) then 
    1852           0 :         WRITE(*,*) 'info',info
    1853           0 :         stop 'inverse_z: zpotrf failed.' 
    1854             :       endif
    1855           0 :       call zpotri('U',n,mati,n,info) 
    1856           0 :       if(info.ne.0) stop 'inverse_z: zpotri failed.'
    1857           0 :       do i = 1,n 
    1858           0 :         do j = 1,i 
    1859           0 :           mati(i,j) = conjg(mati(j,i)) 
    1860             :         enddo 
    1861             :       enddo
    1862           0 :       end subroutine inverse_z
    1863             : 
    1864           0 :       subroutine inverse_zp(mati,mat)
    1865             :       implicit none
    1866             :       complex, intent(out) :: mati(:)
    1867             :       complex, intent(in)  :: mat(:)
    1868             :       integer                 :: n,nn,info
    1869           0 :       nn = size(mat,1) ; n = nint(sqrt(0.25+2*nn)-0.5)
    1870           0 :       if(size(mati).ne.nn)
    1871           0 :      +  stop 'inverse_zp: output array has wrong size.'
    1872           0 :       mati = mat
    1873           0 :       call zpptrf('U',n,mati,n,info) 
    1874           0 :       if(info.ne.0) stop 'inverse_zp: zpptrf failed.'
    1875           0 :       call zpptri('U',n,mati,n,info) 
    1876           0 :       if(info.ne.0) stop 'inverse_zp: zpptri failed.'
    1877           0 :       end subroutine inverse_zp
    1878             : 
    1879           0 :       subroutine inverse_d1(mat)
    1880             :       implicit none
    1881             :       real, intent(inout) :: mat(:,:)
    1882           0 :       real, allocatable   :: work(:)
    1883           0 :       integer, allocatable   :: ipiv(:)
    1884             :       integer                :: n,info,i,j
    1885           0 :       n = size(mat,1)
    1886           0 :       if(size(mat,2).ne.n) stop 'inverse_d1: array dimensions differ.'
    1887           0 :       allocate ( ipiv(n),work(n) )
    1888           0 :       call dgetrf(n,n,mat,n,ipiv,info)
    1889           0 :       if(info.ne.0) stop 'inverse_d1: dpotrf failed.'
    1890           0 :       call dgetri(n,mat,n,ipiv,work,n,info)
    1891           0 :       if(info.ne.0) stop 'inverse_d1: dpotri failed.'
    1892           0 :       end subroutine inverse_d1
    1893             : 
    1894           0 :       subroutine inverse_dp1(mat)
    1895             :       implicit none
    1896             :       real, intent(inout) :: mat(:)
    1897             :       integer                :: n,nn,info
    1898           0 :       nn = size(mat,1) ; n = nint(sqrt(0.25+2*nn)-0.5)
    1899           0 :       call dpptrf('U',n,mat,info)
    1900           0 :       if(info.ne.0) stop 'inverse_dp1: dpptrf failed.'
    1901           0 :       call dpptri('U',n,mat,info) 
    1902           0 :       if(info.ne.0) stop 'inverse_dp1: dpptri failed.'
    1903           0 :       end subroutine inverse_dp1
    1904             : 
    1905           0 :       subroutine inverse_z1(mat)
    1906             :       implicit none
    1907             :       complex, intent(inout) :: mat(:,:)
    1908           0 :       complex, allocatable   :: work(:)
    1909           0 :       integer,    allocatable   :: ipiv(:)
    1910             :       integer                   :: n,info
    1911           0 :       n = size(mat,1)
    1912           0 :       if(size(mat,2).ne.n) stop 'inverse_z1: array dimensions differ.'
    1913           0 :       allocate(ipiv(n),work(n))
    1914           0 :       call zgetrf(n,n,mat,n,ipiv,info)
    1915           0 :       if(info.ne.0) stop 'inverse_z1: zgetrf failed.'
    1916           0 :       call zgetri(n,mat,n,ipiv,work,n,info)
    1917           0 :       if(info.ne.0) stop 'inverse_z1: zgetri failed.'
    1918           0 :       end subroutine inverse_z1
    1919             : 
    1920           0 :       subroutine inverse_zp1(mat)
    1921             :       implicit none
    1922             :       complex, intent(inout) :: mat(:)
    1923           0 :       complex, allocatable   :: work(:)
    1924           0 :       integer,    allocatable   :: ipiv(:)
    1925             :       integer                   :: n,nn,info
    1926           0 :       nn = size(mat,1) ; n = nint(sqrt(0.25+2*nn)-0.5) 
    1927           0 :       allocate(ipiv(n),work(n))
    1928           0 :       call zsptrf('U',n,mat,ipiv,info)
    1929           0 :       if(info.ne.0) stop 'inverse_zp1: zpptrf failed.'
    1930           0 :       call zsptri('U',n,mat,ipiv,work,info)
    1931           0 :       if(info.ne.0) stop 'inverse_zp1: zpptri failed.'
    1932           0 :       end subroutine inverse_zp1
    1933             : 
    1934             : c     --------
    1935             : 
    1936           0 :       subroutine sqrtmat_d(matout,matin)
    1937             :       implicit none
    1938             :       real, intent(out) :: matout(:,:)
    1939             :       real, intent(in)  :: matin(:,:)
    1940           0 :       real, allocatable :: eval(:),evec(:,:)
    1941             :       integer              :: n,i
    1942           0 :       n = size(matin,1)
    1943           0 :       if(size(matin,2).ne.n) 
    1944           0 :      +  stop 'sqrtmat_d: dimensions of input array differ.'
    1945           0 :       allocate ( evec(n,n),eval(n) )
    1946           0 :       call diagonalize(evec,eval,matin) 
    1947           0 :       if(any(eval.lt.0.0)) stop 'sqrtmat_d: negative eigenvalue.'
    1948           0 :       do i = 1,n
    1949           0 :         evec(:,i) = sqrt(sqrt(eval(i))) * evec(:,i)
    1950             :       enddo
    1951           0 :       matout = matmul(evec,transpose(evec))
    1952           0 :       deallocate ( evec,eval )
    1953           0 :       end subroutine sqrtmat_d
    1954             : 
    1955           0 :       subroutine sqrtmat_dp(matout,matin)
    1956             :       implicit none
    1957             :       real, intent(out) :: matout(:)
    1958             :       real, intent(in)  :: matin(:)
    1959           0 :       real, allocatable :: eval(:),evec(:,:)
    1960             :       integer              :: nn,n,i,j
    1961           0 :       nn = size(matin,1) ; n = nint(sqrt(0.25+2*nn)-0.5)
    1962           0 :       if(size(matout).ne.nn)
    1963           0 :      +  stop 'sqrtmat_dp: output array has wrong size.'
    1964           0 :       allocate ( evec(n,n),eval(n) )
    1965           0 :       call diagonalize(evec,eval,matin) 
    1966           0 :       if(any(eval.lt.0.0)) stop 'sqrtmat_dp: negative eigenvalue.'
    1967           0 :       do i = 1,n
    1968           0 :         evec(:,i) = sqrt(sqrt(eval(i))) * evec(:,i)
    1969             :       enddo
    1970           0 :       matout = packmat ( matmul(evec,transpose(evec)) )
    1971           0 :       deallocate ( evec,eval )
    1972           0 :       end subroutine sqrtmat_dp
    1973             : 
    1974           0 :       subroutine sqrtmat_z(matout,matin)
    1975             :       implicit none
    1976             :       complex, intent(out) :: matout(:,:)
    1977             :       complex, intent(in)  :: matin(:,:)
    1978           0 :       complex, allocatable :: evec(:,:)
    1979           0 :       real,    allocatable :: eval(:)
    1980             :       integer                 :: n,i
    1981           0 :       n = size(matin,1)
    1982           0 :       if(size(matin,2).ne.n) 
    1983           0 :      +  stop 'sqrtmat_z: dimensions of input array differ.'
    1984           0 :       allocate ( evec(n,n),eval(n) )
    1985           0 :       call diagonalize(evec,eval,matin) 
    1986           0 :       if(any(eval.lt.0.0)) stop 'sqrtmat_z: negative eigenvalue.'
    1987           0 :       do i = 1,n
    1988           0 :         evec(:,i) = sqrt(sqrt(eval(i))) * evec(:,i)
    1989             :       enddo
    1990           0 :       matout = matmul(evec,conjg(transpose(evec)))
    1991           0 :       deallocate ( evec,eval )
    1992           0 :       end subroutine sqrtmat_z
    1993             : 
    1994           0 :       subroutine sqrtmat_zp(matout,matin)
    1995             :       implicit none
    1996             :       complex, intent(out) :: matout(:)
    1997             :       complex, intent(in)  :: matin(:)
    1998           0 :       complex, allocatable :: evec(:,:)
    1999           0 :       real,    allocatable :: eval(:)
    2000             :       integer                 :: nn,n,i,j
    2001           0 :       nn = size(matin,1) ; n = nint(sqrt(0.25+2*nn)-0.5)
    2002           0 :       if(size(matout).ne.nn)
    2003           0 :      +  stop 'sqrtmat_zp: output array has wrong size.'
    2004           0 :       allocate ( eval(n),evec(n,n) )
    2005           0 :       call diagonalize(evec,eval,matin) 
    2006           0 :       if(any(eval.lt.0.0)) stop 'sqrtmat_zp: negative eigenvalue.'
    2007           0 :       do i = 1,n
    2008           0 :         evec(:,i) = sqrt(sqrt(eval(i))) * evec(:,i) !sqrt(sqrt(eval(i))) * evec(:,i)
    2009             :       enddo
    2010           0 :       matout = packmat(matmul(evec,conjg(transpose(evec))))
    2011           0 :       deallocate ( evec,eval )
    2012           0 :       end subroutine sqrtmat_zp
    2013             : 
    2014           0 :       subroutine sqrtmat_d1(mat)
    2015             :       implicit none
    2016             :       real, intent(inout) :: mat(:,:)
    2017           0 :       real, allocatable   :: eval(:),evec(:,:)
    2018             :       integer                :: n,i
    2019           0 :       n = size(mat,1)
    2020           0 :       if(size(mat,2).ne.n) stop 'sqrtmat_d1: array dimensions differ.'
    2021           0 :       allocate ( evec(n,n),eval(n) )
    2022           0 :       call diagonalize(evec,eval,mat) 
    2023           0 :       if(any(eval.lt.0.0)) stop 'sqrtmat_d1: negative eigenvalue.'
    2024           0 :       do i = 1,n
    2025           0 :         evec(:,i) = sqrt(sqrt(eval(i))) * evec(:,i)
    2026             :       enddo
    2027           0 :       mat = matmul(evec,transpose(evec))
    2028           0 :       deallocate ( evec,eval )
    2029           0 :       end subroutine sqrtmat_d1
    2030             : 
    2031           0 :       subroutine sqrtmat_dp1(mat)
    2032             :       implicit none
    2033             :       real, intent(inout) :: mat(:)
    2034           0 :       real, allocatable   :: eval(:),evec(:,:)
    2035             :       integer                :: nn,n,i,j
    2036           0 :       nn = size(mat,1) ; n = nint(sqrt(0.25+2*nn)-0.5)
    2037           0 :       allocate ( evec(n,n),eval(n) )
    2038           0 :       call diagonalize(evec,eval,mat) 
    2039           0 :       if(any(eval.lt.0.0)) stop 'sqrtmat_dp1: negative eigenvalue.'
    2040           0 :       do i = 1,n
    2041           0 :         evec(:,i) = sqrt(sqrt(eval(i))) * evec(:,i)
    2042             :       enddo
    2043           0 :       mat = packmat ( matmul(evec,transpose(evec)) )
    2044           0 :       deallocate ( eval,evec )
    2045           0 :       end subroutine sqrtmat_dp1
    2046             : 
    2047           0 :       subroutine sqrtmat_z1(mat)
    2048             :       implicit none
    2049             :       complex, intent(inout) :: mat(:,:)
    2050           0 :       complex, allocatable   :: evec(:,:)
    2051           0 :       real,    allocatable   :: eval(:)
    2052             :       integer                   :: n,i
    2053           0 :       n = size(mat,1) 
    2054           0 :       if(size(mat,2).ne.n) stop 'sqrtmat_z1: array dimensions differ.'
    2055           0 :       allocate ( evec(n,n),eval(n) )
    2056           0 :       call diagonalize(evec,eval,mat) 
    2057           0 :       if(any(eval.lt.0.0)) stop 'sqrtmat_z1: negative eigenvalue.'
    2058           0 :       do i = 1,n
    2059           0 :         evec(:,i) = sqrt(sqrt(eval(i))) * evec(:,i)
    2060             :       enddo
    2061           0 :       mat = matmul(evec,conjg(transpose(evec)))
    2062           0 :       deallocate ( eval,evec )
    2063           0 :       end subroutine sqrtmat_z1
    2064             : 
    2065           0 :       subroutine sqrtmat_zp1(mat)
    2066             :       implicit none
    2067             :       complex, intent(inout) :: mat(:)
    2068           0 :       complex, allocatable   :: evec(:,:)
    2069             :       real,    allocatable   :: eval(:)
    2070             :       integer                   :: nn,n,i,j
    2071           0 :       nn = size(mat,1) ; n = nint(sqrt(0.25+2*nn)-0.5) 
    2072           0 :       allocate ( eval(n),evec(n,n) )
    2073           0 :       call diagonalize(evec,eval,mat) 
    2074           0 :       if(any(eval.lt.0.0)) stop 'sqrtmat_zp1: negative eigenvalue.'
    2075           0 :       do i = 1,n
    2076           0 :         evec(:,i) = sqrt(sqrt(eval(i))) * evec(:,i)
    2077             :       enddo
    2078           0 :       mat = packmat(matmat(evec,conjg(transpose(evec))))
    2079           0 :       deallocate ( evec,eval )
    2080           0 :       end subroutine sqrtmat_zp1
    2081             :       
    2082             : c     --------
    2083             : 
    2084             :       end module m_wrapper
    2085             :       

Generated by: LCOV version 1.13