|
5164
|
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 |
|
|
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 } |
