<p>The reason why the problem persists might be simply because your <code>lib3j6j9j.a</code> does not include necessary files (such as d1mach). Actually, we can compile necessary files rather directly, so I will summarize the procedure below:</p>

<p>1) Download <code>drc3jm.f</code> (which calculates 3j-symbols) and dependencies from the Netlib/Slatec page (<a href="http://www.netlib.org/slatec/src/" rel="nofollow">here</a> or <a href="http://www.netlib.org/cgi-bin/netlibfiles.pl?filename=/slatec/src/drc3jm.f" rel="nofollow">here</a>). Unpack the archive file to get Fortran files (*.f).</p>

<pre><code>tar xvf netlibfiles.tgz
</code></pre>

<p>2) Remove <code>d1mach.f</code>, <code>i1mach.f</code>, and <code>r1mach.f</code> (if any). Instead, download their alternative versions from Netlib/blas (*):</p>

<pre><code>rm -f i1mach.f r1mach.f d1mach.f
wget http://www.netlib.org/blas/i1mach.f 
wget http://www.netlib.org/blas/r1mach.f
wget http://www.netlib.org/blas/d1mach.f 
</code></pre>

<p>3) Compile all *.f files</p>

<pre><code>gfortran testf.f90 *.f
</code></pre>

<p>together with a main program testf.f90 (in free-format), e.g.,</p>

<pre><code>program main
implicit none
integer, parameter :: N = 1000
double precision coef( N ), M2min, M2max, M2
integer ier
ier = 0 ; coef(:) = 0.0d0

call DRC3JM( 15.0d0, 30.0d0, 40.d0, 2.0d0,  M2min, M2max, coef, N, ier )
print *, "M2min, M2max, ier = ", M2min, M2max, ier

M2 = 2.0d0
print "(a, f20.15)", "coef = ", coef( nint(M2 - M2min+1) )  !! -0.019081579799192
end
</code></pre>

<p>Then running the executable gives the desired result.</p>

<hr>

<p>3-a) We can also make these *.f as a library and link with C++ codes, e.g., as follows:</p>

<pre><code>gfortran -c *.f
ar rv mylib.a *.o
g++ testc.cpp mylib.a -lgfortran
</code></pre>

<p>with a main program (testc.cpp)</p>

<pre><code>#include &lt;cstdio&gt;
extern "C"
double drc3jm_ (double*, double*, double*, 
                double*, double*, double*, double*, int*, int*);

int main()
{
    double* coef;
    double L1, L2, L3, M1, M2min, M2max, M2;
    int ier, k, N = 1000;

    coef = new double [ N ];
    L1 = 15.0; L2 = 30.0; L3 = 40.0; M1 = 2.0;

    drc3jm_ ( &amp;L1, &amp;L2,    &amp;L3,
              &amp;M1, &amp;M2min, &amp;M2max, coef, &amp;N, &amp;ier );  
    printf( "M2min, M2max, ierr = %10.5f%10.5f%d\n", M2min, M2max, ier );

    M2 = 2.0;   
    k = (int)(M2 - M2min + 1.0e-3);
    printf( "coef = %20.15f\n", coef[ k ] );  // -0.019081579799192
    return 0;
}
</code></pre>

<p>We can see that the two programs give the same coefficient (-0.019081579799192) for</p>

<pre><code>j1=15, j2=30, j3=40, m1=2, m2=2, m3=-4
</code></pre>

<p>You can also get the same result with an online tool, e.g., <a href="http://www-stone.ch.cam.ac.uk/cgi-bin/wigner.cgi?symbol=3j&amp;j1=15&amp;j2=30&amp;j3=40&amp;m1=2&amp;m2=2&amp;m3=-4" rel="nofollow">here</a>.</p>

<hr>

<p>But depending on cases, it may be simpler to use other libraries. One approach is to use the corresponding GSL routines (<a href="https://www.gnu.org/software/gsl/manual/html_node/3_002dj-Symbols.html" rel="nofollow">here</a>) as</p>

<pre><code>#include &lt;cstdio&gt;
extern "C"
double gsl_sf_coupling_3j (int two_ja, int two_jb, int two_jc,
                           int two_ma, int two_mb, int two_mc);  
int main()
{
    double coef;
    coef = gsl_sf_coupling_3j( 30, 60, 80, 4, 4, -8 );  // -0.019081579799205
    // NOTE: all j's and m's need to be doubled.

    printf( "coef = %20.15f\n", coef );
    return 0;
}
</code></pre>

<p>Here you need to link necessary GSL libraries (e.g., <code>g++ test.cpp -lgsl</code> or <code>g++ test.cpp /usr/lib64/libgsl.so.0 /usr/lib64/libgslcblas.so.0</code> etc). </p>

<hr>

<p>Yet another approach is to use a latest program <a href="http://fy.chalmers.se/subatom/wigxjpf/" rel="nofollow">WIGXJPF</a> (the related paper is <a href="http://arxiv.org/abs/1504.08329" rel="nofollow">here</a>). I tried this a bit and it seems extremely easy to install (only one <code>make</code>) and use. For example, enter the <code>example/</code> directory and try <code>gcc -I../inc csimple.c ../lib/libwigxjpf.a</code>.
According to the above paper, this program may offer some accuracy and performance advantage.</p>

<hr>

<p>(*) For more details, please see <a href="http://www.netlib.org/misc/faq.html#2.17" rel="nofollow">the Netlib/FAQ page</a> (thanks to @VladimirF in the comment). We could utilize the original d1mach.f etc in Slatec, but we need to modify them so as to obtain correct machine-dependent constants. The above BLAS versions of d1mach.f etc handle this automatically, so they are more convenient. </p>
undefined reference to 'd1mach_'
<p>I'm trying to link a fortran subroutine with c++, but can't quite figure out what exactly is wrong here:
The fortran subroutine calls some functions eg. d1mach or xermsg, which aren't defined in the fortran subroutine but called externally. When compiling, the error is "undefined reference to d1mach_"(or xermsg). 
I tried linking a library I think might contain the said functions (There seems to be a file called d1mach.o and xermsg.o inside the library) but the same error still persists. What might I be doing wrong?</p>

<pre><code>extern"C" {
void drc3jm_(double *L1,double *L2,double *L3,double *M1,double *M2MIN,
    double *M2MAX,double *THRCOF,int *NDIM,int *IER);
}
</code></pre>

<p>This is the function I use to call the subroutine, and haven't used any new headers beside iostream</p>

<pre><code>*DECK DRC3JM
      SUBROUTINE DRC3JM (L1, L2, L3, M1, M2MIN, M2MAX, THRCOF, NDIM,
     +   IER)
CALL XERMSG('SLATEC','DRC3JM','L1-ABS(M1) less than zero or '//
     +      'L1+ABS(M1) not integer.',IER,1)
</code></pre>

<p>This is the declaration of the fortran subroutine which calls the undeclared function xermsg.</p>

<p>I link the library using the -L/<em>path</em>/lib instruction but to no avail.
The subroutine is to calculate a mathematical function and is part of the slatec codes.</p>

<p>Please let me know what other information you might need.</p>