1 //------------------------------------------------------------------------------
2 // GB_mex_complex: convert a real matrix into a complex one
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 // If A is real, C has an all-zero imaginary part.
11 // If A is complex, then C = A.
12
13 // This is a sparse version of the MATLAB 'complex' function, which does not
14 // work for sparse matrices. This is self-contained and does not use GraphBLAS
15 // at all.
16
17 #include "mex.h"
18 #include "matrix.h"
19 #include <stdlib.h>
20 #include <string.h>
21 #include <stdint.h>
22
23 #define USAGE "C = GB_mex_complex (A)"
24
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])25 void mexFunction
26 (
27 int nargout,
28 mxArray *pargout [ ],
29 int nargin,
30 const mxArray *pargin [ ]
31 )
32 {
33
34 // check inputs
35 if (nargout > 1 || nargin != 1)
36 {
37 mexErrMsgTxt ("Usage: " USAGE) ;
38 }
39
40 // get the input matrix
41 const mxArray *A = pargin [0] ;
42 bool A_is_sparse = mxIsSparse (A) ;
43 if (mxIsLogical (A))
44 {
45 mexErrMsgTxt ("A must be double or double complex") ;
46 }
47
48 int64_t *Ap, *Ai ;
49
50 if (A_is_sparse)
51 {
52 Ap = (int64_t *) mxGetJc (A) ;
53 Ai = (int64_t *) mxGetIr (A) ;
54 }
55 else
56 {
57 Ap = NULL ;
58 Ai = NULL ;
59 }
60
61 double *Ax = NULL ;
62 if (mxIsComplex (A))
63 {
64 Ax = (double *) mxGetComplexDoubles (pargin [0]) ;
65 }
66 else
67 {
68 Ax = (double *) mxGetDoubles (pargin [0]) ;
69 }
70
71 int64_t m = mxGetM (A) ;
72 int64_t n = mxGetN (A) ;
73 int64_t anz = (A_is_sparse) ? Ap [n] : (m*n) ;
74
75 // create the output matrix
76 if (A_is_sparse)
77 {
78 // A and C are sparse
79 pargout [0] = mxCreateSparse (m, n, anz+1, mxCOMPLEX) ;
80 int64_t *Cp = (int64_t *) mxGetJc (pargout [0]) ;
81 int64_t *Ci = (int64_t *) mxGetIr (pargout [0]) ;
82 // copy the pattern of A into C
83 memcpy (Cp, Ap, (n+1) * sizeof (int64_t)) ;
84 memcpy (Ci, Ai, anz * sizeof (int64_t)) ;
85 }
86 else
87 {
88 // A and C are full
89 pargout [0] = mxCreateDoubleMatrix (m, n, mxCOMPLEX) ;
90 }
91
92 // copy the values of A into C
93 double *Cx = (double *) mxGetComplexDoubles (pargout [0]) ;
94 if (mxIsComplex (A))
95 {
96 memcpy (Cx, Ax, anz * 2 * sizeof (double)) ;
97 }
98 else
99 {
100 for (int64_t k = 0 ; k < anz ; k++)
101 {
102 Cx [2*k ] = Ax [k] ;
103 Cx [2*k+1] = 0 ;
104 }
105 }
106 }
107
108