1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8! This code segment has been fully created by: 9! Nick Papior Andersen, 2013, nickpapior@gmail.com 10! Please conctact the author, prior to re-using this code. 11 12! This particular solution method relies on solving the GF 13! with the tri-diagonalization routine. 14! This will leverage memory usage and also the execution time. 15! It is customized to deal with calculating the inverted matrix 16! only in a certain region. 17 18module m_ts_trimat_invert 19 20 ! We use the general inversion module to obtain the 21 ! same routines etc. 22 23 use precision, only : dp 24 use class_zTrimat 25 use m_trimat_invert, only : calc_Mnn_inv 26 use m_trimat_invert, only : calc_Xn_div_Cn_p1, calc_Yn_div_Bn_m1 27 use m_trimat_invert, only : Xn_div_Cn_p1, Yn_div_Bn_m1 28 29 use m_pivot_array, only : Npiv, ipiv, init_pivot 30 31 implicit none 32 33 private 34 35 ! Used for BLAS calls (local variables) 36 complex(dp), parameter :: z0 = dcmplx( 0._dp, 0._dp) 37 complex(dp), parameter :: z1 = dcmplx( 1._dp, 0._dp) 38 complex(dp), parameter :: zm1 = dcmplx(-1._dp, 0._dp) 39 40 public :: invert_BiasTriMat_prep 41#ifdef TBTRANS 42 public :: BiasTriMat_prep 43#endif 44 public :: invert_BiasTriMat_col 45 public :: invert_BiasTriMat_rgn 46 public :: TriMat_Bias_idxs 47 48contains 49 50 subroutine invert_BiasTriMat_prep(M,Minv, part_cols, all_nn) 51 52 type(zTriMat), intent(inout) :: M, Minv 53 ! This is a 3x... list of 54 ! 1 == part 55 ! 2 == smallest column in the part 56 ! 3 == largest column in the part 57 integer, intent(in), optional :: part_cols(:,:) 58 logical, intent(in), optional :: all_nn 59 60 complex(dp), pointer :: Mpinv(:) 61 62 integer :: n, np, sNm1, sNp1 63 integer :: sCol, eCol, i 64 logical :: piv_initialized 65 66 np = parts(M) 67#ifndef TS_NOCHECKS 68 if ( np /= parts(Minv) ) then 69 call die('Could not calculate the inverse on non equal sized & 70 &matrices') 71 end if 72 if ( np == 1 ) then 73 call die('This matrix is not tri-diagonal') 74 end if 75#endif 76 piv_initialized = .true. 77#ifndef TS_NOCHECKS 78 do n = 1 , np 79 if ( Npiv < nrows_g(M,n) ) piv_initialized = .false. 80 end do 81 if ( .not. piv_initialized ) then 82 call die('Pivoting array for inverting matrix not set.') 83 end if 84#endif 85 86 call timer('V_TM_Pinv',1) 87 88 ! Calculate all Xn/Cn+1 89 do n = np - 1 , 1 , -1 90 Mpinv => val(Minv,n+1,n+1) 91 sNp1 = nrows_g(M,n+1) 92 call calc_Xn_div_Cn_p1(M,Minv, n, Mpinv, sNp1**2 ) 93 end do 94 ! Calculate all Yn/Bn-1 95 do n = 2 , np 96 Mpinv => val(Minv,n-1,n-1) 97 sNm1 = nrows_g(M,n-1) 98 call calc_Yn_div_Bn_m1(M,Minv, n, Mpinv, sNm1**2 ) 99 end do 100 101 ! At this point the status is: 102 ! - M contains the original matrix 103 ! - Minv contains all Xn/Cn+1 and Yn/Bn-1 in all parts 104 105 piv_initialized = .false. 106 if ( present(all_nn) ) piv_initialized = all_nn 107 if ( .not. present(part_cols) ) piv_initialized = .true. 108 109 if ( piv_initialized ) then 110 111 ! We need to calculate all Mnn 112 do n = 1 , np 113 call calc_Mnn_inv(M,Minv,n) 114 end do 115 116 else if ( present(part_cols) ) then 117 118 ! This loops an array of parts to invert 119 do i = 1 , size(part_cols,2) 120 121 ! Note that these *must* be separate parts 122 n = part_cols(1,i) 123 sCol = part_cols(2,i) 124 eCol = part_cols(3,i) 125 126 call calc_Mnn_inv_cols(M,Minv,n,sCol,eCol) 127 128 end do 129 130 else 131 call die('inv-BTM-prep: Wrong options.') 132 end if 133 134 ! At this point we have calculated the 135 ! Mnn matrices for the overlapping regions for the 136 ! electrodes. 137 ! This is now saved in Minv. 138 ! Minv contains all information to calculate 139 ! all Gf columns 140 141 call timer('V_TM_Pinv',2) 142 143 end subroutine invert_BiasTriMat_prep 144 145#ifdef TBTRANS 146 subroutine BiasTrimat_prep(M,N_Elec,Elecs,has_El,part_cols) 147 148 use m_ts_electype 149 150 type(zTriMat), intent(inout) :: M 151 integer, intent(in) :: N_Elec 152 type(Elec), intent(in) :: Elecs(N_Elec) 153 logical, intent(in) :: has_El(N_Elec) 154 integer, intent(inout), allocatable :: part_cols(:,:) 155 156 integer :: iEl, off, sCol, eCol 157 integer :: n, np, sN, Ne, sNm1, sNp1 158 ! We need to allow the electrodes to be split 159 ! across two blocks. 160 integer :: tmp(3,N_Elec*2) 161 162 np = parts(M) 163 164 Ne = 0 165 off = 0 166 do n = 1 , np 167 ! Get part information 168 sN = nrows_g(M,n) 169 170 sCol = huge(1) 171 eCol = -1 172 173 do iEl = 1 , N_Elec 174 175 if ( .not. has_El(iEl) ) cycle 176 177 ! We know that any diagonal region with the 178 ! electrodes will only occupy at most 2 parts 179 ! Hence, it makes no sense to loop them individually 180 ! Also inDpvt is a sorted array with device indices 181 sNm1 = Elecs(iEl)%inDpvt%r(1) 182 sNp1 = Elecs(iEl)%inDpvt%r(Elecs(iEl)%inDpvt%n) 183 if ( which_part(M,sNm1) <= n .and. & 184 n <= which_part(M,sNp1) ) then 185 186 ! get number of columns that belongs to 187 ! the electrode in the 'n' diagonal part 188 ! this means we only calculate "what is needed" 189 sCol = min(sCol, max(sNm1 - off , 1) ) 190 eCol = max(eCol, min(sNp1 - off , sN) ) 191 192 end if 193 194 end do 195 196 if ( eCol /= -1 ) then 197 Ne = Ne + 1 198 tmp(1,Ne) = n 199 tmp(2,Ne) = sCol 200 tmp(3,Ne) = eCol 201 end if 202 203 off = off + sN 204 205 end do 206 207 allocate(part_cols(3,Ne)) 208 part_cols(:,:) = tmp(:,1:Ne) 209 210 end subroutine BiasTrimat_prep 211#endif 212 213 subroutine invert_BiasTriMat_col(M,Minv,El,calc_parts) 214 215 use m_ts_electype 216 use m_ts_method, only : orb_offset 217 218 type(zTriMat), intent(inout) :: M, Minv 219 type(Elec), intent(in) :: El 220 logical, intent(in) :: calc_parts(:) 221 222 complex(dp), pointer :: Mpinv(:), Mp(:) 223 complex(dp), pointer :: Xn(:), Yn(:) 224 complex(dp), pointer :: z(:) 225 226 integer :: nr, np, no 227 integer :: idx_o 228 integer :: sPart, ePart, lsPart, lePart 229 integer :: sColF, eColF, sIdxF, eIdxF 230 integer :: sColT, eColT, sIdxT, eIdxT 231 integer :: sN, sNc, sNm1, sNp1, n, s 232 integer :: off 233 logical :: piv_initialized 234 235 ! In this routine M should have been processed through invert_PrepTriMat 236 ! So M contains all *needed* inv(Mnn) and all Xn/Cn+1 and Yn/Bn-1. 237 ! So we will save the result in Minv 238#ifndef TS_NOCHECKS 239 if ( parts(M) /= parts(Minv) ) then 240 call die('Could not calculate the inverse on non equal sized & 241 &matrices') 242 end if 243 if ( parts(M) == 1 ) then 244 call die('This matrix is not tri-diagonal') 245 end if 246 piv_initialized = .true. 247 do n = 1 , parts(M) 248 if ( Npiv < nrows_g(M,n) ) piv_initialized = .false. 249 end do 250 if ( .not. piv_initialized ) then 251 call die('Pivoting array for inverting matrix not set.') 252 end if 253 254 if ( parts(M) /= size(calc_parts) ) then 255 call die('Error in code, calc_parts, not consistent') 256 end if 257#endif 258 259 call timer('V_TM_inv',1) 260 261 nr = nrows_g(M) 262 np = parts(M) 263 264 idx_o = El%idx_o - orb_offset(El%idx_o) 265 no = TotUsedOrbs(El) 266 267 sPart = which_part(M,idx_o) 268 ePart = which_part(M,idx_o+no-1) 269#ifndef TS_NOCHECKS 270 if ( sPart < 1 ) call die('Error in the Bias inversion, sPart') 271 if ( ePart - sPart + 1 > 2 ) call die('Error in trimat partition') 272 if ( ePart > parts(M) ) call die('Error in the Bias inversion, ePart') 273#endif 274 275 ! Point to the matrices 276 z => val(Minv,all=.true.) 277 278 ! First we need to copy over the Mnn with the electrode part! 279 280 ! with this offset we can calculate the column offset for 281 ! the current part 282 off = 0 283 do n = 1 , sPart - 1 284 off = off + nrows_g(M,n) 285 end do 286 idx_o = idx_o - off 287 if ( idx_o <= 0 ) call die('Error in electrode setup') 288 do n = sPart , ePart 289 290 ! current count of orbitals in the tri-diagonal segment 291 if ( n > 1 ) sNm1 = nrows_g(M,n-1) 292 sN = nrows_g(M,n ) 293 if ( n < np ) sNp1 = nrows_g(M,n+1) 294 295 ! placement of the already inverted matrix 296 Mp => val(M,n,n) 297 298 ! get number of columns that belongs to 299 ! the electrode in the 'n' diagonal part 300 sColF = max(idx_o , 1) 301 eColF = min(idx_o + no - 1 , sN) 302#ifndef TS_NOCHECKS 303 if ( eColF < sColF ) & 304 call die('Here: Something went wrong') 305#endif 306 sIdxF = (sColF-1) * sN + 1 307 eIdxF = eColF * sN 308 309 ! get placement of the diagonal block in the column 310 call TriMat_Bias_idxs(Minv,no,n,sIdxT,eIdxT) 311 Mpinv => z(sIdxT:eIdxT) 312 313 ! get the placement in the inversed column 314 if ( 1 <= idx_o ) then 315 ! we are taking the first part of the inversed matrix 316 sColT = 1 317 eColT = min(eColF-sColF+1,no) 318 else 319 sColT = -idx_o + 2 ! we have to pass zero 320 eColT = no 321 end if 322 sIdxT = (sColT-1) * sN + 1 323 eIdxT = eColT * sN 324 325#ifndef TS_NOCHECKS 326 if ( eIdxT - sIdxT /= eIdxF - sIdxF ) & 327 call die('Error in determining column') 328#endif 329 call zcopy(eIdxT-sIdxT+1,Mp(sIdxF),1,Mpinv(sIdxT),1) 330 331 if ( sPart == ePart ) cycle ! we have everything! :) 332 333 ! We need to calculate the remaining 334 ! inverted matrix (they currently reside 335 ! with the neighbouring cells) 336 337 ! first calculate the missing number of columns 338 sNc = no - (eColF - sColF + 1) 339 340 if ( n == sPart ) then 341 ! we miss the right part of Mnn 342 ! we thus need the 343 ! Mmn = -Ym+1/Bm * Mm+1n 344 345 Yn => Yn_div_Bn_m1(M,n+1) 346 347 ! placement of the inverted matrix 348 Mp => val(M,n+1,n+1) 349 350#ifdef USE_GEMM3M 351 call zgemm3m( & 352#else 353 call zgemm( & 354#endif 355 'N','N',sN,sNc,sNp1, & 356 zm1, Yn, sN, Mp(1),sNp1,z0, Mpinv(eIdxT+1),sN) 357 358 else if ( n == ePart ) then 359 ! we miss the left part of Mnn 360 ! we thus need the 361 ! Mmn = -Xm-1/Cm * Mm-1n 362 363 Xn => Xn_div_Cn_p1(M,n-1) 364 365 ! placement of the inverted matrix 366 Mp => val(M,n-1,n-1) 367 s = sNm1 ** 2 368 369#ifdef USE_GEMM3M 370 call zgemm3m( & 371#else 372 call zgemm( & 373#endif 374 'N','N',sN,sNc,sNm1, & 375 zm1, Xn, sN, Mp(s-sNm1*sNc+1),sNm1,z0, Mpinv(1),sN) 376 377 end if 378 379 ! update offset on rows 380 off = off + sN 381 idx_o = idx_o - sN 382 383 end do 384 385 ! We now have inv(Mnn) in the correct place. 386 ! figure out what we should calculate 387 ! We can not save a lot of computations here 388 ! as we can not be sure of the full extend of 389 ! the "end"-blocks. Only if two consecutive 390 ! blocks are fully encompassed in one electrode 391 ! can we reduce the computational work. 392 do n = 1 , np 393 if ( calc_parts(n) ) then 394 lsPart = max(1,n-1) 395 exit 396 end if 397 end do 398 do n = np , 1 , -1 399 if ( calc_parts(n) ) then 400 lePart = min(n+1,np) 401 exit 402 end if 403 end do 404 405 ! We now calculate: 406 ! Mmn = -Ym+1/Bm * Mm+1n, for m<n 407 do n = sPart - 1 , lsPart , - 1 408 409 sN = nrows_g(M,n) 410 sNp1 = nrows_g(M,n+1) 411 412 ! get Ym+1/Bm 413 Yn => Yn_div_Bn_m1(M,n+1) 414 415 ! Get Mm+1n 416 call TriMat_Bias_idxs(Minv,no,n+1,sIdxF,eIdxF) 417 Mp => z(sIdxF:eIdxF) 418 ! Get Mmn 419 call TriMat_Bias_idxs(Minv,no,n,sIdxF,eIdxF) 420 Mpinv => z(sIdxF:eIdxF) 421 422#ifdef USE_GEMM3M 423 call zgemm3m( & 424#else 425 call zgemm( & 426#endif 427 'N','N',sN,no,sNp1, & 428 zm1, Yn, sN, Mp(1),sNp1,z0, Mpinv(1),sN) 429 430 end do 431 432 ! We now calculate: 433 ! Mmn = -Xm-1/Cm * Mm-1n, for m>n 434 do n = ePart + 1 , lePart 435 436 sNm1 = nrows_g(M,n-1) 437 sN = nrows_g(M,n) 438 439 ! get Xm-1/Cm 440 Xn => Xn_div_Cn_p1(M,n-1) 441 442 ! Get Mm-1n 443 call TriMat_Bias_idxs(Minv,no,n-1,sIdxF,eIdxF) 444 Mp => z(sIdxF:eIdxF) 445 ! Get Mmn 446 call TriMat_Bias_idxs(Minv,no,n,sIdxF,eIdxF) 447 Mpinv => z(sIdxF:eIdxF) 448 449#ifdef USE_GEMM3M 450 call zgemm3m( & 451#else 452 call zgemm( & 453#endif 454 'N','N',sN,no,sNm1, & 455 zm1, Xn, sN, Mp(1),sNm1,z0, Mpinv(1),sN) 456 457 end do 458 459 ! At this point the total 460 ! inverted column is placed at the end of 461 ! the tri-mat inversion. 462 463 call timer('V_TM_inv',2) 464 465 end subroutine invert_BiasTriMat_col 466 467 subroutine invert_BiasTriMat_rgn(M,Minv,r,pvt,r_col,only_diag) 468 469 use m_region 470 471 type(zTriMat), intent(inout) :: M, Minv 472 ! The pivoting table for the tri-diagonal matrix 473 ! pvt is the pivoting back to original index (to remove 474 ! rgn_pivot calls, i.e. pvt%r(r%r) == (/1,...,no_u/) 475 type(tRgn), intent(in) :: r, pvt 476 ! The siesta indices for the column we wish to calculate 477 type(tRgn), intent(in) :: r_col 478 ! Whether we should only calculate the diagonal Green function 479 logical, intent(in), optional :: only_diag 480 481 complex(dp), pointer :: Mpinv(:), Mp(:) 482 complex(dp), pointer :: XY(:) 483 complex(dp), pointer :: z(:) 484 485 ! Used for BLAS calls (local variables) 486 complex(dp), parameter :: z0 = dcmplx( 0._dp, 0._dp) 487 complex(dp), parameter :: zm1 = dcmplx(-1._dp, 0._dp) 488 489 integer :: nr, np 490 integer :: sPart, ePart 491 integer :: i, i_Elec, idx_Elec 492 integer :: sIdxF, eIdxF, sIdxT, eIdxT 493 integer :: sN, n, in, s, sNo, nb 494 integer, pointer :: crows(:) 495 496 ! In this routine M should have been processed through invert_PrepTriMat 497 ! So M contains all *needed* inv(Mnn) and all Xn/Cn+1 and Yn/Bn-1. 498 ! So we will save the result in Minv 499#ifndef TS_NOCHECKS 500 if ( parts(M) /= parts(Minv) ) then 501 call die('Could not calculate the inverse on non equal sized & 502 &matrices') 503 end if 504 if ( parts(M) == 1 ) then 505 call die('This matrix is not tri-diagonal') 506 end if 507#endif 508 call timer('V_TM_inv',1) 509 510 nr = nrows_g(M) 511 np = parts(M) 512 crows => cum_rows(M) 513 514 ! This code is based on the down-folded self-energies 515 ! which are determined by the col region 516 517 sPart = huge(1) 518 ePart = 0 519 do n = 1 , r_col%n 520 s = pvt%r(r_col%r(n)) 521 sPart = min(sPart,which_part(M,s)) 522 ePart = max(ePart,which_part(M,s)) 523 end do 524#ifndef TS_NOCHECKS 525 if ( sPart < 1 ) call die('Error in the Bias inversion, sPart') 526 if ( ePart - sPart + 1 > 2 ) call die('Error in trimat partition') 527 if ( ePart > np ) call die('Error in the Bias inversion, ePart') 528#endif 529 530 ! Point to the matrices 531 z => val(Minv,all=.true.) 532 533 ! CHECK 534 ! This requires that the o_inD is sorted 535 ! according to the device region. 536 ! Check m_tbt_regions to assert this! 537 538 i_Elec = 1 539 do while ( i_Elec <= r_col%n ) 540 541 idx_Elec = pvt%r(r_col%r(i_Elec)) 542 543 ! We start by copying over the Mnn in blocks 544 545 ! We start by creating a region of consecutive 546 ! memory. 547 n = which_part(M,idx_Elec) 548 sN = nrows_g(M,n) 549 nb = 1 550 do while ( i_Elec + nb <= r_col%n ) 551 i = pvt%r(r_col%r(i_Elec+nb)) 552 ! In case it is not consecutive 553 if ( i - idx_Elec /= nb ) exit 554 ! In case the block changes, then 555 ! we cut the block size here. 556 if ( n /= which_part(M,i) ) exit 557 nb = nb + 1 558 end do 559 560 ! Copy over this portion of the Mnn array 561 562 ! Figure out which part we have Mnn in 563 i = pvt%r(r_col%r(i_Elec)) 564 n = which_part(M,i) 565 566 ! get placement of the diagonal block in the column 567 call TriMat_Bias_idxs(Minv,r_col%n,n,sIdxT,eIdxT) 568 Mpinv => z(sIdxT:eIdxT) 569 570 ! Correct the copied elements 571 ! Figure out the placement in the copied to array 572 ! First we calculate the starting index of the block 573 sIdxT = ( i_Elec - 1 ) * sN + 1 574 eIdxT = ( i_Elec + nb - 1) * sN 575 576 ! *** Now we have the matrix that we can save the 577 ! block Mnn in 578 579 ! We now need to figure out the placement of the 580 ! Mnn part that we have already calculated. 581 Mp => val(M,n,n) 582 i = pvt%r(r_col%r(i_Elec)) 583 sIdxF = (i-(crows(n)-sN)-1) * sN + 1 584 i = pvt%r(r_col%r(i_Elec+nb-1)) 585 eIdxF = (i-(crows(n)-sN)) * sN 586 587 ! Check that we have something correct... 588!print *,trim(El%name),sIdxT,eIdxT,sIdxF,eIdxF 589!print *,trim(El%name),eIdxT-sIdxT,eIdxF-sIdxF 590 591#ifndef TS_NOCHECKS 592 if ( eIdxT - sIdxT /= eIdxF - sIdxF ) & 593 call die('Error in determining column') 594#endif 595 596 ! Prepare for next segment 597 Mp => Mp(sIdxF:) 598 599 ! Copy over diagonal block 600 call zcopy(eIdxT-sIdxT+1,Mp(1),1,Mpinv(sIdxT),1) 601 602 ! Calculate the off-diagonal Green function in the regions 603 ! of interest 604 do in = n - 1 , sPart , - 1 605 606!print *,'calculating Mn-1n',in,n, i_Elec 607 ! Number of orbitals in the other segment 608 sNo = nrows_g(M,in) 609 ! grab the Yn matrix to perform the 610 ! Gf off-diagonal calculation. 611 XY => Yn_div_Bn_m1(M,in+1) 612 613 ! placement of the inverted matrix 614 call TriMat_Bias_idxs(Minv,r_col%n,in,sIdxT,eIdxT) 615 Mpinv => z(sIdxT:eIdxT) 616 sIdxT = ( i_Elec - 1 ) * sNo + 1 617 618#ifdef USE_GEMM3M 619 call zgemm3m( & 620#else 621 call zgemm( & 622#endif 623 'N','N',sNo,nb,sN, & 624 zm1, XY, sNo, Mp(1),sN,z0, Mpinv(sIdxT),sNo) 625 626 Mp => Mpinv(sIdxT:) 627 sN = sNo 628 629 end do 630 631 ! Reset arrays to just before off-diagonal 632 sN = nrows_g(M,n) 633 Mp => val(M,n,n) 634 Mp => Mp(sIdxF:) 635 636 ! Calculate the off-diagonal Green function in the regions 637 ! of interest 638 do in = n + 1 , ePart 639 640!print *,'calculating Mn+1n',in,n, i_Elec 641 ! Number of orbitals in the other segment 642 sNo = nrows_g(M,in) 643 ! grab the Xn matrix to perform the 644 ! Gf off-diagonal calculation. 645 XY => Xn_div_Cn_p1(M,in-1) 646 647 ! placement of the inverted matrix 648 call TriMat_Bias_idxs(Minv,r_col%n,in,sIdxT,eIdxT) 649 Mpinv => z(sIdxT:eIdxT) 650 sIdxT = ( i_Elec - 1 ) * sNo + 1 651 652#ifdef USE_GEMM3M 653 call zgemm3m( & 654#else 655 call zgemm( & 656#endif 657 'N','N',sNo,nb,sN, & 658 zm1, XY, sNo, Mp(1),sN,z0, Mpinv(sIdxT),sNo) 659 660 Mp => Mpinv(sIdxT:) 661 sN = sNo 662 663 end do 664 665 ! Update current segment of the electrode copied entries. 666 i_Elec = i_Elec + nb 667 668 end do 669 670 if ( present(only_diag) ) then 671 if ( only_diag ) then 672 call timer('V_TM_inv',2) 673 return 674 end if 675 end if 676 677 ! Now we need to calculate the remaining column 678 679 ! We now calculate: 680 ! Mmn = -Ym+1/Bm * Mm+1n, for m<n 681 do n = sPart - 1 , 1 , - 1 682 683 sN = nrows_g(M,n) 684 sNo = nrows_g(M,n+1) 685 686 ! get Ym+1/Bm 687 XY => Yn_div_Bn_m1(M,n+1) 688 689 ! Get Mm+1n 690 call TriMat_Bias_idxs(Minv,r_col%n,n+1,sIdxF,eIdxF) 691 Mp => z(sIdxF:eIdxF) 692 ! Get Mmn 693 call TriMat_Bias_idxs(Minv,r_col%n,n,sIdxF,eIdxF) 694 Mpinv => z(sIdxF:eIdxF) 695 696#ifdef USE_GEMM3M 697 call zgemm3m( & 698#else 699 call zgemm( & 700#endif 701 'N','N',sN,r_col%n,sNo, & 702 zm1, XY, sN, Mp(1),sNo,z0, Mpinv(1),sN) 703 704 end do 705 706 ! We now calculate: 707 ! Mmn = -Xm-1/Cm * Mm-1n, for m>n 708 do n = ePart + 1 , np 709 710 sNo = nrows_g(M,n-1) 711 sN = nrows_g(M,n) 712 713 ! get Xm-1/Cm 714 XY => Xn_div_Cn_p1(M,n-1) 715 716 ! Get Mm-1n 717 call TriMat_Bias_idxs(Minv,r_col%n,n-1,sIdxF,eIdxF) 718 Mp => z(sIdxF:eIdxF) 719 ! Get Mmn 720 call TriMat_Bias_idxs(Minv,r_col%n,n,sIdxF,eIdxF) 721 Mpinv => z(sIdxF:eIdxF) 722 723#ifdef USE_GEMM3M 724 call zgemm3m( & 725#else 726 call zgemm( & 727#endif 728 'N','N',sN,r_col%n,sNo, & 729 zm1, XY, sN, Mp(1),sNo,z0, Mpinv(1),sN) 730 731 end do 732 733 ! At this point the total 734 ! inverted column is placed at the end of 735 ! the tri-mat inversion. 736 737 call timer('V_TM_inv',2) 738 739 end subroutine invert_BiasTriMat_rgn 740 741 742 ! We will partition the system by: 743 ! 1. nrows_g(tri,1) x no 744 ! 2. nrows_g(tri,2) x no 745 ! ... 746 subroutine TriMat_Bias_idxs(M,no,p,sIdx,eIdx) 747 type(zTriMat), intent(in) :: M 748 ! no is the number of orbitals we wish to take out 749 ! p is the part that we wish to point to 750 integer, intent(in) :: no, p 751 integer, intent(out) :: sIdx, eIdx 752 integer, pointer :: crows(:) 753 754 integer :: cum 755 756 crows => cum_rows(M) 757 758 cum = nrows_g(M,p) 759 eIdx = no * cum - 1 760 cum = nrows_g(M) - crows(p) + cum 761 762 ! This is the number of elements already occupied 763 sIdx = elements(M, all=.true.) - no * cum + 1 764 eIdx = sIdx + eIdx 765 766 end subroutine TriMat_Bias_idxs 767 768 769 subroutine calc_Mnn_inv_cols(M,Minv,n,sCol,eCol) 770 type(zTriMat), intent(inout) :: M, Minv 771 integer, intent(in) :: n, sCol, eCol 772 ! Local variables 773 complex(dp), pointer :: Mp(:), Mpinv(:) 774 complex(dp), pointer :: Xn(:), Yn(:), Cn(:), Bn(:) 775 integer :: sNm1, sN, sNp1, ierr, i 776 777 if ( 1 < n ) sNm1 = nrows_g(M,n-1) 778 sN = nrows_g(M,n) 779 if ( n < parts(M) ) sNp1 = nrows_g(M,n+1) 780 781 if ( sCol == 1 .and. eCol == sN ) then 782 call calc_Mnn_inv(M,Minv,n) 783 return 784 end if 785 786 ! Retrieve Ann 787 Mp => val(M,n,n) 788 if ( n == 1 ) then 789 ! First we calculate M11^-1 790 ! Retrieve the X1/C2 array 791 Xn => Xn_div_Cn_p1(Minv,n) 792 ! The C2 array 793 Cn => val(M,n,n+1) 794 ! Calculate: A1 - X1 795#ifdef USE_GEMM3M 796 call zgemm3m( & 797#else 798 call zgemm( & 799#endif 800 'N','N',sN,sN,sNp1, & 801 zm1, Cn,sN, Xn,sNp1,z1, Mp,sN) 802 803 else if ( n == parts(M) ) then 804 805 ! Retrieve the Yn/Bn-1 array 806 Yn => Yn_div_Bn_m1(Minv,n) 807 ! The Bn-1 array 808 Bn => val(M,n,n-1) 809 ! Calculate: An - Yn 810#ifdef USE_GEMM3M 811 call zgemm3m( & 812#else 813 call zgemm( & 814#endif 815 'N','N',sN,sN,sNm1, & 816 zm1, Bn,sN, Yn,sNm1,z1, Mp,sN) 817 818 else 819 ! Retrieve the Xn/Cn+1 array 820 Xn => Xn_div_Cn_p1(Minv,n) 821 ! The Cn+1 array 822 Cn => val(M,n,n+1) 823 ! Calculate: An - Xn 824#ifdef USE_GEMM3M 825 call zgemm3m( & 826#else 827 call zgemm( & 828#endif 829 'N','N',sN,sN,sNp1, & 830 zm1, Cn,sN, Xn,sNp1,z1, Mp,sN) 831 ! Retrieve the Yn/Bn-1 array 832 Yn => Yn_div_Bn_m1(Minv,n) 833 ! The Bn-1 array 834 Bn => val(M,n,n-1) 835 ! Calculate: An - Xn - Yn 836#ifdef USE_GEMM3M 837 call zgemm3m( & 838#else 839 call zgemm( & 840#endif 841 'N','N',sN,sN,sNm1, & 842 zm1, Bn,sN, Yn,sNm1,z1, Mp,sN) 843 844 end if 845 846 ! Retrieve the position in the inverted matrix 847 Mpinv => val(Minv,n,n) 848 Mpinv((sCol-1)*sN+1:eCol*sN) = dcmplx(0._dp,0._dp) 849 do i = sCol - 1 , eCol - 1 850 Mpinv(i * sN + i + 1) = dcmplx(1._dp,0._dp) 851 end do 852 853 i = eCol - sCol + 1 854 call zgesv(sN,i,Mp,sN,ipiv,Mpinv((sCol-1)*sN+1),sN,ierr) 855 if ( ierr /= 0 ) then 856 call die('Error in inverting the partial bias block') 857 end if 858 859 end subroutine calc_Mnn_inv_cols 860 861end module m_ts_trimat_invert 862 863