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 nonzero 
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=1e12 
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, j1 
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 + (j1)*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 + (j1)*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,N1 
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,N1 
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,IP1 
207 
G=A(J,IP) 
208 
H=A(J,IQ) 
209 
A(J,IP)=GS*(H+G*TAU) 
210 
A(J,IQ)=H+S*(GH*TAU) 
211 
END DO 
212 
DO J=IP+1,IQ1 
213 
G=A(IP,J) 
214 
H=A(J,IQ) 
215 
A(IP,J)=GS*(H+G*TAU) 
216 
A(J,IQ)=H+S*(GH*TAU) 
217 
END DO 
218 
DO J=IQ+1,N 
219 
G=A(IP,J) 
220 
H=A(IQ,J) 
221 
A(IP,J)=GS*(H+G*TAU) 
222 
A(IQ,J)=H+S*(GH*TAU) 
223 
END DO 
224 
DO J=1,N 
225 
G=V(J,IP) 
226 
H=V(J,IQ) 
227 
V(J,IP)=GS*(H+G*TAU) 
228 
V(J,IQ)=H+S*(GH*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 