C++ Array1::Resize方法代码示例

//----------------------------------------------------------------------
//Keeping this separate from calcLoc for efficiency reasons..
//It's most likely faster to add to a double than fill matrices..
//We also don't need to calculate the ion-ion interaction for this
void Molecular_system::separatedLocal(Sample_point * sample,Array1 <doublevar> & onebody,
    Array2 <doublevar> & twobody)
{
  int nions=sample->ionSize();
  int nelectrons=sample->electronSize();
  onebody.Resize(nelectrons);
  twobody.Resize(nelectrons,nelectrons);
  onebody=0.0; twobody=0.0;
  Array1 <doublevar> R(5);
  sample->updateEIDist();
  sample->updateEEDist();

  for(int e=0; e< nelectrons; e++) {
    for(int i=0; i < nions; i++) {
      sample->getEIDist(e,i, R);
      onebody(e)-=sample->getIonCharge(i)/R(0);
    }
  }
  Array1 <doublevar> R2(5);
  for(int i=0; i< nelectrons; i++) {
    for(int j=i+1; j<nelectrons; j++) {
      sample->getEEDist(i,j,R2);
      twobody(i,j)+= 1/R2(0);
    }
  }

  Array1 <doublevar> pos(3);
  for(int e=0; e< nelectrons; e++) {
    sample->getElectronPos(e,pos);
    for(int d=0; d< 3; d++) 
      onebody(e)-=electric_field(d)*pos(d);
  }
}

int binary_read_checksum(Array1 <doublevar> & a, FILE * f,int nhint) { 
  int ntry=0;
  while(true) { 
  //First try to read in zeros
    ntry++;
    if(ntry>1) cout << "WARNING: detected corruption" << endl;
    doublevar zero1,zero2;
    int ntry_zero=0;
    while(true) { 
      ntry_zero++;
      if(!fread(&zero1,sizeof(doublevar),1,f)) return 0;
      //cout << "zero1 " << zero1 << endl;
      if(zero1==0) {
          if(!fread(&zero2,sizeof(doublevar),1,f)) return 0;
          if(zero2==0) break;
      }
    }
    if(ntry_zero> 1) cout << "WARNING: possible corruption: had to search " << ntry_zero << endl;
    doublevar n_d;
    fread(&n_d,sizeof(doublevar),1,f);
    int n=int(n_d);
    if(n_d-n == 0 and ( nhint<0 or n==nhint))  { 
      a.Resize(n);
      fread(a.v,sizeof(doublevar),n,f);
      doublevar check_array=checksum(a);
      doublevar check_file=1e99;
      fread(&check_file,sizeof(doublevar),1,f);
      if(fabs(check_array-check_file) < 1e-10) return ntry;
      else { 
        cout << "WARNING: detected corruption" << endl;
      }
      if(ferror(f)) return 0;
    }
  }
}

void Force_fitter::fit_force(Array1 <doublevar> & r, 
                             Array1 <doublevar> & bare_force,
                             Array1 <doublevar> & fit_force) {
  int ndim=bare_force.GetDim(0);

  fit_force=bare_force; return;

  fit_force.Resize(ndim);
  fit_force=0;

  if(r(0) < cut) {
    for(int d=0; d< ndim; d++) {
      doublevar basis=0;
      for(int i=0; i< nexp; i++) {
        basis+=coeff(i)*pow(r(0),i+m);
      }
      fit_force(d)=bare_force(d)*basis;
    }
  }
  else {
    fit_force=bare_force;
  }


}

void Jastrow_threebody_piece_diffspin::getParms(Array1 <doublevar> & parms) {
  if(freeze) {
    parms.Resize(0);
  }
  else {
    parms.Resize(nparms());
    int counter=0;
    int natoms=parm_centers.GetDim(0);
    for(int s=0; s< unique_parameters_spin.GetSize();s++){
      for(int i=0; i< unique_parameters_spin(s).GetDim(0); i++) {
	for(int j=0; j< _nparms(i); j++) {
	  parms(counter++) = unique_parameters_spin(s)(i,j);//*natoms;
	}
      }
    }
  }

}

void getBlocks(vector <string> & words, vector <string> & labels, 
               Array1 < Array1 < Properties_block> > & allblocks, 
               Array1 < Array1 < Average_generator * > > & avg_gen) {

  vector < vector < string> > blocktext;
  vector <string> tmp;
  unsigned int pos=0;
  while(readsection(words, pos, tmp, "BLOCK"))
    blocktext.push_back(tmp);

  int nblock=blocktext.size();
  //cout << blocktext.size() << " blocks \n";

  //vector <string> labels;
  vector <int> belong_to;
  vector <int> lab_nb; //number of blocks that belong to each label
  string label;
  for(int b=0; b< nblock; b++) {
    if(!readvalue(blocktext[b], pos=0, label, "LABEL"))
      error("didn't find label in the file");
    vector<string>::iterator place=find(labels.begin(), labels.end(), label);
    belong_to.push_back(place-labels.begin());
    if(place==labels.end()) {
      labels.push_back(label);
      lab_nb.push_back(1);
    }
    else lab_nb[belong_to[b]]++;
  }
  
  
  
  
  int nlabels=labels.size();
  //cout << nlabels << " labels" << endl;
  //for(int i=0; i< nlabels; i++) {
  //  cout << labels[i] << "  " << lab_nb[i] << "  " << lab_ninit[i] << endl;
  //}

  allblocks.Resize(nlabels);
  for(int i=0; i< nlabels; i++) {
    allblocks(i).Resize(lab_nb[i]);
  }
  
  
  Array1 <int> nb_read(nlabels, 0);
  for(int b=0; b< nblock; b++) {
    int lab=belong_to[b];
    int i=nb_read(lab);
    nb_read(lab)++;
    allblocks(lab)(i).restoreFromLog(blocktext[b]);
  }
  
  getInit(words, labels, avg_gen);
  
}

void Rgaussian_function::getVarParms(Array1 <doublevar> & parms) {
  //cout << "getVarParms " << endl;
  parms.Resize(nparms());

  int count=0;
  for(int l=0; l < nmax; l++) { 
    for(int e=0; e< numExpansion(l); e++) { 
      parms(count++)=log(gaussexp(l,e));
      //parms(count++)=gausscoeff(l,e);
    }
  }
}

int HEG_system::enforcePbc(Array1 <doublevar> & pos, Array1 <int> & nshifted) {
  assert(pos.GetDim(0) >=3);
  int shifted=0;
  nshifted.Resize(3);
  nshifted=0;

  int nshift=0;
  for(int i=0; i< 3; i++) {
    int shouldcontinue=1;
    while(shouldcontinue) {

      //Whether we're past the origin side
      doublevar tooshort=0;
      for(int j=0; j<3;j++) tooshort+=normVec(i,j)*(pos(j)-origin(j));

      //Whether we're past the lattice vector
      doublevar toofar=0;
      for(int j=0; j< 3; j++) toofar+=normVec(i,j)*(pos(j)-corners(i,j));

      //the 1e-12 seems to help avoid numerical problems, esp 
      //when integrating over a grid(which tends to hit the edges)
      if(tooshort < -1e-12) {
	//cout <<"tooshort " << tooshort << endl;
        for(int j=0; j< 3; j++) pos(j)+=latVec(i,j);
        shifted=1;
        nshifted(i)+=1;
        nshift++;
      }
      else if(toofar > 1e-12) {
	//cout << "toofar " << toofar << endl;
        for(int j=0; j< 3; j++) pos(j)-=latVec(i,j);
        shifted=1;
	// JK: here was +1, which works fine for real k-points that
	// correspond to standing waves. For general (complex) k-points the
	// wavefunction is "directional", however.
        nshifted(i)-=1;

        nshift++;
      }
      else {
        shouldcontinue=0;
      }
      if(nshift > 1000)
	
        error("Did over 1000 shifts and we're still out of the simulation cell."
	            "  There's probably something wrong.  Position : ", pos(i));
    }
  }

  return shifted;

}

//We write different versions in case we want to use BLAS routines, which don't template well..
dcomplex InverseUpdateRow(Array2 <dcomplex> & a1, const Array1 <dcomplex> & newRow,
                           const int lRow, const int n)
{
  Array1 <dcomplex>  tmpColL;
  tmpColL.Resize(n);
  Array1 <dcomplex>  prod;
  prod.Resize(n);

  dcomplex f=0.0;
  for(int i=0;i<n;++i) {
    f += newRow[i]*a1(i,lRow);
  }
  f =-1.0/f;

  for(int j=0;j<n;++j){
    prod[j]   =0.0;
    tmpColL[j]=a1(j,lRow);
    for(int i=0;i<n;++i) {
      prod[j] += newRow[i]*a1(i,j);
    }
    prod[j] *= f;
  }

  for(int i=0;i<n;++i) {
    dcomplex t(tmpColL[i]);
    for(int j=0;j<n;++j) {
      a1(i,j) += t*prod[j];
    }
  }

  f = -f;
  for(int i=0;i<n;++i) {
    a1(i,lRow) = f*tmpColL[i];
  }
  return f;
}

void getInit(vector <string> & words, vector <string> & labels, Array1 < Array1 <Average_generator*> > & avg_gen) {
  unsigned int pos=0;
  pos=0;
  vector<string> tmp;
  vector <vector<string> > inittext;
  while(readsection(words, pos, tmp, "INIT")) 
    inittext.push_back(tmp);
  int ninit=inittext.size();
  int nlabels=labels.size();
  //cout << "ninit " << ninit << endl;
  //if we have multiple init sections, we should check to see that they're 
  //all identical.  For now, let's take the first one.
  avg_gen.Resize(nlabels);
  
  if(ninit==0) { 
    for(int i=0; i< nlabels; i++) avg_gen(i).Resize(0); 
    return; 
  }

  for(int i=0; i< nlabels; i++) { 
    int ini=0;
    string inilabel;
    readvalue(inittext[ini],pos=0,inilabel, "LABEL");
    while(labels[i]!=inilabel) {
      ini++;
      readvalue(inittext[ini],pos=0,inilabel, "LABEL");
      if(ini>= ninit)  {
        cout << "Couldn't find init section for label " << labels[i] 
        << ".  AVERAGE { } sections will not be reported.\n";
        break;
      }
    }
    if(ini < ninit) { 
      vector <vector <string> > gen_secs;
      pos=0;
      while(readsection(inittext[ini],pos,tmp, "AVERAGE_GENERATOR"))
        gen_secs.push_back(tmp);
      avg_gen(i).Resize(gen_secs.size());
      for(int j=0; j< gen_secs.size(); j++) { 
        allocate(gen_secs[j],avg_gen(i)(j));
      }
    }
  }
  
}

void GeneralizedEigenSystemSolverRealGeneralMatrices(Array2 < doublevar > & Ain, 
         Array1 <dcomplex> & W, Array2 <doublevar> & VL, Array2 <doublevar> & VR) { 
#ifdef USE_LAPACK //if LAPACK
  int N=Ain.dim[0];
  
  Array2 <doublevar> A_temp=Ain; //,VL(N,N),VR(N,N);
  Array1 <doublevar> WORK,RWORK(2*N),WI(N),WR(N);
  WI.Resize(N);
  VL.Resize(N,N);
  VR.Resize(N,N);

  int info;
  int NB=64;
  int NMAX=N;
  int lda=NMAX;
  int ldb=NMAX;
  int LWORK=5*NMAX;
  WORK.Resize(LWORK);


  info=  dgeev('V','V',N,A_temp.v, lda,WR.v,WI.v,VL.v,lda,VR.v,lda,WORK.v,LWORK);

  if(info>0)
    error("Internal error in the LAPACK routine dgeev",info);
  if(info<0)
    error("Problem with the input parameter of LAPACK routine dgeev in position ",-info);
  W.Resize(N);
  for(int i=0; i< N; i++) { 
    W(i)=dcomplex(WR(i),WI(i));
  }

//  for (int i=0; i<N; i++)
//    evals(i)=W[N-1-i];

 // for (int i=0; i<N; i++) {
 //   for (int j=0; j<N; j++) {
 //     evecs(j,i)=A_temp(N-1-i,j);
 //   }
 // }
 //END OF LAPACK 
#else //IF NO LAPACK
  error("need LAPACK for eigensystem solver for general matrices");
#endif //END OF NO LAPACK
}

void Slat_Jastrow::evalTestPos(Array1 <doublevar> & pos, Sample_point * sample,Array1 <Wf_return> & wf) {
  Array1 <Wf_return> slat_val,jast_val;
  slater_wf->evalTestPos(pos,sample,slat_val);
  jastrow_wf->evalTestPos(pos,sample,jast_val);
  assert(slat_val.GetDim(0)==jast_val.GetDim(0));
  int n=slat_val.GetDim(0);
  wf.Resize(n);
  for(int i=0; i< n; i++) { 
    wf(i).Resize(nfunc_,2);
    if ( slat_val(i).is_complex==1 || jast_val(i).is_complex==1 )
      wf(i).is_complex=1;
    
    for(int j=0; j< nfunc_; j++) { 
      wf(i).phase(j,0)=slat_val(i).phase(j,0)+jast_val(i).phase(j,0); 
      wf(i).amp(j,0)=slat_val(i).amp(j,0)+jast_val(i).amp(j,0);  //add the logarithm
    }
  }

}

void Molecular_system::locDerivative(int ion, Sample_point * sample,
                                     Force_fitter & fitter,
                                     Array1 <doublevar> & der) {

  sample->updateEIDist();
  der.Resize(3);
  der=0;
  Array1 <doublevar> R(5);
  int ne=sample->electronSize();
  Array1 <doublevar> tmpder(3), tmpder_fit(3);
  for(int e=0; e< ne; e++) {
    sample->getEIDist(e,ion, R);
    for(int d=0; d< 3; d++) 
      //der(d)-=sample->getIonCharge(ion)*R(d+2)/(R(1)*R(0));
      tmpder(d)=-sample->getIonCharge(ion)*R(d+2)/(R(1)*R(0));
    fitter.fit_force(R,tmpder, tmpder_fit);
    for(int d=0; d< 3; d++)
      der(d)+=tmpder_fit(d);
  }

  Array1 <doublevar> vec(3);
  Array1 <doublevar> r_ion(3);
  sample->getIonPos(ion, r_ion);

  int nIon=sample->ionSize();

  for(int j=0; j< nIon; j++) {
    if(j!=ion) {
    sample->getIonPos(j,R);
    doublevar r2=0;
    for(int d=0; d <3; d++) 
      vec(d)=r_ion(d)-R(d);

    for(int d=0; d< 3; d++)
      r2+=vec(d)*vec(d);
    doublevar r=sqrt(r2);

    for(int d=0; d< 3; d++) 
      der(d)-=sample->getIonCharge(ion)*sample->getIonCharge(j)*vec(d)/(r2*r);
    }
  }

}

void Molecular_system::calcLocWithTestPos(Sample_point * sample,
                                          Array1 <doublevar> & tpos,
                                          Array1 <doublevar> & Vtest) {

  int nions=sample->ionSize();
  int nelectrons=sample->electronSize();
  Vtest.Resize(nelectrons + 1);
  Vtest = 0; 
  Array1 <doublevar> oldpos(3);
  //cout << "Calculating local energy\n";
  sample->getElectronPos(0, oldpos);
  sample->setElectronPosNoNotify(0, tpos);
  
  Array1 <doublevar> R(5);

  sample->updateEIDist();
  sample->updateEEDist();

  for(int i=0; i < nions; i++) {
    sample->getEIDist(0, i, R);
    Vtest(nelectrons)+=-sample->getIonCharge(i)/R(0);
  }
  doublevar dist = 0.0; 
  for(int d=0; d<3; d++)  {
    dist += (oldpos(d) - tpos(d))*(oldpos(d) - tpos(d)); 
  }
  dist = sqrt(dist); 
  Vtest(0) = 1/dist; 
  Array1 <doublevar> R2(5);
  for(int i=1; i< nelectrons; i++) {
    sample->getEEDist(0,i,R2);
    Vtest(i)= 1/R2(0);
  }
  sample->setElectronPosNoNotify(0, oldpos);
  sample->updateEIDist();
  sample->updateEEDist();
  //cout << "elec-elec: " << elecElec << endl;
  //cout << "pot " << pot << endl;

  for(int d=0; d< 3; d++) 
    Vtest(nelectrons) -=electric_field(d)*tpos(d);
}

int reblock_average(Array1 <Properties_block> & orig_block, int reblock, int equil,
            Properties_final_average & avg) {
  int nblock=orig_block.GetDim(0);
  if(reblock > nblock) 
    return 0;
  
  int neffblock=0;
  Array1<Properties_block> reblocks;

  reblocks.Resize(nblock/reblock);
  int start=nblock%reblock;
  int count=0;
  for(int i=start; i< nblock; i+=reblock) {
    int top=min(i+reblock, nblock);
    reblocks(count++).reduceBlocks(orig_block, i, top);
  }
  neffblock=nblock/reblock;
  avg.blockReduce(reblocks, 0, nblock/reblock, equil);

  return neffblock;
}

void sort_abs_values_descending(const Array1 <doublevar> & vals, Array1 <doublevar> & newvals, Array1 <int> & list){
  int n=vals.GetSize();
  list.Resize(n); 
  for (int i=0; i < n; i++) {
    list(i) = i;
    newvals(i)=vals(i);
  }
  
  for (int i=1; i < n; i++) {
    const double temp = newvals[i];
    const double abstemp =  fabs(newvals[i]);
    int j;
    for (j=i-1; j>=0 && fabs(newvals[j])<abstemp; j--) {
      newvals[j+1] = newvals[j];
      list[j+1] = list[j];
    }
    newvals[j+1] = temp;
    list[j+1] = i;
  }
  //for (int i=0; i < n; i++) 
  //cout << list[i]<<endl;
}