1 //------------------------------------------------------------------------------
2 // GraphBLAS/Demo/Program/complex_demo.c: demo for user-defined complex type
3 //------------------------------------------------------------------------------
4
5 // SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
6 // SPDX-License-Identifier: Apache-2.0
7
8 //------------------------------------------------------------------------------
9
10 // This demo illustrates the creation and use of a user-defined type for double
11 // complex matrices and vectors. Run the output of this program in MATLAB
12 // to check the results.
13
14 #include "graphblas_demos.h"
15
16 //------------------------------------------------------------------------------
17 // print a complex matrix
18 //------------------------------------------------------------------------------
19
20 // when printed, 1 is added to all row indices so the results can be
21 // checked in MATLAB
22
print_complex_matrix(GrB_Matrix A,char * name)23 void print_complex_matrix (GrB_Matrix A, char *name)
24 {
25 GrB_Index nrows, ncols, nentries ;
26
27 GrB_Matrix_nvals (&nentries, A) ;
28 GrB_Matrix_nrows (&nrows, A) ;
29 GrB_Matrix_ncols (&ncols, A) ;
30
31 printf (
32 "\n%% GraphBLAS matrix %s: nrows: %.16g ncols %.16g entries: %.16g\n",
33 name, (double) nrows, (double) ncols, (double) nentries) ;
34
35 GrB_Index *I = (GrB_Index *) malloc (MAX (nentries,1) * sizeof (GrB_Index));
36 GrB_Index *J = (GrB_Index *) malloc (MAX (nentries,1) * sizeof (GrB_Index));
37 GxB_FC64_t *X = (GxB_FC64_t *)
38 malloc (MAX (nentries,1) * sizeof (GxB_FC64_t)) ;
39
40 if (Complex == GxB_FC64)
41 {
42 GxB_Matrix_extractTuples_FC64 (I, J, X, &nentries, A) ;
43 }
44 else
45 {
46 GrB_Matrix_extractTuples_UDT (I, J, X, &nentries, A) ;
47 }
48
49 printf ("%s = sparse (%.16g,%.16g) ;\n", name,
50 (double) nrows, (double) ncols) ;
51 for (int64_t k = 0 ; k < nentries ; k++)
52 {
53 printf (" %s (%.16g,%.16g) = (%20.16g) + (%20.16g)*1i ;\n",
54 name, (double) (1 + I [k]), (double) (1 + J [k]),
55 creal (X [k]), cimag (X [k])) ;
56 }
57 printf ("%s\n", name) ;
58
59 free (I) ;
60 free (J) ;
61 free (X) ;
62 }
63
64 //------------------------------------------------------------------------------
65 // C = A*B for complex matrices
66 //------------------------------------------------------------------------------
67
main(int argc,char ** argv)68 int main (int argc, char **argv)
69 {
70 GrB_Index m = 3, k = 5, n = 4 ;
71 GrB_Matrix A, B, C ;
72
73 // initialize GraphBLAS and create the user-defined Complex type
74 GrB_Info info ;
75 GrB_init (GrB_NONBLOCKING) ;
76 int nthreads ;
77 GxB_Global_Option_get (GxB_GLOBAL_NTHREADS, &nthreads) ;
78 fprintf (stderr, "complex_demo: nthreads: %d\n", nthreads) ;
79
80 bool predefined = (argc > 1) ;
81 if (predefined)
82 {
83 fprintf (stderr, "Using pre-defined GxB_FC64 complex type\n") ;
84 }
85 else
86 {
87 fprintf (stderr, "Using user-defined Complex type\n") ;
88 }
89
90 info = Complex_init (predefined) ;
91 if (info != GrB_SUCCESS)
92 {
93 fprintf (stderr, "Complex init failed\n") ;
94 abort ( ) ;
95 }
96
97 // generate random matrices A and B
98 simple_rand_seed (1) ;
99 random_matrix (&A, false, false, m, k, 6, 0, true) ;
100 random_matrix (&B, false, false, k, n, 8, 0, true) ;
101
102 GxB_Matrix_fprint (A, "C", GxB_SHORT, stderr) ;
103 GxB_Matrix_fprint (B, "C", GxB_SHORT, stderr) ;
104
105 // C = A*B
106 GrB_Matrix_new (&C, Complex, m, n) ;
107 GrB_mxm (C, NULL, NULL, Complex_plus_times, A, B, NULL) ;
108
109 GxB_Matrix_fprint (C, "C", GxB_SHORT, stderr) ;
110
111 // print the results
112 printf ("\n%% run this output of this program as a script in MATLAB:\n") ;
113 print_complex_matrix (A, "A") ;
114 print_complex_matrix (B, "B") ;
115 print_complex_matrix (C, "C") ;
116
117 printf ("E = A*B\n") ;
118 printf ("err = norm (C-E,1)\n") ;
119 printf ("assert (err < 1e-12)\n") ;
120
121 // free all matrices
122 GrB_Matrix_free (&A) ;
123 GrB_Matrix_free (&B) ;
124 GrB_Matrix_free (&C) ;
125
126 // free the Complex types, operators, monoids, and semiring
127 Complex_finalize ( ) ;
128
129 // finalize GraphBLAS
130 GrB_finalize ( ) ;
131 }
132
133 //------------------------------------------------------------------------------
134
135 #if 0
136
137 int main ( )
138 {
139 printf ("complex data type not available (ANSI C11 or higher required)\n") ;
140 }
141
142 #endif
143
144