+static void assign_constraint(struct gmx_lincsdata *li,
+ int constraint_index,
+ int a1, int a2,
+ real lenA, real lenB,
+ const t_blocka *at2con)
+{
+ int con;
+
+ con = li->nc;
+
+ /* Make an mapping of local topology constraint index to LINCS index */
+ li->con_index[constraint_index] = con;
+
+ li->bllen0[con] = lenA;
+ li->ddist[con] = lenB - lenA;
+ /* Set the length to the topology A length */
+ li->bllen[con] = lenA;
+ li->bla[2*con] = a1;
+ li->bla[2*con+1] = a2;
+
+ /* Make space in the constraint connection matrix for constraints
+ * connected to both end of the current constraint.
+ */
+ li->ncc +=
+ at2con->index[a1 + 1] - at2con->index[a1] - 1 +
+ at2con->index[a2 + 1] - at2con->index[a2] - 1;
+
+ li->blnr[con + 1] = li->ncc;
+
+ /* Increase the constraint count */
+ li->nc++;
+}
+
+/* Check if constraint with topology index constraint_index is connected
+ * to other constraints, and if so add those connected constraints to our task.
+ */
+static void check_assign_connected(struct gmx_lincsdata *li,
+ const t_iatom *iatom,
+ const t_idef *idef,
+ int bDynamics,
+ int a1, int a2,
+ const t_blocka *at2con)
+{
+ /* Currently this function only supports constraint groups
+ * in which all constraints share at least one atom
+ * (e.g. H-bond constraints).
+ * Check both ends of the current constraint for
+ * connected constraints. We need to assign those
+ * to the same task.
+ */
+ int end;
+
+ for (end = 0; end < 2; end++)
+ {
+ int a, k;
+
+ a = (end == 0 ? a1 : a2);
+
+ for (k = at2con->index[a]; k < at2con->index[a + 1]; k++)
+ {
+ int cc;
+
+ cc = at2con->a[k];
+ /* Check if constraint cc has not yet been assigned */
+ if (li->con_index[cc] == -1)
+ {
+ int type;
+ real lenA, lenB;
+
+ type = iatom[cc*3];
+ lenA = idef->iparams[type].constr.dA;
+ lenB = idef->iparams[type].constr.dB;
+
+ if (bDynamics || lenA != 0 || lenB != 0)
+ {
+ assign_constraint(li, cc, iatom[3*cc + 1], iatom[3*cc + 2], lenA, lenB, at2con);
+ }
+ }
+ }
+ }
+}
+
+/* Check if constraint with topology index constraint_index is involved
+ * in a constraint triangle, and if so add the other two constraints
+ * in the triangle to our task.
+ */
+static void check_assign_triangle(struct gmx_lincsdata *li,
+ const t_iatom *iatom,
+ const t_idef *idef,
+ int bDynamics,
+ int constraint_index,
+ int a1, int a2,
+ const t_blocka *at2con)
+{
+ int nca, cc[32], ca[32], k;
+ int c_triangle[2] = { -1, -1 };
+
+ nca = 0;
+ for (k = at2con->index[a1]; k < at2con->index[a1 + 1]; k++)
+ {
+ int c;
+
+ c = at2con->a[k];
+ if (c != constraint_index)
+ {
+ int aa1, aa2;
+
+ aa1 = iatom[c*3 + 1];
+ aa2 = iatom[c*3 + 2];
+ if (aa1 != a1)
+ {
+ cc[nca] = c;
+ ca[nca] = aa1;
+ nca++;
+ }
+ if (aa2 != a1)
+ {
+ cc[nca] = c;
+ ca[nca] = aa2;
+ nca++;
+ }
+ }
+ }
+
+ for (k = at2con->index[a2]; k < at2con->index[a2 + 1]; k++)
+ {
+ int c;
+
+ c = at2con->a[k];
+ if (c != constraint_index)
+ {
+ int aa1, aa2, i;
+
+ aa1 = iatom[c*3 + 1];
+ aa2 = iatom[c*3 + 2];
+ if (aa1 != a2)
+ {
+ for (i = 0; i < nca; i++)
+ {
+ if (aa1 == ca[i])
+ {
+ c_triangle[0] = cc[i];
+ c_triangle[1] = c;
+ }
+ }
+ }
+ if (aa2 != a2)
+ {
+ for (i = 0; i < nca; i++)
+ {
+ if (aa2 == ca[i])
+ {
+ c_triangle[0] = cc[i];
+ c_triangle[1] = c;
+ }
+ }
+ }
+ }
+ }
+
+ if (c_triangle[0] >= 0)
+ {
+ int end;
+
+ for (end = 0; end < 2; end++)
+ {
+ /* Check if constraint c_triangle[end] has not yet been assigned */
+ if (li->con_index[c_triangle[end]] == -1)
+ {
+ int i, type;
+ real lenA, lenB;
+
+ i = c_triangle[end]*3;
+ type = iatom[i];
+ lenA = idef->iparams[type].constr.dA;
+ lenB = idef->iparams[type].constr.dB;
+
+ if (bDynamics || lenA != 0 || lenB != 0)
+ {
+ assign_constraint(li, c_triangle[end], iatom[i + 1], iatom[i + 2], lenA, lenB, at2con);
+ }
+ }
+ }
+ }
+}
+
+static void set_matrix_indices(struct gmx_lincsdata *li,
+ const lincs_task_t *li_task,
+ const t_blocka *at2con,
+ gmx_bool bSortMatrix)
+{
+ int b;
+
+ for (b = li_task->b0; b < li_task->b1; b++)
+ {
+ int a1, a2, i, k;
+
+ a1 = li->bla[b*2];
+ a2 = li->bla[b*2 + 1];
+
+ i = li->blnr[b];
+ for (k = at2con->index[a1]; k < at2con->index[a1 + 1]; k++)
+ {
+ int concon;
+
+ concon = li->con_index[at2con->a[k]];
+ if (concon != b)
+ {
+ li->blbnb[i++] = concon;
+ }
+ }
+ for (k = at2con->index[a2]; k < at2con->index[a2 + 1]; k++)
+ {
+ int concon;
+
+ concon = li->con_index[at2con->a[k]];
+ if (concon != b)
+ {
+ li->blbnb[i++] = concon;
+ }
+ }
+
+ if (bSortMatrix)
+ {
+ /* Order the blbnb matrix to optimize memory access */
+ qsort(&(li->blbnb[li->blnr[b]]), li->blnr[b + 1] - li->blnr[b],
+ sizeof(li->blbnb[0]), int_comp);
+ }
+ }
+}
+