Statistics
| Revision:

root / src / Mat_util.f90 @ 7

History | View | Annotate | Download (8.2 kB)

1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2
!     
3
SUBROUTINE GenInv(N,A,InvA,NReal)
4
!!!!!!!!!!!!!!!!
5
  !     
6
  !     This subroutines calculates the generalized inverse of a matrix
7
  !     It first diagonalize the matrix A, then inverse all non-zero
8
  !     eigenvalues, and forms the  InvA matrix using these new eigenvalues
9
  !     
10
  !     Input:
11
  !     N : dimension of A
12
  !     NReal :Actual dimension of A
13
  !     A(N,N) : Initial Matrix, stored in A(Nreal,Nreal)
14
  !     
15
  !     Output:
16
  !     InvA(N,N) : Inversed Matrix, stored in a (Nreal,NReal) matrix
17
  !     
18
!!!!!!!!!!!!!!!!!!!!!!!
19

    
20
  Use Vartypes
21
  IMPLICIT NONE
22

    
23
  INTEGER(KINT), INTENT(IN) :: N,Nreal
24
  REAL(KREAL), INTENT(IN) :: A(NReal,NReal)
25
  REAL(KREAL), INTENT(OUT) :: InvA(NReal,NReal)
26
  !     
27

    
28
  INTEGER(KINT) :: I,J,K
29
  REAL(KREAL), ALLOCATABLE :: EigVec(:,:) ! (Nreal,Nreal)
30
  REAL(KREAL), ALLOCATABLE :: EigVal(:) ! (Nreal)
31
  REAL(KREAL), ALLOCATABLE :: ATmp(:,:) ! (NReal,Nreal)
32
  REAL(KREAL) :: ss
33
  !     
34
  REAL(KREAL), PARAMETER :: eps=1e-12
35

    
36
  ALLOCATE(EigVec(Nreal,Nreal), EigVal(Nreal),ATmp(NReal,NReal))
37
! A will be destroyed in Jacobi so we save it
38
  ATmp=A
39
  CALL JAcobi(ATmp,N,EigVal,EigVec,NReal)
40

    
41
  DO I=1,N
42
     IF (abs(EigVal(I)).GT.eps) EigVal(I)=1.d0/EigVal(I)
43
  END DO
44

    
45
  InvA=0.d0
46
  do k = 1, n
47
     do j = 1, n
48
        ss = eigval(k) * eigvec(j,k)
49
        do i = 1, j
50
            InvA(i,j) = InvA(i,j) + ss * eigvec(i,k)
51
         end do
52
      end do
53
   end do
54
   do j = 1, n
55
      do i = 1, j-1
56
         InvA(j,i) = InvA(i,j)
57
      end do
58
   end do
59

    
60

    
61
  DEALLOCATE(EigVec, EigVal,ATmp)
62
END SUBROUTINE GenInv
63

    
64
!============================================================
65
!
66
!     ++  Matrix diagonalization Using jacobi
67
!  Works only for symmetric matrices
68
!     EigvenVectors  : V(i,i)
69
!     Eigenvalues : D(i)
70
! PFL 30/05/03
71
! This versioin uses packed matrices.
72
! it unpacks them before calling Jacobi !
73
! we have  AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
74
!
75
!============================================================
76
!
77
      SUBROUTINE JacPacked(N,AP,D,V,nreal)
78

    
79
      Use Vartypes
80
      IMPLICIT NONE
81

    
82
      INTEGER(KINT), INTENT(IN) :: N,NREAL
83
      REAL(KREAL) :: AP(N*(N+1)/2)
84
      REAL(KREAL), ALLOCATABLE ::  A(:,:)
85
      REAL(KREAL) :: V(Nreal,Nreal),D(Nreal)
86
      INTEGER(KINT) :: i,j,nn
87

    
88
      allocate(A(nreal,nreal))
89
      nn=n*(n+1)/2
90
!      WRITE(*,*) 'Jacpa 0'
91
!      WRITE(*,'(I3,10(1X,F15.6))') n,(AP(i),i=1,min(nn,5))
92
 !     WRITE(*,*) 'Jacpa 0'
93
      do j=1,n
94
         do i=1,j
95
!            WRITE(*,*) i,j
96
            A(i,J)=AP(i + (j-1)*j/2)
97
            A(j,I)=A(i,J)
98
         end do
99
      end do
100
!      do j=1,n
101
!         WRITE(*,'(10(1X,F15.6))') (A(i,J),i=1,min(n,5))
102
!      end do
103

    
104
!      WRITE(*,*) 'Jacpa 1'
105
      call Jacobi(A,n,D,V,Nreal)
106
!      WRITE(*,*) 'Jacpa 2'
107
!      DO i=1,n
108
!         WRITE(*,'(1X,I3,10(1X,F15.6))') i,D(i),(V(j,i),j=1,min(n,5))
109
!      end do
110
      deallocate(a)
111

    
112
    end SUBROUTINE JacPacked
113

    
114

    
115

    
116
!     
117
!============================================================
118
!     
119
!     ++  Matrix diagonalization Using jacobi
120
!  Works only for symmetric matrices
121
!     EigvenVectors  : V
122
!     Eigenvalues : D
123
!     
124
!============================================================
125
!     
126
SUBROUTINE JACOBI(A,N,D,V,Nreal)
127

    
128
!!!!!!!!!!!!!!!!
129
  !     
130
  !     Input:
131
  !     N      :  Dimension of A
132
  !     NReal  : Actual dimensions of A, D and V.
133
  !
134
  !     Input/output:
135
  !     A(N,N) : Matrix to be diagonalized, store in a (Nreal,Nreal) matrix
136
  !              Destroyed in output.
137
  !     Output:
138
  !     V(N,N) : Eigenvectors, stored in V(NReal, NReal)
139
  !     D(N)   : Eigenvalues, stored in D(NReal)
140
  !     
141

    
142
  Use Vartypes
143

    
144
  IMPLICIT NONE
145
  INTEGER(KINT), parameter :: max_it=500
146
  REAL(KREAL), ALLOCATABLE ::  B(:),Z(:)
147

    
148
  INTEGER(KINT) :: N,NReal
149
  REAL(KREAL) :: A(NReal,NReal)
150
  REAL(KREAL) :: V(Nreal,Nreal),D(Nreal)
151

    
152
  INTEGER(KINT) :: I, J,IP, IQ, NROT
153
  REAL(KREAL) :: SM, H, Tresh, G, T, Theta, C, S, Tau
154

    
155
  allocate(B(N),Z(N))
156

    
157
  DO  IP=1,N
158
     DO IQ=1,N
159
        V(IP,IQ)=0.
160
     END DO
161
     V(IP,IP)=1.
162
  END DO
163
  DO  IP=1,N
164
     B(IP)=A(IP,IP)
165
     D(IP)=B(IP)
166
     Z(IP)=0.
167
  END DO
168
  NROT=0
169
  DO  I=1,max_it
170
     SM=0.
171
     DO  IP=1,N-1
172
        DO IQ=IP+1,N
173
           SM=SM+ABS(A(IP,IQ))
174
        END DO
175
     END DO
176
     IF(SM.EQ.0.) GOTO 100
177
     IF(I.LT.4)THEN
178
        TRESH=0.2*SM/N**2
179
     ELSE
180
        TRESH=0.
181
     ENDIF
182
     DO  IP=1,N-1
183
        DO  IQ=IP+1,N
184
           G=100.*ABS(A(IP,IQ))
185
           IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) &
186
                .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
187
              A(IP,IQ)=0.
188
           ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN
189
              H=D(IQ)-D(IP)
190
              IF(ABS(H)+G.EQ.ABS(H))THEN
191
                 T=A(IP,IQ)/H
192
              ELSE
193
                 THETA=0.5*H/A(IP,IQ)
194
                 T=1./(ABS(THETA)+SQRT(1.+THETA**2))
195
                 IF(THETA.LT.0.)T=-T
196
              ENDIF
197
              C=1./SQRT(1+T**2)
198
              S=T*C
199
              TAU=S/(1.+C)
200
              H=T*A(IP,IQ)
201
              Z(IP)=Z(IP)-H
202
              Z(IQ)=Z(IQ)+H
203
              D(IP)=D(IP)-H
204
              D(IQ)=D(IQ)+H
205
              A(IP,IQ)=0.
206
              DO J=1,IP-1
207
                 G=A(J,IP)
208
                 H=A(J,IQ)
209
                 A(J,IP)=G-S*(H+G*TAU)
210
                 A(J,IQ)=H+S*(G-H*TAU)
211
              END DO
212
              DO  J=IP+1,IQ-1
213
                 G=A(IP,J)
214
                 H=A(J,IQ)
215
                 A(IP,J)=G-S*(H+G*TAU)
216
                 A(J,IQ)=H+S*(G-H*TAU)
217
              END DO
218
              DO  J=IQ+1,N
219
                 G=A(IP,J)
220
                 H=A(IQ,J)
221
                 A(IP,J)=G-S*(H+G*TAU)
222
                 A(IQ,J)=H+S*(G-H*TAU)
223
              END DO
224
              DO  J=1,N
225
                 G=V(J,IP)
226
                 H=V(J,IQ)
227
                 V(J,IP)=G-S*(H+G*TAU)
228
                 V(J,IQ)=H+S*(G-H*TAU)
229
              END DO
230
              NROT=NROT+1
231
           ENDIF
232
        END DO
233
     END DO
234
     DO  IP=1,N
235
        B(IP)=B(IP)+Z(IP)
236
        D(IP)=B(IP)
237
        Z(IP)=0.
238
     END DO
239
  END DO
240
  write(6,*) max_it,' iterations should never happen'
241
  STOP
242
100 DEALLOCATE(B,Z)
243
  RETURN
244
END SUBROUTINE JACOBI
245
!
246
!============================================================
247
!c
248
!c ++  Sort Eigenvectors in ascending eigenvalues order
249
!c
250
!c============================================================
251
!c
252
SUBROUTINE trie(nb_niv,ene,psi,nreal)
253

    
254
  !
255
  ! Input:
256
  !   Nb_niv      :  Dimension of Psi
257
  !   NReal  : Actual dimensions of Psi, Ene
258
  ! Input/output:
259
  !   Psi(N,N) : Eigvenvectors, changed in output. Stored in a (Nreal,Nreal) matrix
260
  !   Ene(N)   : Eigenvalues, changed in output. Stored in Ene(NReal)
261
  !   
262
!!!!!!!!!!!!!!!!
263

    
264
  Use VarTypes
265
  IMPLICIT NONE
266

    
267
  integer(KINT) :: i, j, k, nb_niv, nreal
268
  real(KREAL) :: ene(nreal),psi(nreal,nreal)
269
  real(KREAL) :: a
270

    
271

    
272
!!!!!!!!!!!!!!!!
273

    
274

    
275
  DO i=1,nb_niv
276
     DO j=i+1,nb_niv
277
        IF (ene(i) .GT. ene(j)) THEN
278
           !              permutation
279
           a=ene(i)
280
           ene(i)=ene(j)
281
           ene(j)=a
282

    
283
           DO k=1,nb_niv
284
              a=psi(k,i)
285
              psi(k,i)=psi(k,j)
286
              psi(k,j)=a
287
           END DO
288
        END IF
289
     END DO
290
  END DO
291

    
292
END SUBROUTINE trie
293

    
294
!============================================================
295
!c
296
!c ++  Sort Eigenvectors in ascending eigenvalues order
297
!c
298
!c============================================================
299
!c
300
SUBROUTINE SortEigenSys(N,Eigval,Eigvec)
301

    
302

    
303
  !
304
  ! Input/output:
305
  !   N : dimension of the system
306
  !   EigVec(N,N) : Eigvenvectors, changed in output. Stored in a (N,N) matrix
307
  !   EigVal(N)   : Eigenvalues, changed in output. Stored in a (n) vector
308
  !   
309
 ! Process: 
310
 ! We use first a ranking, then a working array to reorder the eigenvalues and eigenvectors
311
 
312
!!!!!!!!!!!!!!!!
313

    
314
  Use VarTypes
315
  use m_mrgrnk
316
  IMPLICIT NONE
317

    
318

    
319
  INTEGER(KINT), INTENT(IN) :: N
320
  REAL(KREAL), INTENT(OUT) :: EigVal(N), Eigvec(N,N)
321

    
322
  integer(KINT) :: i
323

    
324
  INTEGER(KINT), ALLOCATABLE :: Rank(:) !N
325
  REAL(KREAL), ALLOCATABLE :: EigValT(:) !n
326
  REAL(KREAL), ALLOCATABLE :: EigVecT(:,:) !n,n
327

    
328
!!!!!!!!!!!!!!!!
329
  ALLOCATE(Rank(n),EigValT(n),EigVecT(n,n))
330
  CALL MrgRnk(EigVal,Rank)
331

    
332
  DO I=1,N
333
     EigValT(I)=EigVal(Rank(I))
334
     EigVecT(I,1:N)=EigVec(Rank(I),1:N)
335
  END DO
336
  EigVal=EigValT
337
  EigVec=EigVecT
338

    
339
  DEALLOCATE(Rank,EigValT,EigVecT)
340

    
341

    
342
END SUBROUTINE SortEigenSys
343