Changeset 487 for trunk/matml

Timestamp:

10/12/2009 02:31:51 PM (3 years ago)

Author:

powell

Message:

New function adjust_gaussians, needs lapack, bump libfreenergy version.

Location:

trunk/matml/src/ternary

Files:

• ChangeLog (modified) (1 diff)
• configure.in (modified) (1 diff)
• freenergy.c (modified) (1 diff)
• freenergy.h (modified) (1 diff)
• Makefile.am (modified) (1 diff)
• ternary.tex.in (modified) (1 diff)

trunk/matml/src/ternary/ChangeLog

@@ -3 +3 @@
   * New arbitrary polynomial terms in the free energy functions.
   * Fixed and considerably improved spinodal caculation.
+  * New adjust_gaussians function to aid in setting Gaussian parameters.
 
  --

trunk/matml/src/ternary/configure.in

@@ -28 +28 @@
   AC_MSG_ERROR([Cannot find qhull library which ternary requires; you can download it at http://www.qhull.org/ or install the libqhull-dev .deb package.]))
 
+MATH_BLASLAPACK_CHECKS
+
 dnl Documentation tools
 CXREF_LATEX_CHECKS

trunk/matml/src/ternary/freenergy.c

@@ -592 +592 @@
   return 0;
 }
+
+
+#if HAVE_LAPACK == yes
+#include <stdlib.h> /*+ For malloc and free +*/
+#include <stdio.h> /*+ For printf +*/
+void dgetrf_ (int *, int *, double *, int *, int *, int *);
+void dlaswp_ (int *, double *, int *, int *, int *, int *, int *);
+void dtrsm_ (char *, char *, char *, char *, int *, int *, double *,
+             double *, int *, double *, int *);
+
+/*++++++++++++++++++++++++++++++++++++++
+  Adjust the Gaussian energy parameters (G) such that on return, the free
+  energy at each Gaussian concentration is equal to the incoming value of G0.
+
+  int adjust_gaussians Retunrs 0 or an error code.
+
+  double T Temperature.
+
+  double P Pressure (ignored for now).
+
+  double *G0 Intended free energy values at Gaussian center compositions.
+
+  energy_params *eparams Free energy function parameters.
+
+  int efunc Index indicating which set of energy parameters to use
+  ++++++++++++++++++++++++++++++++++++++*/
+
+int adjust_gaussians
+(double T, double P, double *G0, energy_params *eparams, int efunc)
+{
+  double *G_scratch, *G_matrix, one=1.;
+  int i,j, gs = eparams [efunc].n_gauss, piv[gs];
+
+  // Do one malloc for both arrays
+  if (!(G_scratch = malloc (gs*(gs+1) * sizeof (double))))
+    return 1;
+  G_matrix = G_scratch + gs;
+
+  // Calculate the free energy at each point without any Gaussians.
+  for (i=0; i<gs; i++)
+    (eparams[efunc].gauss)[i].G = 0.;
+  for (i=0; i<gs; i++)
+    G_scratch [i] = free_energy
+      ((eparams[efunc].gauss)[i].C2, (eparams[efunc].gauss)[i].C3, T, P,
+       eparams+efunc);
+
+  // Calculate the effect of each Gaussian on the energy at all of the other
+  // points.
+  for (i=0; i<gs; i++)
+    {
+      (eparams[efunc].gauss)[i].G = 1.;
+      for (j=0; j<gs; j++)
+        G_matrix [i*gs+j] = free_energy
+          ((eparams[efunc].gauss)[j].C2, (eparams[efunc].gauss)[j].C3, T, P,
+           eparams+efunc) - G_scratch [j];
+      (eparams[efunc].gauss)[i].G = 0.;
+    }
+
+  // Set G_scratch values to G0 for solving below
+  for (i=0; i<gs; i++)
+    G_scratch [i] = G0 [i] - G_scratch [i];
+
+  // Solve for the new G values values using LAPACK
+  dgetrf_ (&gs, &gs, G_matrix, &gs, piv, &i);
+  if (i != 0)
+    return i;
+  i=1;
+  dlaswp_ (&i, G_scratch, &gs, &i, &gs, piv, &i);
+  dtrsm_("Left", "Lower", "No transpose", "Unit", &gs, &i, &one, G_matrix, &gs,
+         G_scratch, &gs);
+  dtrsm_("Left", "Upper", "No transpose", "Non-Unit", &gs, &i, &one, G_matrix,
+         &gs, G_scratch, &gs);
+
+  // Put them into the Gaussian structures
+  for (i=0; i<gs; i++)
+    (eparams[efunc].gauss)[i].G = G_scratch [i];
+
+#ifdef DEBUG
+  // Test the results
+  for (i=0; i<gs; i++)
+    printf ("Point %d at %g,%g: G=%g, energy=%g, G0=%g\n", i,
+            (eparams[efunc].gauss)[i].C2, (eparams[efunc].gauss)[i].C3,
+            (eparams[efunc].gauss)[i].G,
+            free_energy ((eparams[efunc].gauss)[i].C2,
+                         (eparams[efunc].gauss)[i].C3, T, P, eparams+efunc),
+            G0 [i]);
+#endif
+
+  // Clean up and go home
+  free(G_scratch);
+  return 0;
+}
+#endif

trunk/matml/src/ternary/freenergy.h

@@ -111 +111 @@
 (energy_point *points, int n, double T, double P, energy_params *eparams,
  int efunc, free_energy_flags flags);
+#if HAVE_LAPACK == yes
+int adjust_gaussians
+(double T, double P, double *G0, energy_params *eparams, int efunc);
+#endif
 
 #endif // FREENERGY_H

trunk/matml/src/ternary/Makefile.am

@@ -12 +12 @@
 
 libfreenergy_la_SOURCES = freenergy.c
-libfreenergy_la_LIBADD  = -lm
-libfreenergy_la_CFLAGS  = -std=c99 #-DDEBUG
-libfreenergy_la_LDFLAGS = -version-info 0:0:0
+libfreenergy_la_LIBADD  = -lm $(BLASLAPACK_LIBS)
+libfreenergy_la_CFLAGS  = -std=c99 -DHAVE_LAPACK=$(HAVE_LAPACK) #-DDEBUG
+libfreenergy_la_LDFLAGS = -version-info 1:0:0
 
 libgibbs_la_SOURCES = geomview.c qhull.c book.c spinodal.c

trunk/matml/src/ternary/ternary.tex.in

@@ -49 +49 @@
 
 To compile {\tt Ternary}, you need the qhull library and geomview program, as
-well as a C compiler.  On Ubuntu, installing the {\tt build-essential}, {\tt
-  libqhull-dev} and {\tt geomview} packages satisfy this.  If you want to build
-the documentation as well, you will also need {\tt cxref} along with \Latex (in
-the {\tt texlive-base-bin} package) for the PDF, and either hevea or latex2html
-for HTML documentation.
+well as a C compiler.  On Debian and Ubuntu, installing the {\tt
+  build-essential}, {\tt liblapack-dev}, {\tt libqhull-dev} and {\tt geomview}
+packages satisfy this.  If you want to build the documentation as well, you
+will also need {\tt cxref} along with \LaTeX (in the {\tt texlive-latex-base}
+package) for the PDF, and either {\tt hevea} or {\tt latex2html} for HTML
+documentation.
 
-Next, unpack the .tar.gz source archive by typing (at the command line):
+Next, unpack the {\tt .tar.gz} source archive by typing (at the command line):
 \begin{quote}
   \tt tar xzf Ternary-[version].tar.gz