1 /* ========================================================================= */
2 /* === AMD mexFunction ===================================================== */
3 /* ========================================================================= */
4
5 /* ------------------------------------------------------------------------- */
6 /* AMD, Copyright (c) Timothy A. Davis, */
7 /* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
8 /* email: DrTimothyAldenDavis@gmail.com */
9 /* ------------------------------------------------------------------------- */
10
11 /*
12 * Usage:
13 * p = amd (A)
14 * p = amd (A, Control)
15 * [p, Info] = amd (A)
16 * [p, Info] = amd (A, Control)
17 * Control = amd ; % return the default Control settings for AMD
18 * amd ; % print the default Control settings for AMD
19 *
20 * Given a square matrix A, compute a permutation P suitable for a Cholesky
21 * factorization of the matrix B (P,P), where B = spones (A) + spones (A').
22 * The method used is the approximate minimum degree ordering method. See
23 * amd.m and amd.h for more information.
24 *
25 * The input matrix need not have sorted columns, and can have duplicate
26 * entries.
27 */
28
29 #include "amd.h"
30 #include "mex.h"
31 #include "matrix.h"
32 #define Long SuiteSparse_long
33
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])34 void mexFunction
35 (
36 int nargout,
37 mxArray *pargout [ ],
38 int nargin,
39 const mxArray *pargin [ ]
40 )
41 {
42 Long i, m, n, *Ap, *Ai, *P, nc, result, spumoni, full ;
43 double *Pout, *InfoOut, Control [AMD_CONTROL], Info [AMD_INFO], *ControlIn ;
44 mxArray *A ;
45
46 /* --------------------------------------------------------------------- */
47 /* get control parameters */
48 /* --------------------------------------------------------------------- */
49
50 spumoni = 0 ;
51 if (nargin == 0)
52 {
53 /* get the default control parameters, and return */
54 pargout [0] = mxCreateDoubleMatrix (AMD_CONTROL, 1, mxREAL) ;
55 amd_l_defaults (mxGetPr (pargout [0])) ;
56 if (nargout == 0)
57 {
58 amd_l_control (mxGetPr (pargout [0])) ;
59 }
60 return ;
61 }
62
63 amd_l_defaults (Control) ;
64 if (nargin > 1)
65 {
66 ControlIn = mxGetPr (pargin [1]) ;
67 nc = mxGetM (pargin [1]) * mxGetN (pargin [1]) ;
68 Control [AMD_DENSE]
69 = (nc > 0) ? ControlIn [AMD_DENSE] : AMD_DEFAULT_DENSE ;
70 Control [AMD_AGGRESSIVE]
71 = (nc > 1) ? ControlIn [AMD_AGGRESSIVE] : AMD_DEFAULT_AGGRESSIVE ;
72 spumoni = (nc > 2) ? (ControlIn [2] != 0) : 0 ;
73 }
74
75 if (spumoni > 0)
76 {
77 amd_l_control (Control) ;
78 }
79
80 /* --------------------------------------------------------------------- */
81 /* get inputs */
82 /* --------------------------------------------------------------------- */
83
84 if (nargout > 2 || nargin > 2)
85 {
86 mexErrMsgTxt ("Usage: p = amd (A)\nor [p, Info] = amd (A, Control)") ;
87 }
88
89 A = (mxArray *) pargin [0] ;
90 n = mxGetN (A) ;
91 m = mxGetM (A) ;
92 if (spumoni > 0)
93 {
94 mexPrintf (" input matrix A is %d-by-%d\n", m, n) ;
95 }
96 if (mxGetNumberOfDimensions (A) != 2)
97 {
98 mexErrMsgTxt ("amd: A must be 2-dimensional") ;
99 }
100 if (m != n)
101 {
102 mexErrMsgTxt ("amd: A must be square") ;
103 }
104
105 /* --------------------------------------------------------------------- */
106 /* allocate workspace for output permutation */
107 /* --------------------------------------------------------------------- */
108
109 P = mxMalloc ((n+1) * sizeof (Long)) ;
110
111 /* --------------------------------------------------------------------- */
112 /* if A is full, convert to a sparse matrix */
113 /* --------------------------------------------------------------------- */
114
115 full = !mxIsSparse (A) ;
116 if (full)
117 {
118 if (spumoni > 0)
119 {
120 mexPrintf (
121 " input matrix A is full (sparse copy of A will be created)\n");
122 }
123 mexCallMATLAB (1, &A, 1, (mxArray **) pargin, "sparse") ;
124 }
125 Ap = (Long *) mxGetJc (A) ;
126 Ai = (Long *) mxGetIr (A) ;
127 if (spumoni > 0)
128 {
129 mexPrintf (" input matrix A has %d nonzero entries\n", Ap [n]) ;
130 }
131
132 /* --------------------------------------------------------------------- */
133 /* order the matrix */
134 /* --------------------------------------------------------------------- */
135
136 result = amd_l_order (n, Ap, Ai, P, Control, Info) ;
137
138 /* --------------------------------------------------------------------- */
139 /* if A is full, free the sparse copy of A */
140 /* --------------------------------------------------------------------- */
141
142 if (full)
143 {
144 mxDestroyArray (A) ;
145 }
146
147 /* --------------------------------------------------------------------- */
148 /* print results (including return value) */
149 /* --------------------------------------------------------------------- */
150
151 if (spumoni > 0)
152 {
153 amd_l_info (Info) ;
154 }
155
156 /* --------------------------------------------------------------------- */
157 /* check error conditions */
158 /* --------------------------------------------------------------------- */
159
160 if (result == AMD_OUT_OF_MEMORY)
161 {
162 mexErrMsgTxt ("amd: out of memory") ;
163 }
164 else if (result == AMD_INVALID)
165 {
166 mexErrMsgTxt ("amd: input matrix A is corrupted") ;
167 }
168
169 /* --------------------------------------------------------------------- */
170 /* copy the outputs to MATLAB */
171 /* --------------------------------------------------------------------- */
172
173 /* output permutation, P */
174 pargout [0] = mxCreateDoubleMatrix (1, n, mxREAL) ;
175 Pout = mxGetPr (pargout [0]) ;
176 for (i = 0 ; i < n ; i++)
177 {
178 Pout [i] = P [i] + 1 ; /* change to 1-based indexing for MATLAB */
179 }
180 mxFree (P) ;
181
182 /* Info */
183 if (nargout > 1)
184 {
185 pargout [1] = mxCreateDoubleMatrix (AMD_INFO, 1, mxREAL) ;
186 InfoOut = mxGetPr (pargout [1]) ;
187 for (i = 0 ; i < AMD_INFO ; i++)
188 {
189 InfoOut [i] = Info [i] ;
190 }
191 }
192 }
193