1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
|
/* spdgemm.c
*
* Copyright (C) 2014 Patrick Alken
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_spblas.h>
#include <gsl/gsl_errno.h>
/*
gsl_spblas_dgemm()
Multiply two sparse matrices
Inputs: alpha - scalar factor
A - sparse matrix
B - sparse matrix
C - (output) C = alpha * A * B
Return: success or error
Notes:
1) based on CSparse routine cs_multiply
*/
int
gsl_spblas_dgemm(const double alpha, const gsl_spmatrix *A,
const gsl_spmatrix *B, gsl_spmatrix *C)
{
if (A->size2 != B->size1 || A->size1 != C->size1 || B->size2 != C->size2)
{
GSL_ERROR("matrix dimensions do not match", GSL_EBADLEN);
}
else if (A->sptype != B->sptype || A->sptype != C->sptype)
{
GSL_ERROR("matrix storage formats do not match", GSL_EINVAL);
}
else if (!GSL_SPMATRIX_ISCCS(A))
{
GSL_ERROR("compressed column format required", GSL_EINVAL);
}
else
{
int status = GSL_SUCCESS;
const size_t M = A->size1;
const size_t N = B->size2;
size_t *Bi = B->i;
size_t *Bp = B->p;
double *Bd = B->data;
size_t *w = (size_t *) A->work; /* workspace of length M */
double *x = (double *) C->work; /* workspace of length M */
size_t *Cp, *Ci;
double *Cd;
size_t j, p;
size_t nz = 0;
if (C->nzmax < A->nz + B->nz)
{
status = gsl_spmatrix_realloc(A->nz + B->nz, C);
if (status)
{
GSL_ERROR("unable to realloc matrix C", status);
}
}
/* initialize workspace to 0 */
for (j = 0; j < M; ++j)
w[j] = 0;
Cp = C->p;
Ci = C->i;
Cd = C->data;
for (j = 0; j < N; ++j)
{
if (nz + M > C->nzmax)
{
status = gsl_spmatrix_realloc(2 * C->nzmax + M, C);
if (status)
{
GSL_ERROR("unable to realloc matrix C", status);
}
/* these pointers could have changed due to reallocation */
Ci = C->i;
Cd = C->data;
}
Cp[j] = nz; /* column j of C starts here */
for (p = Bp[j]; p < Bp[j + 1]; ++p)
{
nz = gsl_spblas_scatter(A, Bi[p], Bd[p], w, x, j + 1, C, nz);
}
for (p = Cp[j]; p < nz; ++p)
Cd[p] = x[Ci[p]];
}
Cp[N] = nz;
C->nz = nz;
/* scale by alpha */
gsl_spmatrix_scale(C, alpha);
return status;
}
} /* gsl_spblas_dgemm() */
/*
gsl_spblas_scatter()
Keep a running total x -> x + alpha*A(:,j) for adding matrices together in CCS,
which will eventually be stored in C(:,j)
When a new non-zero element with row index i is found, update C->i with
the row index. C->data is updated only by the calling function after all
matrices have been added via this function.
Inputs: A - sparse matrix m-by-n
j - column index
alpha - scalar factor
w - keeps track which rows of column j have been added to C;
initialize to 0 prior to first call
x - column vector of length m
mark -
C - output matrix whose jth column will be added to A(:,j)
nz - (input/output) number of non-zeros in matrix C
Notes:
1) This function is designed to be called successively when adding multiple
matrices together. Column j of C is stored contiguously as per CCS but not
necessarily in order - ie: the row indices C->i may not be in ascending order.
2) based on CSparse routine cs_scatter
*/
size_t
gsl_spblas_scatter(const gsl_spmatrix *A, const size_t j, const double alpha,
size_t *w, double *x, const size_t mark, gsl_spmatrix *C,
size_t nz)
{
size_t p;
size_t *Ai = A->i;
size_t *Ap = A->p;
double *Ad = A->data;
size_t *Ci = C->i;
for (p = Ap[j]; p < Ap[j + 1]; ++p)
{
size_t i = Ai[p]; /* A(i,j) is nonzero */
if (w[i] < mark) /* check if row i has been stored in column j yet */
{
w[i] = mark; /* i is new entry in column j */
Ci[nz++] = i; /* add i to pattern of C(:,j) */
x[i] = alpha * Ad[p]; /* x(i) = alpha * A(i,j) */
}
else /* this (i,j) exists in C from a previous call */
{
x[i] += alpha * Ad[p]; /* add alpha*A(i,j) to C(i,j) */
}
}
return (nz) ;
} /* gsl_spblas_scatter() */
|
