1 /* -*- mode: C -*- */
2 /*
3 IGraph library.
4 Copyright (C) 2009-2012 Gabor Csardi <csardi.gabor@gmail.com>
5 334 Harvard st, Cambridge MA, 02139 USA
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20 02110-1301 USA
21
22 */
23
24 #define NCOMPLEX /* to make it compile with MSVC on Windows */
25
26 #include <cs/cs.h>
27 #include <igraph.h>
28 #include "linalg/blas_internal.h"
29 #include "linalg/arpack_internal.h"
30
31 #include "test_utilities.inc"
32
igraph_matrix_dgemv(const igraph_matrix_t * m,const igraph_vector_t * v,igraph_vector_t * res,igraph_real_t alpha,igraph_real_t beta,igraph_bool_t transpose_m)33 int igraph_matrix_dgemv(const igraph_matrix_t *m,
34 const igraph_vector_t *v,
35 igraph_vector_t *res,
36 igraph_real_t alpha,
37 igraph_real_t beta,
38 igraph_bool_t transpose_m) {
39
40 int nrow = igraph_matrix_nrow(m);
41 int ncol = igraph_matrix_ncol(m);
42 long int vlen = igraph_vector_size(v);
43 int one = 1;
44 char t = transpose_m ? 't' : 'n';
45 long int input_len = transpose_m ? nrow : ncol;
46 long int output_len = transpose_m ? ncol : nrow;
47
48 if (vlen != input_len) {
49 IGRAPH_ERROR("Matrix and vector sizes are incompatible", IGRAPH_EINVAL);
50 }
51
52 if (beta != 0 && igraph_vector_size(res) != output_len) {
53 IGRAPH_ERROR("Non-zero beta and bad `res' vector size, possible mistake",
54 IGRAPH_EINVAL);
55 }
56
57 IGRAPH_CHECK(igraph_vector_resize(res, output_len));
58
59 igraphdgemv_(&t, &nrow, &ncol, &alpha, &MATRIX(*m, 0, 0),
60 &nrow, VECTOR(*v), &one, &beta, VECTOR(*res), &one);
61
62 return 0;
63 }
64
igraph_matrix_vector_prod(const igraph_matrix_t * m,const igraph_vector_t * v,igraph_vector_t * res)65 int igraph_matrix_vector_prod(const igraph_matrix_t *m,
66 const igraph_vector_t *v,
67 igraph_vector_t *res) {
68 return igraph_matrix_dgemv(m, v, res, 1.0, 0.0, /*transpose=*/ 0);
69 }
70
my_dgemv(const igraph_matrix_t * m,const igraph_vector_t * v,igraph_vector_t * res,igraph_real_t alpha,igraph_real_t beta,igraph_bool_t transpose_m)71 int my_dgemv(const igraph_matrix_t *m,
72 const igraph_vector_t *v,
73 igraph_vector_t *res,
74 igraph_real_t alpha,
75 igraph_real_t beta,
76 igraph_bool_t transpose_m) {
77
78 int nrow = igraph_matrix_nrow(m);
79 int ncol = igraph_matrix_ncol(m);
80 long int vlen = igraph_vector_size(v);
81 int one = 1;
82 char t = transpose_m ? 't' : 'n';
83 long int input_len = transpose_m ? nrow : ncol;
84 long int output_len = transpose_m ? ncol : nrow;
85
86 if (vlen != input_len) {
87 IGRAPH_ERROR("Matrix and vector sizes are incompatible", IGRAPH_EINVAL);
88 }
89
90 if (beta != 0 && igraph_vector_size(res) != output_len) {
91 IGRAPH_ERROR("Non-zero beta and bad `res' vector size, possible mistake",
92 IGRAPH_EINVAL);
93 }
94
95 IGRAPH_CHECK(igraph_vector_resize(res, output_len));
96
97 igraphdgemv_(&t, &nrow, &ncol, &alpha, &MATRIX(*m, 0, 0),
98 &nrow, VECTOR(*v), &one, &beta, VECTOR(*res), &one);
99
100 return 0;
101 }
102
my_gaxpy(const igraph_matrix_t * m,const igraph_vector_t * v,igraph_vector_t * res)103 int my_gaxpy(const igraph_matrix_t *m,
104 const igraph_vector_t *v,
105 igraph_vector_t *res) {
106 return my_dgemv(m, v, res, 1.0, 0.0, /*transpose=*/ 0);
107 }
108
my_dgemm(const igraph_matrix_t * m1,const igraph_matrix_t * m2,igraph_matrix_t * res)109 int my_dgemm(const igraph_matrix_t *m1,
110 const igraph_matrix_t *m2,
111 igraph_matrix_t *res) {
112
113 long int m1_r = igraph_matrix_nrow(m1);
114 long int m1_c = igraph_matrix_ncol(m1);
115 long int m2_r = igraph_matrix_nrow(m2);
116 long int m2_c = igraph_matrix_ncol(m2);
117 long int i, j, k;
118
119 if (m1_c != m2_r) {
120 IGRAPH_ERROR("Cannot multiply matrices, invalid dimensions", IGRAPH_EINVAL);
121 }
122
123 IGRAPH_CHECK(igraph_matrix_resize(res, m1_r, m2_c));
124 igraph_matrix_null(res);
125
126 for (i = 0; i < m1_r; i++) {
127 for (j = 0; j < m2_c; j++) {
128 for (k = 0; k < m1_c /* which is also m2_r*/; k++) {
129 MATRIX(*res, i, j) += MATRIX(*m1, i, k) * MATRIX(*m2, k, j);
130 }
131 }
132 }
133
134 return 0;
135 }
136
check_same(const igraph_sparsemat_t * A,const igraph_matrix_t * M)137 igraph_bool_t check_same(const igraph_sparsemat_t *A,
138 const igraph_matrix_t *M) {
139
140 long int nrow = igraph_sparsemat_nrow(A);
141 long int ncol = igraph_sparsemat_ncol(A);
142 long int j, p, nzero = 0;
143
144 if (nrow != igraph_matrix_nrow(M) ||
145 ncol != igraph_matrix_ncol(M)) {
146 return 0;
147 }
148
149 for (j = 0; j < A->cs->n; j++) {
150 for (p = A->cs->p[j]; p < A->cs->p[j + 1]; p++) {
151 long int to = A->cs->i[p];
152 igraph_real_t value = A->cs->x[p];
153 if (value != MATRIX(*M, to, j)) {
154 return 0;
155 }
156 nzero += 1;
157 }
158 }
159
160 for (j = 0; j < nrow; j++) {
161 for (p = 0; p < ncol; p++) {
162 if (MATRIX(*M, j, p) != 0) {
163 nzero -= 1;
164 }
165 }
166 }
167
168 return nzero == 0;
169 }
170
main()171 int main() {
172
173 igraph_sparsemat_t A, B, C, D;
174 igraph_vector_t v, w, x, y;
175 igraph_matrix_t M, N, O;
176 long int i;
177
178 srand(1);
179
180 /* Matrix-vector product */
181 #define NROW 10
182 #define NCOL 5
183 #define EDGES NROW*NCOL/3
184 igraph_matrix_init(&M, NROW, NCOL);
185 igraph_sparsemat_init(&A, NROW, NCOL, EDGES);
186 for (i = 0; i < EDGES; i++) {
187 long int r = RNG_INTEGER(0, NROW - 1);
188 long int c = RNG_INTEGER(0, NCOL - 1);
189 igraph_real_t value = RNG_INTEGER(1, 5);
190 MATRIX(M, r, c) = MATRIX(M, r, c) + value;
191 igraph_sparsemat_entry(&A, r, c, value);
192 }
193 igraph_sparsemat_compress(&A, &B);
194 igraph_sparsemat_destroy(&A);
195
196 igraph_vector_init(&v, NCOL);
197 igraph_vector_init(&w, NCOL);
198 for (i = 0; i < NCOL; i++) {
199 VECTOR(v)[i] = VECTOR(w)[i] = RNG_INTEGER(1, 5);
200 }
201
202 igraph_vector_init(&x, NROW);
203 igraph_vector_init(&y, NROW);
204 my_gaxpy(&M, &v, &x);
205 igraph_vector_null(&y);
206 igraph_sparsemat_gaxpy(&B, &w, &y);
207
208 if (!igraph_vector_all_e(&x, &y)) {
209 return 1;
210 }
211
212 igraph_vector_destroy(&x);
213 igraph_vector_destroy(&y);
214 igraph_vector_destroy(&v);
215 igraph_vector_destroy(&w);
216 igraph_sparsemat_destroy(&B);
217 igraph_matrix_destroy(&M);
218
219 #undef NROW
220 #undef NCOL
221 #undef EDGES
222
223 /* Matrix-matrix product */
224 #define NROW_A 10
225 #define NCOL_A 7
226 #define EDGES_A NROW_A*NCOL_A/3
227 #define NROW_B 7
228 #define NCOL_B 9
229 #define EDGES_B NROW_B*NCOL_B/3
230 igraph_matrix_init(&M, NROW_A, NCOL_A);
231 igraph_sparsemat_init(&A, NROW_A, NCOL_A, EDGES_A);
232 for (i = 0; i < EDGES_A; i++) {
233 long int r = RNG_INTEGER(0, NROW_A - 1);
234 long int c = RNG_INTEGER(0, NCOL_A - 1);
235 igraph_real_t value = RNG_INTEGER(1, 5);
236 MATRIX(M, r, c) = MATRIX(M, r, c) + value;
237 igraph_sparsemat_entry(&A, r, c, value);
238 }
239 igraph_sparsemat_compress(&A, &C);
240 igraph_sparsemat_destroy(&A);
241
242 igraph_matrix_init(&N, NROW_B, NCOL_B);
243 igraph_sparsemat_init(&B, NROW_B, NCOL_B, EDGES_B);
244 for (i = 0; i < EDGES_B; i++) {
245 long int r = RNG_INTEGER(0, NROW_B - 1);
246 long int c = RNG_INTEGER(0, NCOL_B - 1);
247 igraph_real_t value = RNG_INTEGER(1, 5);
248 MATRIX(N, r, c) = MATRIX(N, r, c) + value;
249 igraph_sparsemat_entry(&B, r, c, value);
250 }
251 igraph_sparsemat_compress(&B, &D);
252 igraph_sparsemat_destroy(&B);
253
254 igraph_matrix_init(&O, 0, 0);
255 my_dgemm(&M, &N, &O);
256 igraph_sparsemat_multiply(&C, &D, &A);
257
258 if (! check_same(&A, &O)) {
259 return 2;
260 }
261
262 igraph_sparsemat_destroy(&C);
263 igraph_sparsemat_destroy(&D);
264 igraph_sparsemat_destroy(&A);
265 igraph_matrix_destroy(&M);
266 igraph_matrix_destroy(&N);
267 igraph_matrix_destroy(&O);
268
269 VERIFY_FINALLY_STACK();
270
271 return 0;
272 }
273