Bisection.C 22.2 KB
 Francois Gygi committed Apr 08, 2012 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 //////////////////////////////////////////////////////////////////////////////// // // Copyright (c) 2011 The Regents of the University of California // // This file is part of Qbox // // Qbox is distributed under the terms of the GNU General Public License // as published by the Free Software Foundation, either version 2 of // the License, or (at your option) any later version. // See the file COPYING in the root directory of this distribution // or . // //////////////////////////////////////////////////////////////////////////////// // // Bisection.C // //////////////////////////////////////////////////////////////////////////////// #include "Bisection.h" #include #include  Francois Gygi committed May 26, 2012 21 22 23 24 #include "jade.h" #include "FourierTransform.h" using namespace std;  Francois Gygi committed Apr 08, 2012 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76  //////////////////////////////////////////////////////////////////////////////// int walsh(int l, int n, int i) { // walsh function at level l at position i in a grid of n points assert(i>=0); assert(i= n/2 ) return 0; else return 1; } else if ( l == 1 ) { assert(n%4==0); if ( i >= n/4 && i < (3*n)/4 ) return 0; else return 1; } else if ( l == 2 ) { assert(n%8==0); if ( (i >= n/8 && i < 3*n/8) || (i >= 5*n/8 && i < 7*n/8) ) return 0; else return 1; } else if ( l == 3 ) { assert(n%16==0); if ( (i >= n/16 && i < 3*n/16) || (i >= 5*n/16 && i < 7*n/16) || (i >= 9*n/16 && i < 11*n/16) || (i >= 13*n/16 && i < 15*n/16) ) return 0; else return 1; } else if ( l == 4 ) { assert(n%32==0); if ( (i >= n/32 && i < 3*n/32) || (i >= 5*n/32 && i < 7*n/32) || (i >= 9*n/32 && i < 11*n/32) || (i >= 13*n/32 && i < 15*n/32) || (i >= 17*n/32 && i < 19*n/32) || (i >= 21*n/32 && i < 23*n/32) || (i >= 25*n/32 && i < 27*n/32) || (i >= 29*n/32 && i < 31*n/32) ) return 0; else return 1; } else assert(false); } ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed Sep 17, 2019 77 78 Bisection::Bisection(const SlaterDet& sd, const int nlevels[3]) : ctxt_(sd.context())  Francois Gygi committed Apr 08, 2012 79 80 { // localization vectors are long int  Francois Gygi committed Apr 19, 2012 81  assert(sizeof(long int) >= 4);  Francois Gygi committed Apr 08, 2012 82   Francois Gygi committed May 26, 2012 83 84  nst_ = sd.nst(); nstloc_ = sd.nstloc();  Francois Gygi committed Apr 08, 2012 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100  // number of bisection levels in each direction nlevels_[0]=nlevels[0]; nlevels_[1]=nlevels[1]; nlevels_[2]=nlevels[2]; // number of subdivisions required in each direction // ndiv = 2^nlevel ndiv_[0] = 1 << nlevels[0]; ndiv_[1] = 1 << nlevels[1]; ndiv_[2] = 1 << nlevels[2]; // largest number of levels nlevelsmax_=max(nlevels[0],max(nlevels[1],nlevels[2])); // real-space grid size for wave functions  Francois Gygi committed May 26, 2012 101 102 103  np_[0] = sd.basis().np(0); np_[1] = sd.basis().np(1); np_[2] = sd.basis().np(2);  Francois Gygi committed Apr 08, 2012 104 105 106 107 108 109 110 111 112  // adapt the grid dimensions to the levels of bisection for ( int idim = 0; idim < 3; idim++ ) { for ( int j=1; j<=nlevels[idim]; j++ ) { int base = 1 << j; if ( np_[idim] % base != 0 ) np_[idim] += base/2; } }  Francois Gygi committed Sep 24, 2014 113 114 115  while (!sd.basis().factorizable(np_[0])) np_[0] += (1<np2_loc(); np012loc_ = ft_->np012loc(); // xy projector index: xy_proj[i+j*np0] is the xy projector // associated with a point located at i + j*np0 + k*np0*np1 xy_proj_.resize(np01_); for ( int i = 0; i < np_[0]; i++ ) for ( int j = 0; j < np_[1]; j++ ) { int i_slice_x= ( i * ndiv_[0] ) / np_[0]; int i_slice_y= ( j * ndiv_[1] ) / np_[1]; int xy_proj_index = i_slice_x + ndiv_[0] * i_slice_y; xy_proj_[i+np_[0]*j] = xy_proj_index; } // nmat_: number of A matrices nmat_ = nlevels[0] + nlevels[1] + nlevels[2];  Francois Gygi committed May 26, 2012 139 140 141 142 143  // each projector uses two bits of the localization vector // check that sizeof(long int) is sufficient // number of bits in a long int is 8*sizeof(long int) assert(2*nmat_ <= 8*sizeof(long int));  Francois Gygi committed Apr 08, 2012 144 145  // allocate A matrices amat_.resize(nmat_);  Francois Gygi committed May 26, 2012 146  const ComplexMatrix &c = sd.c();  Francois Gygi committed Apr 08, 2012 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162  for ( int i = 0; i < nmat_; i++ ) amat_[i] = new DoubleMatrix(c.context(),c.n(),c.n(),c.nb(),c.nb()); // allocate rotation matrix u_ = new DoubleMatrix(c.context(),c.n(),c.n(),c.nb(),c.nb()); // matrices of real space wave functions in subdomains rmat_.resize( ndiv_[0]*ndiv_[1] ); { // xyproj_rsize: real-space size of xy projectors vector xyproj_rsize( ndiv_[0]*ndiv_[1] , 0 ); for ( int ixy = 0; ixy < np01_; ixy++ ) xyproj_rsize[xy_proj_[ixy]] += np2loc_; // max_xyproj_rsize: maximum real-space size of xy projectors vector max_xyproj_rsize( ndiv_[0]*ndiv_[1] , 0 ); MPI_Allreduce((void*)&xyproj_rsize[0] , (void*)&max_xyproj_rsize[0],  Francois Gygi committed Mar 04, 2016 163  (int)xyproj_rsize.size(), MPI_INT ,MPI_MAX , ctxt_.comm());  Francois Gygi committed Apr 08, 2012 164 165 166 167  // allocate matrices rmat_[i] for ( int i = 0; i < rmat_.size(); i++ ) {  Francois Gygi committed May 26, 2012 168 169 170  int n = c.n(); int nb = c.nb(); int m = max_xyproj_rsize[i] * c.context().nprow();  Francois Gygi committed Apr 08, 2012 171  int mb = max_xyproj_rsize[i];  Francois Gygi committed May 26, 2012 172  rmat_[i] = new DoubleMatrix(c.context(),m,n,mb,nb);  Francois Gygi committed Apr 08, 2012 173 174 175 176  rmat_[i]->clear(); } }  Francois Gygi committed May 26, 2012 177  localization_.resize(nst_);  Francois Gygi committed Apr 08, 2012 178 179 180 181 182 183 184 185 186 187 188 189 } //////////////////////////////////////////////////////////////////////////////// Bisection::~Bisection(void) { delete ft_; for ( int i = 0; i < nmat_; i++ ) delete amat_[i]; delete u_; for ( int i = 0; i < rmat_.size(); i++ ) delete rmat_[i]; } ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed May 26, 2012 190 void Bisection::compute_transform(const SlaterDet& sd, int maxsweep, double tol)  Francois Gygi committed Apr 08, 2012 191 {  Francois Gygi committed May 26, 2012 192 193  // compute the transformation with tolerance tol  Francois Gygi committed Apr 08, 2012 194 195  map tmap;  Francois Gygi committed May 26, 2012 196 197 198  // check that basis is real // jade is not implemented for complex matrices assert( sd.basis().real() );  Francois Gygi committed Apr 08, 2012 199   Francois Gygi committed May 26, 2012 200 201 202  const ComplexMatrix& c = sd.c(); const vector& occ = sd.occ(); assert(occ.size() == nst_);  Francois Gygi committed Apr 08, 2012 203 204  #ifdef TIMING  Francois Gygi committed May 26, 2012 205  tmap["wf FFT"].start();  Francois Gygi committed Apr 08, 2012 206 #endif  Francois Gygi committed May 26, 2012 207 208 209 210 211 212 213 214 215  // Fourier transform states and save real-space values in // rmat_ matrices vector > wftmp(np012loc_,0.0); for ( int n = 0; n < nstloc_; n++ ) { ft_->backward(c.cvalptr(c.mloc()*n),&wftmp[0]); // pointers to rmat vector p_rmat( ndiv_[0]*ndiv_[1] ); for (int iproj=0; iprojmloc(); p_rmat[iproj] = rmat_[iproj]->valptr(index); } // copy wf to rmat arrays for ( int iz = 0; iz < np2loc_; iz++ ) { for ( int ixy = 0; ixy < np01_; ixy++ )  Francois Gygi committed Apr 08, 2012 224  {  Francois Gygi committed May 26, 2012 225 226 227  int xy_proj = xy_proj_[ixy]; (*p_rmat[xy_proj]) = wftmp[ixy + iz*np01_].real(); p_rmat[xy_proj]++;  Francois Gygi committed Apr 08, 2012 228 229  } }  Francois Gygi committed May 26, 2012 230  }  Francois Gygi committed Apr 08, 2012 231 #ifdef TIMING  Francois Gygi committed May 26, 2012 232  tmap["wf FFT"].stop();  Francois Gygi committed Apr 08, 2012 233 234 #endif  Francois Gygi committed May 26, 2012 235  // compute the x/y A matrices  Francois Gygi committed Apr 08, 2012 236 #ifdef TIMING  Francois Gygi committed May 26, 2012 237  tmap["xy products"].start();  Francois Gygi committed Apr 08, 2012 238 #endif  Francois Gygi committed May 26, 2012 239 240 241  // clear a matrices for ( int i = 0; i < nmat_; i++ ) amat_[i]->clear();  Francois Gygi committed Apr 08, 2012 242   Francois Gygi committed May 26, 2012 243  int size=amat_[0]->size();  Francois Gygi committed Apr 08, 2012 244   Francois Gygi committed May 26, 2012 245 246  // allocate matrix for products of projected parts DoubleMatrix products(c.context(),c.n(),c.n(),c.nb(),c.nb());  Francois Gygi committed Apr 08, 2012 247   Francois Gygi committed May 26, 2012 248 249 250 251 252 253  for ( int i_proj=0; i_projvalptr(0); for ( int i=0; ivalptr(0); double fac = 1.0 / (np_[0]*np_[1]*np_[2]); for ( int i = 0; i < size; i++ ) coeff[i] *= fac; imat++; } // else: matrix A is a projector in the z direction else if ( ilevelreal matrix product DoubleMatrix c_proxy(c); DoubleMatrix c_pz_proxy(c_pz); for ( int ilevel=0; ilevel p_rmat( ndiv_[0]*ndiv_[1] ); for ( int iproj=0; iprojmloc(); p_rmat[iproj] = rmat_[iproj]->valptr(index);  Francois Gygi committed Apr 08, 2012 344  }  Francois Gygi committed May 26, 2012 345 346  // save values of wf*walsh in rmat arrays for ( int iz = 0; iz < np2loc_; iz++ )  Francois Gygi committed Apr 08, 2012 347  {  Francois Gygi committed May 26, 2012 348 349 350  int izglobal = ft_->np2_first() + iz; double walsh_z = walsh(ilevel, np_[2], izglobal); for ( int ixy = 0; ixy < np01_; ixy++ )  Francois Gygi committed Apr 08, 2012 351  {  Francois Gygi committed May 26, 2012 352 353 354 355  int i = ixy + iz * np01_; int xy_proj = xy_proj_[ixy]; wftmp[i] = (*p_rmat[xy_proj]) * walsh_z; p_rmat[xy_proj]++;  Francois Gygi committed Apr 08, 2012 356 357  } }  Francois Gygi committed May 26, 2012 358  ft_->forward(&wftmp[0],c_pz.valptr(c_pz.mloc()*n));  Francois Gygi committed Apr 08, 2012 359  }  Francois Gygi committed May 26, 2012 360 361 362 363 364 365  // compute the product // factor -2.0 in next line: G and -G amat_[imat]->gemm('t','n',2.0,c_pz_proxy,c_proxy,0.0); // rank-1 update using first row of cd_proxy() and c_proxy amat_[imat]->ger(-1.0,c_pz_proxy,0,c_proxy,0); imat++;  Francois Gygi committed Apr 08, 2012 366 367  } }  Francois Gygi committed May 26, 2012 368  }  Francois Gygi committed Apr 08, 2012 369 #ifdef TIMING  Francois Gygi committed May 26, 2012 370  tmap["z products"].stop();  Francois Gygi committed Apr 08, 2012 371 372 373 #endif #ifdef DEBUG  Francois Gygi committed May 26, 2012 374  // check the values of the amat matrices  Francois Gygi committed Apr 08, 2012 375 #ifdef TIMING  Francois Gygi committed May 26, 2012 376  tmap["check_amat"].start();  Francois Gygi committed Apr 08, 2012 377 #endif  Francois Gygi committed May 26, 2012 378  check_amat(c);  Francois Gygi committed Apr 08, 2012 379 #ifdef TIMING  Francois Gygi committed May 26, 2012 380  tmap["check_amat"].stop();  Francois Gygi committed Apr 08, 2012 381 382 383 #endif #endif  Francois Gygi committed May 26, 2012 384 385 386  // set to zero matrix elements of the matrices amat_[i] if they couple // states with differing occupation numbers trim_amat(occ);  Francois Gygi committed Apr 08, 2012 387 388  #ifdef DEBUG_PRINT_MAT  Francois Gygi committed May 26, 2012 389 390 391 392 393  for ( int k = 0; k < amat_.size(); k++ ) { cout << "A(k=" << k << "):" << endl; cout << *amat_[k]; }  Francois Gygi committed Apr 08, 2012 394 395 #endif  Francois Gygi committed May 26, 2012 396  // joint approximate diagonalization of the matrices amat (jade)  Francois Gygi committed Apr 08, 2012 397 #ifdef TIMING  Francois Gygi committed May 26, 2012 398  tmap["jade"].start();  Francois Gygi committed Apr 08, 2012 399 400 #endif  Francois Gygi committed May 26, 2012 401 402  // diagonal values adiag_[k][i] // adiag_ is resized by jade  Francois Gygi committed Apr 08, 2012 403   Francois Gygi committed May 26, 2012 404  // diagonalize projectors  Francois Gygi committed Mar 27, 2018 405 406  // int nsweep = jade(maxsweep,tol,amat_,*u_,adiag_); jade(maxsweep,tol,amat_,*u_,adiag_);  Francois Gygi committed Feb 23, 2015 407 #ifdef TIMING  Francois Gygi committed Mar 04, 2016 408  if ( ctxt_.onpe0() )  Francois Gygi committed Feb 23, 2015 409 410 411  cout << "Bisection::compute_transform: nsweep=" << nsweep << " maxsweep=" << maxsweep << " tol=" << tol << endl; #endif  Francois Gygi committed Apr 08, 2012 412 413  #ifdef TIMING  Francois Gygi committed May 26, 2012 414  tmap["jade"].stop();  Francois Gygi committed Apr 08, 2012 415 416 #endif  Francois Gygi committed May 26, 2012 417 418 419 #ifdef DEBUG_PRINT_MAT cout << "U:" << endl; cout << *u_;  Francois Gygi committed Apr 08, 2012 420 421 #endif  Francois Gygi committed May 26, 2012 422 423 424 425  for ( int imat=0; imat::iterator i = tmap.begin(); i != tmap.end(); i++ ) { double time = (*i).second.real(); double tmin = time; double tmax = time;  Francois Gygi committed Mar 04, 2016 435 436 437  ctxt_.dmin(1,1,&tmin,1); ctxt_.dmax(1,1,&tmax,1); if ( ctxt_.myproc()==0 )  Francois Gygi committed Apr 08, 2012 438  {  Francois Gygi committed Sep 23, 2016 439 440  string s = "name=\"" + (*i).first + "\""; cout << ""  Francois Gygi committed Apr 08, 2012 443 444 445 446 447 448 449  << endl; } } #endif } ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed May 26, 2012 450 void Bisection::compute_localization(double epsilon)  Francois Gygi committed Apr 08, 2012 451 {  Francois Gygi committed May 26, 2012 452 453  // compute localization vector for a threshold epsilon for ( int n = 0; n < nst_; n++ )  Francois Gygi committed Apr 08, 2012 454  {  Francois Gygi committed May 26, 2012 455 456  localization_[n] = 0; for ( int imat = 0; imat < nmat_; imat++ )  Francois Gygi committed Apr 08, 2012 457  {  Francois Gygi committed May 26, 2012 458 459 460 461 462 463  if ( adiag_[imat][n] < epsilon ) localization_[n] += 1<<(2*imat); else if ( adiag_[imat][n] > 1.0-epsilon) localization_[n] += 1<<(2*imat+1); else localization_[n] += (1<<(2*imat)) + (1<<(2*imat+1));  Francois Gygi committed Apr 08, 2012 464 465 466  } }  Francois Gygi committed May 26, 2012 467 468 469 #ifdef DEBUG // print localization vector and number of overlaps (including self) // for each state  Francois Gygi committed Mar 04, 2016 470  if ( ctxt_.onpe0() )  Francois Gygi committed Apr 08, 2012 471  {  Francois Gygi committed May 26, 2012 472 473  int sum = 0; for ( int i = 0; i < nst_; i++ )  Francois Gygi committed Apr 08, 2012 474  {  Francois Gygi committed May 26, 2012 475 476 477 478 479 480 481 482 483 484 485  int count = 0; for ( int j = 0; j < nst_; j++ ) { if ( overlap(i,j) ) count++; } cout << "localization[" << i << "]: " << localization_[i] << " " << bitset<30>(localization_[i]) << " overlaps: " << count << endl; sum += count;  Francois Gygi committed Apr 08, 2012 486  }  Francois Gygi committed May 26, 2012 487 488  cout << "total overlaps: " << sum << " / " << nst_*nst_ << " = " << ((double) sum)/(nst_*nst_) << endl;  Francois Gygi committed Apr 08, 2012 489  }  Francois Gygi committed May 26, 2012 490 #endif  Francois Gygi committed Apr 08, 2012 491   Francois Gygi committed May 26, 2012 492 493  // broadcast localization to all tasks to ensure consistency MPI_Bcast( (void *) &localization_[0], localization_.size(),  Francois Gygi committed Mar 04, 2016 494  MPI_LONG, 0, ctxt_.comm() );  Francois Gygi committed Apr 08, 2012 495 496 497 498  } ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed May 26, 2012 499 void Bisection::forward(SlaterDet& sd)  Francois Gygi committed Apr 08, 2012 500 {  Francois Gygi committed May 26, 2012 501 502  forward(*u_,sd); }  Francois Gygi committed Apr 08, 2012 503   Francois Gygi committed May 26, 2012 504 505 506 507 508 509 510 511 512 513 514 //////////////////////////////////////////////////////////////////////////////// void Bisection::forward(DoubleMatrix& u, SlaterDet& sd) { // apply the bisection transformation to the SlaterDet sd // apply the rotation u to sd.c() ComplexMatrix& c = sd.c(); ComplexMatrix cp(c); DoubleMatrix cp_proxy(cp); DoubleMatrix c_proxy(c); c_proxy.gemm('n','n',1.0,cp_proxy,u,0.0); }  Francois Gygi committed Apr 08, 2012 515   Francois Gygi committed May 26, 2012 516 517 518 519 520 //////////////////////////////////////////////////////////////////////////////// void Bisection::backward(SlaterDet& sd) { backward(*u_,sd); }  Francois Gygi committed Apr 08, 2012 521   Francois Gygi committed May 26, 2012 522 523 524 525 526 527 528 529 530 531 //////////////////////////////////////////////////////////////////////////////// void Bisection::backward(DoubleMatrix& u, SlaterDet& sd) { // apply the inverse bisection transformation to SlaterDet sd // apply rotation u^T to sd ComplexMatrix& c = sd.c(); ComplexMatrix cp(c); DoubleMatrix cp_proxy(cp); DoubleMatrix c_proxy(c); c_proxy.gemm('n','t',1.0,cp_proxy,u,0.0);  Francois Gygi committed Apr 08, 2012 532 533 534 } ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed May 26, 2012 535 bool Bisection::check_amat(const ComplexMatrix &c)  Francois Gygi committed Apr 08, 2012 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 { // allocate memory for copy of the wave functions ComplexMatrix cd(c.context(),c.m(),c.n(),c.mb(),c.nb()); // precompute values of Walsh functions in three directions // wx[l][i] = value of level l Walsh function at position i in direction x vector > > w(3); for ( int idir=0; idir<3; idir++ ) { w[idir].resize(nlevels_[idir]); // for ( int l = 0; l < nlevels_[idir]; l++ ) { w[idir][l].resize(np_[idir]); for ( int i = 0; i < np_[idir]; i++ ) { w[idir][l][i] = walsh(l,np_[idir],i); } } } // compute matrices B_k = vector bmat(nmat_); for ( int k = 0; k < bmat.size(); k++ ) { bmat[k] = new DoubleMatrix(c.context(),c.n(),c.n(),c.nb(),c.nb()); } DoubleMatrix c_proxy(c); DoubleMatrix cd_proxy(cd); vector > wftmp(ft_->np012loc()); complex *f = &wftmp[0]; // compute matrices A at all levels for ( int l=0 , imat=0; l < nlevelsmax_; l++ ) { // x direction if ( lbackward(cd.cvalptr(cd.mloc()*n),&wftmp[0]); for ( int i = 0; i < np_[0]; i++ ) { if ( w[0][l][i] == 0 ) { for ( int j = 0; j < np_[1]; j++ ) for ( int k = 0; k < ft_->np2_loc(); k++ ) f[i + np_[0] * ( j + np_[1] * k )] = 0.0; } } ft_->forward(&wftmp[0],cd.valptr(cd.mloc()*n)); } // factor -2.0 in next line: G and -G bmat[imat]->gemm('t','n',2.0,cd_proxy,c_proxy,0.0); // rank-1 update using first row of cd_proxy() and c_proxy bmat[imat]->ger(-1.0,cd_proxy,0,c_proxy,0); imat++; } // y direction if ( lbackward(cd.cvalptr(cd.mloc()*n),&wftmp[0]); for ( int j = 0; j < np_[1]; j++ ) { if ( w[1][l][j] == 0 ) { for ( int i = 0; i < np_[0]; i++ ) for ( int k = 0; k < ft_->np2_loc(); k++ ) f[i + np_[0] * ( j + np_[1] * k )] = 0.0; } } ft_->forward(&wftmp[0],cd.valptr(cd.mloc()*n)); } // factor -2.0 in next line: G and -G bmat[imat]->gemm('t','n',2.0,cd_proxy,c_proxy,0.0); // rank-1 update using first row of cd_proxy() and c_proxy bmat[imat]->ger(-1.0,cd_proxy,0,c_proxy,0); imat++; } // z direction if ( lbackward(cd.cvalptr(cd.mloc()*n),&wftmp[0]); for ( int k = 0; k < ft_->np2_loc(); k++ ) { int kglobal = ft_->np2_first() + k; const int istart = k * np_[0] * np_[1]; if ( w[2][l][kglobal] == 0 ) { for ( int ij = 0; ij < np_[0]*np_[1]; ij++ ) f[istart+ij] = 0.0; } } ft_->forward(&wftmp[0],cd.valptr(cd.mloc()*n)); } // factor -2.0 in next line: G and -G bmat[imat]->gemm('t','n',2.0,cd_proxy,c_proxy,0.0); // rank-1 update using first row of cd_proxy() and c_proxy bmat[imat]->ger(-1.0,cd_proxy,0,c_proxy,0); imat++; } } // for l // testing the matrices for ( int imat=0; imatvalptr(0); double *b=bmat[imat]->valptr(0); int ncoeff=amat_[imat]->mloc()*amat_[imat]->nloc(); for ( int i=0; i1e-10) { cout << "error > 1.e-10 for matrix " << imat << endl; } } } return true; } //////////////////////////////////////////////////////////////////////////////// void Bisection::trim_amat(const vector& occ) { // set to zero the matrix elements of the matrices amat_[k] if they couple // states with differing occupation numbers  Francois Gygi committed Aug 09, 2014 679  const double trim_tol = 1.e-6;  Francois Gygi committed Apr 08, 2012 680 681 682 683 684 685 686 687  // check if all occupation numbers are the same double occ_max = 0.0, occ_min = 2.0; for ( int i = 0; i < occ.size(); i++ ) { occ_max = max(occ_max, occ[i]); occ_min = min(occ_min, occ[i]); } // return if all occupation numbers are equal  Francois Gygi committed Aug 09, 2014 688  if ( fabs(occ_max-occ_min) < trim_tol )  Francois Gygi committed Apr 08, 2012 689 690 691 692 693 694 695 696 697 698 699 700 701  return; const int mloc = amat_[0]->mloc(); const int nloc = amat_[0]->nloc(); // loop over elements local to this task for ( int i = 0; i < mloc; i++ ) { const int iglobal = amat_[0]->iglobal(i); for ( int j = 0; j < nloc; j++ ) { const int jglobal = amat_[0]->jglobal(j); const int ival = i + mloc * j;  Francois Gygi committed Aug 09, 2014 702  if ( fabs(occ[iglobal] - occ[jglobal]) > trim_tol )  Francois Gygi committed Apr 08, 2012 703 704 705 706 707 708 709  for ( int k = 0; k < amat_.size(); k++ ) (*amat_[k])[ival] = 0.0; } } } ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed Feb 12, 2015 710 bool Bisection::overlap(int i, int j) const  Francois Gygi committed May 26, 2012 711 712 713 714 715 { return overlap(localization_,i,j); } ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed Feb 12, 2015 716 bool Bisection::overlap(const vector& loc_, int i, int j) const  Francois Gygi committed Apr 08, 2012 717 {  Francois Gygi committed May 26, 2012 718 719 720 721 722  // overlap: return true if the functions i and j overlap according // to the localization vector loc_ long int loc_i = loc_[i]; long int loc_j = loc_[j]; while ( loc_i!=0 && loc_j!=0 )  Francois Gygi committed Apr 08, 2012 723 724  { // get the weight of projections for each state  Francois Gygi committed May 26, 2012 725 726 727 728  bool p_right_i = loc_i & 1; bool p_left_i = loc_i & 2; bool p_right_j = loc_j & 1; bool p_left_j = loc_j & 2;  Francois Gygi committed Apr 08, 2012 729 730 731 732 733  // return false as soon as the states are found to be separated if ( !( ( p_right_i && p_right_j ) || ( p_left_i && p_left_j ) ) ) return false;  Francois Gygi committed May 26, 2012 734 735  loc_i >>= 2; loc_j >>= 2;  Francois Gygi committed Apr 08, 2012 736 737 738 739 740 741 742  } // return true if the states overlap return true; } ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed Feb 12, 2015 743 double Bisection::pair_fraction(void) const  Francois Gygi committed Apr 08, 2012 744 745 { // pair_fraction: return fraction of pairs having non-zero overlap  Francois Gygi committed Feb 12, 2015 746  // count pairs (i,j) having non-zero overlap for i != j only  Francois Gygi committed Apr 08, 2012 747  int sum = 0;  Francois Gygi committed Feb 12, 2015 748  for ( int i = 1; i < nst_; i++ )  Francois Gygi committed Apr 08, 2012 749 750  { int count = 0;  Francois Gygi committed Feb 12, 2015 751  for ( int j = i+1; j < nst_; j++ )  Francois Gygi committed Apr 08, 2012 752  {  Francois Gygi committed May 26, 2012 753  if ( overlap(i,j) )  Francois Gygi committed Apr 08, 2012 754 755 756 757  count++; } sum += count; }  Francois Gygi committed Feb 12, 2015 758 759 760  // add overlap with self: (i,i) sum += nst_; return ((double) sum)/((nst_*(nst_+1))/2);  Francois Gygi committed Apr 08, 2012 761 762 763 } ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed Feb 12, 2015 764 double Bisection::size(int i) const  Francois Gygi committed Apr 08, 2012 765 766 { // size: return fraction of the domain on which state i is non-zero  Francois Gygi committed May 26, 2012 767  long int loc = localization_[i];  Francois Gygi committed Apr 08, 2012 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786  double size = 1.0; // process all projectors while ( loc != 0 ) { // weight of projections bool p_right = loc & 1; bool p_left = loc & 2; if ( (p_right && !p_left) || (!p_right && p_left) ) size *= 0.5; loc >>= 2; } // return true if the states overlap return size; } ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed Feb 12, 2015 787 double Bisection::total_size(void) const  Francois Gygi committed Apr 08, 2012 788 789 790 { // total_size: return normalized sum of sizes double sum = 0.0;  Francois Gygi committed May 26, 2012 791 792 793  for ( int i = 0; i < nst_; i++ ) sum += size(i); return sum / nst_;  Francois Gygi committed Apr 08, 2012 794 }