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
sparse-dmsolve.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2006-2013 David Bateman
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include <vector>
28 
29 #include "MArray.h"
30 #include "MSparse.h"
31 #include "SparseQR.h"
32 #include "SparseCmplxQR.h"
33 #include "MatrixType.h"
34 #include "oct-sort.h"
35 #include "oct-locbuf.h"
36 #include "oct-inttypes.h"
37 
38 template <class T>
39 static MSparse<T>
40 dmsolve_extract (const MSparse<T> &A, const octave_idx_type *Pinv,
41  const octave_idx_type *Q, octave_idx_type rst,
42  octave_idx_type rend, octave_idx_type cst,
43  octave_idx_type cend, octave_idx_type maxnz = -1,
44  bool lazy = false)
45 {
46  octave_idx_type nr = rend - rst, nc = cend - cst;
47  maxnz = (maxnz < 0 ? A.nnz () : maxnz);
48  octave_idx_type nz;
49 
50  // Cast to uint64 to handle overflow in this multiplication
51  if (octave_uint64 (nr)*octave_uint64 (nc) < octave_uint64 (maxnz))
52  nz = nr*nc;
53  else
54  nz = maxnz;
55 
56  MSparse<T> B (nr, nc, (nz < maxnz ? nz : maxnz));
57  // Some sparse functions can support lazy indexing (where elements
58  // in the row are in no particular order), even though octave in
59  // general can't. For those functions that can using it is a big
60  // win here in terms of speed.
61  if (lazy)
62  {
63  nz = 0;
64  for (octave_idx_type j = cst ; j < cend ; j++)
65  {
66  octave_idx_type qq = (Q ? Q[j] : j);
67  B.xcidx (j - cst) = nz;
68  for (octave_idx_type p = A.cidx (qq) ; p < A.cidx (qq+1) ; p++)
69  {
70  octave_quit ();
71  octave_idx_type r = (Pinv ? Pinv[A.ridx (p)] : A.ridx (p));
72  if (r >= rst && r < rend)
73  {
74  B.xdata (nz) = A.data (p);
75  B.xridx (nz++) = r - rst ;
76  }
77  }
78  }
79  B.xcidx (cend - cst) = nz ;
80  }
81  else
82  {
83  OCTAVE_LOCAL_BUFFER (T, X, rend - rst);
84  octave_sort<octave_idx_type> sort;
85  octave_idx_type *ri = B.xridx ();
86  nz = 0;
87  for (octave_idx_type j = cst ; j < cend ; j++)
88  {
89  octave_idx_type qq = (Q ? Q[j] : j);
90  B.xcidx (j - cst) = nz;
91  for (octave_idx_type p = A.cidx (qq) ; p < A.cidx (qq+1) ; p++)
92  {
93  octave_quit ();
94  octave_idx_type r = (Pinv ? Pinv[A.ridx (p)] : A.ridx (p));
95  if (r >= rst && r < rend)
96  {
97  X[r-rst] = A.data (p);
98  B.xridx (nz++) = r - rst ;
99  }
100  }
101  sort.sort (ri + B.xcidx (j - cst), nz - B.xcidx (j - cst));
102  for (octave_idx_type p = B.cidx (j - cst); p < nz; p++)
103  B.xdata (p) = X[B.xridx (p)];
104  }
105  B.xcidx (cend - cst) = nz ;
106  }
107 
108  return B;
109 }
110 
111 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
112 static MSparse<double>
113 dmsolve_extract (const MSparse<double> &A, const octave_idx_type *Pinv,
114  const octave_idx_type *Q, octave_idx_type rst,
115  octave_idx_type rend, octave_idx_type cst,
116  octave_idx_type cend, octave_idx_type maxnz,
117  bool lazy);
118 
119 static MSparse<Complex>
120 dmsolve_extract (const MSparse<Complex> &A, const octave_idx_type *Pinv,
121  const octave_idx_type *Q, octave_idx_type rst,
122  octave_idx_type rend, octave_idx_type cst,
123  octave_idx_type cend, octave_idx_type maxnz,
124  bool lazy);
125 #endif
126 
127 template <class T>
128 static MArray<T>
129 dmsolve_extract (const MArray<T> &m, const octave_idx_type *,
130  const octave_idx_type *, octave_idx_type r1,
131  octave_idx_type r2, octave_idx_type c1,
132  octave_idx_type c2)
133 {
134  r2 -= 1;
135  c2 -= 1;
136  if (r1 > r2) { std::swap (r1, r2); }
137  if (c1 > c2) { std::swap (c1, c2); }
138 
139  octave_idx_type new_r = r2 - r1 + 1;
140  octave_idx_type new_c = c2 - c1 + 1;
141 
142  MArray<T> result (dim_vector (new_r, new_c));
143 
144  for (octave_idx_type j = 0; j < new_c; j++)
145  for (octave_idx_type i = 0; i < new_r; i++)
146  result.xelem (i, j) = m.elem (r1+i, c1+j);
147 
148  return result;
149 }
150 
151 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
152 static MArray<double>
153 dmsolve_extract (const MArray<double> &m, const octave_idx_type *,
154  const octave_idx_type *, octave_idx_type r1,
155  octave_idx_type r2, octave_idx_type c1,
156  octave_idx_type c2)
157 
158 static MArray<Complex>
159 dmsolve_extract (const MArray<Complex> &m, const octave_idx_type *,
160  const octave_idx_type *, octave_idx_type r1,
161  octave_idx_type r2, octave_idx_type c1,
162  octave_idx_type c2)
163 #endif
164 
165 template <class T>
166 static void
167 dmsolve_insert (MArray<T> &a, const MArray<T> &b, const octave_idx_type *Q,
168  octave_idx_type r, octave_idx_type c)
169 {
170  T *ax = a.fortran_vec ();
171  const T *bx = b.fortran_vec ();
172  octave_idx_type anr = a.rows ();
173  octave_idx_type nr = b.rows ();
174  octave_idx_type nc = b.cols ();
175  for (octave_idx_type j = 0; j < nc; j++)
176  {
177  octave_idx_type aoff = (c + j) * anr;
178  octave_idx_type boff = j * nr;
179  for (octave_idx_type i = 0; i < nr; i++)
180  {
181  octave_quit ();
182  ax[Q[r + i] + aoff] = bx[i + boff];
183  }
184  }
185 }
186 
187 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
188 static void
189 dmsolve_insert (MArray<double> &a, const MArray<double> &b,
190  const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
191 
192 static void
193 dmsolve_insert (MArray<Complex> &a, const MArray<Complex> &b,
194  const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
195 #endif
196 
197 template <class T>
198 static void
199 dmsolve_insert (MSparse<T> &a, const MSparse<T> &b, const octave_idx_type *Q,
200  octave_idx_type r, octave_idx_type c)
201 {
202  octave_idx_type b_rows = b.rows ();
203  octave_idx_type b_cols = b.cols ();
204  octave_idx_type nr = a.rows ();
205  octave_idx_type nc = a.cols ();
206 
207  OCTAVE_LOCAL_BUFFER (octave_idx_type, Qinv, nr);
208  for (octave_idx_type i = 0; i < nr; i++)
209  Qinv[Q[i]] = i;
210 
211  // First count the number of elements in the final array
212  octave_idx_type nel = a.xcidx (c) + b.nnz ();
213 
214  if (c + b_cols < nc)
215  nel += a.xcidx (nc) - a.xcidx (c + b_cols);
216 
217  for (octave_idx_type i = c; i < c + b_cols; i++)
218  for (octave_idx_type j = a.xcidx (i); j < a.xcidx (i+1); j++)
219  if (Qinv[a.xridx (j)] < r || Qinv[a.xridx (j)] >= r + b_rows)
220  nel++;
221 
222  OCTAVE_LOCAL_BUFFER (T, X, nr);
223  octave_sort<octave_idx_type> sort;
224  MSparse<T> tmp (a);
225  a = MSparse<T> (nr, nc, nel);
226  octave_idx_type *ri = a.xridx ();
227 
228  for (octave_idx_type i = 0; i < tmp.cidx (c); i++)
229  {
230  a.xdata (i) = tmp.xdata (i);
231  a.xridx (i) = tmp.xridx (i);
232  }
233  for (octave_idx_type i = 0; i < c + 1; i++)
234  a.xcidx (i) = tmp.xcidx (i);
235 
236  octave_idx_type ii = a.xcidx (c);
237 
238  for (octave_idx_type i = c; i < c + b_cols; i++)
239  {
240  octave_quit ();
241 
242  for (octave_idx_type j = tmp.xcidx (i); j < tmp.xcidx (i+1); j++)
243  if (Qinv[tmp.xridx (j)] < r || Qinv[tmp.xridx (j)] >= r + b_rows)
244  {
245  X[tmp.xridx (j)] = tmp.xdata (j);
246  a.xridx (ii++) = tmp.xridx (j);
247  }
248 
249  octave_quit ();
250 
251  for (octave_idx_type j = b.cidx (i-c); j < b.cidx (i-c+1); j++)
252  {
253  X[Q[r + b.ridx (j)]] = b.data (j);
254  a.xridx (ii++) = Q[r + b.ridx (j)];
255  }
256 
257  sort.sort (ri + a.xcidx (i), ii - a.xcidx (i));
258  for (octave_idx_type p = a.xcidx (i); p < ii; p++)
259  a.xdata (p) = X[a.xridx (p)];
260  a.xcidx (i+1) = ii;
261  }
262 
263  for (octave_idx_type i = c + b_cols; i < nc; i++)
264  {
265  for (octave_idx_type j = tmp.xcidx (i); j < tmp.cidx (i+1); j++)
266  {
267  a.xdata (ii) = tmp.xdata (j);
268  a.xridx (ii++) = tmp.xridx (j);
269  }
270  a.xcidx (i+1) = ii;
271  }
272 }
273 
274 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
275 static void
276 dmsolve_insert (MSparse<double> &a, const SparseMatrix &b,
277  const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
278 
279 static void
280 dmsolve_insert (MSparse<Complex> &a, const MSparse<Complex> &b,
281  const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
282 #endif
283 
284 template <class T, class RT>
285 static void
286 dmsolve_permute (MArray<RT> &a, const MArray<T>& b, const octave_idx_type *p)
287 {
288  octave_idx_type b_nr = b.rows ();
289  octave_idx_type b_nc = b.cols ();
290  const T *Bx = b.fortran_vec ();
291  a.resize (dim_vector (b_nr, b_nc));
292  RT *Btx = a.fortran_vec ();
293  for (octave_idx_type j = 0; j < b_nc; j++)
294  {
295  octave_idx_type off = j * b_nr;
296  for (octave_idx_type i = 0; i < b_nr; i++)
297  {
298  octave_quit ();
299  Btx[p[i] + off] = Bx[ i + off];
300  }
301  }
302 }
303 
304 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
305 static void
306 dmsolve_permute (MArray<double> &a, const MArray<double>& b,
307  const octave_idx_type *p);
308 
309 static void
310 dmsolve_permute (MArray<Complex> &a, const MArray<double>& b,
311  const octave_idx_type *p);
312 
313 static void
314 dmsolve_permute (MArray<Complex> &a, const MArray<Complex>& b,
315  const octave_idx_type *p);
316 #endif
317 
318 template <class T, class RT>
319 static void
320 dmsolve_permute (MSparse<RT> &a, const MSparse<T>& b, const octave_idx_type *p)
321 {
322  octave_idx_type b_nr = b.rows ();
323  octave_idx_type b_nc = b.cols ();
324  octave_idx_type b_nz = b.nnz ();
325  octave_idx_type nz = 0;
326  a = MSparse<RT> (b_nr, b_nc, b_nz);
327  octave_sort<octave_idx_type> sort;
328  octave_idx_type *ri = a.xridx ();
329  OCTAVE_LOCAL_BUFFER (RT, X, b_nr);
330  a.xcidx (0) = 0;
331  for (octave_idx_type j = 0; j < b_nc; j++)
332  {
333  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
334  {
335  octave_quit ();
336  octave_idx_type r = p[b.ridx (i)];
337  X[r] = b.data (i);
338  a.xridx (nz++) = p[b.ridx (i)];
339  }
340  sort.sort (ri + a.xcidx (j), nz - a.xcidx (j));
341  for (octave_idx_type i = a.cidx (j); i < nz; i++)
342  {
343  octave_quit ();
344  a.xdata (i) = X[a.xridx (i)];
345  }
346  a.xcidx (j+1) = nz;
347  }
348 }
349 
350 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
351 static void
352 dmsolve_permute (MSparse<double> &a, const MSparse<double>& b,
353  const octave_idx_type *p);
354 
355 static void
356 dmsolve_permute (MSparse<Complex> &a, const MSparse<double>& b,
357  const octave_idx_type *p);
358 
359 static void
360 dmsolve_permute (MSparse<Complex> &a, const MSparse<Complex>& b,
361  const octave_idx_type *p);
362 #endif
363 
364 static void
365 solve_singularity_warning (double)
366 {
367  // Dummy singularity handler so that LU solver doesn't flag
368  // an error for numerically rank defficient matrices
369 }
370 
371 template <class RT, class ST, class T>
372 RT
373 dmsolve (const ST &a, const T &b, octave_idx_type &info)
374 {
375 #ifdef HAVE_CXSPARSE
376  octave_idx_type nr = a.rows ();
377  octave_idx_type nc = a.cols ();
378  octave_idx_type b_nr = b.rows ();
379  octave_idx_type b_nc = b.cols ();
380  RT retval;
381 
382  if (nr < 0 || nc < 0 || nr != b_nr)
383  (*current_liboctave_error_handler)
384  ("matrix dimension mismatch in solution of minimum norm problem");
385  else if (nr == 0 || nc == 0 || b_nc == 0)
386  retval = RT (nc, b_nc, 0.0);
387  else
388  {
389  octave_idx_type nnz_remaining = a.nnz ();
390  CXSPARSE_DNAME () csm;
391  csm.m = nr;
392  csm.n = nc;
393  csm.x = 0;
394  csm.nz = -1;
395  csm.nzmax = a.nnz ();
396  // Cast away const on A, with full knowledge that CSparse won't touch it.
397  // Prevents the methods below making a copy of the data.
398  csm.p = const_cast<octave_idx_type *>(a.cidx ());
399  csm.i = const_cast<octave_idx_type *>(a.ridx ());
400 
401 #if defined (CS_VER) && (CS_VER >= 2)
402  CXSPARSE_DNAME (d) *dm = CXSPARSE_DNAME(_dmperm) (&csm, 0);
403  octave_idx_type *p = dm->p;
404  octave_idx_type *q = dm->q;
405 #else
406  CXSPARSE_DNAME (d) *dm = CXSPARSE_DNAME(_dmperm) (&csm);
407  octave_idx_type *p = dm->P;
408  octave_idx_type *q = dm->Q;
409 #endif
410  OCTAVE_LOCAL_BUFFER (octave_idx_type, pinv, nr);
411  for (octave_idx_type i = 0; i < nr; i++)
412  pinv[p[i]] = i;
413  RT btmp;
414  dmsolve_permute (btmp, b, pinv);
415  info = 0;
416  retval.resize (nc, b_nc);
417 
418  // Leading over-determined block
419  if (dm->rr[2] < nr && dm->cc[3] < nc)
420  {
421  ST m = dmsolve_extract (a, pinv, q, dm->rr[2], nr, dm->cc[3], nc,
422  nnz_remaining, true);
423  nnz_remaining -= m.nnz ();
424  RT mtmp =
425  qrsolve (m, dmsolve_extract (btmp, 0, 0, dm->rr[2], b_nr, 0,
426  b_nc), info);
427  dmsolve_insert (retval, mtmp, q, dm->cc[3], 0);
428  if (dm->rr[2] > 0 && !info)
429  {
430  m = dmsolve_extract (a, pinv, q, 0, dm->rr[2],
431  dm->cc[3], nc, nnz_remaining, true);
432  nnz_remaining -= m.nnz ();
433  RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
434  dm->rr[2], 0, b_nc);
435  btmp.insert (ctmp - m * mtmp, 0, 0);
436  }
437  }
438 
439  // Structurally non-singular blocks
440  // FIXME: Should use fine Dulmange-Mendelsohn decomposition here.
441  if (dm->rr[1] < dm->rr[2] && dm->cc[2] < dm->cc[3] && !info)
442  {
443  ST m = dmsolve_extract (a, pinv, q, dm->rr[1], dm->rr[2],
444  dm->cc[2], dm->cc[3], nnz_remaining, false);
445  nnz_remaining -= m.nnz ();
446  RT btmp2 = dmsolve_extract (btmp, 0, 0, dm->rr[1], dm->rr[2],
447  0, b_nc);
448  double rcond = 0.0;
449  MatrixType mtyp (MatrixType::Full);
450  RT mtmp = m.solve (mtyp, btmp2, info, rcond,
451  solve_singularity_warning, false);
452  if (info != 0)
453  {
454  info = 0;
455  mtmp = qrsolve (m, btmp2, info);
456  }
457 
458  dmsolve_insert (retval, mtmp, q, dm->cc[2], 0);
459  if (dm->rr[1] > 0 && !info)
460  {
461  m = dmsolve_extract (a, pinv, q, 0, dm->rr[1], dm->cc[2],
462  dm->cc[3], nnz_remaining, true);
463  nnz_remaining -= m.nnz ();
464  RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
465  dm->rr[1], 0, b_nc);
466  btmp.insert (ctmp - m * mtmp, 0, 0);
467  }
468  }
469 
470  // Trailing under-determined block
471  if (dm->rr[1] > 0 && dm->cc[2] > 0 && !info)
472  {
473  ST m = dmsolve_extract (a, pinv, q, 0, dm->rr[1], 0,
474  dm->cc[2], nnz_remaining, true);
475  RT mtmp =
476  qrsolve (m, dmsolve_extract (btmp, 0, 0, 0, dm->rr[1] , 0,
477  b_nc), info);
478  dmsolve_insert (retval, mtmp, q, 0, 0);
479  }
480 
481  CXSPARSE_DNAME (_dfree) (dm);
482  }
483  return retval;
484 #else
485  (*current_liboctave_error_handler)
486  ("CXSPARSE unavailable; cannot solve minimum norm problem");
487  return RT ();
488 #endif
489 }
490 
491 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
492 extern Matrix
493 dmsolve (const SparseMatrix &a, const Matrix &b,
494  octave_idx_type &info);
495 
496 extern ComplexMatrix
497 dmsolve (const SparseMatrix &a, const ComplexMatrix &b,
498  octave_idx_type &info);
499 
500 extern ComplexMatrix
501 dmsolve (const SparseComplexMatrix &a, const Matrix &b,
502  octave_idx_type &info);
503 
504 extern ComplexMatrix
505 dmsolve (const SparseComplexMatrix &a, const ComplexMatrix &b,
506  octave_idx_type &info);
507 
508 extern SparseMatrix
509 dmsolve (const SparseMatrix &a, const SparseMatrix &b,
510  octave_idx_type &info);
511 
512 extern SparseComplexMatrix
513 dmsolve (const SparseMatrix &a, const SparseComplexMatrix &b,
514  octave_idx_type &info);
515 
516 extern SparseComplexMatrix
517 dmsolve (const SparseComplexMatrix &a, const SparseMatrix &b,
518  octave_idx_type &info);
519 
520 extern SparseComplexMatrix
521 dmsolve (const SparseComplexMatrix &a, const SparseComplexMatrix &b,
522  octave_idx_type &info);
523 #endif
Generated on Mon Dec 30 2013 03:04:50 for GNU Octave by   doxygen 1.8.1.2