GNU Octave  4.0.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-2015 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 @nospell{Dulmage-Mendelsohn} decomposition\n\
140 Perform a @nospell{Dulmage-Mendelsohn} permutation of the sparse matrix\n\
141 @var{S}.\n\
142 \n\
143 With a single output argument @code{dmperm} performs the row permutations\n\
144 @var{p} such that @code{@var{S}(@var{p},:)} has no zero elements on the\n\
145 diagonal.\n\
146 \n\
147 Called with two or more output arguments, returns the row and column\n\
148 permutations, such that @code{@var{S}(@var{p}, @var{q})} is in block\n\
149 triangular form. The values of @var{r} and @var{S} define the boundaries\n\
150 of the blocks. If @var{S} is square then @code{@var{r} == @var{S}}.\n\
151 \n\
152 The method used is described in: @nospell{A. Pothen & C.-J. Fan.}\n\
153 @cite{Computing the Block Triangular Form of a Sparse Matrix}.\n\
154 ACM Trans. Math. Software, 16(4):303-324, 1990.\n\
155 @seealso{colamd, ccolamd}\n\
156 @end deftypefn")
157 {
158  int nargin = args.length ();
159  octave_value_list retval;
160 
161  if (nargin != 1)
162  {
163  print_usage ();
164  return retval;
165  }
166 
167 #if HAVE_CXSPARSE
168  retval = dmperm_internal (false, args(0), nargout);
169 #else
170  error ("dmperm: not available in this version of Octave");
171 #endif
172 
173  return retval;
174 }
175 
176 /*
177 %!testif HAVE_CXSPARSE
178 %! n = 20;
179 %! a = speye (n,n);
180 %! a = a(randperm (n),:);
181 %! assert (a(dmperm (a),:), speye (n));
182 
183 %!testif HAVE_CXSPARSE
184 %! n = 20;
185 %! d = 0.2;
186 %! a = tril (sprandn (n,n,d), -1) + speye (n,n);
187 %! a = a(randperm (n), randperm (n));
188 %! [p,q,r,s] = dmperm (a);
189 %! assert (tril (a(p,q), -1), sparse (n, n));
190 */
191 
192 DEFUN_DLD (sprank, args, nargout,
193  "-*- texinfo -*-\n\
194 @deftypefn {Loadable Function} {@var{p} =} sprank (@var{S})\n\
195 @cindex structural rank\n\
196 \n\
197 Calculate the structural rank of the sparse matrix @var{S}.\n\
198 \n\
199 Note that only the structure of the matrix is used in this calculation based\n\
200 on a @nospell{Dulmage-Mendelsohn} permutation to block triangular form. As\n\
201 such the numerical rank of the matrix @var{S} is bounded by\n\
202 @code{sprank (@var{S}) >= rank (@var{S})}. Ignoring floating point errors\n\
203 @code{sprank (@var{S}) == rank (@var{S})}.\n\
204 @seealso{dmperm}\n\
205 @end deftypefn")
206 {
207  int nargin = args.length ();
208  octave_value_list retval;
209 
210  if (nargin != 1)
211  {
212  print_usage ();
213  return retval;
214  }
215 
216 #if HAVE_CXSPARSE
217  retval = dmperm_internal (true, args(0), nargout);
218 #else
219  error ("sprank: not available in this version of Octave");
220 #endif
221 
222  return retval;
223 }
224 
225 /*
226 %!testif HAVE_CXSPARSE
227 %! assert (sprank (speye (20)), 20)
228 %!testif HAVE_CXSPARSE
229 %! assert (sprank ([1,0,2,0;2,0,4,0]), 2)
230 
231 %!error sprank (1,2)
232 */
octave_idx_type * xridx(void)
Definition: Sparse.h:524
bool is_real_type(void) const
Definition: ov.h:651
octave_idx_type rows(void) const
Definition: ov.h:473
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type * xcidx(void)
Definition: Sparse.h:537
void error(const char *fmt,...)
Definition: error.cc:476
F77_RET_T const double const double double * d
#define CXSPARSE_NAME(name)
Definition: dmperm.cc:43
octave_idx_type nnz(void) const
Definition: Sparse.h:248
octave_idx_type nzmax(void) const
Definition: Sparse.h:246
octave_idx_type columns(void) const
Definition: ov.h:475
int error_state
Definition: error.cc:101
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:824
double arg(double x)
Definition: lo-mappers.h:37
T & xelem(octave_idx_type n)
Definition: Array.h:353
#define DEFUN_DLD(name, args_name, nargout_name, doc)
Definition: defun-dld.h:59
static RowVector put_int(octave_idx_type *p, octave_idx_type n)
Definition: dmperm.cc:47
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:820
static octave_value_list dmperm_internal(bool rank, const octave_value arg, int nargout)
Definition: dmperm.cc:57