1* 2* $Id$ 3* 4 5* *********************************** 6* * * 7* * geodesic2_init * 8* * * 9* *********************************** 10* 11* Uses - geodesic common block 12* 13 14 subroutine geodesic2_init() 15 implicit none 16#include "errquit.fh" 17 18#include "bafdecls.fh" 19#include "geodesic2.fh" 20 21* **** local variables **** 22 logical value 23 integer npack1 24 25* **** external functions **** 26 integer psi_ne,psi_neq 27 external psi_ne,psi_neq 28 logical Dneall_m_allocate,Dneall_4m_allocate 29 external Dneall_m_allocate,Dneall_4m_allocate 30 31 call Pack_npack(1,npack1) 32c nemax = psi_neq(1)+psi_neq(2) 33c ne(1) = psi_neq(1) 34c ne(2) = psi_neq(2) 35 36 value = BA_alloc_get(mt_dcpl,npack1*(psi_neq(1)+psi_neq(2)), 37 > 'Hold',Hold(2),Hold(1)) 38 value = value.and. 39 > BA_alloc_get(mt_dcpl,npack1*(psi_neq(1)+psi_neq(2)), 40 > 'Q',Q(2),Q(1)) 41 value = value.and.Dneall_m_allocate(0,R) 42 value = value.and.Dneall_m_allocate(0,A) 43 value = value.and.Dneall_4m_allocate(0,V) 44 value = value.and.Dneall_4m_allocate(0,W) 45 46 value = value.and. 47 > BA_alloc_get(mt_dbl,2*(psi_ne(1)+psi_ne(2)), 48 > 'S',S(2),S(1)) 49 if (.not. value) call errquit('geodesic2_init:allocating heap',0, 50 & MA_ERR) 51 52 return 53 end 54 55* *********************************** 56* * * 57* * geodesic2_finalize * 58* * * 59* *********************************** 60* 61* Uses - geodesic common block 62* 63 64 subroutine geodesic2_finalize() 65 implicit none 66#include "errquit.fh" 67 68#include "bafdecls.fh" 69#include "geodesic2.fh" 70 71* **** local variables **** 72 logical value 73 74 logical Dneall_m_free 75 external Dneall_m_free 76 77 value = BA_free_heap(S(2)) 78 value = value.and.Dneall_m_free(W) 79 value = value.and.Dneall_m_free(V) 80 value = value.and.Dneall_m_free(A) 81 value = value.and.Dneall_m_free(R) 82 value = value.and.BA_free_heap(Q(2)) 83 value = value.and.BA_free_heap(Hold(2)) 84 if (.not. value) 85 > call errquit('geodesic2_finalize:freeing heap memory',0, 86 & MA_ERR) 87 88 return 89 end 90 91 92 93* *********************************** 94* * * 95* * geodesic2_start * 96* * * 97* *********************************** 98* 99* This routine determines the pxp matrices R and YA, and 100* the orthogonal nxp matrix Q. Q and R are determined from 101* the QR decomposition of the projected direction (I-YY^t)H, and 102* YH is defined as the Lagrange Multiplier pxp matrix Y^tH. 103* 104* Uses - geodesic2 common block 105* 106 subroutine geodesic2_start(Y,H,max_sigma,dE) 107 implicit none 108 complex*16 Y(*) 109 complex*16 H(*) 110 real*8 max_sigma,dE 111 112#include "bafdecls.fh" 113#include "errquit.fh" 114#include "geodesic2.fh" 115 116* **** local variables **** 117 logical value 118 integer npack1,nemax 119 integer ms,n,ispin,ne(2) 120 integer shift,shift2,i,j 121 integer T(2) 122 123* **** external functions **** 124 logical Dneall_4m_allocate,Dneall_m_free 125 integer psi_ispin,psi_neq 126 real*8 electron_eorbit_noocc 127 external Dneall_4m_allocate,Dneall_m_free 128 external psi_ispin,psi_neq 129 external electron_eorbit_noocc 130 131 call nwpw_timing_start(10) 132 call Pack_npack(1,npack1) 133 nemax = psi_neq(1) + psi_neq(2) 134 ispin = psi_ispin() 135 ne(1) = psi_neq(1) 136 ne(2) = psi_neq(2) 137 138* **** allocate tmp space **** 139 value = Dneall_4m_allocate(0,T) 140 if (.not. value) 141 > call errquit('geodesic2_start: pushing stack',0, MA_ERR) 142 143* **** Hold <-- H **** 144 !call dcopy(2*npack1*nemax,H,1,dcpl_mb(Hold(1)),1) 145 call Parallel_shared_vector_copy(.true.,2*npack1*nemax, 146 > H,dcpl_mb(Hold(1))) 147 148* **** calculate A=<Y|H> **** 149 call Dneall_ffm_Multiply(0,Y,H,npack1,dbl_mb(A(1))) 150 151* **** calculate Q=(I-YYt)H - should not be necessary but just in case **** 152 call Dneall_fmf_Multiply(0,Y,npack1, 153 > dbl_mb(A(1)),1.0d0, 154 > dcpl_mb(Q(1)),0.0d0) 155 call DAXPY_OMP(2*npack1*nemax,(-1.0d0),H,1,dcpl_mb(Q(1)),1) 156 call DSCAL_OMP(2*npack1*nemax,(-1.0d0),dcpl_mb(Q(1)),1) 157 158* **** calculate QR using Modified Gram-Schmidt **** 159 call Dneall_fm_QR(0,dcpl_mb(Q(1)),npack1,dbl_mb(R(1))) 160 161* **** generate T from A and R **** 162* - - 163* T = |A, -R^t| 164* |R, 0 | 165* - - 166 call Dneall_AR_to_4m(0,dbl_mb(A(1)),dbl_mb(R(1)),dbl_mb(T(1))) 167 168* **** Factor T--> V,W,and S **** 169 call Dneall_4m_FactorSkew(0, dbl_mb(T(1)), 170 > dbl_mb(V(1)), 171 > dbl_mb(W(1)), 172 > dbl_mb(S(1))) 173 174* **** calculate dE **** 175 dE = 2.0d0*electron_eorbit_noocc(H) 176 177 178* **** deallocate tmp space **** 179 value = Dneall_m_free(T) 180 if (.not. value) 181 > call errquit('geodesic2_start:popping stack',1, MA_ERR) 182 183 call nwpw_timing_end(10) 184 return 185 end 186 187* *********************************** 188* * * 189* * geodesic2_generate_T * 190* * * 191* *********************************** 192* 193* This routine determines T. T is defined 194* to be a 2nx2n skew symmetric matrix. 195* 196* - - 197* T = |A, -R^t| 198* |R, 0 | 199* - - 200* 201* Entry - n: dimension of matrices A and R 202* A: an nxn skew symmetric matrix 203* R: an nxn matrix 204* Exit - T: a 2nx2n skew symmetric matrix 205 206 subroutine geodesic2_generate_T(n,A,R,T) 207 implicit none 208 integer n 209 real*8 A(n,n) 210 real*8 R(n,n) 211 real*8 T(2*n,2*n) 212 213* **** local variables **** 214 integer i,j 215 216 call dcopy(4*n*n,0.0d0,0,T,1) 217 218* **** copy A to upper-left of T **** 219 do j=1,n 220 do i=1,n 221 T(i,j) = A(i,j) 222 end do 223 end do 224 225* **** copy R to lower-left of T **** 226 do j=1,n 227 do i=1,n 228 T(i+n,j) = R(i,j) 229 end do 230 end do 231 232* **** copy -R^t to upper-right of T **** 233 do j=1,n 234 do i=1,n 235 T(i,j+n) = -R(j,i) 236 end do 237 end do 238 239 return 240 end 241 242 243* *********************************** 244* * * 245* * geodesic2_get * 246* * * 247* *********************************** 248* 249* This routine calculates 250* 251* Ynew = Yold*M(t) + Q*N(t) 252* 253* where 254* - - - - 255* | M(t) | = Exp(t*T)* | I | 256* | N(t) | | 0 | 257* - - - - 258* 259 subroutine geodesic2_get(t,Yold,Ynew) 260 implicit none 261 real*8 t 262 complex*16 Yold(*) 263 complex*16 Ynew(*) 264 265#include "bafdecls.fh" 266#include "errquit.fh" 267#include "geodesic2.fh" 268 269* **** local variables **** 270 logical value 271 integer npack1,ispin,ne(2) 272 integer ms,n,j,k,shift1,shift2 273 integer MM(2),NN(2) 274 275* **** external functions **** 276 logical Dneall_m_push_get,Dneall_m_pop_stack 277 external Dneall_m_push_get,Dneall_m_pop_stack 278 279 280 call Pack_npack(1,npack1) 281 282* **** allocate tmp space **** 283 value = Dneall_m_push_get(0,MM) 284 value = value.and.Dneall_m_push_get(0,NN) 285 if (.not. value) 286 > call errquit('geodesic2_get: pushing stack',0, MA_ERR) 287 288 call geodesic2_get_MandN(t,dbl_mb(MM(1)),dbl_mb(NN(1))) 289 290 call Dneall_fmf_Multiply(0, 291 > Yold,npack1, 292 > dbl_mb(MM(1)), 293 > 1.0d0, 294 > Ynew, 295 > 0.0d0) 296 call Dneall_fmf_Multiply(0, 297 > dcpl_mb(Q(1)),npack1, 298 > dbl_mb(NN(1)), 299 > 1.0d0, 300 > Ynew, 301 > 1.0d0) 302 303* **** deallocate tmp space **** 304 value = Dneall_m_pop_stack(NN) 305 value = value.and.Dneall_m_pop_stack(MM) 306 if (.not. value) 307 > call errquit('geodesic2_get:popping stack',1, MA_ERR) 308 309 return 310 end 311 312 313* *********************************** 314* * * 315* * geodesic2_transport * 316* * * 317* *********************************** 318* 319* This routine calculates 320* 321* Hnew = Hold*M(t) + Yold*R^t*N(t) 322* 323* where 324* - - - - 325* | M(t) | = Exp(t*T)* | I | 326* | N(t) | | 0 | 327* - - - - 328* 329 subroutine geodesic2_transport(t,Yold,Hnew) 330 implicit none 331#include "errquit.fh" 332 real*8 t 333 complex*16 Yold(*) 334 complex*16 Hnew(*) 335 336#include "bafdecls.fh" 337#include "geodesic2.fh" 338 339 340* **** local variables **** 341 logical value 342 integer npack1 343 integer MM(2),NN(2),TT(2) 344 345* **** external functions **** 346 logical Dneall_m_push_get,Dneall_m_pop_stack 347 external Dneall_m_push_get,Dneall_m_pop_stack 348 349 350 call Pack_npack(1,npack1) 351 352* **** allocate tmp space **** 353 value = Dneall_m_push_get(0,MM) 354 value = value.and.Dneall_m_push_get(0,NN) 355 value = value.and.Dneall_m_push_get(0,TT) 356 if (.not. value) 357 > call errquit('geodesic2_transport: pushing stack',0, MA_ERR) 358 359 call geodesic2_get_MandN(t,dbl_mb(MM(1)),dbl_mb(NN(1))) 360 361* **** TT(t) = -R^t*NN(t) **** 362 call Dneall_mmm_Multiply2ab(0, 363 > dbl_mb(R(1)), 364 > dbl_mb(NN(1)), 365 > (-1.0d0), 366 > dbl_mb(TT(1)), 367 > 0.0d0) 368 369 370* *** Hnew <-- Hold*M(t) + Yold*TT(t) *** 371 call Dneall_fmf_Multiply(0, 372 > dcpl_mb(Hold(1)),npack1, 373 > dbl_mb(MM(1)), 374 > 1.0d0, 375 > Hnew, 376 > 0.0d0) 377 call Dneall_fmf_Multiply(0, 378 > Yold,npack1, 379 > dbl_mb(TT(1)), 380 > 1.0d0, 381 > Hnew, 382 > 1.0d0) 383 384* **** deallocate tmp space **** 385 value = Dneall_m_pop_stack(TT) 386 value = value.and.Dneall_m_pop_stack(NN) 387 value = value.and.Dneall_m_pop_stack(MM) 388 if (.not. value) 389 > call errquit('geodesic2_transport:popping stack',1, MA_ERR) 390 391 392 return 393 end 394 395 subroutine geodesic2_checkMN(name,n,MM,NN) 396 implicit none 397 character*(*) name 398 integer n 399 real*8 MM(n,n),NN(n,n) 400 401* **** local variables **** 402 integer i,j,k 403 real*8 tmp(n,n),sum 404 405 do j=1,n 406 do i=1,n 407 sum = 0.0d0 408 do k=1,n 409 sum = sum + MM(k,i)*NN(k,j) 410 end do 411 tmp(i,j) = sum 412 end do 413 end do 414 415 write(*,*) 416 write(*,*) "checkMN:",name 417 do i=1,n 418 write(*,'(22F8.3)') (tmp(i,j),j=1,n) 419 end do 420 write(*,*) 421 return 422 end 423 424* *********************************** 425* * * 426* * geodesic2_get_MandN * 427* * * 428* *********************************** 429* 430* This routine returns 431* - - - - 432* | M(t) | = Exp(t*T)* | I | 433* | N(t) | | 0 | 434* - - - - 435* where 436* 437* T = U*Sigma*U^H, with U=(V+iW) 438* 439* is a skew matrix that is decomposed into V,W,and Sigma 440* 441 subroutine geodesic2_get_MandN(t,M,N) 442 implicit none 443 real*8 t 444 real*8 M(*),N(*) 445 446#include "bafdecls.fh" 447#include "errquit.fh" 448#include "geodesic2.fh" 449 450 451* **** local variables **** 452 logical value 453 integer ispin,ne(2) 454 integer ms,i,j,nn,shift,shift2,shift3 455 integer RR(2) 456 457* **** external functions **** 458 logical Dneall_4m_push_get,Dneall_m_pop_stack 459 external Dneall_4m_push_get,Dneall_m_pop_stack 460 461* **** allocate tmp space **** 462 value = Dneall_4m_push_get(0,RR) 463 if (.not. value) 464 > call errquit('geodesic2_get_MandN: pushing stack',0, MA_ERR) 465 466 call Dneall_4m_RotateSkew(0,t, 467 > dbl_mb(V(1)), 468 > dbl_mb(W(1)), 469 > dbl_mb(S(1)), 470 > dbl_mb(RR(1)) ) 471 472 call Dneall_4m_to_MN(0,dbl_mb(RR(1)),M,N) 473 474* **** deallocate tmp space **** 475 value = Dneall_m_pop_stack(RR) 476 if (.not. value) 477 > call errquit('geodesic2_get_MandN:popping stack',1, MA_ERR) 478 479 return 480 end 481