root / src / Mat_util.f90 @ 4
Historique | Voir | Annoter | Télécharger (8,25 ko)
| 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 |
INTEGER(KINT), PARAMETER :: Max_it=500 |
| 84 |
REAL(KREAL) :: AP(N*(N+1)/2) |
| 85 |
REAL(KREAL), ALLOCATABLE :: A(:,:) |
| 86 |
REAL(KREAL) :: V(Nreal,Nreal),D(Nreal) |
| 87 |
INTEGER(KINT) :: i,j,nn |
| 88 |
|
| 89 |
allocate(A(nreal,nreal)) |
| 90 |
nn=n*(n+1)/2 |
| 91 |
! WRITE(*,*) 'Jacpa 0' |
| 92 |
! WRITE(*,'(I3,10(1X,F15.6))') n,(AP(i),i=1,min(nn,5)) |
| 93 |
! WRITE(*,*) 'Jacpa 0' |
| 94 |
do j=1,n |
| 95 |
do i=1,j |
| 96 |
! WRITE(*,*) i,j |
| 97 |
A(i,J)=AP(i + (j-1)*j/2) |
| 98 |
A(j,I)=A(i,J) |
| 99 |
end do |
| 100 |
end do |
| 101 |
! do j=1,n |
| 102 |
! WRITE(*,'(10(1X,F15.6))') (A(i,J),i=1,min(n,5)) |
| 103 |
! end do |
| 104 |
|
| 105 |
! WRITE(*,*) 'Jacpa 1' |
| 106 |
call Jacobi(A,n,D,V,Nreal) |
| 107 |
! WRITE(*,*) 'Jacpa 2' |
| 108 |
! DO i=1,n |
| 109 |
! WRITE(*,'(1X,I3,10(1X,F15.6))') i,D(i),(V(j,i),j=1,min(n,5)) |
| 110 |
! end do |
| 111 |
deallocate(a) |
| 112 |
|
| 113 |
end SUBROUTINE JacPacked |
| 114 |
|
| 115 |
|
| 116 |
|
| 117 |
! |
| 118 |
!============================================================ |
| 119 |
! |
| 120 |
! ++ Matrix diagonalization Using jacobi |
| 121 |
! Works only for symmetric matrices |
| 122 |
! EigvenVectors : V |
| 123 |
! Eigenvalues : D |
| 124 |
! |
| 125 |
!============================================================ |
| 126 |
! |
| 127 |
SUBROUTINE JACOBI(A,N,D,V,Nreal) |
| 128 |
|
| 129 |
!!!!!!!!!!!!!!!! |
| 130 |
! |
| 131 |
! Input: |
| 132 |
! N : Dimension of A |
| 133 |
! NReal : Actual dimensions of A, D and V. |
| 134 |
! |
| 135 |
! Input/output: |
| 136 |
! A(N,N) : Matrix to be diagonalized, store in a (Nreal,Nreal) matrix |
| 137 |
! Destroyed in output. |
| 138 |
! Output: |
| 139 |
! V(N,N) : Eigenvectors, stored in V(NReal, NReal) |
| 140 |
! D(N) : Eigenvalues, stored in D(NReal) |
| 141 |
! |
| 142 |
|
| 143 |
Use Vartypes |
| 144 |
|
| 145 |
IMPLICIT NONE |
| 146 |
INTEGER(KINT), parameter :: max_it=500 |
| 147 |
REAL(KREAL), ALLOCATABLE :: B(:),Z(:) |
| 148 |
|
| 149 |
INTEGER(KINT) :: N,NReal |
| 150 |
REAL(KREAL) :: A(NReal,NReal) |
| 151 |
REAL(KREAL) :: V(Nreal,Nreal),D(Nreal) |
| 152 |
|
| 153 |
INTEGER(KINT) :: I, J,IP, IQ, NROT |
| 154 |
REAL(KREAL) :: SM, H, Tresh, G, T, Theta, C, S, Tau |
| 155 |
|
| 156 |
allocate(B(N),Z(N)) |
| 157 |
|
| 158 |
DO IP=1,N |
| 159 |
DO IQ=1,N |
| 160 |
V(IP,IQ)=0. |
| 161 |
END DO |
| 162 |
V(IP,IP)=1. |
| 163 |
END DO |
| 164 |
DO IP=1,N |
| 165 |
B(IP)=A(IP,IP) |
| 166 |
D(IP)=B(IP) |
| 167 |
Z(IP)=0. |
| 168 |
END DO |
| 169 |
NROT=0 |
| 170 |
DO I=1,max_it |
| 171 |
SM=0. |
| 172 |
DO IP=1,N-1 |
| 173 |
DO IQ=IP+1,N |
| 174 |
SM=SM+ABS(A(IP,IQ)) |
| 175 |
END DO |
| 176 |
END DO |
| 177 |
IF(SM.EQ.0.) GOTO 100 |
| 178 |
IF(I.LT.4)THEN |
| 179 |
TRESH=0.2*SM/N**2 |
| 180 |
ELSE |
| 181 |
TRESH=0. |
| 182 |
ENDIF |
| 183 |
DO IP=1,N-1 |
| 184 |
DO IQ=IP+1,N |
| 185 |
G=100.*ABS(A(IP,IQ)) |
| 186 |
IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) & |
| 187 |
.AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN |
| 188 |
A(IP,IQ)=0. |
| 189 |
ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN |
| 190 |
H=D(IQ)-D(IP) |
| 191 |
IF(ABS(H)+G.EQ.ABS(H))THEN |
| 192 |
T=A(IP,IQ)/H |
| 193 |
ELSE |
| 194 |
THETA=0.5*H/A(IP,IQ) |
| 195 |
T=1./(ABS(THETA)+SQRT(1.+THETA**2)) |
| 196 |
IF(THETA.LT.0.)T=-T |
| 197 |
ENDIF |
| 198 |
C=1./SQRT(1+T**2) |
| 199 |
S=T*C |
| 200 |
TAU=S/(1.+C) |
| 201 |
H=T*A(IP,IQ) |
| 202 |
Z(IP)=Z(IP)-H |
| 203 |
Z(IQ)=Z(IQ)+H |
| 204 |
D(IP)=D(IP)-H |
| 205 |
D(IQ)=D(IQ)+H |
| 206 |
A(IP,IQ)=0. |
| 207 |
DO J=1,IP-1 |
| 208 |
G=A(J,IP) |
| 209 |
H=A(J,IQ) |
| 210 |
A(J,IP)=G-S*(H+G*TAU) |
| 211 |
A(J,IQ)=H+S*(G-H*TAU) |
| 212 |
END DO |
| 213 |
DO J=IP+1,IQ-1 |
| 214 |
G=A(IP,J) |
| 215 |
H=A(J,IQ) |
| 216 |
A(IP,J)=G-S*(H+G*TAU) |
| 217 |
A(J,IQ)=H+S*(G-H*TAU) |
| 218 |
END DO |
| 219 |
DO J=IQ+1,N |
| 220 |
G=A(IP,J) |
| 221 |
H=A(IQ,J) |
| 222 |
A(IP,J)=G-S*(H+G*TAU) |
| 223 |
A(IQ,J)=H+S*(G-H*TAU) |
| 224 |
END DO |
| 225 |
DO J=1,N |
| 226 |
G=V(J,IP) |
| 227 |
H=V(J,IQ) |
| 228 |
V(J,IP)=G-S*(H+G*TAU) |
| 229 |
V(J,IQ)=H+S*(G-H*TAU) |
| 230 |
END DO |
| 231 |
NROT=NROT+1 |
| 232 |
ENDIF |
| 233 |
END DO |
| 234 |
END DO |
| 235 |
DO IP=1,N |
| 236 |
B(IP)=B(IP)+Z(IP) |
| 237 |
D(IP)=B(IP) |
| 238 |
Z(IP)=0. |
| 239 |
END DO |
| 240 |
END DO |
| 241 |
write(6,*) max_it,' iterations should never happen' |
| 242 |
STOP |
| 243 |
100 DEALLOCATE(B,Z) |
| 244 |
RETURN |
| 245 |
END SUBROUTINE JACOBI |
| 246 |
! |
| 247 |
!============================================================ |
| 248 |
!c |
| 249 |
!c ++ Sort Eigenvectors in ascending eigenvalues order |
| 250 |
!c |
| 251 |
!c============================================================ |
| 252 |
!c |
| 253 |
SUBROUTINE trie(nb_niv,ene,psi,nreal) |
| 254 |
|
| 255 |
! |
| 256 |
! Input: |
| 257 |
! Nb_niv : Dimension of Psi |
| 258 |
! NReal : Actual dimensions of Psi, Ene |
| 259 |
! Input/output: |
| 260 |
! Psi(N,N) : Eigvenvectors, changed in output. Stored in a (Nreal,Nreal) matrix |
| 261 |
! Ene(N) : Eigenvalues, changed in output. Stored in Ene(NReal) |
| 262 |
! |
| 263 |
!!!!!!!!!!!!!!!! |
| 264 |
|
| 265 |
Use VarTypes |
| 266 |
IMPLICIT NONE |
| 267 |
|
| 268 |
integer(KINT) :: i,j,k,nb_niv,max_niv, nreal |
| 269 |
real(KREAL) :: ene(nreal),psi(nreal,nreal) |
| 270 |
real(KREAL) :: a |
| 271 |
|
| 272 |
|
| 273 |
!!!!!!!!!!!!!!!! |
| 274 |
|
| 275 |
|
| 276 |
DO i=1,nb_niv |
| 277 |
DO j=i+1,nb_niv |
| 278 |
IF (ene(i) .GT. ene(j)) THEN |
| 279 |
! permutation |
| 280 |
a=ene(i) |
| 281 |
ene(i)=ene(j) |
| 282 |
ene(j)=a |
| 283 |
|
| 284 |
DO k=1,nb_niv |
| 285 |
a=psi(k,i) |
| 286 |
psi(k,i)=psi(k,j) |
| 287 |
psi(k,j)=a |
| 288 |
END DO |
| 289 |
END IF |
| 290 |
END DO |
| 291 |
END DO |
| 292 |
|
| 293 |
END SUBROUTINE trie |
| 294 |
|
| 295 |
!============================================================ |
| 296 |
!c |
| 297 |
!c ++ Sort Eigenvectors in ascending eigenvalues order |
| 298 |
!c |
| 299 |
!c============================================================ |
| 300 |
!c |
| 301 |
SUBROUTINE SortEigenSys(N,Eigval,Eigvec) |
| 302 |
|
| 303 |
|
| 304 |
! |
| 305 |
! Input/output: |
| 306 |
! N : dimension of the system |
| 307 |
! EigVec(N,N) : Eigvenvectors, changed in output. Stored in a (N,N) matrix |
| 308 |
! EigVal(N) : Eigenvalues, changed in output. Stored in a (n) vector |
| 309 |
! |
| 310 |
! Process: |
| 311 |
! We use first a ranking, then a working array to reorder the eigenvalues and eigenvectors |
| 312 |
|
| 313 |
!!!!!!!!!!!!!!!! |
| 314 |
|
| 315 |
Use VarTypes |
| 316 |
use m_mrgrnk |
| 317 |
IMPLICIT NONE |
| 318 |
|
| 319 |
|
| 320 |
INTEGER(KINT), INTENT(IN) :: N |
| 321 |
REAL(KREAL), INTENT(OUT) :: EigVal(N), Eigvec(N,N) |
| 322 |
|
| 323 |
integer(KINT) :: i,j,k |
| 324 |
|
| 325 |
INTEGER(KINT), ALLOCATABLE :: Rank(:) !N |
| 326 |
REAL(KREAL), ALLOCATABLE :: EigValT(:) !n |
| 327 |
REAL(KREAL), ALLOCATABLE :: EigVecT(:,:) !n,n |
| 328 |
|
| 329 |
!!!!!!!!!!!!!!!! |
| 330 |
ALLOCATE(Rank(n),EigValT(n),EigVecT(n,n)) |
| 331 |
CALL MrgRnk(EigVal,Rank) |
| 332 |
|
| 333 |
DO I=1,N |
| 334 |
EigValT(I)=EigVal(Rank(I)) |
| 335 |
EigVecT(I,1:N)=EigVec(Rank(I),1:N) |
| 336 |
END DO |
| 337 |
EigVal=EigValT |
| 338 |
EigVec=EigVecT |
| 339 |
|
| 340 |
DEALLOCATE(Rank,EigValT,EigVecT) |
| 341 |
|
| 342 |
|
| 343 |
END SUBROUTINE SortEigenSys |
| 344 |
|