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