← Back to branch summary

~ubuntu-branches/debian/squeeze/gmsh/squeeze

« back to all changes in this revision

Viewing changes to contrib/Chaco/eigen/splarax.c

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2009-01-07 16:02:08 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20090107160208-vklhtj69br5yw5bh
Tags: 2.2.6.dfsg-2
* debian/control: fixed lintian warning "debhelper-but-no-misc-depends"
* debian/watch: fixed lintian warning
  "debian-watch-file-should-mangle-version"

Show diffs side-by-side

added added

removed removed

Lines of Context:
contrib/Chaco/eigen/splarax.c
 
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 "defs.h"
 
7
#include "structs.h"
 
8
 
 
9
 
 
10
/* Sparse linked A(matrix) times x(vector), double precision. */
 
11
void      splarax(result, mat, n, vec, vwsqrt, work)
 
12
double   *result;               /* result of matrix vector multiplication */
 
13
struct    vtx_data **mat;       /* graph data structure */
 
14
int       n;                    /* number of rows/columns in matrix */
 
15
double   *vec;                  /* vector being multiplied by matrix */
 
16
double   *vwsqrt;               /* square roots of vertex weights */
 
17
double   *work;                 /* work vector from 1-n */
 
18
{
 
19
    extern int PERTURB;         /* perturb matrix? */
 
20
    extern int NPERTURB;        /* if so, number of edges to perturb */
 
21
    extern double PERTURB_MAX;  /* maximum value of perturbation */
 
22
    struct vtx_data *mat_i;     /* an entry in "mat" */
 
23
    double    sum;              /* sums inner product of matrix-row & vector */
 
24
    int      *colpntr;          /* loops through indices of nonzeros in a row */
 
25
    float    *wgtpntr;          /* loops through values of nonzeros */
 
26
    int       i, j;             /* loop counters */
 
27
    double   *wrkpntr;          /* loops through indices of work vector */
 
28
    double   *vwsqpntr;         /* loops through indices of vwsqrt */
 
29
    double   *vecpntr;          /* loops through indices of vec */
 
30
    double   *respntr;          /* loops through indices of result */
 
31
    int       last_edge;        /* last edge in edge list */
 
32
    void      perturb();
 
33
 
 
34
    if (vwsqrt == NULL) {               /* No vertex weights */
 
35
        if (mat[1]->ewgts == NULL) {    /* No edge weights */
 
36
            respntr = result;
 
37
            for (i = 1; i <= n; i++) {
 
38
                mat_i = mat[i];
 
39
                colpntr = mat_i->edges;
 
40
                last_edge = mat_i->nedges - 1;
 
41
                sum = last_edge * vec[*colpntr++];
 
42
                for (j = last_edge; j; j--) {
 
43
                    sum -= vec[*colpntr++];
 
44
                }
 
45
                *(++respntr) = sum;
 
46
            }
 
47
        }
 
48
        else {                          /* Edge weights */
 
49
            respntr = result;
 
50
            for (i = 1; i <= n; i++) {
 
51
                mat_i = mat[i];
 
52
                colpntr = mat_i->edges;
 
53
                wgtpntr = mat_i->ewgts;
 
54
                sum = 0.0;
 
55
                for (j = mat_i->nedges; j; j--) {
 
56
                    sum -= *wgtpntr++ * vec[*colpntr++];
 
57
                }
 
58
                *(++respntr) = sum;     /* -sum if want -Ax */
 
59
            }
 
60
        }
 
61
        if (PERTURB && NPERTURB > 0 && PERTURB_MAX > 0.0)
 
62
            perturb(result, vec);
 
63
    }
 
64
    else {                              /* Vertex weights */
 
65
        if (vwsqrt != NULL) {
 
66
            wrkpntr = work;
 
67
            vecpntr = vec;
 
68
            vwsqpntr = vwsqrt;
 
69
            for (i = n; i; i--) {
 
70
                *(++wrkpntr) = *(++vecpntr) / *(++vwsqpntr);
 
71
            }
 
72
        }
 
73
 
 
74
        if (mat[1]->ewgts == NULL) {    /* No edge weights. */
 
75
            respntr = result;
 
76
            for (i = 1; i <= n; i++) {
 
77
                mat_i = mat[i];
 
78
                colpntr = mat_i->edges;
 
79
                last_edge = mat_i->nedges - 1;
 
80
                sum = (last_edge) * work[*colpntr++];
 
81
                for (j = last_edge; j; j--) {
 
82
                    sum -= work[*colpntr++];
 
83
                }
 
84
                *(++respntr) = sum;
 
85
            }
 
86
        }
 
87
        else {                          /* Edge weights. */
 
88
            respntr = result;
 
89
            for (i = 1; i <= n; i++) {
 
90
                mat_i = mat[i];
 
91
                colpntr = mat_i->edges;
 
92
                wgtpntr = mat_i->ewgts;
 
93
                sum = 0.0;
 
94
                for (j = mat_i->nedges; j; j--) {
 
95
                    sum -= *wgtpntr++ * work[*colpntr++];
 
96
                }
 
97
                *(++respntr) = sum;     /* -sum if want -Ax */
 
98
            }
 
99
        }
 
100
        if (PERTURB && NPERTURB > 0 && PERTURB_MAX > 0.0)
 
101
            perturb(result, work);
 
102
 
 
103
        if (vwsqrt != NULL) {
 
104
            respntr = result;
 
105
            vwsqpntr = vwsqrt;
 
106
            for (i = n; i; i--) {
 
107
                *(++respntr) /= *(++vwsqpntr);
 
108
            }
 
109
        }
 
110
    }
 
111
}
 
112
 
 
113
/* Sparse linked A(matrix) times x(vector)i, float version. */
 
114
void      splarax_float(result, mat, n, vec, vwsqrt, work)
 
115
float    *result;               /* result of matrix vector multiplication */
 
116
struct    vtx_data **mat;       /* graph data structure */
 
117
int       n;                    /* number of rows/columns in matrix */
 
118
float    *vec;                  /* vector being multiplied by matrix */
 
119
float    *vwsqrt;               /* square roots of vertex weights */
 
120
float    *work;                 /* work vector from 1-n */
 
121
{
 
122
    extern int PERTURB;         /* perturb matrix? */
 
123
    extern int NPERTURB;        /* if so, number of edges to perturb */
 
124
    extern double PERTURB_MAX;  /* maximum value of perturbation */
 
125
    struct vtx_data *mat_i;     /* an entry in "mat" */
 
126
    double    sum;              /* sums inner product of matrix-row & vector */
 
127
    int      *colpntr;          /* loops through indices of nonzeros in a row */
 
128
    float    *wgtpntr;          /* loops through values of nonzeros */
 
129
    int       i, j;             /* loop counters */
 
130
    float    *wrkpntr;          /* loops through indices of work vector */
 
131
    float    *vwsqpntr;         /* loops through indices of vwsqrt */
 
132
    float    *vecpntr;          /* loops through indices of vec */
 
133
    float    *respntr;          /* loops through indices of result */
 
134
    int       last_edge;        /* last edge in edge list */
 
135
    void      perturb_float();
 
136
 
 
137
    if (vwsqrt == NULL) {               /* No vertex weights */
 
138
        if (mat[1]->ewgts == NULL) {    /* No edge weights */
 
139
            respntr = result;
 
140
            for (i = 1; i <= n; i++) {
 
141
                mat_i = mat[i];
 
142
                colpntr = mat_i->edges;
 
143
                last_edge = mat_i->nedges - 1;
 
144
                sum = (last_edge) * vec[*colpntr++];
 
145
                for (j = last_edge; j; j--) {
 
146
                    sum -= vec[*colpntr++];
 
147
                }
 
148
                *(++respntr) = sum;
 
149
            }
 
150
        }
 
151
        else {                          /* Edge weights */
 
152
            respntr = result;
 
153
            for (i = 1; i <= n; i++) {
 
154
                mat_i = mat[i];
 
155
                colpntr = mat_i->edges;
 
156
                wgtpntr = mat_i->ewgts;
 
157
                sum = 0.0;
 
158
                for (j = mat_i->nedges; j; j--) {
 
159
                    sum -= *wgtpntr++ * vec[*colpntr++];
 
160
                }
 
161
                *(++respntr) = sum;     /* -sum if want -Ax */
 
162
            }
 
163
        }
 
164
        if (PERTURB && NPERTURB > 0 && PERTURB_MAX > 0.0)
 
165
            perturb_float(result, vec);
 
166
    }
 
167
    else {                              /* Vertex weights */
 
168
        if (vwsqrt != NULL) {
 
169
            wrkpntr = work;
 
170
            vecpntr = vec;
 
171
            vwsqpntr = vwsqrt;
 
172
            for (i = n; i; i--) {
 
173
                *(++wrkpntr) = *(++vecpntr) / *(++vwsqpntr);
 
174
            }
 
175
        }
 
176
 
 
177
        if (mat[1]->ewgts == NULL) {    /* No edge weights. */
 
178
            respntr = result;
 
179
            for (i = 1; i <= n; i++) {
 
180
                mat_i = mat[i];
 
181
                colpntr = mat_i->edges;
 
182
                last_edge = mat_i->nedges - 1;
 
183
                sum = (last_edge) * work[*colpntr++];
 
184
                for (j = last_edge; j; j--) {
 
185
                    sum -= work[*colpntr++];
 
186
                }
 
187
                *(++respntr) = sum;
 
188
            }
 
189
        }
 
190
        else {                          /* Edge weights. */
 
191
            respntr = result;
 
192
            for (i = 1; i <= n; i++) {
 
193
                mat_i = mat[i];
 
194
                colpntr = mat_i->edges;
 
195
                wgtpntr = mat_i->ewgts;
 
196
                sum = 0.0;
 
197
                for (j = mat_i->nedges; j; j--) {
 
198
                    sum -= *wgtpntr++ * work[*colpntr++];
 
199
                }
 
200
                *(++respntr) = sum;     /* -sum if want -Ax */
 
201
            }
 
202
        }
 
203
        if (PERTURB && NPERTURB > 0 && PERTURB_MAX > 0.0)
 
204
            perturb_float(result, work);
 
205
 
 
206
        if (vwsqrt != NULL) {
 
207
            respntr = result;
 
208
            vwsqpntr = vwsqrt;
 
209
            for (i = n; i; i--) {
 
210
                *(++respntr) /= *(++vwsqpntr);
 
211
            }
 
212
        }
 
213
    }
 
214
}