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