C++ Array2类代码示例

void Wf_return::setVals(Array2 <dcomplex> & vals ,
                        Array1 <doublevar> &  p) {
  
  // here we extract amplitude and phase, and their gradients and laplacians,
  // from the complex value and its gradient and derivative
  
  is_complex=1;
  cvals=vals;
  int ntype=vals.GetDim(1);
  int nwf=vals.GetDim(0);
  for (int w=0; w< nwf; w++) {
    amp(w,0)=vals(w,0).real();
    phase(w,0)=p(w);
    
    doublevar sum_ii=0;
    doublevar sum_ri=0;
    if (ntype>=4) {
      for (int i=1; i<4; i++) {
        amp(w,i)=vals(w,i).real();
        phase(w,i)=vals(w,i).imag();
        sum_ii+=phase(w,i)*phase(w,i);
        sum_ri+=amp(w,i)*phase(w,i);
      }
    }
    
    if (ntype>=5) {
      amp(w,4)=vals(w,4).real()+sum_ii;
      phase(w,4)=vals(w,4).imag()-2*sum_ri;
    }
  }
  
}

int main(int argc, char **argv)
{
  Triolet_init(&argc, &argv);
  List<TriInput> arr = make_inputs();
  Array2<Float> result = f(arr);

  // Check output
  int y, x;
  bool ok = true;
  for (y = 0; y < 2; y++) {
    for (x = 0; x < 3; x++) {
      float result_val = (Float)result.at(y, x);
      if (result_val != outputs[y][x]) ok = false;
    }
  }

  if (ok) {
    fputs("ok", stdout);
    return 0;
  } else {
    // Result does not match; print all outputs
    for (y = 0; y < 2; y++) {
      for (x = 0; x < 3; x++) {
        printf("%6f\t%6f\n", outputs[y][x], (float)(Float)result.at(y, x));
      }
    }
    return 1;
  }
}

void recenter(Array2 <doublevar> & pos, Array1 <doublevar> & mass) {
  int natoms=pos.GetDim(0);
  assert(pos.GetDim(0)==mass.GetDim(0));
  assert(pos.GetDim(1)==3);

  Array1 <doublevar> center(3,0.0); //center of mass
  doublevar tot=0.0;
  for(int at=0; at< natoms; at++) {
    for(int d=0; d< 3; d++) {
      center(d)+=mass(at)*pos(at,d);
    }
    tot+=mass(at);
  }

  for(int d=0; d< 3; d++) {
    center(d)/=tot;
  }

  for(int at=0; at< natoms; at++) {
    for(int d=0; d< 3; d++) {
     pos(at,d)-=center(d);
     }
  }

}

void copy_array (Array1& source, Array2& dest) {
  assert(std::equal(source.shape(),source.shape()+source.num_dimensions(),
                    dest.shape()));
  // Dispatch to the proper function
  typedef typename Array1::element element_type;
  copy_dispatch<element_type>::
    copy_array(source.begin(),source.end(),dest.begin());
}

inline void output_array(Array2 <doublevar> & arr) {
  for(int i=0; i< arr.GetDim(0); i++) {
    for(int j=0; j < arr.GetDim(1); j++) {
      cout << arr(i,j) << "  ";
    }
    cout << endl;
  }
}

/*
   Return the molecule orbital values
   */
void Average_ekt::calc_mos(Sample_point * sample, int e, Array2 <dcomplex> & movals) { 
  if(complex_orbitals) { 
    movals.Resize(nmo,1);
    cmomat->updateVal(sample,e,0,movals); 
  }
  else { 
    Array2 <doublevar> movals_d(nmo,1);
    momat->updateVal(sample,e,0,movals_d);
    movals.Resize(nmo,1);
    for(int i=0; i< nmo; i++) { 
      movals(i,0)=movals_d(i,0);
    }
  }
}

/*
   Calculate the values, gradient and laplacians of molecule orbitals, 
   returned as: [val, grad, lap] 
   */
void Average_ekt::calc_mosLap(Sample_point * sample, int e, Array2 <dcomplex> & molaps) { 
  if(complex_orbitals) { 
    molaps.Resize(nmo,5);
    cmomat->updateLap(sample, e, 0, molaps);
  }
  else { 
    Array2 <doublevar> molaps_d(nmo,5);
    momat->updateLap(sample,e,0,molaps_d);
    molaps.Resize(nmo,5);
    for(int i=0; i< nmo; i++) { 
      for (int d=0; d<5; d++) 
        molaps(i,d)=molaps_d(i,d);
    }
  }
}

// Extracts the inverted covariance matrix for a given branch from the full covariance matrix.
void branch_inverted_covariance_matrix(Model &Mod, Array2 &Covfull, long b, Array2 &Covbri) {
  long i,j;

  ublas::matrix<double> CC(Mod.df, Mod.df);
  ublas::matrix<double> II(Mod.df, Mod.df);

  for(i=0; i < Mod.df; i++) {
    for(j=0; j< Mod.df; j++) {
      CC(i,j) = Covfull[Mod.df*b + i][Mod.df*b + j];
    }
  }

  bool res = matrix_inverse(CC, II);
  if (!res) {
    throw std::out_of_range( "Could not invert the Covariance matrix." );
  }

  Covbri.resize(Mod.df);
  for(i=0; i < Mod.df; i++) {
    Covbri[i].resize(Mod.df);
    for(j=0; j< Mod.df; j++) {
      Covbri[i][j] = II(i,j);
    }
  }
}

// Computes the full covariance matrix for the MLE, using observed Fisher info.
void full_MLE_observed_covariance_matrix(Tree &T, Model &Mod, Parameters &Par, Counts &data, Array2 &Cov) {
  long i, j;

  long npars = Mod.df*T.nedges + Mod.rdf;
  ublas::matrix<double> CC(npars, npars);
  ublas::matrix<double> II(npars, npars);
  Array2 I;

  Observed_Fisher_information(T, Mod, Par, data, I);

  // Need to convert from Matrix to ublas::matrix<double>
  for(i=0; i < npars; i++) {
    for(j=0; j < npars; j++) {
      II(i, j) = I[i][j];
    }
  }

  bool res = matrix_inverse(II, CC);
  if (!res) {
    throw std::out_of_range("Could not invert the Fisher information matrix.");
  }

  Cov.resize(npars);
  for(i=0; i < npars; i++) {
    Cov[i].resize(npars);
    for(j=0; j < npars; j++) {
      Cov[i][j] = CC(i,j);
    }
  }
}

void MO_matrix_blas::writeorb(ostream & os, Array2 <doublevar> & rotation, Array1 <int>  &moList) {


  int nmo_write=moList.GetDim(0);
  assert(rotation.GetDim(0)==nmo_write);
  assert(rotation.GetDim(1)==nmo_write);
  os.precision(15);
  int counter=0;
  for(int m=0; m < nmo_write; m++) {
    int mo=moList(m);
    for(int ion=0; ion<centers.size(); ion++)
    {
      int f=0;

      for(int n=0; n< centers.nbasis(ion); n++)
      {
        int fnum=centers.basis(ion,n);
        int imax=basis(fnum)->nfunc();

        for(int i=0; i<imax; i++)
        {
          os << mo+1 << "  "   << f+1 << "   " << ion+1 << "   " << counter+1 << endl;
          f++;  //keep a total of functions on center
          counter++;
        } //i
      } //n
    }  //ion
  }
  os << "COEFFICIENTS\n";
  ifstream orbin(orbfile.c_str());
  rotate_orb(orbin, os, rotation, moList, totbasis);
  orbin.close();


}

void Cutoff_cusp::calcLap(
  const Array1 <doublevar> & r,
  Array2 <doublevar> & symvals,
  const int startfill
)
{
  assert(r.GetDim(0) >=5);
  assert(symvals.GetDim(0) >= 1+startfill);
  assert(symvals.GetDim(1) >= 5);
  if(r(0) < rcut) {
    doublevar zz=r(0)*rcutinv;
    doublevar zz2=zz*zz;
    doublevar pp=zz-zz2+zz*zz2/3;
    doublevar pade=1./(1+gamma*pp);
    symvals(startfill,0)=cusp*rcut*(pp*pade-pade0);
    doublevar pade2=pade*pade;
    doublevar ppd=1.-2.*zz+zz2;
    doublevar ppdd=-2.+2.*zz;
    doublevar dadr=ppd*pade2/r(0);
    doublevar dadr2=ppdd*pade2*rcutinv
                    -2.*gamma*ppd*ppd*pade2*pade*rcutinv
                    +2.*dadr;

    for(int i=1; i< 4; i++) {
      symvals(startfill, i)=cusp*dadr*r(i+1);
    }
    symvals(startfill, 4)=cusp*dadr2;
  }
  else {
    for(int i=0; i< 5; i++)
      symvals(startfill, i)=0;
  }
}

//----------------------------------------------------------------------
//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 Molecular_system::getBounds(Array2 <doublevar> & latticevec,Array1<doublevar> & origin) { 
  cout << "getBounds()" << endl;
  latticevec.Resize(3,3);
  latticevec=0.0;
  Array1 <doublevar> minmax(6);
  for(int d=0; d< 3; d++) { 
    minmax(2*d)=minmax(2*d+1)=0.0;
  }
  
  int nions=nIons();
  Array1 <doublevar> ionpos(3);
  for(int i=0; i< nions; i++) {
    getIonPos(i,ionpos);
    for(int d=0; d< 3; d++) {
      if(ionpos(d) < minmax(2*d)) minmax(2*d) = ionpos(d);
      if(ionpos(d) > minmax(2*d+1)) minmax(2*d+1)=ionpos(d);
    }
  }

  for(int d=0; d< 3; d++) {
   // minmax(2*d)-=4.0;
   //minmax(2*d+1)+=4.0;
    origin(d)=minmax(2*d)-4.0;
    latticevec(d,d)=minmax(2*d+1)-origin(d)+4.0;
  }
  cout << "getBounds done" << endl;
  return 1;
}

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
}

doublevar Wannier_method::eval_tstep(Array3 <dcomplex> & eikr, Array2 <doublevar> & Rgen,
                                     Array2 <doublevar> & Rgen_save, Array2 <doublevar> & deriv, doublevar tstep,
                                     Array2 <doublevar> & R) {
    int norb=Rgen.GetDim(0);
    for(int ii=0; ii< norb; ii++) {
        for(int jj=ii+1; jj < norb; jj++) {
            Rgen(ii,jj)=Rgen_save(ii,jj)-tstep*deriv(ii,jj);
        }
    }
    return evaluate_local(eikr,Rgen,R);
}