in tools/ADOLC/main.cpp [406:506]
double compute_ba_J(int nruns, double time_limit, int n, int m, int p,
double *cams, double *X, double *w, int *obs, double *feats,
double *reproj_err, double *w_err, BASparseMat *J, double* t_sparsity)
{
if (nruns == 0)
return 0.;
bool doRowCompression = false;
high_resolution_clock::time_point start, end;
int tapeTag = 1;
start = high_resolution_clock::now();
adouble *acams, *aX, *aw, *areproj_err,
*af_prior_err, *aw_err;
areproj_err = new adouble[2 * p];
af_prior_err = new adouble[n - 2];
aw_err = new adouble[p];
// Record on a tape
trace_on(tapeTag);
acams = new adouble[BA_NCAMPARAMS*n];
for (int i = 0; i < BA_NCAMPARAMS*n; i++)
acams[i] <<= cams[i];
aX = new adouble[3 * m];
for (int i = 0; i < 3 * m; i++)
aX[i] <<= X[i];
aw = new adouble[p];
for (int i = 0; i < p; i++)
aw[i] <<= w[i];
ba_objective(n, m, p, acams, aX, aw, obs, feats,
areproj_err, aw_err);
for (int i = 0; i < 2 * p; i++)
areproj_err[i] >>= reproj_err[i];
for (int i = 0; i < p; i++)
aw_err[i] >>= w_err[i];
trace_off();
delete[] acams;
delete[] aX;
delete[] aw;
delete[] areproj_err;
delete[] af_prior_err;
delete[] aw_err;
end = high_resolution_clock::now();
double t_tape = duration_cast<duration<double>>(end - start).count();
//////// Compute J and compute sparsity always again and again
vector<double> in(J->ncols);
memcpy(&in[0], cams, BA_NCAMPARAMS*n*sizeof(double));
int off = BA_NCAMPARAMS*n;
memcpy(&in[off], X, 3 * m*sizeof(double));
off += 3 * m;
memcpy(&in[off], w, p*sizeof(double));
int opt[4];
opt[0] = 0; // default
opt[1] = 0; // default
opt[2] = 0; // 0=auto 1=F 2=R
opt[3] = doRowCompression ? 1 : 0;
int nnz;
unsigned int *ridxs = nullptr, *cidxs = nullptr;
double *nzvals = nullptr;
int samePattern = 0;
start = high_resolution_clock::now();
sparse_jac(tapeTag, J->nrows, J->ncols, samePattern,
in.data(), &nnz, &ridxs, &cidxs, &nzvals, opt);
end = high_resolution_clock::now();
*t_sparsity = duration_cast<duration<double>>(end - start).count();
samePattern = 1;
start = high_resolution_clock::now();
for (int i = 0; i < nruns; i++)
{
sparse_jac(tapeTag, J->nrows, J->ncols, samePattern,
in.data(), &nnz, &ridxs, &cidxs, &nzvals, opt);
}
end = high_resolution_clock::now();
double t_J = duration_cast<duration<double>>(end - start).count() / nruns;
convert_J(nnz, ridxs, cidxs, nzvals, J);
delete[] ridxs;
delete[] cidxs;
delete[] nzvals;
cout << "t_tape: " << t_tape << endl;
cout << "t_sparsity: " << *t_sparsity << endl;
cout << "t_J:" << t_J << endl;
return t_J;
}