Statistics
| Revision:

## root / src / Mat_util.f90 @ 7

 1 2 3 1 pfleura2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  1 pfleura2 !  1 pfleura2 SUBROUTINE GenInv(N,A,InvA,NReal)  1 pfleura2 !!!!!!!!!!!!!!!!  1 pfleura2  !  1 pfleura2  ! This subroutines calculates the generalized inverse of a matrix  1 pfleura2  ! It first diagonalize the matrix A, then inverse all non-zero  1 pfleura2  ! eigenvalues, and forms the InvA matrix using these new eigenvalues  1 pfleura2  !  1 pfleura2  ! Input:  1 pfleura2  ! N : dimension of A  1 pfleura2  ! NReal :Actual dimension of A  1 pfleura2  ! A(N,N) : Initial Matrix, stored in A(Nreal,Nreal)  1 pfleura2  !  1 pfleura2  ! Output:  1 pfleura2  ! InvA(N,N) : Inversed Matrix, stored in a (Nreal,NReal) matrix  1 pfleura2  !  1 pfleura2 !!!!!!!!!!!!!!!!!!!!!!!  1 pfleura2 1 pfleura2  Use Vartypes  1 pfleura2  IMPLICIT NONE  1 pfleura2 1 pfleura2  INTEGER(KINT), INTENT(IN) :: N,Nreal  1 pfleura2  REAL(KREAL), INTENT(IN) :: A(NReal,NReal)  1 pfleura2  REAL(KREAL), INTENT(OUT) :: InvA(NReal,NReal)  1 pfleura2  !  1 pfleura2 1 pfleura2  INTEGER(KINT) :: I,J,K  1 pfleura2  REAL(KREAL), ALLOCATABLE :: EigVec(:,:) ! (Nreal,Nreal)  1 pfleura2  REAL(KREAL), ALLOCATABLE :: EigVal(:) ! (Nreal)  1 pfleura2  REAL(KREAL), ALLOCATABLE :: ATmp(:,:) ! (NReal,Nreal)  1 pfleura2  REAL(KREAL) :: ss  1 pfleura2  !  1 pfleura2  REAL(KREAL), PARAMETER :: eps=1e-12  1 pfleura2 1 pfleura2  ALLOCATE(EigVec(Nreal,Nreal), EigVal(Nreal),ATmp(NReal,NReal))  1 pfleura2 ! A will be destroyed in Jacobi so we save it  1 pfleura2  ATmp=A  1 pfleura2  CALL JAcobi(ATmp,N,EigVal,EigVec,NReal)  1 pfleura2 1 pfleura2  DO I=1,N  1 pfleura2  IF (abs(EigVal(I)).GT.eps) EigVal(I)=1.d0/EigVal(I)  1 pfleura2  END DO  1 pfleura2 1 pfleura2  InvA=0.d0  1 pfleura2  do k = 1, n  1 pfleura2  do j = 1, n  1 pfleura2  ss = eigval(k) * eigvec(j,k)  1 pfleura2  do i = 1, j  1 pfleura2  InvA(i,j) = InvA(i,j) + ss * eigvec(i,k)  1 pfleura2  end do  1 pfleura2  end do  1 pfleura2  end do  1 pfleura2  do j = 1, n  1 pfleura2  do i = 1, j-1  1 pfleura2  InvA(j,i) = InvA(i,j)  1 pfleura2  end do  1 pfleura2  end do  1 pfleura2 1 pfleura2 1 pfleura2  DEALLOCATE(EigVec, EigVal,ATmp)  1 pfleura2 END SUBROUTINE GenInv  1 pfleura2 1 pfleura2 !============================================================  1 pfleura2 !  1 pfleura2 ! ++ Matrix diagonalization Using jacobi  1 pfleura2 ! Works only for symmetric matrices  1 pfleura2 ! EigvenVectors : V(i,i)  1 pfleura2 ! Eigenvalues : D(i)  1 pfleura2 ! PFL 30/05/03  1 pfleura2 ! This versioin uses packed matrices.  1 pfleura2 ! it unpacks them before calling Jacobi !  1 pfleura2 ! we have AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;  1 pfleura2 !  1 pfleura2 !============================================================  1 pfleura2 !  1 pfleura2  SUBROUTINE JacPacked(N,AP,D,V,nreal)  1 pfleura2 1 pfleura2  Use Vartypes  1 pfleura2  IMPLICIT NONE  1 pfleura2 1 pfleura2  INTEGER(KINT), INTENT(IN) :: N,NREAL  1 pfleura2  REAL(KREAL) :: AP(N*(N+1)/2)  1 pfleura2  REAL(KREAL), ALLOCATABLE :: A(:,:)  1 pfleura2  REAL(KREAL) :: V(Nreal,Nreal),D(Nreal)  1 pfleura2  INTEGER(KINT) :: i,j,nn  1 pfleura2 1 pfleura2  allocate(A(nreal,nreal))  1 pfleura2  nn=n*(n+1)/2  1 pfleura2 ! WRITE(*,*) 'Jacpa 0'  1 pfleura2 ! WRITE(*,'(I3,10(1X,F15.6))') n,(AP(i),i=1,min(nn,5))  1 pfleura2  ! WRITE(*,*) 'Jacpa 0'  1 pfleura2  do j=1,n  1 pfleura2  do i=1,j  1 pfleura2 ! WRITE(*,*) i,j  1 pfleura2  A(i,J)=AP(i + (j-1)*j/2)  1 pfleura2  A(j,I)=A(i,J)  1 pfleura2  end do  1 pfleura2  end do  1 pfleura2 ! do j=1,n  1 pfleura2 ! WRITE(*,'(10(1X,F15.6))') (A(i,J),i=1,min(n,5))  1 pfleura2 ! end do  1 pfleura2 1 pfleura2 ! WRITE(*,*) 'Jacpa 1'  1 pfleura2  call Jacobi(A,n,D,V,Nreal)  1 pfleura2 ! WRITE(*,*) 'Jacpa 2'  1 pfleura2 ! DO i=1,n  1 pfleura2 ! WRITE(*,'(1X,I3,10(1X,F15.6))') i,D(i),(V(j,i),j=1,min(n,5))  1 pfleura2 ! end do  1 pfleura2  deallocate(a)  1 pfleura2 1 pfleura2  end SUBROUTINE JacPacked  1 pfleura2 1 pfleura2 1 pfleura2 1 pfleura2 !  1 pfleura2 !============================================================  1 pfleura2 !  1 pfleura2 ! ++ Matrix diagonalization Using jacobi  1 pfleura2 ! Works only for symmetric matrices  1 pfleura2 ! EigvenVectors : V  1 pfleura2 ! Eigenvalues : D  1 pfleura2 !  1 pfleura2 !============================================================  1 pfleura2 !  1 pfleura2 SUBROUTINE JACOBI(A,N,D,V,Nreal)  1 pfleura2 1 pfleura2 !!!!!!!!!!!!!!!!  1 pfleura2  !  1 pfleura2  ! Input:  1 pfleura2  ! N : Dimension of A  1 pfleura2  ! NReal : Actual dimensions of A, D and V.  1 pfleura2  !  1 pfleura2  ! Input/output:  1 pfleura2  ! A(N,N) : Matrix to be diagonalized, store in a (Nreal,Nreal) matrix  1 pfleura2  ! Destroyed in output.  1 pfleura2  ! Output:  1 pfleura2  ! V(N,N) : Eigenvectors, stored in V(NReal, NReal)  1 pfleura2  ! D(N) : Eigenvalues, stored in D(NReal)  1 pfleura2  !  1 pfleura2 1 pfleura2  Use Vartypes  1 pfleura2 1 pfleura2  IMPLICIT NONE  1 pfleura2  INTEGER(KINT), parameter :: max_it=500  1 pfleura2  REAL(KREAL), ALLOCATABLE :: B(:),Z(:)  1 pfleura2 1 pfleura2  INTEGER(KINT) :: N,NReal  1 pfleura2  REAL(KREAL) :: A(NReal,NReal)  1 pfleura2  REAL(KREAL) :: V(Nreal,Nreal),D(Nreal)  1 pfleura2 1 pfleura2  INTEGER(KINT) :: I, J,IP, IQ, NROT  1 pfleura2  REAL(KREAL) :: SM, H, Tresh, G, T, Theta, C, S, Tau  1 pfleura2 1 pfleura2  allocate(B(N),Z(N))  1 pfleura2 1 pfleura2  DO IP=1,N  1 pfleura2  DO IQ=1,N  1 pfleura2  V(IP,IQ)=0.  1 pfleura2  END DO  1 pfleura2  V(IP,IP)=1.  1 pfleura2  END DO  1 pfleura2  DO IP=1,N  1 pfleura2  B(IP)=A(IP,IP)  1 pfleura2  D(IP)=B(IP)  1 pfleura2  Z(IP)=0.  1 pfleura2  END DO  1 pfleura2  NROT=0  1 pfleura2  DO I=1,max_it  1 pfleura2  SM=0.  1 pfleura2  DO IP=1,N-1  1 pfleura2  DO IQ=IP+1,N  1 pfleura2  SM=SM+ABS(A(IP,IQ))  1 pfleura2  END DO  1 pfleura2  END DO  1 pfleura2  IF(SM.EQ.0.) GOTO 100  1 pfleura2  IF(I.LT.4)THEN  1 pfleura2  TRESH=0.2*SM/N**2  1 pfleura2  ELSE  1 pfleura2  TRESH=0.  1 pfleura2  ENDIF  1 pfleura2  DO IP=1,N-1  1 pfleura2  DO IQ=IP+1,N  1 pfleura2  G=100.*ABS(A(IP,IQ))  1 pfleura2  IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) &  1 pfleura2  .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN  1 pfleura2  A(IP,IQ)=0.  1 pfleura2  ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN  1 pfleura2  H=D(IQ)-D(IP)  1 pfleura2  IF(ABS(H)+G.EQ.ABS(H))THEN  1 pfleura2  T=A(IP,IQ)/H  1 pfleura2  ELSE  1 pfleura2  THETA=0.5*H/A(IP,IQ)  1 pfleura2  T=1./(ABS(THETA)+SQRT(1.+THETA**2))  1 pfleura2  IF(THETA.LT.0.)T=-T  1 pfleura2  ENDIF  1 pfleura2  C=1./SQRT(1+T**2)  1 pfleura2  S=T*C  1 pfleura2  TAU=S/(1.+C)  1 pfleura2  H=T*A(IP,IQ)  1 pfleura2  Z(IP)=Z(IP)-H  1 pfleura2  Z(IQ)=Z(IQ)+H  1 pfleura2  D(IP)=D(IP)-H  1 pfleura2  D(IQ)=D(IQ)+H  1 pfleura2  A(IP,IQ)=0.  1 pfleura2  DO J=1,IP-1  1 pfleura2  G=A(J,IP)  1 pfleura2  H=A(J,IQ)  1 pfleura2  A(J,IP)=G-S*(H+G*TAU)  1 pfleura2  A(J,IQ)=H+S*(G-H*TAU)  1 pfleura2  END DO  1 pfleura2  DO J=IP+1,IQ-1  1 pfleura2  G=A(IP,J)  1 pfleura2  H=A(J,IQ)  1 pfleura2  A(IP,J)=G-S*(H+G*TAU)  1 pfleura2  A(J,IQ)=H+S*(G-H*TAU)  1 pfleura2  END DO  1 pfleura2  DO J=IQ+1,N  1 pfleura2  G=A(IP,J)  1 pfleura2  H=A(IQ,J)  1 pfleura2  A(IP,J)=G-S*(H+G*TAU)  1 pfleura2  A(IQ,J)=H+S*(G-H*TAU)  1 pfleura2  END DO  1 pfleura2  DO J=1,N  1 pfleura2  G=V(J,IP)  1 pfleura2  H=V(J,IQ)  1 pfleura2  V(J,IP)=G-S*(H+G*TAU)  1 pfleura2  V(J,IQ)=H+S*(G-H*TAU)  1 pfleura2  END DO  1 pfleura2  NROT=NROT+1  1 pfleura2  ENDIF  1 pfleura2  END DO  1 pfleura2  END DO  1 pfleura2  DO IP=1,N  1 pfleura2  B(IP)=B(IP)+Z(IP)  1 pfleura2  D(IP)=B(IP)  1 pfleura2  Z(IP)=0.  1 pfleura2  END DO  1 pfleura2  END DO  1 pfleura2  write(6,*) max_it,' iterations should never happen'  1 pfleura2  STOP  1 pfleura2 100 DEALLOCATE(B,Z)  1 pfleura2  RETURN  1 pfleura2 END SUBROUTINE JACOBI  1 pfleura2 !  1 pfleura2 !============================================================  1 pfleura2 !c  1 pfleura2 !c ++ Sort Eigenvectors in ascending eigenvalues order  1 pfleura2 !c  1 pfleura2 !c============================================================  1 pfleura2 !c  1 pfleura2 SUBROUTINE trie(nb_niv,ene,psi,nreal)  1 pfleura2 1 pfleura2  !  1 pfleura2  ! Input:  1 pfleura2  ! Nb_niv : Dimension of Psi  1 pfleura2  ! NReal : Actual dimensions of Psi, Ene  1 pfleura2  ! Input/output:  1 pfleura2  ! Psi(N,N) : Eigvenvectors, changed in output. Stored in a (Nreal,Nreal) matrix  1 pfleura2  ! Ene(N) : Eigenvalues, changed in output. Stored in Ene(NReal)  1 pfleura2  !  1 pfleura2 !!!!!!!!!!!!!!!!  1 pfleura2 1 pfleura2  Use VarTypes  1 pfleura2  IMPLICIT NONE  1 pfleura2 2 pfleura2  integer(KINT) :: i, j, k, nb_niv, nreal  1 pfleura2  real(KREAL) :: ene(nreal),psi(nreal,nreal)  1 pfleura2  real(KREAL) :: a  1 pfleura2 1 pfleura2 1 pfleura2 !!!!!!!!!!!!!!!!  1 pfleura2 1 pfleura2 1 pfleura2  DO i=1,nb_niv  1 pfleura2  DO j=i+1,nb_niv  1 pfleura2  IF (ene(i) .GT. ene(j)) THEN  1 pfleura2  ! permutation  1 pfleura2  a=ene(i)  1 pfleura2  ene(i)=ene(j)  1 pfleura2  ene(j)=a  1 pfleura2 1 pfleura2  DO k=1,nb_niv  1 pfleura2  a=psi(k,i)  1 pfleura2  psi(k,i)=psi(k,j)  1 pfleura2  psi(k,j)=a  1 pfleura2  END DO  1 pfleura2  END IF  1 pfleura2  END DO  1 pfleura2  END DO  1 pfleura2 1 pfleura2 END SUBROUTINE trie  1 pfleura2 1 pfleura2 !============================================================  1 pfleura2 !c  1 pfleura2 !c ++ Sort Eigenvectors in ascending eigenvalues order  1 pfleura2 !c  1 pfleura2 !c============================================================  1 pfleura2 !c  1 pfleura2 SUBROUTINE SortEigenSys(N,Eigval,Eigvec)  1 pfleura2 1 pfleura2 1 pfleura2  !  1 pfleura2  ! Input/output:  1 pfleura2  ! N : dimension of the system  1 pfleura2  ! EigVec(N,N) : Eigvenvectors, changed in output. Stored in a (N,N) matrix  1 pfleura2  ! EigVal(N) : Eigenvalues, changed in output. Stored in a (n) vector  1 pfleura2  !  1 pfleura2  ! Process:  1 pfleura2  ! We use first a ranking, then a working array to reorder the eigenvalues and eigenvectors  1 pfleura2 1 pfleura2 !!!!!!!!!!!!!!!!  1 pfleura2 1 pfleura2  Use VarTypes  1 pfleura2  use m_mrgrnk  1 pfleura2  IMPLICIT NONE  1 pfleura2 1 pfleura2 1 pfleura2  INTEGER(KINT), INTENT(IN) :: N  1 pfleura2  REAL(KREAL), INTENT(OUT) :: EigVal(N), Eigvec(N,N)  1 pfleura2 2 pfleura2  integer(KINT) :: i  1 pfleura2 1 pfleura2  INTEGER(KINT), ALLOCATABLE :: Rank(:) !N  1 pfleura2  REAL(KREAL), ALLOCATABLE :: EigValT(:) !n  1 pfleura2  REAL(KREAL), ALLOCATABLE :: EigVecT(:,:) !n,n  1 pfleura2 1 pfleura2 !!!!!!!!!!!!!!!!  1 pfleura2  ALLOCATE(Rank(n),EigValT(n),EigVecT(n,n))  1 pfleura2  CALL MrgRnk(EigVal,Rank)  1 pfleura2 1 pfleura2  DO I=1,N  1 pfleura2  EigValT(I)=EigVal(Rank(I))  1 pfleura2  EigVecT(I,1:N)=EigVec(Rank(I),1:N)  1 pfleura2  END DO  1 pfleura2  EigVal=EigValT  1 pfleura2  EigVec=EigVecT  1 pfleura2 1 pfleura2  DEALLOCATE(Rank,EigValT,EigVecT)  1 pfleura2 1 pfleura2 1 pfleura2 END SUBROUTINE SortEigenSys