1 /* This software was developed by Bruce Hendrickson and Robert Leland   *
2  * at Sandia National Laboratories under US Department of Energy        *
3  * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
4 
5 #include <stdio.h>
6 #include "structs.h"
7 #include "defs.h"
8 
9 static struct ipairs *pedges;	/* perturbed edges */
10 static double *pvals;		/* perturbed values */
11 
12 /* Inititialize the perturbation */
perturb_init(n)13 void      perturb_init(n)
14 int       n;			/* graph size at this level */
15 {
16     extern int NPERTURB;	/* number of edges to perturb */
17     extern double PERTURB_MAX;	/* maximum perturbation */
18     int       i, j;		/* loop counter */
19     double   *smalloc();
20     double    drandom();
21 
22     /* Initialize the diagonal perturbation weights */
23     pedges = (struct ipairs *) smalloc((unsigned) NPERTURB * sizeof(struct ipairs));
24     pvals = (double *) smalloc((unsigned) NPERTURB * sizeof(double));
25 
26     if (n <= 1) {
27 	for (i = 0; i < NPERTURB; i++) {
28 	    pedges[i].val1 = pedges[i].val2 = 0;
29 	    pvals[i] = 0;
30 	}
31 	return;
32     }
33 
34     for (i = 0; i < NPERTURB; i++) {
35 	pedges[i].val1 = 1 + (n * drandom());
36 
37 	/* Find another vertex to define an edge. */
38 	j = 1 + (n * drandom());
39 	while (j == i)
40 	    j = 1 + (n * drandom());
41 	pedges[i].val2 = 1 + (n * drandom());
42 
43 	pvals[i] = PERTURB_MAX * drandom();
44     }
45 }
46 
perturb_clear()47 void      perturb_clear()
48 {
49     int       sfree();
50 
51     sfree((char *) pedges);
52     sfree((char *) pvals);
53     pedges = NULL;
54     pvals = NULL;
55 }
56 
57 
58 /* Modify the result of splarax to break any graph symmetry */
perturb(result,vec)59 void      perturb(result, vec)
60 double   *result;		/* result of matrix-vector multiply */
61 double   *vec;			/* vector matrix multiplies */
62 {
63     extern int NPERTURB;	/* number of edges to perturb */
64     int       i;		/* loop counter */
65 
66     for (i = 0; i < NPERTURB; i++) {
67 	result[pedges[i].val1] +=
68 	   pvals[i] * (vec[pedges[i].val2] - vec[pedges[i].val1]);
69 	result[pedges[i].val2] +=
70 	   pvals[i] * (vec[pedges[i].val1] - vec[pedges[i].val2]);
71     }
72 }
73 
74 
75 /* Modify the result of splarax to break any graph symmetry, float version */
perturb_float(result,vec)76 void      perturb_float(result, vec)
77 float    *result;		/* result of matrix-vector multiply */
78 float    *vec;			/* vector matrix multiplies */
79 {
80     extern int NPERTURB;	/* number of edges to perturb */
81     int       i;		/* loop counter */
82 
83     for (i = 0; i < NPERTURB; i++) {
84 	result[pedges[i].val1] +=
85 	   pvals[i] * (vec[pedges[i].val2] - vec[pedges[i].val1]);
86 	result[pedges[i].val2] +=
87 	   pvals[i] * (vec[pedges[i].val1] - vec[pedges[i].val2]);
88     }
89 }
90