GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
dmperm.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2005-2013 David Bateman
4 Copyright (C) 1998-2005 Andy Adler
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include "defun-dld.h"
29 #include "error.h"
30 #include "gripes.h"
31 #include "oct-obj.h"
32 #include "utils.h"
33 
34 #include "oct-sparse.h"
35 #include "ov-re-sparse.h"
36 #include "ov-cx-sparse.h"
37 #include "SparseQR.h"
38 #include "SparseCmplxQR.h"
39 
40 #ifdef USE_64_BIT_IDX_T
41 #define CXSPARSE_NAME(name) cs_dl ## name
42 #else
43 #define CXSPARSE_NAME(name) cs_di ## name
44 #endif
45 
46 static RowVector
48 {
49  RowVector ret (n);
50  for (octave_idx_type i = 0; i < n; i++)
51  ret.xelem (i) = p[i] + 1;
52  return ret;
53 }
54 
55 #if HAVE_CXSPARSE
56 static octave_value_list
57 dmperm_internal (bool rank, const octave_value arg, int nargout)
58 {
59  octave_value_list retval;
60  octave_idx_type nr = arg.rows ();
61  octave_idx_type nc = arg.columns ();
62  SparseMatrix m;
64  CXSPARSE_NAME () csm;
65  csm.m = nr;
66  csm.n = nc;
67  csm.x = 0;
68  csm.nz = -1;
69 
70  if (arg.is_real_type ())
71  {
72  m = arg.sparse_matrix_value ();
73  csm.nzmax = m.nnz ();
74  csm.p = m.xcidx ();
75  csm.i = m.xridx ();
76  }
77  else
78  {
79  cm = arg.sparse_complex_matrix_value ();
80  csm.nzmax = cm.nnz ();
81  csm.p = cm.xcidx ();
82  csm.i = cm.xridx ();
83  }
84 
85  if (!error_state)
86  {
87  if (nargout <= 1 || rank)
88  {
89 #if defined(CS_VER) && (CS_VER >= 2)
90  octave_idx_type *jmatch = CXSPARSE_NAME (_maxtrans) (&csm, 0);
91 #else
92  octave_idx_type *jmatch = CXSPARSE_NAME (_maxtrans) (&csm);
93 #endif
94  if (rank)
95  {
96  octave_idx_type r = 0;
97  for (octave_idx_type i = 0; i < nc; i++)
98  if (jmatch[nr+i] >= 0)
99  r++;
100  retval(0) = static_cast<double>(r);
101  }
102  else
103  retval(0) = put_int (jmatch + nr, nc);
104  CXSPARSE_NAME (_free) (jmatch);
105  }
106  else
107  {
108 #if defined(CS_VER) && (CS_VER >= 2)
109  CXSPARSE_NAME (d) *dm = CXSPARSE_NAME(_dmperm) (&csm, 0);
110 #else
111  CXSPARSE_NAME (d) *dm = CXSPARSE_NAME(_dmperm) (&csm);
112 #endif
113 
114  //retval(5) = put_int (dm->rr, 5);
115  //retval(4) = put_int (dm->cc, 5);
116 #if defined(CS_VER) && (CS_VER >= 2)
117  retval(3) = put_int (dm->s, dm->nb+1);
118  retval(2) = put_int (dm->r, dm->nb+1);
119  retval(1) = put_int (dm->q, nc);
120  retval(0) = put_int (dm->p, nr);
121 #else
122  retval(3) = put_int (dm->S, dm->nb+1);
123  retval(2) = put_int (dm->R, dm->nb+1);
124  retval(1) = put_int (dm->Q, nc);
125  retval(0) = put_int (dm->P, nr);
126 #endif
127  CXSPARSE_NAME (_dfree) (dm);
128  }
129  }
130  return retval;
131 }
132 #endif
133 
134 DEFUN_DLD (dmperm, args, nargout,
135  "-*- texinfo -*-\n\
136 @deftypefn {Loadable Function} {@var{p} =} dmperm (@var{S})\n\
137 @deftypefnx {Loadable Function} {[@var{p}, @var{q}, @var{r}, @var{S}] =} dmperm (@var{S})\n\
138 \n\
139 @cindex Dulmage-Mendelsohn decomposition\n\
140 Perform a Dulmage-Mendelsohn permutation of the sparse matrix @var{S}.\n\
141 With a single output argument @code{dmperm} performs the row permutations\n\
142 @var{p} such that @code{@var{S}(@var{p},:)} has no zero elements on the\n\
143 diagonal.\n\
144 \n\
145 Called with two or more output arguments, returns the row and column\n\
146 permutations, such that @code{@var{S}(@var{p}, @var{q})} is in block\n\
147 triangular form. The values of @var{r} and @var{S} define the boundaries\n\
148 of the blocks. If @var{S} is square then @code{@var{r} == @var{S}}.\n\
149 \n\
150 The method used is described in: A. Pothen & C.-J. Fan. @cite{Computing the\n\
151 Block Triangular Form of a Sparse Matrix}. ACM Trans. Math. Software,\n\
152 16(4):303-324, 1990.\n\
153 @seealso{colamd, ccolamd}\n\
154 @end deftypefn")
155 {
156  int nargin = args.length ();
157  octave_value_list retval;
158 
159  if (nargin != 1)
160  {
161  print_usage ();
162  return retval;
163  }
164 
165 #if HAVE_CXSPARSE
166  retval = dmperm_internal (false, args(0), nargout);
167 #else
168  error ("dmperm: not available in this version of Octave");
169 #endif
170 
171  return retval;
172 }
173 
174 /*
175 %!testif HAVE_CXSPARSE
176 %! n = 20;
177 %! a = speye (n,n);
178 %! a = a(randperm (n),:);
179 %! assert (a(dmperm (a),:), speye (n));
180 
181 %!testif HAVE_CXSPARSE
182 %! n = 20;
183 %! d = 0.2;
184 %! a = tril (sprandn (n,n,d), -1) + speye (n,n);
185 %! a = a(randperm (n), randperm (n));
186 %! [p,q,r,s] = dmperm (a);
187 %! assert (tril (a(p,q), -1), sparse (n, n));
188 */
189 
190 DEFUN_DLD (sprank, args, nargout,
191  "-*- texinfo -*-\n\
192 @deftypefn {Loadable Function} {@var{p} =} sprank (@var{S})\n\
193 @cindex structural rank\n\
194 \n\
195 Calculate the structural rank of the sparse matrix @var{S}. Note that\n\
196 only the structure of the matrix is used in this calculation based on\n\
197 a Dulmage-Mendelsohn permutation to block triangular form. As such the\n\
198 numerical rank of the matrix @var{S} is bounded by\n\
199 @code{sprank (@var{S}) >= rank (@var{S})}. Ignoring floating point errors\n\
200 @code{sprank (@var{S}) == rank (@var{S})}.\n\
201 @seealso{dmperm}\n\
202 @end deftypefn")
203 {
204  int nargin = args.length ();
205  octave_value_list retval;
206 
207  if (nargin != 1)
208  {
209  print_usage ();
210  return retval;
211  }
212 
213 #if HAVE_CXSPARSE
214  retval = dmperm_internal (true, args(0), nargout);
215 #else
216  error ("sprank: not available in this version of Octave");
217 #endif
218 
219  return retval;
220 }
221 
222 /*
223 %!testif HAVE_CXSPARSE
224 %! assert (sprank (speye (20)), 20)
225 %!testif HAVE_CXSPARSE
226 %! assert (sprank ([1,0,2,0;2,0,4,0]), 2)
227 
228 %!error sprank (1,2)
229 */