Chimie Théorique » Opt'n Path

!----------------------------------------------------------------------

!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon,

!  Centre National de la Recherche Scientifique,

!  Université Claude Bernard Lyon 1. All rights reserved.

!

!  This work is registered with the Agency for the Protection of Programs

!  as IDDN.FR.001.100009.000.S.P.2014.000.30625

!

!  Authors: P. Fleurat-Lessard, P. Dayal

!  Contact: optnpath@gmail.com

!

! This file is part of "Opt'n Path".

!

!  "Opt'n Path" is free software: you can redistribute it and/or modify

!  it under the terms of the GNU Affero General Public License as

!  published by the Free Software Foundation, either version 3 of the License,

!  or (at your option) any later version.

!

!  "Opt'n Path" is distributed in the hope that it will be useful,

!  but WITHOUT ANY WARRANTY; without even the implied warranty of

!

!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

!  GNU Affero General Public License for more details.

!

!  You should have received a copy of the GNU Affero General Public License

!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.

!

! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,

! for commercial licensing opportunities.

!----------------------------------------------------------------------

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

SUBROUTINE GenInv(N,A,InvA,NReal)

!!!!!!!!!!!!!!!!

  !

  !     This subroutines calculates the generalized inverse of a matrix

  !     It first diagonalize the matrix A, then inverse all non-zero

  !     eigenvalues, and forms the  InvA matrix using these new eigenvalues

  !

  !     Input:

  !     N : dimension of A

  !     NReal :Actual dimension of A

  !     A(N,N) : Initial Matrix, stored in A(Nreal,Nreal)

  !

  !     Output:

  !     InvA(N,N) : Inversed Matrix, stored in a (Nreal,NReal) matrix

  !

!!!!!!!!!!!!!!!!!!!!!!!

  Use Vartypes

  IMPLICIT NONE

  INTEGER(KINT), INTENT(IN) :: N,Nreal

  REAL(KREAL), INTENT(IN) :: A(NReal,NReal)

  REAL(KREAL), INTENT(OUT) :: InvA(NReal,NReal)

  !

  INTEGER(KINT) :: I,J,K

  REAL(KREAL), ALLOCATABLE :: EigVec(:,:) ! (Nreal,Nreal)

  REAL(KREAL), ALLOCATABLE :: EigVal(:) ! (Nreal)

  REAL(KREAL), ALLOCATABLE :: ATmp(:,:) ! (NReal,Nreal)

  REAL(KREAL) :: ss

  !

  REAL(KREAL), PARAMETER :: eps=1e-12

  ALLOCATE(EigVec(Nreal,Nreal), EigVal(Nreal),ATmp(NReal,NReal))

! A will be destroyed in Jacobi so we save it

  ATmp=A

  CALL JAcobi(ATmp,N,EigVal,EigVec,NReal)

  DO I=1,N

     IF (abs(EigVal(I)).GT.eps) EigVal(I)=1.d0/EigVal(I)

  END DO

  InvA=0.d0

  do k = 1, n

     do j = 1, n

        ss = eigval(k) * eigvec(j,k)

        do i = 1, j

            InvA(i,j) = InvA(i,j) + ss * eigvec(i,k)

         end do

      end do

   end do

   do j = 1, n

      do i = 1, j-1

         InvA(j,i) = InvA(i,j)

      end do

   end do

  DEALLOCATE(EigVec, EigVal,ATmp)

END SUBROUTINE GenInv

!============================================================

!

!     ++  Matrix diagonalization Using jacobi

!  Works only for symmetric matrices

!     EigvenVectors  : V(i,i)

!     Eigenvalues : D(i)

! PFL 30/05/03

! This version uses packed matrices.

! it unpacks them before calling Jacobi !

! we have  AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;

!

!============================================================

!

      SUBROUTINE JacPacked(N,AP,D,V,nreal)

      Use Vartypes

      IMPLICIT NONE

      INTEGER(KINT), INTENT(IN) :: N,NREAL

      REAL(KREAL) :: AP(N*(N+1)/2)

      REAL(KREAL), ALLOCATABLE ::  A(:,:)

      REAL(KREAL) :: V(Nreal,Nreal),D(Nreal)

      INTEGER(KINT) :: i,j,nn

      allocate(A(nreal,nreal))

      nn=n*(n+1)/2

!      WRITE(*,*) 'Jacpa 0'

!      WRITE(*,'(I3,10(1X,F15.6))') n,(AP(i),i=1,min(nn,5))

 !     WRITE(*,*) 'Jacpa 0'

      do j=1,n

         do i=1,j

!            WRITE(*,*) i,j

            A(i,J)=AP(i + (j-1)*j/2)

            A(j,I)=A(i,J)

         end do

      end do

!      do j=1,n

!         WRITE(*,'(10(1X,F15.6))') (A(i,J),i=1,min(n,5))

!      end do

!      WRITE(*,*) 'Jacpa 1'

      call Jacobi(A,n,D,V,Nreal)

!      WRITE(*,*) 'Jacpa 2'

!      DO i=1,n

!         WRITE(*,'(1X,I3,10(1X,F15.6))') i,D(i),(V(j,i),j=1,min(n,5))

!      end do

      deallocate(a)

    end SUBROUTINE JacPacked

!

!============================================================

!

!     ++  Matrix diagonalization Using jacobi

!  Works only for symmetric matrices

!     EigvenVectors  : V

!     Eigenvalues : D

!

!============================================================

!

SUBROUTINE JACOBI(A,N,D,V,Nreal)

!!!!!!!!!!!!!!!!

  !

  !     Input:

  !     N      :  Dimension of A

  !     NReal  : Actual dimensions of A, D and V.

  !

  !     Input/output:

  !     A(N,N) : Matrix to be diagonalized, store in a (Nreal,Nreal) matrix

  !              Destroyed in output.

  !     Output:

  !     V(N,N) : Eigenvectors, stored in V(NReal, NReal)

  !     D(N)   : Eigenvalues, stored in D(NReal)

  !

  Use Vartypes

  IMPLICIT NONE

  INTEGER(KINT), parameter :: max_it=500

  REAL(KREAL), ALLOCATABLE ::  B(:),Z(:)

  INTEGER(KINT) :: N,NReal

  REAL(KREAL) :: A(NReal,NReal)

  REAL(KREAL) :: V(Nreal,Nreal),D(Nreal)

  INTEGER(KINT) :: I, J,IP, IQ, NROT

  REAL(KREAL) :: SM, H, Tresh, G, T, Theta, C, S, Tau

  allocate(B(N),Z(N))

  DO  IP=1,N

     DO IQ=1,N

        V(IP,IQ)=0.

     END DO

     V(IP,IP)=1.

  END DO

  DO  IP=1,N

     B(IP)=A(IP,IP)

     D(IP)=B(IP)

     Z(IP)=0.

  END DO

  NROT=0

  DO  I=1,max_it

     SM=0.

     DO  IP=1,N-1

        DO IQ=IP+1,N

           SM=SM+ABS(A(IP,IQ))

        END DO

     END DO

     IF(SM.EQ.0.) GOTO 100

     IF(I.LT.4)THEN

        TRESH=0.2*SM/N**2

     ELSE

        TRESH=0.

     ENDIF

     DO  IP=1,N-1

        DO  IQ=IP+1,N

           G=100.*ABS(A(IP,IQ))

           IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) &

                .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN

              A(IP,IQ)=0.

           ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN

              H=D(IQ)-D(IP)

              IF(ABS(H)+G.EQ.ABS(H))THEN

                 T=A(IP,IQ)/H

              ELSE

                 THETA=0.5*H/A(IP,IQ)

                 T=1./(ABS(THETA)+SQRT(1.+THETA**2))

                 IF(THETA.LT.0.)T=-T

              ENDIF

              C=1./SQRT(1+T**2)

              S=T*C

              TAU=S/(1.+C)

              H=T*A(IP,IQ)

              Z(IP)=Z(IP)-H

              Z(IQ)=Z(IQ)+H

              D(IP)=D(IP)-H

              D(IQ)=D(IQ)+H

              A(IP,IQ)=0.

              DO J=1,IP-1

                 G=A(J,IP)

                 H=A(J,IQ)

                 A(J,IP)=G-S*(H+G*TAU)

                 A(J,IQ)=H+S*(G-H*TAU)

              END DO

              DO  J=IP+1,IQ-1

                 G=A(IP,J)

                 H=A(J,IQ)

                 A(IP,J)=G-S*(H+G*TAU)

                 A(J,IQ)=H+S*(G-H*TAU)

              END DO

              DO  J=IQ+1,N

                 G=A(IP,J)

                 H=A(IQ,J)

                 A(IP,J)=G-S*(H+G*TAU)

                 A(IQ,J)=H+S*(G-H*TAU)

              END DO

              DO  J=1,N

                 G=V(J,IP)

                 H=V(J,IQ)

                 V(J,IP)=G-S*(H+G*TAU)

                 V(J,IQ)=H+S*(G-H*TAU)

              END DO

              NROT=NROT+1

           ENDIF

        END DO

     END DO

     DO  IP=1,N

        B(IP)=B(IP)+Z(IP)

        D(IP)=B(IP)

        Z(IP)=0.

     END DO

  END DO

  write(6,*) max_it,' iterations should never happen'

  STOP

100 DEALLOCATE(B,Z)

  RETURN

END SUBROUTINE JACOBI

!

!============================================================

!c

!c ++  Sort Eigenvectors in ascending eigenvalues order

!c

!c============================================================

!c

SUBROUTINE trie(nb_niv,ene,psi,nreal)

  !

  ! Input:

  !   Nb_niv      :  Dimension of Psi

  !   NReal  : Actual dimensions of Psi, Ene

  ! Input/output:

  !   Psi(N,N) : Eigvenvectors, changed in output. Stored in a (Nreal,Nreal) matrix

  !   Ene(N)   : Eigenvalues, changed in output. Stored in Ene(NReal)

  !

!!!!!!!!!!!!!!!!

  Use VarTypes

  IMPLICIT NONE

  integer(KINT) :: i, j, k, nb_niv, nreal

  real(KREAL) :: ene(nreal),psi(nreal,nreal)

  real(KREAL) :: a

!!!!!!!!!!!!!!!!

  DO i=1,nb_niv

     DO j=i+1,nb_niv

        IF (ene(i) .GT. ene(j)) THEN

           !              permutation

           a=ene(i)

           ene(i)=ene(j)

           ene(j)=a

           DO k=1,nb_niv

              a=psi(k,i)

              psi(k,i)=psi(k,j)

              psi(k,j)=a

           END DO

        END IF

     END DO

  END DO

END SUBROUTINE trie