Statistics
| Revision:

root / src / Mat_util.f90 @ 7

History | View | Annotate | Download (8.2 kB)

1 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 1 pfleura2
!
3 1 pfleura2
SUBROUTINE GenInv(N,A,InvA,NReal)
4 1 pfleura2
!!!!!!!!!!!!!!!!
5 1 pfleura2
  !
6 1 pfleura2
  !     This subroutines calculates the generalized inverse of a matrix
7 1 pfleura2
  !     It first diagonalize the matrix A, then inverse all non-zero
8 1 pfleura2
  !     eigenvalues, and forms the  InvA matrix using these new eigenvalues
9 1 pfleura2
  !
10 1 pfleura2
  !     Input:
11 1 pfleura2
  !     N : dimension of A
12 1 pfleura2
  !     NReal :Actual dimension of A
13 1 pfleura2
  !     A(N,N) : Initial Matrix, stored in A(Nreal,Nreal)
14 1 pfleura2
  !
15 1 pfleura2
  !     Output:
16 1 pfleura2
  !     InvA(N,N) : Inversed Matrix, stored in a (Nreal,NReal) matrix
17 1 pfleura2
  !
18 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!
19 1 pfleura2
20 1 pfleura2
  Use Vartypes
21 1 pfleura2
  IMPLICIT NONE
22 1 pfleura2
23 1 pfleura2
  INTEGER(KINT), INTENT(IN) :: N,Nreal
24 1 pfleura2
  REAL(KREAL), INTENT(IN) :: A(NReal,NReal)
25 1 pfleura2
  REAL(KREAL), INTENT(OUT) :: InvA(NReal,NReal)
26 1 pfleura2
  !
27 1 pfleura2
28 1 pfleura2
  INTEGER(KINT) :: I,J,K
29 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: EigVec(:,:) ! (Nreal,Nreal)
30 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: EigVal(:) ! (Nreal)
31 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: ATmp(:,:) ! (NReal,Nreal)
32 1 pfleura2
  REAL(KREAL) :: ss
33 1 pfleura2
  !
34 1 pfleura2
  REAL(KREAL), PARAMETER :: eps=1e-12
35 1 pfleura2
36 1 pfleura2
  ALLOCATE(EigVec(Nreal,Nreal), EigVal(Nreal),ATmp(NReal,NReal))
37 1 pfleura2
! A will be destroyed in Jacobi so we save it
38 1 pfleura2
  ATmp=A
39 1 pfleura2
  CALL JAcobi(ATmp,N,EigVal,EigVec,NReal)
40 1 pfleura2
41 1 pfleura2
  DO I=1,N
42 1 pfleura2
     IF (abs(EigVal(I)).GT.eps) EigVal(I)=1.d0/EigVal(I)
43 1 pfleura2
  END DO
44 1 pfleura2
45 1 pfleura2
  InvA=0.d0
46 1 pfleura2
  do k = 1, n
47 1 pfleura2
     do j = 1, n
48 1 pfleura2
        ss = eigval(k) * eigvec(j,k)
49 1 pfleura2
        do i = 1, j
50 1 pfleura2
            InvA(i,j) = InvA(i,j) + ss * eigvec(i,k)
51 1 pfleura2
         end do
52 1 pfleura2
      end do
53 1 pfleura2
   end do
54 1 pfleura2
   do j = 1, n
55 1 pfleura2
      do i = 1, j-1
56 1 pfleura2
         InvA(j,i) = InvA(i,j)
57 1 pfleura2
      end do
58 1 pfleura2
   end do
59 1 pfleura2
60 1 pfleura2
61 1 pfleura2
  DEALLOCATE(EigVec, EigVal,ATmp)
62 1 pfleura2
END SUBROUTINE GenInv
63 1 pfleura2
64 1 pfleura2
!============================================================
65 1 pfleura2
!
66 1 pfleura2
!     ++  Matrix diagonalization Using jacobi
67 1 pfleura2
!  Works only for symmetric matrices
68 1 pfleura2
!     EigvenVectors  : V(i,i)
69 1 pfleura2
!     Eigenvalues : D(i)
70 1 pfleura2
! PFL 30/05/03
71 1 pfleura2
! This versioin uses packed matrices.
72 1 pfleura2
! it unpacks them before calling Jacobi !
73 1 pfleura2
! we have  AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
74 1 pfleura2
!
75 1 pfleura2
!============================================================
76 1 pfleura2
!
77 1 pfleura2
      SUBROUTINE JacPacked(N,AP,D,V,nreal)
78 1 pfleura2
79 1 pfleura2
      Use Vartypes
80 1 pfleura2
      IMPLICIT NONE
81 1 pfleura2
82 1 pfleura2
      INTEGER(KINT), INTENT(IN) :: N,NREAL
83 1 pfleura2
      REAL(KREAL) :: AP(N*(N+1)/2)
84 1 pfleura2
      REAL(KREAL), ALLOCATABLE ::  A(:,:)
85 1 pfleura2
      REAL(KREAL) :: V(Nreal,Nreal),D(Nreal)
86 1 pfleura2
      INTEGER(KINT) :: i,j,nn
87 1 pfleura2
88 1 pfleura2
      allocate(A(nreal,nreal))
89 1 pfleura2
      nn=n*(n+1)/2
90 1 pfleura2
!      WRITE(*,*) 'Jacpa 0'
91 1 pfleura2
!      WRITE(*,'(I3,10(1X,F15.6))') n,(AP(i),i=1,min(nn,5))
92 1 pfleura2
 !     WRITE(*,*) 'Jacpa 0'
93 1 pfleura2
      do j=1,n
94 1 pfleura2
         do i=1,j
95 1 pfleura2
!            WRITE(*,*) i,j
96 1 pfleura2
            A(i,J)=AP(i + (j-1)*j/2)
97 1 pfleura2
            A(j,I)=A(i,J)
98 1 pfleura2
         end do
99 1 pfleura2
      end do
100 1 pfleura2
!      do j=1,n
101 1 pfleura2
!         WRITE(*,'(10(1X,F15.6))') (A(i,J),i=1,min(n,5))
102 1 pfleura2
!      end do
103 1 pfleura2
104 1 pfleura2
!      WRITE(*,*) 'Jacpa 1'
105 1 pfleura2
      call Jacobi(A,n,D,V,Nreal)
106 1 pfleura2
!      WRITE(*,*) 'Jacpa 2'
107 1 pfleura2
!      DO i=1,n
108 1 pfleura2
!         WRITE(*,'(1X,I3,10(1X,F15.6))') i,D(i),(V(j,i),j=1,min(n,5))
109 1 pfleura2
!      end do
110 1 pfleura2
      deallocate(a)
111 1 pfleura2
112 1 pfleura2
    end SUBROUTINE JacPacked
113 1 pfleura2
114 1 pfleura2
115 1 pfleura2
116 1 pfleura2
!
117 1 pfleura2
!============================================================
118 1 pfleura2
!
119 1 pfleura2
!     ++  Matrix diagonalization Using jacobi
120 1 pfleura2
!  Works only for symmetric matrices
121 1 pfleura2
!     EigvenVectors  : V
122 1 pfleura2
!     Eigenvalues : D
123 1 pfleura2
!
124 1 pfleura2
!============================================================
125 1 pfleura2
!
126 1 pfleura2
SUBROUTINE JACOBI(A,N,D,V,Nreal)
127 1 pfleura2
128 1 pfleura2
!!!!!!!!!!!!!!!!
129 1 pfleura2
  !
130 1 pfleura2
  !     Input:
131 1 pfleura2
  !     N      :  Dimension of A
132 1 pfleura2
  !     NReal  : Actual dimensions of A, D and V.
133 1 pfleura2
  !
134 1 pfleura2
  !     Input/output:
135 1 pfleura2
  !     A(N,N) : Matrix to be diagonalized, store in a (Nreal,Nreal) matrix
136 1 pfleura2
  !              Destroyed in output.
137 1 pfleura2
  !     Output:
138 1 pfleura2
  !     V(N,N) : Eigenvectors, stored in V(NReal, NReal)
139 1 pfleura2
  !     D(N)   : Eigenvalues, stored in D(NReal)
140 1 pfleura2
  !
141 1 pfleura2
142 1 pfleura2
  Use Vartypes
143 1 pfleura2
144 1 pfleura2
  IMPLICIT NONE
145 1 pfleura2
  INTEGER(KINT), parameter :: max_it=500
146 1 pfleura2
  REAL(KREAL), ALLOCATABLE ::  B(:),Z(:)
147 1 pfleura2
148 1 pfleura2
  INTEGER(KINT) :: N,NReal
149 1 pfleura2
  REAL(KREAL) :: A(NReal,NReal)
150 1 pfleura2
  REAL(KREAL) :: V(Nreal,Nreal),D(Nreal)
151 1 pfleura2
152 1 pfleura2
  INTEGER(KINT) :: I, J,IP, IQ, NROT
153 1 pfleura2
  REAL(KREAL) :: SM, H, Tresh, G, T, Theta, C, S, Tau
154 1 pfleura2
155 1 pfleura2
  allocate(B(N),Z(N))
156 1 pfleura2
157 1 pfleura2
  DO  IP=1,N
158 1 pfleura2
     DO IQ=1,N
159 1 pfleura2
        V(IP,IQ)=0.
160 1 pfleura2
     END DO
161 1 pfleura2
     V(IP,IP)=1.
162 1 pfleura2
  END DO
163 1 pfleura2
  DO  IP=1,N
164 1 pfleura2
     B(IP)=A(IP,IP)
165 1 pfleura2
     D(IP)=B(IP)
166 1 pfleura2
     Z(IP)=0.
167 1 pfleura2
  END DO
168 1 pfleura2
  NROT=0
169 1 pfleura2
  DO  I=1,max_it
170 1 pfleura2
     SM=0.
171 1 pfleura2
     DO  IP=1,N-1
172 1 pfleura2
        DO IQ=IP+1,N
173 1 pfleura2
           SM=SM+ABS(A(IP,IQ))
174 1 pfleura2
        END DO
175 1 pfleura2
     END DO
176 1 pfleura2
     IF(SM.EQ.0.) GOTO 100
177 1 pfleura2
     IF(I.LT.4)THEN
178 1 pfleura2
        TRESH=0.2*SM/N**2
179 1 pfleura2
     ELSE
180 1 pfleura2
        TRESH=0.
181 1 pfleura2
     ENDIF
182 1 pfleura2
     DO  IP=1,N-1
183 1 pfleura2
        DO  IQ=IP+1,N
184 1 pfleura2
           G=100.*ABS(A(IP,IQ))
185 1 pfleura2
           IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) &
186 1 pfleura2
                .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
187 1 pfleura2
              A(IP,IQ)=0.
188 1 pfleura2
           ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN
189 1 pfleura2
              H=D(IQ)-D(IP)
190 1 pfleura2
              IF(ABS(H)+G.EQ.ABS(H))THEN
191 1 pfleura2
                 T=A(IP,IQ)/H
192 1 pfleura2
              ELSE
193 1 pfleura2
                 THETA=0.5*H/A(IP,IQ)
194 1 pfleura2
                 T=1./(ABS(THETA)+SQRT(1.+THETA**2))
195 1 pfleura2
                 IF(THETA.LT.0.)T=-T
196 1 pfleura2
              ENDIF
197 1 pfleura2
              C=1./SQRT(1+T**2)
198 1 pfleura2
              S=T*C
199 1 pfleura2
              TAU=S/(1.+C)
200 1 pfleura2
              H=T*A(IP,IQ)
201 1 pfleura2
              Z(IP)=Z(IP)-H
202 1 pfleura2
              Z(IQ)=Z(IQ)+H
203 1 pfleura2
              D(IP)=D(IP)-H
204 1 pfleura2
              D(IQ)=D(IQ)+H
205 1 pfleura2
              A(IP,IQ)=0.
206 1 pfleura2
              DO J=1,IP-1
207 1 pfleura2
                 G=A(J,IP)
208 1 pfleura2
                 H=A(J,IQ)
209 1 pfleura2
                 A(J,IP)=G-S*(H+G*TAU)
210 1 pfleura2
                 A(J,IQ)=H+S*(G-H*TAU)
211 1 pfleura2
              END DO
212 1 pfleura2
              DO  J=IP+1,IQ-1
213 1 pfleura2
                 G=A(IP,J)
214 1 pfleura2
                 H=A(J,IQ)
215 1 pfleura2
                 A(IP,J)=G-S*(H+G*TAU)
216 1 pfleura2
                 A(J,IQ)=H+S*(G-H*TAU)
217 1 pfleura2
              END DO
218 1 pfleura2
              DO  J=IQ+1,N
219 1 pfleura2
                 G=A(IP,J)
220 1 pfleura2
                 H=A(IQ,J)
221 1 pfleura2
                 A(IP,J)=G-S*(H+G*TAU)
222 1 pfleura2
                 A(IQ,J)=H+S*(G-H*TAU)
223 1 pfleura2
              END DO
224 1 pfleura2
              DO  J=1,N
225 1 pfleura2
                 G=V(J,IP)
226 1 pfleura2
                 H=V(J,IQ)
227 1 pfleura2
                 V(J,IP)=G-S*(H+G*TAU)
228 1 pfleura2
                 V(J,IQ)=H+S*(G-H*TAU)
229 1 pfleura2
              END DO
230 1 pfleura2
              NROT=NROT+1
231 1 pfleura2
           ENDIF
232 1 pfleura2
        END DO
233 1 pfleura2
     END DO
234 1 pfleura2
     DO  IP=1,N
235 1 pfleura2
        B(IP)=B(IP)+Z(IP)
236 1 pfleura2
        D(IP)=B(IP)
237 1 pfleura2
        Z(IP)=0.
238 1 pfleura2
     END DO
239 1 pfleura2
  END DO
240 1 pfleura2
  write(6,*) max_it,' iterations should never happen'
241 1 pfleura2
  STOP
242 1 pfleura2
100 DEALLOCATE(B,Z)
243 1 pfleura2
  RETURN
244 1 pfleura2
END SUBROUTINE JACOBI
245 1 pfleura2
!
246 1 pfleura2
!============================================================
247 1 pfleura2
!c
248 1 pfleura2
!c ++  Sort Eigenvectors in ascending eigenvalues order
249 1 pfleura2
!c
250 1 pfleura2
!c============================================================
251 1 pfleura2
!c
252 1 pfleura2
SUBROUTINE trie(nb_niv,ene,psi,nreal)
253 1 pfleura2
254 1 pfleura2
  !
255 1 pfleura2
  ! Input:
256 1 pfleura2
  !   Nb_niv      :  Dimension of Psi
257 1 pfleura2
  !   NReal  : Actual dimensions of Psi, Ene
258 1 pfleura2
  ! Input/output:
259 1 pfleura2
  !   Psi(N,N) : Eigvenvectors, changed in output. Stored in a (Nreal,Nreal) matrix
260 1 pfleura2
  !   Ene(N)   : Eigenvalues, changed in output. Stored in Ene(NReal)
261 1 pfleura2
  !
262 1 pfleura2
!!!!!!!!!!!!!!!!
263 1 pfleura2
264 1 pfleura2
  Use VarTypes
265 1 pfleura2
  IMPLICIT NONE
266 1 pfleura2
267 2 pfleura2
  integer(KINT) :: i, j, k, nb_niv, nreal
268 1 pfleura2
  real(KREAL) :: ene(nreal),psi(nreal,nreal)
269 1 pfleura2
  real(KREAL) :: a
270 1 pfleura2
271 1 pfleura2
272 1 pfleura2
!!!!!!!!!!!!!!!!
273 1 pfleura2
274 1 pfleura2
275 1 pfleura2
  DO i=1,nb_niv
276 1 pfleura2
     DO j=i+1,nb_niv
277 1 pfleura2
        IF (ene(i) .GT. ene(j)) THEN
278 1 pfleura2
           !              permutation
279 1 pfleura2
           a=ene(i)
280 1 pfleura2
           ene(i)=ene(j)
281 1 pfleura2
           ene(j)=a
282 1 pfleura2
283 1 pfleura2
           DO k=1,nb_niv
284 1 pfleura2
              a=psi(k,i)
285 1 pfleura2
              psi(k,i)=psi(k,j)
286 1 pfleura2
              psi(k,j)=a
287 1 pfleura2
           END DO
288 1 pfleura2
        END IF
289 1 pfleura2
     END DO
290 1 pfleura2
  END DO
291 1 pfleura2
292 1 pfleura2
END SUBROUTINE trie
293 1 pfleura2
294 1 pfleura2
!============================================================
295 1 pfleura2
!c
296 1 pfleura2
!c ++  Sort Eigenvectors in ascending eigenvalues order
297 1 pfleura2
!c
298 1 pfleura2
!c============================================================
299 1 pfleura2
!c
300 1 pfleura2
SUBROUTINE SortEigenSys(N,Eigval,Eigvec)
301 1 pfleura2
302 1 pfleura2
303 1 pfleura2
  !
304 1 pfleura2
  ! Input/output:
305 1 pfleura2
  !   N : dimension of the system
306 1 pfleura2
  !   EigVec(N,N) : Eigvenvectors, changed in output. Stored in a (N,N) matrix
307 1 pfleura2
  !   EigVal(N)   : Eigenvalues, changed in output. Stored in a (n) vector
308 1 pfleura2
  !
309 1 pfleura2
 ! Process:
310 1 pfleura2
 ! We use first a ranking, then a working array to reorder the eigenvalues and eigenvectors
311 1 pfleura2
312 1 pfleura2
!!!!!!!!!!!!!!!!
313 1 pfleura2
314 1 pfleura2
  Use VarTypes
315 1 pfleura2
  use m_mrgrnk
316 1 pfleura2
  IMPLICIT NONE
317 1 pfleura2
318 1 pfleura2
319 1 pfleura2
  INTEGER(KINT), INTENT(IN) :: N
320 1 pfleura2
  REAL(KREAL), INTENT(OUT) :: EigVal(N), Eigvec(N,N)
321 1 pfleura2
322 2 pfleura2
  integer(KINT) :: i
323 1 pfleura2
324 1 pfleura2
  INTEGER(KINT), ALLOCATABLE :: Rank(:) !N
325 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: EigValT(:) !n
326 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: EigVecT(:,:) !n,n
327 1 pfleura2
328 1 pfleura2
!!!!!!!!!!!!!!!!
329 1 pfleura2
  ALLOCATE(Rank(n),EigValT(n),EigVecT(n,n))
330 1 pfleura2
  CALL MrgRnk(EigVal,Rank)
331 1 pfleura2
332 1 pfleura2
  DO I=1,N
333 1 pfleura2
     EigValT(I)=EigVal(Rank(I))
334 1 pfleura2
     EigVecT(I,1:N)=EigVec(Rank(I),1:N)
335 1 pfleura2
  END DO
336 1 pfleura2
  EigVal=EigValT
337 1 pfleura2
  EigVec=EigVecT
338 1 pfleura2
339 1 pfleura2
  DEALLOCATE(Rank,EigValT,EigVecT)
340 1 pfleura2
341 1 pfleura2
342 1 pfleura2
END SUBROUTINE SortEigenSys