static double work_volume0;
static M3P work_sout, work_U;

/* How much work is done by changing U x dc */
static double work_integrand (double c)
{
    M3 M, MI, FSUM;
    double volumec;
    M3MULADDDIAG (c, work_U, 1, M);
    M3InV (M, MI, volumec);
    volumec *= work_volume0;
    M3MUL (work_sout, MI, FSUM);
    return (volumec*M3TRPROD(FSUM,work_U));
} /* end work_integrand() */


/* Calculate the work done by external stress in symmetric */
/* deformation space by straight path Romberg integration: */
double work (M3 H0, M3 sout, M3 U)
{
    work_volume0 = M3VOLUME(H0);
    work_sout = sout;
    work_U = U;
    return(qromb(&work_integrand,0.,1.));
} /* end work() */