GNU Octave  4.2.1
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
sparse-dmsolve.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2016-2017 John W. Eaton
4 Copyright (C) 2006-2016 David Bateman
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 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include <vector>
29 
30 #include "MArray.h"
31 #include "MSparse.h"
32 #include "MatrixType.h"
33 #include "oct-inttypes.h"
34 #include "oct-locbuf.h"
35 #include "oct-sort.h"
36 #include "oct-sparse.h"
37 #include "sparse-qr.h"
38 
39 template <typename T>
40 static MSparse<T>
42  const octave_idx_type *Q, octave_idx_type rst,
44  octave_idx_type cend, octave_idx_type maxnz = -1,
45  bool lazy = false)
46 {
47  octave_idx_type nr = rend - rst;
48  octave_idx_type nc = cend - cst;
49 
50  maxnz = (maxnz < 0 ? A.nnz () : maxnz);
51 
52  octave_idx_type nz;
53 
54  // Cast to uint64 to handle overflow in this multiplication
55  if (octave_uint64 (nr)*octave_uint64 (nc) < octave_uint64 (maxnz))
56  nz = nr*nc;
57  else
58  nz = maxnz;
59 
60  MSparse<T> B (nr, nc, (nz < maxnz ? nz : maxnz));
61 
62  // Some sparse functions can support lazy indexing (where elements
63  // in the row are in no particular order), even though octave in
64  // general can't. For those functions that can using it is a big
65  // win here in terms of speed.
66 
67  if (lazy)
68  {
69  nz = 0;
70 
71  for (octave_idx_type j = cst ; j < cend ; j++)
72  {
73  octave_idx_type qq = (Q ? Q[j] : j);
74 
75  B.xcidx (j - cst) = nz;
76 
77  for (octave_idx_type p = A.cidx (qq) ; p < A.cidx (qq+1) ; p++)
78  {
79  octave_quit ();
80 
81  octave_idx_type r = (Pinv ? Pinv[A.ridx (p)] : A.ridx (p));
82 
83  if (r >= rst && r < rend)
84  {
85  B.xdata (nz) = A.data (p);
86  B.xridx (nz++) = r - rst;
87  }
88  }
89  }
90 
91  B.xcidx (cend - cst) = nz;
92  }
93  else
94  {
95  OCTAVE_LOCAL_BUFFER (T, X, rend - rst);
96 
98  octave_idx_type *ri = B.xridx ();
99 
100  nz = 0;
101 
102  for (octave_idx_type j = cst ; j < cend ; j++)
103  {
104  octave_idx_type qq = (Q ? Q[j] : j);
105 
106  B.xcidx (j - cst) = nz;
107 
108  for (octave_idx_type p = A.cidx (qq) ; p < A.cidx (qq+1) ; p++)
109  {
110  octave_quit ();
111 
112  octave_idx_type r = (Pinv ? Pinv[A.ridx (p)] : A.ridx (p));
113 
114  if (r >= rst && r < rend)
115  {
116  X[r-rst] = A.data (p);
117  B.xridx (nz++) = r - rst;
118  }
119  }
120 
121  sort.sort (ri + B.xcidx (j - cst), nz - B.xcidx (j - cst));
122 
123  for (octave_idx_type p = B.cidx (j - cst); p < nz; p++)
124  B.xdata (p) = X[B.xridx (p)];
125  }
126 
127  B.xcidx (cend - cst) = nz;
128  }
129 
130  return B;
131 }
132 
133 template <typename T>
134 static MArray<T>
136  const octave_idx_type *, octave_idx_type r1,
138  octave_idx_type c2)
139 {
140  r2 -= 1;
141  c2 -= 1;
142 
143  if (r1 > r2)
144  std::swap (r1, r2);
145 
146  if (c1 > c2)
147  std::swap (c1, c2);
148 
149  octave_idx_type new_r = r2 - r1 + 1;
150  octave_idx_type new_c = c2 - c1 + 1;
151 
152  MArray<T> result (dim_vector (new_r, new_c));
153 
154  for (octave_idx_type j = 0; j < new_c; j++)
155  {
156  for (octave_idx_type i = 0; i < new_r; i++)
157  result.xelem (i, j) = m.elem (r1+i, c1+j);
158  }
159 
160  return result;
161 }
162 
163 template <typename T>
164 static void
167 {
168  T *ax = a.fortran_vec ();
169 
170  const T *bx = b.fortran_vec ();
171 
172  octave_idx_type anr = a.rows ();
173 
174  octave_idx_type nr = b.rows ();
175  octave_idx_type nc = b.cols ();
176 
177  for (octave_idx_type j = 0; j < nc; j++)
178  {
179  octave_idx_type aoff = (c + j) * anr;
180  octave_idx_type boff = j * nr;
181 
182  for (octave_idx_type i = 0; i < nr; i++)
183  {
184  octave_quit ();
185  ax[Q[r + i] + aoff] = bx[i + boff];
186  }
187  }
188 }
189 
190 template <typename T>
191 static void
194 {
195  octave_idx_type b_rows = b.rows ();
196  octave_idx_type b_cols = b.cols ();
197 
198  octave_idx_type nr = a.rows ();
199  octave_idx_type nc = a.cols ();
200 
202 
203  for (octave_idx_type i = 0; i < nr; i++)
204  Qinv[Q[i]] = i;
205 
206  // First count the number of elements in the final array
207  octave_idx_type nel = a.xcidx (c) + b.nnz ();
208 
209  if (c + b_cols < nc)
210  nel += a.xcidx (nc) - a.xcidx (c + b_cols);
211 
212  for (octave_idx_type i = c; i < c + b_cols; i++)
213  {
214  for (octave_idx_type j = a.xcidx (i); j < a.xcidx (i+1); j++)
215  {
216  if (Qinv[a.xridx (j)] < r || Qinv[a.xridx (j)] >= r + b_rows)
217  nel++;
218  }
219  }
220 
221  OCTAVE_LOCAL_BUFFER (T, X, nr);
222 
224 
225  MSparse<T> tmp (a);
226 
227  a = MSparse<T> (nr, nc, nel);
228 
229  octave_idx_type *ri = a.xridx ();
230 
231  for (octave_idx_type i = 0; i < tmp.cidx (c); i++)
232  {
233  a.xdata (i) = tmp.xdata (i);
234  a.xridx (i) = tmp.xridx (i);
235  }
236 
237  for (octave_idx_type i = 0; i < c + 1; i++)
238  a.xcidx (i) = tmp.xcidx (i);
239 
240  octave_idx_type ii = a.xcidx (c);
241 
242  for (octave_idx_type i = c; i < c + b_cols; i++)
243  {
244  octave_quit ();
245 
246  for (octave_idx_type j = tmp.xcidx (i); j < tmp.xcidx (i+1); j++)
247  {
248  if (Qinv[tmp.xridx (j)] < r || Qinv[tmp.xridx (j)] >= r + b_rows)
249  {
250  X[tmp.xridx (j)] = tmp.xdata (j);
251  a.xridx (ii++) = tmp.xridx (j);
252  }
253  }
254 
255  octave_quit ();
256 
257  for (octave_idx_type j = b.cidx (i-c); j < b.cidx (i-c+1); j++)
258  {
259  X[Q[r + b.ridx (j)]] = b.data (j);
260  a.xridx (ii++) = Q[r + b.ridx (j)];
261  }
262 
263  sort.sort (ri + a.xcidx (i), ii - a.xcidx (i));
264 
265  for (octave_idx_type p = a.xcidx (i); p < ii; p++)
266  a.xdata (p) = X[a.xridx (p)];
267 
268  a.xcidx (i+1) = ii;
269  }
270 
271  for (octave_idx_type i = c + b_cols; i < nc; i++)
272  {
273  for (octave_idx_type j = tmp.xcidx (i); j < tmp.cidx (i+1); j++)
274  {
275  a.xdata (ii) = tmp.xdata (j);
276  a.xridx (ii++) = tmp.xridx (j);
277  }
278 
279  a.xcidx (i+1) = ii;
280  }
281 }
282 
283 template <typename T, typename RT>
284 static void
286 {
287  octave_idx_type b_nr = b.rows ();
288  octave_idx_type b_nc = b.cols ();
289 
290  const T *Bx = b.fortran_vec ();
291 
292  a.resize (dim_vector (b_nr, b_nc));
293 
294  RT *Btx = a.fortran_vec ();
295 
296  for (octave_idx_type j = 0; j < b_nc; j++)
297  {
298  octave_idx_type off = j * b_nr;
299  for (octave_idx_type i = 0; i < b_nr; i++)
300  {
301  octave_quit ();
302  Btx[p[i] + off] = Bx[ i + off];
303  }
304  }
305 }
306 
307 template <typename T, typename RT>
308 static void
310 {
311  octave_idx_type b_nr = b.rows ();
312  octave_idx_type b_nc = b.cols ();
313  octave_idx_type b_nz = b.nnz ();
314 
315  octave_idx_type nz = 0;
316 
317  a = MSparse<RT> (b_nr, b_nc, b_nz);
319  octave_idx_type *ri = a.xridx ();
320 
321  OCTAVE_LOCAL_BUFFER (RT, X, b_nr);
322 
323  a.xcidx (0) = 0;
324 
325  for (octave_idx_type j = 0; j < b_nc; j++)
326  {
327  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
328  {
329  octave_quit ();
330  octave_idx_type r = p[b.ridx (i)];
331  X[r] = b.data (i);
332  a.xridx (nz++) = p[b.ridx (i)];
333  }
334 
335  sort.sort (ri + a.xcidx (j), nz - a.xcidx (j));
336 
337  for (octave_idx_type i = a.cidx (j); i < nz; i++)
338  {
339  octave_quit ();
340  a.xdata (i) = X[a.xridx (i)];
341  }
342 
343  a.xcidx (j+1) = nz;
344  }
345 }
346 
347 #if defined (HAVE_CXSPARSE)
348 
349 static void
351 {
352  // Dummy singularity handler so that LU solver doesn't flag
353  // an error for numerically rank defficient matrices
354 }
355 
356 #endif
357 
358 template <typename RT, typename ST, typename T>
359 RT
360 dmsolve (const ST &a, const T &b, octave_idx_type &info)
361 {
362  RT retval;
363 
364 #if defined (HAVE_CXSPARSE)
365 
366  octave_idx_type nr = a.rows ();
367  octave_idx_type nc = a.cols ();
368 
369  octave_idx_type b_nr = b.rows ();
370  octave_idx_type b_nc = b.cols ();
371 
372  if (nr < 0 || nc < 0 || nr != b_nr)
373  (*current_liboctave_error_handler)
374  ("matrix dimension mismatch in solution of minimum norm problem");
375 
376  if (nr == 0 || nc == 0 || b_nc == 0)
377  retval = RT (nc, b_nc, 0.0);
378  else
379  {
380  octave_idx_type nnz_remaining = a.nnz ();
381 
382  CXSPARSE_DNAME () csm;
383 
384  csm.m = nr;
385  csm.n = nc;
386  csm.x = 0;
387  csm.nz = -1;
388  csm.nzmax = a.nnz ();
389 
390  // Cast away const on A, with full knowledge that CSparse won't touch it.
391  // Prevents the methods below making a copy of the data.
392  csm.p = const_cast<octave_idx_type *>(a.cidx ());
393  csm.i = const_cast<octave_idx_type *>(a.ridx ());
394 
395  CXSPARSE_DNAME (d) *dm = CXSPARSE_DNAME(_dmperm) (&csm, 0);
396  octave_idx_type *p = dm->p;
397  octave_idx_type *q = dm->q;
398 
400 
401  for (octave_idx_type i = 0; i < nr; i++)
402  pinv[p[i]] = i;
403 
404  RT btmp;
405  dmsolve_permute (btmp, b, pinv);
406  info = 0;
407 
408  retval.resize (nc, b_nc);
409 
410  // Leading over-determined block
411  if (dm->rr[2] < nr && dm->cc[3] < nc)
412  {
413  ST m = dmsolve_extract (a, pinv, q, dm->rr[2], nr, dm->cc[3], nc,
414  nnz_remaining, true);
415  nnz_remaining -= m.nnz ();
416  RT mtmp = octave::math::qrsolve (m, dmsolve_extract (btmp, 0, 0, dm->rr[2],
417  b_nr, 0, b_nc), info);
418  dmsolve_insert (retval, mtmp, q, dm->cc[3], 0);
419 
420  if (dm->rr[2] > 0 && ! info)
421  {
422  m = dmsolve_extract (a, pinv, q, 0, dm->rr[2],
423  dm->cc[3], nc, nnz_remaining, true);
424  nnz_remaining -= m.nnz ();
425  RT ctmp = dmsolve_extract (btmp, 0, 0, 0, dm->rr[2], 0, b_nc);
426  btmp.insert (ctmp - m * mtmp, 0, 0);
427  }
428  }
429 
430  // Structurally non-singular blocks
431  // FIXME: Should use fine Dulmange-Mendelsohn decomposition here.
432  if (dm->rr[1] < dm->rr[2] && dm->cc[2] < dm->cc[3] && ! info)
433  {
434  ST m = dmsolve_extract (a, pinv, q, dm->rr[1], dm->rr[2],
435  dm->cc[2], dm->cc[3], nnz_remaining, false);
436  nnz_remaining -= m.nnz ();
437  RT btmp2 = dmsolve_extract (btmp, 0, 0, dm->rr[1], dm->rr[2],
438  0, b_nc);
439  double rcond = 0.0;
441  RT mtmp = m.solve (mtyp, btmp2, info, rcond,
443  if (info != 0)
444  {
445  info = 0;
446  mtmp = octave::math::qrsolve (m, btmp2, info);
447  }
448 
449  dmsolve_insert (retval, mtmp, q, dm->cc[2], 0);
450  if (dm->rr[1] > 0 && ! info)
451  {
452  m = dmsolve_extract (a, pinv, q, 0, dm->rr[1], dm->cc[2],
453  dm->cc[3], nnz_remaining, true);
454  nnz_remaining -= m.nnz ();
455  RT ctmp = dmsolve_extract (btmp, 0, 0, 0, dm->rr[1], 0, b_nc);
456  btmp.insert (ctmp - m * mtmp, 0, 0);
457  }
458  }
459 
460  // Trailing under-determined block
461  if (dm->rr[1] > 0 && dm->cc[2] > 0 && ! info)
462  {
463  ST m = dmsolve_extract (a, pinv, q, 0, dm->rr[1], 0,
464  dm->cc[2], nnz_remaining, true);
465  RT mtmp = octave::math::qrsolve (m, dmsolve_extract (btmp, 0, 0, 0, dm->rr[1],
466  0, b_nc), info);
467  dmsolve_insert (retval, mtmp, q, 0, 0);
468  }
469 
470  CXSPARSE_DNAME (_dfree) (dm);
471  }
472 
473 #else
474 
475  octave_unused_parameter (a);
476  octave_unused_parameter (b);
477  octave_unused_parameter (info);
478 
479  (*current_liboctave_error_handler)
480  ("support for CXSparse was unavailable or disabled when liboctave was built");
481 
482 #endif
483 
484  return retval;
485 }
486 
487 // Instantiations we need.
488 
489 template ComplexMatrix
491  (const SparseComplexMatrix&, const Matrix&, octave_idx_type&);
492 
493 template SparseComplexMatrix
496 
497 template ComplexMatrix
500 
501 template SparseComplexMatrix
504 
505 template Matrix
507  (const SparseMatrix&, const Matrix&, octave_idx_type&);
508 
509 template SparseMatrix
511  (const SparseMatrix&, const SparseMatrix&, octave_idx_type&);
512 
513 template ComplexMatrix
515  (const SparseMatrix&, const ComplexMatrix&, octave_idx_type&);
516 
517 template SparseComplexMatrix
template Matrix dmsolve< Matrix, SparseMatrix, Matrix >(const SparseMatrix &, const Matrix &, octave_idx_type &)
octave_idx_type * xridx(void)
Definition: Sparse.h:536
octave_int< uint64_t > octave_uint64
octave_idx_type cols(void) const
Definition: Sparse.h:272
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
Definition: sparse-qr.cc:2288
T * data(void)
Definition: Sparse.h:521
octave_idx_type rows(void) const
Definition: Sparse.h:271
static void solve_singularity_warning(double)
template SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseComplexMatrix, SparseMatrix >(const SparseComplexMatrix &, const SparseMatrix &, octave_idx_type &)
template SparseMatrix dmsolve< SparseMatrix, SparseMatrix, SparseMatrix >(const SparseMatrix &, const SparseMatrix &, octave_idx_type &)
static void dmsolve_insert(MArray< T > &a, const MArray< T > &b, const octave_idx_type *Q, octave_idx_type r, octave_idx_type c)
octave_idx_type * xcidx(void)
Definition: Sparse.h:549
octave_idx_type b_nr
Definition: sylvester.cc:74
template SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseMatrix, SparseComplexMatrix >(const SparseMatrix &, const SparseComplexMatrix &, octave_idx_type &)
octave_idx_type * cidx(void)
Definition: Sparse.h:543
T & elem(octave_idx_type n)
Definition: Array.h:482
Template for N-dimensional array classes with like-type math operators.
Definition: MArray.h:32
octave_idx_type rows(void) const
Definition: Array.h:401
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
octave_idx_type nnz(void) const
Actual number of nonzero terms.
Definition: Sparse.h:253
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
bool swap
Definition: load-save.cc:725
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
template ComplexMatrix dmsolve< ComplexMatrix, SparseComplexMatrix, ComplexMatrix >(const SparseComplexMatrix &, const ComplexMatrix &, octave_idx_type &)
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
Sparse< T > sort(octave_idx_type dim=0, sortmode mode=ASCENDING) const
Definition: Sparse.cc:2233
void resize(const dim_vector &dv, const T &rfv)
Definition: Array.cc:1028
double tmp
Definition: data.cc:6300
octave_value retval
Definition: data.cc:6294
Definition: dMatrix.h:37
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
Definition: lu.cc:138
With real return the complex result
Definition: data.cc:3375
static void dmsolve_permute(MArray< RT > &a, const MArray< T > &b, const octave_idx_type *p)
T & xelem(octave_idx_type n)
Definition: Array.h:455
#define CXSPARSE_DNAME(name)
Definition: oct-sparse.h:117
void sort(T *data, octave_idx_type nel)
Definition: oct-sort.cc:1506
octave_idx_type * ridx(void)
Definition: Sparse.h:530
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
octave_idx_type b_nc
Definition: sylvester.cc:75
p
Definition: lu.cc:138
T * xdata(void)
Definition: Sparse.h:523
b
Definition: cellfun.cc:398
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:200
template ComplexMatrix dmsolve< ComplexMatrix, SparseMatrix, ComplexMatrix >(const SparseMatrix &, const ComplexMatrix &, octave_idx_type &)
static MSparse< T > dmsolve_extract(const MSparse< T > &A, const octave_idx_type *Pinv, const octave_idx_type *Q, octave_idx_type rst, octave_idx_type rend, octave_idx_type cst, octave_idx_type cend, octave_idx_type maxnz=-1, bool lazy=false)
template SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseComplexMatrix, SparseComplexMatrix >(const SparseComplexMatrix &, const SparseComplexMatrix &, octave_idx_type &)
const T * fortran_vec(void) const
Definition: Array.h:584
octave_idx_type cols(void) const
Definition: Array.h:409
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
template ComplexMatrix dmsolve< ComplexMatrix, SparseComplexMatrix, Matrix >(const SparseComplexMatrix &, const Matrix &, octave_idx_type &)
RT dmsolve(const ST &a, const T &b, octave_idx_type &info)
F77_RET_T const F77_INT F77_CMPLX * A
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Q