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
|
//------------------------------------------------------------------------------
// GraphBLAS/Demo/Program/complex_demo.c: demo for user-defined complex type
//------------------------------------------------------------------------------
// SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2020, All Rights Reserved.
// http://suitesparse.com See GraphBLAS/Doc/License.txt for license.
//------------------------------------------------------------------------------
// This demo illustrates the creation and use of a user-defined type for double
// complex matrices and vectors. Run the output of this program in MATLAB
// to check the results.
#include "graphblas_demos.h"
//------------------------------------------------------------------------------
// print a complex matrix
//------------------------------------------------------------------------------
// when printed, 1 is added to all row indices so the results can be
// checked in MATLAB
void print_complex_matrix (GrB_Matrix A, char *name)
{
GrB_Index nrows, ncols, nentries ;
GrB_Matrix_nvals (&nentries, A) ;
GrB_Matrix_nrows (&nrows, A) ;
GrB_Matrix_ncols (&ncols, A) ;
printf (
"\n%% GraphBLAS matrix %s: nrows: %.16g ncols %.16g entries: %.16g\n",
name, (double) nrows, (double) ncols, (double) nentries) ;
GrB_Index *I = (GrB_Index *) malloc (MAX (nentries,1) * sizeof (GrB_Index));
GrB_Index *J = (GrB_Index *) malloc (MAX (nentries,1) * sizeof (GrB_Index));
GxB_FC64_t *X = (GxB_FC64_t *)
malloc (MAX (nentries,1) * sizeof (GxB_FC64_t)) ;
if (Complex == GxB_FC64)
{
GxB_Matrix_extractTuples_FC64 (I, J, X, &nentries, A) ;
}
else
{
GrB_Matrix_extractTuples_UDT (I, J, X, &nentries, A) ;
}
printf ("%s = sparse (%.16g,%.16g) ;\n", name,
(double) nrows, (double) ncols) ;
for (int64_t k = 0 ; k < nentries ; k++)
{
printf (" %s (%.16g,%.16g) = (%20.16g) + (%20.16g)*1i ;\n",
name, (double) (1 + I [k]), (double) (1 + J [k]),
creal (X [k]), cimag (X [k])) ;
}
printf ("%s\n", name) ;
free (I) ;
free (J) ;
free (X) ;
}
//------------------------------------------------------------------------------
// C = A*B for complex matrices
//------------------------------------------------------------------------------
int main (int argc, char **argv)
{
GrB_Index m = 3, k = 5, n = 4 ;
GrB_Matrix A, B, C ;
// initialize GraphBLAS and create the user-defined Complex type
GrB_Info info ;
GrB_init (GrB_NONBLOCKING) ;
int nthreads ;
GxB_Global_Option_get (GxB_GLOBAL_NTHREADS, &nthreads) ;
fprintf (stderr, "complex_demo: nthreads: %d\n", nthreads) ;
bool predefined = (argc > 1) ;
if (predefined)
{
fprintf (stderr, "Using pre-defined GxB_FC64 complex type\n") ;
}
else
{
fprintf (stderr, "Using user-defined Complex type\n") ;
}
info = Complex_init (predefined) ;
if (info != GrB_SUCCESS)
{
fprintf (stderr, "Complex init failed: %s\n", GrB_error ( )) ;
abort ( ) ;
}
// generate random matrices A and B
simple_rand_seed (1) ;
random_matrix (&A, false, false, m, k, 6, 0, true) ;
random_matrix (&B, false, false, k, n, 8, 0, true) ;
GxB_Matrix_fprint (A, "C", GxB_SHORT, stderr) ;
GxB_Matrix_fprint (B, "C", GxB_SHORT, stderr) ;
// C = A*B
GrB_Matrix_new (&C, Complex, m, n) ;
GrB_mxm (C, NULL, NULL, Complex_plus_times, A, B, NULL) ;
GxB_Matrix_fprint (C, "C", GxB_SHORT, stderr) ;
// print the results
printf ("\n%% run this output of this program as a script in MATLAB:\n") ;
print_complex_matrix (A, "A") ;
print_complex_matrix (B, "B") ;
print_complex_matrix (C, "C") ;
printf ("E = A*B\n") ;
printf ("err = norm (C-E,1)\n") ;
printf ("assert (err < 1e-12)\n") ;
// free all matrices
GrB_Matrix_free (&A) ;
GrB_Matrix_free (&B) ;
GrB_Matrix_free (&C) ;
// free the Complex types, operators, monoids, and semiring
Complex_finalize ( ) ;
// finalize GraphBLAS
GrB_finalize ( ) ;
}
//------------------------------------------------------------------------------
#if 0
int main ( )
{
printf ("complex data type not available (ANSI C11 or higher required)\n") ;
}
#endif
|
