Changeset 807e8a

Timestamp:

May 30, 2008, 3:48:16 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

9842d8

Parents:

0fe2ca

Message:

libblas needed for libgsl is now found via AC_SEARCH_LIB and list blas cblas gslblas gslcblas

Files:

• molecuilder/configure.ac (modified) (1 diff)
• pcp/configure.ac (modified) (1 diff)
• pcp/src/perturbed.c (modified) (16 diffs)

molecuilder/configure.ac

@@ -36 +36 @@
 # Checks for library functions.
 # check for GNU Scientific Library
-AC_CHECK_LIB(blas,main)
+AC_SEARCH_LIBS(cblas_sdot, blas cblas gslblas gslcblas)
 AC_CHECK_LIB(gsl,main)
 

pcp/configure.ac

@@ -179 +179 @@
 # check for GNU Scientific Library
 #AC_CHECK_LIB(m,main)
-AC_CHECK_LIB(blas,main)
-AC_CHECK_LIB(gsl,main)
+AC_SEARCH_LIBS(cblas_sdot, blas cblas gslblas gslcblas)
+AC_CHECK_LIB(gsl, main)
 
 dnl Check for "extern inline", using a modified version

pcp/src/perturbed.c

@@ -360 +360 @@
   double lambda, Lambda;
   double RElambda10, RELambda10;
+  //double RElambda01, RELambda01;
   const fftw_complex *source = LevS->LPsi->LocalPsi[l];
   fftw_complex *grad = P->Grad.GradientArray[ActualGradient];
@@ -472 +473 @@
   // "-" due to purely imaginary wave function is on left hand side, thus becomes complex conjugated: i -> -i
   // (-i goes into pert. op., "-" remains when on right hand side)  
-  RElambda10 =  GradSP(P,LevS,varphi_1,TempPsi_0) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
-  //RElambda01 = -GradSP(P,LevS,varphi_0,TempPsi_1) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
+  RElambda10 = GradSP(P,LevS,varphi_1,TempPsi_0) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
+  //RElambda01 = GradSP(P,LevS,varphi_0,TempPsi_1) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
   
   MPI_Allreduce ( &RElambda10, &RELambda10, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
@@ -622 +623 @@
 }
 
-#define stay_above 0.001 //!< value above which the coefficient of the wave function will always remain
+#define stay_above 0.00001 //!< value above which the coefficient of the wave function will always remain
 
 /** Finds the minimum of perturbed energy in regards of actual wave function.
@@ -673 +674 @@
   double LocalSP, PsiSP;
   double dEdt0, ddEddt0, ConDirHConDir, ConDirConDir;//, sourceHsource;
-  double E0, E, delta;
-  //double E0, E, dE, ddE, delta, dcos, dsin;
-  //double EI, dEI, ddEI, deltaI, dcosI, dsinI;
+  //double E0, E, delta;
+  double E0, E, dE, ddE, delta, dcos, dsin;
+  double EI, dEI, ddEI, deltaI, dcosI, dsinI;
   //double HartreeddEddt0, XCddEddt0;
   double d[4],D[4], Diff;
@@ -707 +708 @@
     x = .5*LevS->GArray[g].GSq;
     // FIXME: Good way of accessing reciprocal Lev0 Density coefficients on LevS! (not so trivial)  
-    //x += sqrt(Dens->DensityCArray[HGDensity][g].re*Dens->DensityCArray[HGDensity][g].re+Dens->DensityCArray[HGDensity][g].im*Dens->DensityCArray[HGDensity][g].im);
+    x += sqrt(Dens->DensityCArray[HGDensity][g].re*Dens->DensityCArray[HGDensity][g].re+Dens->DensityCArray[HGDensity][g].im*Dens->DensityCArray[HGDensity][g].im);
     x -= T/(double)Num;
     K = x/(x*x+stay_above*stay_above);
@@ -797 +798 @@
   d[1] = GradSP(P,LevS,ConDir,HConDir);
   d[2] = GradSP(P,LevS,ConDir,ConDir);
-  d[3] = 0.;
+  d[3] = 0;
 
   // gather results  
@@ -808 +809 @@
   ConDirHConDir = D[1];
   ConDirConDir = D[2];
+  //sourceHsource = D[3];
   ddEddt0 = 0.0;
   //if (ConDir != P->Grad.GradientArray[ConDirGradient]) Error(SomeError,"FindPerturbedMinimum: ConDir corrupted");
   //if (H1c_grad != P->Grad.GradientArray[H1cGradient]) Error(SomeError,"FindPerturbedMinimum: H1c_grad corrupted");
+  //fprintf(stderr, "lambda*PsiFactor %lg vs. sourceHSource %lg\n", Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda * Psi->LocalPsiStatus[R->ActualLocalPsiNo].PsiFactor, sourceHsource);
+  // note: they really are exactly the same!
   ddEddt0 = Calculate2ndPerturbedDerivative(P, LevS->LPsi->LocalPsi[p], source, ConDir, Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda * Psi->LocalPsiStatus[R->ActualLocalPsiNo].PsiFactor, ConDirHConDir, ConDirConDir);
+  //ddEddt0 = 1.e+5;
 
   for (i=MAXOLD-1; i > 0; i--)
@@ -822 +827 @@
   //if (isnan(ddEddt0)) { fprintf(stderr,"(%i) WARNING in CalculateLineSearch(): ddEddt0_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
 
-  ////deltaI = CalculateDeltaI(E0, dEdt0, ddEddt0, 
-  ////      &EI, &dEI, &ddEI, &dcosI, &dsinI);
-  ////delta = deltaI; E = EI; dE = dEI; ddE = ddEI; dcos = dcosI; dsin = dsinI;
-  if (ddEddt0 > 0) {
-    delta = - dEdt0/ddEddt0;
-    E = E0 + delta * dEdt0 + delta*delta/2. * ddEddt0;
-  } else {
-    delta = 0.;
-    E = E0;
-    fprintf(stderr,"(%i) Taylor approximation leads not to minimum!\n",P->Par.me);
-  }
+  deltaI = CalculateDeltaI(E0, dEdt0, ddEddt0, 
+        &EI, &dEI, &ddEI, &dcosI, &dsinI);
+  delta = deltaI; E = EI; dE = dEI; ddE = ddEI; dcos = dcosI; dsin = dsinI;
+//  if (ddEddt0 > 0) {
+//    delta = - dEdt0/ddEddt0;
+//    E = E0 + delta * dEdt0 + delta*delta/2. * ddEddt0;
+//  } else {
+//    delta = 0.;
+//    E = E0;
+//    fprintf(stderr,"(%i) Taylor approximation leads not to minimum!\n",P->Par.me);
+//  }
 
   // shift energy delta values
@@ -872 +877 @@
  * As wave functions are stored in the reciprocal basis set, the application is straight-forward,
  * for every G vector, the by \a index specified component is multiplied with the respective
- * coefficient. Afterwards, 1/i is applied by flipping real and imaginary components (and an additional minus sign on the new imaginary term).
+ * coefficient. Afterwards, 1/i is applied by flipping real and imaginary components and an additional minus sign on the new imaginary term.
  * \param *P Problem at hand
  * \param *source complex coefficients of wave function \f$\varphi(G)\f$
@@ -1495 +1500 @@
   result += Elambda;
   //fprintf(stderr,"(%i) 1st: Elambda = \t%lg\n", P->Par.me, Elambda);
-
+ 
   E1 =  2.*GradSP(P,LevS,ConDir,H1c_grad) * sqrt(Psi->AllPsiStatus[ActNum].PsiFactor*Psi->LocalPsiStatus[R->ActualLocalPsiNo].PsiFactor);
   result += E1;
@@ -1532 +1537 @@
   struct RunStruct *R = &P->R;
   struct Psis *Psi = &P->Lat.Psi;
-  //struct Lattice *Lat = &P->Lat;
-  //struct Energy *E = Lat->E;
+  struct Lattice *Lat = &P->Lat;
+  struct Energy *E = Lat->E;
   double result = 0.;
-  double Con0 = 0., Elambda = 0.;//, E0 = 0., E1 = 0.;
-  //int i;
+  double Con0 = 0., Elambda = 0., Elambda2 = 0., E0 = 0., E1 = 0.;
+  int i;
   const int state = R->CurrentMin;
   //const int l_normal = R->ActualLocalPsiNo - Psi->TypeStartIndex[state] + Psi->TypeStartIndex[Occupied]; 
@@ -1543 +1548 @@
   Con0 = 2.*ConDirHConDir;
   result += Con0;
-  ////E0 = -2.*sourceHsource;
-  ////result += E0;
-  ////E1 = -E->PsiEnergy[Perturbed1_0Energy][R->ActualLocalPsiNo] - E->PsiEnergy[Perturbed0_1Energy][R->ActualLocalPsiNo];
-  ////result += E1;
+  E0 = -2.*sourceHsource;
+  result += E0;
+  E1 = -E->PsiEnergy[Perturbed1_0Energy][R->ActualLocalPsiNo] - E->PsiEnergy[Perturbed0_1Energy][R->ActualLocalPsiNo];
+  //result += E1;
   //fprintf(stderr,"(%i) 2nd: E1 = \t%lg\n", P->Par.me, E1);
-
-  ////for (i=0;i<Lat->Psi.NoOfPsis;i++) {
-  ////  if (i != ActNum) Elambda += Psi->lambda[i][ActNum]*Psi->Overlap[i][ActNum]+ Psi->lambda[ActNum][i]*Psi->Overlap[ActNum][i];   // overlap contains PsiFactor
-  ////}
-  ////Elambda = Psi->lambda[ActNum][ActNum]*Psi->Overlap[ActNum][ActNum];
-  Elambda = 2.*Psi->lambda[ActNum][ActNum]*ConDirConDir;
-  result -= Elambda;
-
-  //fprintf(stderr,"(%i) 2ndPerturbedDerivative: Result = Con0 + E0 + E1 + Elambda + dEdt0_XC + ddEddt0_XC + dEdt0_H + ddEddt0_H = %lg + %lg + %lg + %lg + %lg + %lg + %lg + %lg = %lg\n", P->Par.me, Con0, E0, E1, Elambda, VolumeFactorR*dEdt0_XC, VolumeFactorR*ddEddt0_XC, dEdt0_H, ddEddt0_H, result);
+  
+  for (i=0;i<Lat->Psi.NoOfPsis;i++) {
+    if (i != ActNum) 
+      Elambda += Psi->lambda[i][ActNum]*Psi->Overlap[i][ActNum]+ Psi->lambda[ActNum][i]*Psi->Overlap[ActNum][i];   // overlap contains PsiFactor
+  }
+  Elambda = Psi->lambda[ActNum][ActNum]*Psi->Overlap[ActNum][ActNum];
+  result += Elambda;
+  Elambda2 = 2.*Psi->lambda[ActNum][ActNum]*ConDirConDir;
+  result -= Elambda2;
+
+  //fprintf(stderr,"(%i) 2ndPerturbedDerivative: Result = Con0 + E0 + E1 + Elambda + Elambda2 = %lg + %lg + %lg + %lg + %lg = %lg\n", P->Par.me, Con0, E0, E1, Elambda, Elambda2, result);
 
   return (result);
@@ -2206 +2213 @@
   }
   
-  int wished = -1;
+  int CurrentOrbital = -1;
   FILE *file =  fopen(P->Call.MainParameterFile,"r");
-  if (!ParseForParameter(0,file,"Orbital",0,1,1,int_type,&wished, 1, optional)) {
+  if (!ParseForParameter(0,file,"Orbital",0,1,1,int_type,&CurrentOrbital, 1, optional)) {
      if (P->Call.out[ReadOut]) fprintf(stderr,"Desired Orbital missing, using: All!\n");
-    wished = -1;
-  } else if (wished != -1) {
-     if (P->Call.out[ReadOut]) fprintf(stderr,"Desired Orbital is: %i.\n", wished);
+    CurrentOrbital = -1;
+  } else if (CurrentOrbital != -1) {
+     if (P->Call.out[ReadOut]) fprintf(stderr,"Desired Orbital is: %i.\n", CurrentOrbital);
   } else { 
      if (P->Call.out[ReadOut]) fprintf(stderr,"Desired Orbital is: All.\n");
@@ -2220 +2227 @@
    // Commence grid filling
   for (k=Psi->TypeStartIndex[Occupied];k<Psi->TypeStartIndex[Occupied+1];k++)  // every local wave functions adds up its part of the current
-    if ((k + P->Par.me_comm_ST_PsiT*(Psi->TypeStartIndex[UnOccupied]-Psi->TypeStartIndex[Occupied]) == wished) || (wished == -1)) {  // compare with global number
+    if ((k + P->Par.me_comm_ST_PsiT*(Psi->TypeStartIndex[UnOccupied]-Psi->TypeStartIndex[Occupied]) == CurrentOrbital) || (CurrentOrbital == -1)) {  // compare with global number
     if (P->Call.out[StepLeaderOut]) fprintf(stderr,"(%i)Calculating Current Density Summand of type %s for Psi (%i/%i) ... \n", P->Par.me, R->MinimisationName[type], Psi->LocalPsiStatus[k].MyGlobalNo, k);
     //ActNum = k - Psi->TypeStartIndex[Occupied] + Psi->TypeStartIndex[1] * Psi->LocalPsiStatus[k].my_color_comm_ST_Psi; // global number of unperturbed Psi
@@ -2278 +2285 @@
               for (n[2]=0;n[2]<N[2];n[2]++) {
                 i0 = n[2]+N[2]*(n[1]+N[1]*n0);
-                n[0]=n0 + N0*myPE; // global relative coordinate: due to partitoning of x-axis in PEPGamma>1 case
+                n[0]=n0 + N0*myPE; // global relative coordinate: due to partitioning of x-axis in PEPGamma>1 case
                 fac[0] = (double)n[0]/(double)N[0];
                 fac[1] = (double)n[1]/(double)N[1];
@@ -3500 +3507 @@
         for (i=0;i<NDIM;i++)
           fprintf(stderr,"\t%lg",gsl_vector_get(eval,i));
+        fprintf(stderr,"\n");
       }
       if (P->Call.out[ValueOut]) {
         if (P->Call.out[ReadOut])
-          fprintf(stderr,"moment  : %e\n", iso);
+          fprintf(stderr,"moment     : %e\n", iso);
         else
           fprintf(stderr,"%e\n", iso);
       }
-      if (P->Call.out[ReadOut]) {
+      if (P->Call.out[ReadOut]) { 
         fprintf(stderr,"anisotropy : %e\n", delta_moment);
         fprintf(stderr,"asymmetry  : %e\n", eta);