GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
sparse-qr.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2016-2018 John W. Eaton
4 Copyright (C) 2005-2018 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
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License 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 <https://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include "CMatrix.h"
29 #include "CSparse.h"
30 #include "MArray.h"
31 #include "dColVector.h"
32 #include "dMatrix.h"
33 #include "dSparse.h"
34 #include "lo-error.h"
35 #include "oct-locbuf.h"
36 #include "oct-refcount.h"
37 #include "oct-sparse.h"
38 #include "quit.h"
39 #include "sparse-qr.h"
40 
41 namespace octave
42 {
43  namespace math
44  {
45  template <typename SPARSE_T>
46  class
47  cxsparse_types
48  {
49  };
50 
51  template <>
52  class
53  cxsparse_types<SparseMatrix>
54  {
55  public:
56 #if defined (HAVE_CXSPARSE)
57  typedef CXSPARSE_DNAME (s) symbolic_type;
58  typedef CXSPARSE_DNAME (n) numeric_type;
59 #else
60  typedef void symbolic_type;
61  typedef void numeric_type;
62 #endif
63  };
64 
65  template <>
66  class
67  cxsparse_types<SparseComplexMatrix>
68  {
69  public:
70 #if defined (HAVE_CXSPARSE)
71  typedef CXSPARSE_ZNAME (s) symbolic_type;
72  typedef CXSPARSE_ZNAME (n) numeric_type;
73 #else
74  typedef void symbolic_type;
75  typedef void numeric_type;
76 #endif
77  };
78 
79  template <typename SPARSE_T>
80  class sparse_qr<SPARSE_T>::sparse_qr_rep
81  {
82  public:
83 
84  sparse_qr_rep (const SPARSE_T& a, int order);
85 
86  // No copying!
87 
88  sparse_qr_rep (const sparse_qr_rep&) = delete;
89 
90  sparse_qr_rep& operator = (const sparse_qr_rep&) = delete;
91 
92  ~sparse_qr_rep (void);
93 
94  bool ok (void) const
95  {
96 #if defined (HAVE_CXSPARSE)
97  return (N && S);
98 #else
99  return false;
100 #endif
101  }
102 
103  SPARSE_T V (void) const;
104 
105  ColumnVector Pinv (void) const;
106 
107  ColumnVector P (void) const;
108 
109  SPARSE_T R (bool econ) const;
110 
111  typename SPARSE_T::dense_matrix_type
112  C (const typename SPARSE_T::dense_matrix_type& b) const;
113 
114  typename SPARSE_T::dense_matrix_type
115  Q (void) const;
116 
117  refcount<int> count;
118 
119  octave_idx_type nrows;
120  octave_idx_type ncols;
121 
122  typename cxsparse_types<SPARSE_T>::symbolic_type *S;
123  typename cxsparse_types<SPARSE_T>::numeric_type *N;
124 
125  template <typename RHS_T, typename RET_T>
126  RET_T
127  tall_solve (const RHS_T& b, octave_idx_type& info) const;
128 
129  template <typename RHS_T, typename RET_T>
130  RET_T
131  wide_solve (const RHS_T& b, octave_idx_type& info) const;
132  };
133 
134  template <typename SPARSE_T>
135  ColumnVector
136  sparse_qr<SPARSE_T>::sparse_qr_rep::Pinv (void) const
137  {
138 #if defined (HAVE_CXSPARSE)
139 
140  ColumnVector ret (N->L->m);
141 
142  for (octave_idx_type i = 0; i < N->L->m; i++)
143  ret.xelem (i) = S->pinv[i];
144 
145  return ret;
146 
147 #else
148 
149  return ColumnVector ();
150 
151 #endif
152  }
153 
154  template <typename SPARSE_T>
155  ColumnVector
156  sparse_qr<SPARSE_T>::sparse_qr_rep::P (void) const
157  {
158 #if defined (HAVE_CXSPARSE)
159 
160  ColumnVector ret (N->L->m);
161 
162  for (octave_idx_type i = 0; i < N->L->m; i++)
163  ret.xelem (S->pinv[i]) = i;
164 
165  return ret;
166 
167 #else
168 
169  return ColumnVector ();
170 
171 #endif
172  }
173 
174  // Specializations.
175 
176  // Real-valued matrices.
177 
178  template <>
179  sparse_qr<SparseMatrix>::sparse_qr_rep::sparse_qr_rep
180  (const SparseMatrix& a, int order)
181  : count (1), nrows (a.rows ()), ncols (a.columns ())
182 #if defined (HAVE_CXSPARSE)
183  , S (nullptr), N (nullptr)
184 #endif
185  {
186 #if defined (HAVE_CXSPARSE)
187 
188  CXSPARSE_DNAME () A;
189 
190  A.nzmax = a.nnz ();
191  A.m = nrows;
192  A.n = ncols;
193  // Cast away const on A, with full knowledge that CSparse won't touch it
194  // Prevents the methods below making a copy of the data.
195  A.p = const_cast<suitesparse_integer *>
196  (to_suitesparse_intptr (a.cidx ()));
197  A.i = const_cast<suitesparse_integer *>
198  (to_suitesparse_intptr (a.ridx ()));
199  A.x = const_cast<double *> (a.data ());
200  A.nz = -1;
201 
202  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
203  S = CXSPARSE_DNAME (_sqr) (order, &A, 1);
204  N = CXSPARSE_DNAME (_qr) (&A, S);
205  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
206 
207  if (! N)
208  (*current_liboctave_error_handler)
209  ("sparse_qr: sparse matrix QR factorization filled");
210 
211  count = 1;
212 
213 #else
214 
215  octave_unused_parameter (order);
216 
217  (*current_liboctave_error_handler)
218  ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built");
219 
220 #endif
221  }
222 
223  template <>
224  sparse_qr<SparseMatrix>::sparse_qr_rep::~sparse_qr_rep (void)
225  {
226 #if defined (HAVE_CXSPARSE)
227  CXSPARSE_DNAME (_sfree) (S);
228  CXSPARSE_DNAME (_nfree) (N);
229 #endif
230  }
231 
232  template <>
233  SparseMatrix
234  sparse_qr<SparseMatrix>::sparse_qr_rep::V (void) const
235  {
236 #if defined (HAVE_CXSPARSE)
237 
238  // Drop zeros from V and sort
239  // FIXME: Is the double transpose to sort necessary?
240 
241  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
242  CXSPARSE_DNAME (_dropzeros) (N->L);
243  CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->L, 1);
244  CXSPARSE_DNAME (_spfree) (N->L);
245  N->L = CXSPARSE_DNAME (_transpose) (D, 1);
246  CXSPARSE_DNAME (_spfree) (D);
247  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
248 
249  octave_idx_type nc = N->L->n;
250  octave_idx_type nz = N->L->nzmax;
251  SparseMatrix ret (N->L->m, nc, nz);
252 
253  for (octave_idx_type j = 0; j < nc+1; j++)
254  ret.xcidx (j) = N->L->p[j];
255 
256  for (octave_idx_type j = 0; j < nz; j++)
257  {
258  ret.xridx (j) = N->L->i[j];
259  ret.xdata (j) = N->L->x[j];
260  }
261 
262  return ret;
263 
264 #else
265 
266  return SparseMatrix ();
267 
268 #endif
269  }
270 
271  template <>
272  SparseMatrix
273  sparse_qr<SparseMatrix>::sparse_qr_rep::R (bool econ) const
274  {
275 #if defined (HAVE_CXSPARSE)
276 
277  // Drop zeros from R and sort
278  // FIXME: Is the double transpose to sort necessary?
279 
280  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
281  CXSPARSE_DNAME (_dropzeros) (N->U);
282  CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->U, 1);
283  CXSPARSE_DNAME (_spfree) (N->U);
284  N->U = CXSPARSE_DNAME (_transpose) (D, 1);
285  CXSPARSE_DNAME (_spfree) (D);
286  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
287 
288  octave_idx_type nc = N->U->n;
289  octave_idx_type nz = N->U->nzmax;
290 
291  SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
292 
293  for (octave_idx_type j = 0; j < nc+1; j++)
294  ret.xcidx (j) = N->U->p[j];
295 
296  for (octave_idx_type j = 0; j < nz; j++)
297  {
298  ret.xridx (j) = N->U->i[j];
299  ret.xdata (j) = N->U->x[j];
300  }
301 
302  return ret;
303 
304 #else
305 
306  octave_unused_parameter (econ);
307 
308  return SparseMatrix ();
309 
310 #endif
311  }
312 
313  template <>
314  Matrix
315  sparse_qr<SparseMatrix>::sparse_qr_rep::C (const Matrix& b) const
316  {
317 #if defined (HAVE_CXSPARSE)
318 
319  octave_idx_type b_nr = b.rows ();
320  octave_idx_type b_nc = b.cols ();
321 
322  octave_idx_type nc = N->L->n;
323  octave_idx_type nr = nrows;
324 
325  const double *bvec = b.fortran_vec ();
326 
327  Matrix ret (b_nr, b_nc);
328  double *vec = ret.fortran_vec ();
329 
330  if (nr < 0 || nc < 0 || nr != b_nr)
331  (*current_liboctave_error_handler) ("matrix dimension mismatch");
332 
333  if (nr == 0 || nc == 0 || b_nc == 0)
334  ret = Matrix (nc, b_nc, 0.0);
335  else
336  {
337  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
338 
339  for (volatile octave_idx_type j = 0, idx = 0;
340  j < b_nc;
341  j++, idx += b_nr)
342  {
343  octave_quit ();
344 
345  for (octave_idx_type i = nr; i < S->m2; i++)
346  buf[i] = 0.;
347 
348  volatile octave_idx_type nm = (nr < nc ? nr : nc);
349 
350  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
351  CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr);
352  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
353 
354  for (volatile octave_idx_type i = 0; i < nm; i++)
355  {
356  octave_quit ();
357 
358  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
359  CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
360  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
361  }
362 
363  for (octave_idx_type i = 0; i < b_nr; i++)
364  vec[i+idx] = buf[i];
365  }
366  }
367 
368  return ret;
369 
370 #else
371 
372  octave_unused_parameter (b);
373 
374  return Matrix ();
375 
376 #endif
377  }
378 
379  template <>
380  Matrix
381  sparse_qr<SparseMatrix>::sparse_qr_rep::Q (void) const
382  {
383 #if defined (HAVE_CXSPARSE)
384  octave_idx_type nc = N->L->n;
385  octave_idx_type nr = nrows;
386  Matrix ret (nr, nr);
387  double *vec = ret.fortran_vec ();
388 
389  if (nr < 0 || nc < 0)
390  (*current_liboctave_error_handler) ("matrix dimension mismatch");
391 
392  if (nr == 0 || nc == 0)
393  ret = Matrix (nc, nr, 0.0);
394  else
395  {
396  OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1);
397 
398  for (octave_idx_type i = 0; i < nr; i++)
399  bvec[i] = 0.;
400 
401  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
402 
403  for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
404  {
405  octave_quit ();
406 
407  bvec[j] = 1.0;
408  for (octave_idx_type i = nr; i < S->m2; i++)
409  buf[i] = 0.;
410 
411  volatile octave_idx_type nm = (nr < nc ? nr : nc);
412 
413  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
414  CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr);
415  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
416 
417  for (volatile octave_idx_type i = 0; i < nm; i++)
418  {
419  octave_quit ();
420 
421  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
422  CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
423  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
424  }
425 
426  for (octave_idx_type i = 0; i < nr; i++)
427  vec[i+idx] = buf[i];
428 
429  bvec[j] = 0.0;
430  }
431  }
432 
433  return ret.transpose ();
434 
435 #else
436 
437  return Matrix ();
438 
439 #endif
440  }
441 
442  template <>
443  template <>
444  Matrix
445  sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<double>, Matrix>
446  (const MArray<double>& b, octave_idx_type& info) const
447  {
448  info = -1;
449 
450 #if defined (HAVE_CXSPARSE)
451 
452  octave_idx_type nr = nrows;
453  octave_idx_type nc = ncols;
454 
455  octave_idx_type b_nc = b.cols ();
456  octave_idx_type b_nr = b.rows ();
457 
458  const double *bvec = b.data ();
459 
460  Matrix x (nc, b_nc);
461  double *vec = x.fortran_vec ();
462 
463  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
464 
465  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
466  i++, idx+=nc, bidx+=b_nr)
467  {
468  octave_quit ();
469 
470  for (octave_idx_type j = nr; j < S->m2; j++)
471  buf[j] = 0.;
472 
473  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
474  CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + bidx, buf, nr);
475  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
476 
477  for (volatile octave_idx_type j = 0; j < nc; j++)
478  {
479  octave_quit ();
480 
481  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
482  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
483  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
484  }
485 
486  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
487  CXSPARSE_DNAME (_usolve) (N->U, buf);
488  CXSPARSE_DNAME (_ipvec) (S->q, buf, vec + idx, nc);
489  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
490  }
491 
492  info = 0;
493 
494  return x;
495 
496 #else
497 
498  octave_unused_parameter (b);
499 
500  return Matrix ();
501 
502 #endif
503  }
504 
505  template <>
506  template <>
507  Matrix
508  sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<double>, Matrix>
509  (const MArray<double>& b, octave_idx_type& info) const
510  {
511  info = -1;
512 
513 #if defined (HAVE_CXSPARSE)
514 
515  // These are swapped because the original matrix was transposed in
516  // sparse_qr<SparseMatrix>::solve.
517 
518  octave_idx_type nr = ncols;
519  octave_idx_type nc = nrows;
520 
521  octave_idx_type b_nc = b.cols ();
522  octave_idx_type b_nr = b.rows ();
523 
524  const double *bvec = b.data ();
525 
526  Matrix x (nc, b_nc);
527  double *vec = x.fortran_vec ();
528 
529  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
530 
531  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
532 
533  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
534  i++, idx+=nc, bidx+=b_nr)
535  {
536  octave_quit ();
537 
538  for (octave_idx_type j = nr; j < nbuf; j++)
539  buf[j] = 0.;
540 
541  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
542  CXSPARSE_DNAME (_pvec) (S->q, bvec + bidx, buf, nr);
543  CXSPARSE_DNAME (_utsolve) (N->U, buf);
544  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
545 
546  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
547  {
548  octave_quit ();
549 
550  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
551  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
552  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
553  }
554 
555  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
556  CXSPARSE_DNAME (_pvec) (S->pinv, buf, vec + idx, nc);
557  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
558  }
559 
560  info = 0;
561 
562  return x;
563 
564 #else
565 
566  octave_unused_parameter (b);
567 
568  return Matrix ();
569 
570 #endif
571  }
572 
573  template <>
574  template <>
575  SparseMatrix
576  sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseMatrix, SparseMatrix>
577  (const SparseMatrix& b, octave_idx_type& info) const
578  {
579  info = -1;
580 
581 #if defined (HAVE_CXSPARSE)
582 
583  octave_idx_type nr = nrows;
584  octave_idx_type nc = ncols;
585 
586  octave_idx_type b_nr = b.rows ();
587  octave_idx_type b_nc = b.cols ();
588 
589  SparseMatrix x (nc, b_nc, b.nnz ());
590  x.xcidx (0) = 0;
591 
592  volatile octave_idx_type x_nz = b.nnz ();
593  volatile octave_idx_type ii = 0;
594 
595  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
596  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
597 
598  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
599  {
600  octave_quit ();
601 
602  for (octave_idx_type j = 0; j < b_nr; j++)
603  Xx[j] = b.xelem (j,i);
604 
605  for (octave_idx_type j = nr; j < S->m2; j++)
606  buf[j] = 0.;
607 
608  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
609  CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr);
610  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
611 
612  for (volatile octave_idx_type j = 0; j < nc; j++)
613  {
614  octave_quit ();
615 
616  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
617  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
618  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
619  }
620 
621  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
622  CXSPARSE_DNAME (_usolve) (N->U, buf);
623  CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc);
624  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
625 
626  for (octave_idx_type j = 0; j < nc; j++)
627  {
628  double tmp = Xx[j];
629 
630  if (tmp != 0.0)
631  {
632  if (ii == x_nz)
633  {
634  // Resize the sparse matrix
635  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
636  sz = (sz > 10 ? sz : 10) + x_nz;
637  x.change_capacity (sz);
638  x_nz = sz;
639  }
640 
641  x.xdata (ii) = tmp;
642  x.xridx (ii++) = j;
643  }
644  }
645 
646  x.xcidx (i+1) = ii;
647  }
648 
649  info = 0;
650 
651  return x;
652 
653 #else
654 
655  octave_unused_parameter (b);
656 
657  return SparseMatrix ();
658 
659 #endif
660  }
661 
662  template <>
663  template <>
664  SparseMatrix
665  sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<SparseMatrix, SparseMatrix>
666  (const SparseMatrix& b, octave_idx_type& info) const
667  {
668  info = -1;
669 
670 #if defined (HAVE_CXSPARSE)
671 
672  // These are swapped because the original matrix was transposed in
673  // sparse_qr<SparseMatrix>::solve.
674 
675  octave_idx_type nr = ncols;
676  octave_idx_type nc = nrows;
677 
678  octave_idx_type b_nr = b.rows ();
679  octave_idx_type b_nc = b.cols ();
680 
681  SparseMatrix x (nc, b_nc, b.nnz ());
682  x.xcidx (0) = 0;
683 
684  volatile octave_idx_type x_nz = b.nnz ();
685  volatile octave_idx_type ii = 0;
686  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
687 
688  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
689  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
690 
691  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
692  {
693  octave_quit ();
694 
695  for (octave_idx_type j = 0; j < b_nr; j++)
696  Xx[j] = b.xelem (j,i);
697 
698  for (octave_idx_type j = nr; j < nbuf; j++)
699  buf[j] = 0.;
700 
701  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
702  CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr);
703  CXSPARSE_DNAME (_utsolve) (N->U, buf);
704  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
705 
706  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
707  {
708  octave_quit ();
709 
710  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
711  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
712  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
713  }
714 
715  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
716  CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc);
717  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
718 
719  for (octave_idx_type j = 0; j < nc; j++)
720  {
721  double tmp = Xx[j];
722 
723  if (tmp != 0.0)
724  {
725  if (ii == x_nz)
726  {
727  // Resize the sparse matrix
728  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
729  sz = (sz > 10 ? sz : 10) + x_nz;
730  x.change_capacity (sz);
731  x_nz = sz;
732  }
733 
734  x.xdata (ii) = tmp;
735  x.xridx (ii++) = j;
736  }
737  }
738 
739  x.xcidx (i+1) = ii;
740  }
741 
742  info = 0;
743 
744  x.maybe_compress ();
745 
746  return x;
747 
748 #else
749 
750  octave_unused_parameter (b);
751 
752  return SparseMatrix ();
753 
754 #endif
755  }
756 
757  template <>
758  template <>
759  ComplexMatrix
760  sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>, ComplexMatrix>
761  (const MArray<Complex>& b, octave_idx_type& info) const
762  {
763  info = -1;
764 
765 #if defined (HAVE_CXSPARSE)
766 
767  octave_idx_type nr = nrows;
768  octave_idx_type nc = ncols;
769 
770  octave_idx_type b_nc = b.cols ();
771  octave_idx_type b_nr = b.rows ();
772 
773  ComplexMatrix x (nc, b_nc);
774  Complex *vec = x.fortran_vec ();
775 
776  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
777  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
778  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
779 
780  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
781  {
782  octave_quit ();
783 
784  for (octave_idx_type j = 0; j < b_nr; j++)
785  {
786  Complex c = b.xelem (j,i);
787  Xx[j] = c.real ();
788  Xz[j] = c.imag ();
789  }
790 
791  for (octave_idx_type j = nr; j < S->m2; j++)
792  buf[j] = 0.;
793 
794  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
795  CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr);
796  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
797 
798  for (volatile octave_idx_type j = 0; j < nc; j++)
799  {
800  octave_quit ();
801 
802  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
803  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
804  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
805  }
806 
807  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
808  CXSPARSE_DNAME (_usolve) (N->U, buf);
809  CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc);
810 
811  for (octave_idx_type j = nr; j < S->m2; j++)
812  buf[j] = 0.;
813 
814  CXSPARSE_DNAME (_ipvec) (S->pinv, Xz, buf, nr);
815  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
816 
817  for (volatile octave_idx_type j = 0; j < nc; j++)
818  {
819  octave_quit ();
820 
821  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
822  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
823  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
824  }
825 
826  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
827  CXSPARSE_DNAME (_usolve) (N->U, buf);
828  CXSPARSE_DNAME (_ipvec) (S->q, buf, Xz, nc);
829  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
830 
831  for (octave_idx_type j = 0; j < nc; j++)
832  vec[j+idx] = Complex (Xx[j], Xz[j]);
833  }
834 
835  info = 0;
836 
837  return x;
838 
839 #else
840 
841  octave_unused_parameter (b);
842 
843  return ComplexMatrix ();
844 
845 #endif
846  }
847 
848  template <>
849  template <>
850  ComplexMatrix
851  sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>, ComplexMatrix>
852  (const MArray<Complex>& b, octave_idx_type& info) const
853  {
854  info = -1;
855 
856 #if defined (HAVE_CXSPARSE)
857 
858  // These are swapped because the original matrix was transposed in
859  // sparse_qr<SparseMatrix>::solve.
860 
861  octave_idx_type nr = ncols;
862  octave_idx_type nc = nrows;
863 
864  octave_idx_type b_nc = b.cols ();
865  octave_idx_type b_nr = b.rows ();
866 
867  ComplexMatrix x (nc, b_nc);
868  Complex *vec = x.fortran_vec ();
869 
870  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
871 
872  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
873  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
874  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
875 
876  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
877  {
878  octave_quit ();
879 
880  for (octave_idx_type j = 0; j < b_nr; j++)
881  {
882  Complex c = b.xelem (j,i);
883  Xx[j] = c.real ();
884  Xz[j] = c.imag ();
885  }
886 
887  for (octave_idx_type j = nr; j < nbuf; j++)
888  buf[j] = 0.;
889 
890  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
891  CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr);
892  CXSPARSE_DNAME (_utsolve) (N->U, buf);
893  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
894 
895  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
896  {
897  octave_quit ();
898 
899  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
900  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
901  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
902  }
903 
904  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
905  CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc);
906  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
907 
908  for (octave_idx_type j = nr; j < nbuf; j++)
909  buf[j] = 0.;
910 
911  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
912  CXSPARSE_DNAME (_pvec) (S->q, Xz, buf, nr);
913  CXSPARSE_DNAME (_utsolve) (N->U, buf);
914  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
915 
916  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
917  {
918  octave_quit ();
919 
920  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
921  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
922  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
923  }
924 
925  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
926  CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xz, nc);
927  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
928 
929  for (octave_idx_type j = 0; j < nc; j++)
930  vec[j+idx] = Complex (Xx[j], Xz[j]);
931  }
932 
933  info = 0;
934 
935  return x;
936 
937 #else
938 
939  octave_unused_parameter (b);
940 
941  return ComplexMatrix ();
942 
943 #endif
944  }
945 
946  // Complex-valued matrices.
947 
948  template <>
949  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::sparse_qr_rep
950  (const SparseComplexMatrix& a, int order)
951  : count (1), nrows (a.rows ()), ncols (a.columns ())
952 #if defined (HAVE_CXSPARSE)
953  , S (nullptr), N (nullptr)
954 #endif
955  {
956 #if defined (HAVE_CXSPARSE)
957 
958  CXSPARSE_ZNAME () A;
959 
960  A.nzmax = a.nnz ();
961  A.m = nrows;
962  A.n = ncols;
963  // Cast away const on A, with full knowledge that CSparse won't touch it
964  // Prevents the methods below making a copy of the data.
965  A.p = const_cast<suitesparse_integer *>
966  (to_suitesparse_intptr (a.cidx ()));
967  A.i = const_cast<suitesparse_integer *>
968  (to_suitesparse_intptr (a.ridx ()));
969  A.x = const_cast<cs_complex_t *>
970  (reinterpret_cast<const cs_complex_t *> (a.data ()));
971  A.nz = -1;
972 
973  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
974  S = CXSPARSE_ZNAME (_sqr) (order, &A, 1);
975  N = CXSPARSE_ZNAME (_qr) (&A, S);
976  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
977 
978  if (! N)
979  (*current_liboctave_error_handler)
980  ("sparse_qr: sparse matrix QR factorization filled");
981 
982  count = 1;
983 
984 #else
985 
986  octave_unused_parameter (order);
987 
988  (*current_liboctave_error_handler)
989  ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built");
990 
991 #endif
992  }
993 
994  template <>
995  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::~sparse_qr_rep (void)
996  {
997 #if defined (HAVE_CXSPARSE)
998  CXSPARSE_ZNAME (_sfree) (S);
999  CXSPARSE_ZNAME (_nfree) (N);
1000 #endif
1001  }
1002 
1003  template <>
1004  SparseComplexMatrix
1005  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::V (void) const
1006  {
1007 #if defined (HAVE_CXSPARSE)
1008  // Drop zeros from V and sort
1009  // FIXME: Is the double transpose to sort necessary?
1010 
1011  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1012  CXSPARSE_ZNAME (_dropzeros) (N->L);
1013  CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->L, 1);
1014  CXSPARSE_ZNAME (_spfree) (N->L);
1015  N->L = CXSPARSE_ZNAME (_transpose) (D, 1);
1016  CXSPARSE_ZNAME (_spfree) (D);
1017  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1018 
1019  octave_idx_type nc = N->L->n;
1020  octave_idx_type nz = N->L->nzmax;
1021  SparseComplexMatrix ret (N->L->m, nc, nz);
1022 
1023  for (octave_idx_type j = 0; j < nc+1; j++)
1024  ret.xcidx (j) = N->L->p[j];
1025 
1026  for (octave_idx_type j = 0; j < nz; j++)
1027  {
1028  ret.xridx (j) = N->L->i[j];
1029  ret.xdata (j) = reinterpret_cast<Complex *>(N->L->x)[j];
1030  }
1031 
1032  return ret;
1033 
1034 #else
1035 
1036  return SparseComplexMatrix ();
1037 
1038 #endif
1039  }
1040 
1041  template <>
1042  SparseComplexMatrix
1043  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::R (bool econ) const
1044  {
1045 #if defined (HAVE_CXSPARSE)
1046  // Drop zeros from R and sort
1047  // FIXME: Is the double transpose to sort necessary?
1048 
1049  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1050  CXSPARSE_ZNAME (_dropzeros) (N->U);
1051  CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->U, 1);
1052  CXSPARSE_ZNAME (_spfree) (N->U);
1053  N->U = CXSPARSE_ZNAME (_transpose) (D, 1);
1054  CXSPARSE_ZNAME (_spfree) (D);
1055  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1056 
1057  octave_idx_type nc = N->U->n;
1058  octave_idx_type nz = N->U->nzmax;
1059 
1060  SparseComplexMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows),
1061  nc, nz);
1062 
1063  for (octave_idx_type j = 0; j < nc+1; j++)
1064  ret.xcidx (j) = N->U->p[j];
1065 
1066  for (octave_idx_type j = 0; j < nz; j++)
1067  {
1068  ret.xridx (j) = N->U->i[j];
1069  ret.xdata (j) = reinterpret_cast<Complex *>(N->U->x)[j];
1070  }
1071 
1072  return ret;
1073 
1074 #else
1075 
1076  octave_unused_parameter (econ);
1077 
1078  return SparseComplexMatrix ();
1079 
1080 #endif
1081  }
1082 
1083  template <>
1084  ComplexMatrix
1085  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::C (const ComplexMatrix& b) const
1086  {
1087 #if defined (HAVE_CXSPARSE)
1088  octave_idx_type b_nr = b.rows ();
1089  octave_idx_type b_nc = b.cols ();
1090  octave_idx_type nc = N->L->n;
1091  octave_idx_type nr = nrows;
1092  const cs_complex_t *bvec
1093  = reinterpret_cast<const cs_complex_t *> (b.fortran_vec ());
1094  ComplexMatrix ret (b_nr, b_nc);
1095  Complex *vec = ret.fortran_vec ();
1096 
1097  if (nr < 0 || nc < 0 || nr != b_nr)
1098  (*current_liboctave_error_handler) ("matrix dimension mismatch");
1099 
1100  if (nr == 0 || nc == 0 || b_nc == 0)
1101  ret = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
1102  else
1103  {
1104  OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
1105 
1106  for (volatile octave_idx_type j = 0, idx = 0;
1107  j < b_nc;
1108  j++, idx += b_nr)
1109  {
1110  octave_quit ();
1111 
1112  volatile octave_idx_type nm = (nr < nc ? nr : nc);
1113 
1114  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1115  CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + idx,
1116  reinterpret_cast<cs_complex_t *> (buf),
1117  b_nr);
1118  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1119 
1120  for (volatile octave_idx_type i = 0; i < nm; i++)
1121  {
1122  octave_quit ();
1123 
1124  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1125  CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i],
1126  reinterpret_cast<cs_complex_t *> (buf));
1127  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1128  }
1129 
1130  for (octave_idx_type i = 0; i < b_nr; i++)
1131  vec[i+idx] = buf[i];
1132  }
1133  }
1134 
1135  return ret;
1136 
1137 #else
1138 
1139  octave_unused_parameter (b);
1140 
1141  return ComplexMatrix ();
1142 
1143 #endif
1144  }
1145 
1146  template <>
1147  ComplexMatrix
1148  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::Q (void) const
1149  {
1150 #if defined (HAVE_CXSPARSE)
1151  octave_idx_type nc = N->L->n;
1152  octave_idx_type nr = nrows;
1153  ComplexMatrix ret (nr, nr);
1154  Complex *vec = ret.fortran_vec ();
1155 
1156  if (nr < 0 || nc < 0)
1157  (*current_liboctave_error_handler) ("matrix dimension mismatch");
1158 
1159  if (nr == 0 || nc == 0)
1160  ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0));
1161  else
1162  {
1163  OCTAVE_LOCAL_BUFFER (cs_complex_t, bvec, nr);
1164 
1165  for (octave_idx_type i = 0; i < nr; i++)
1166  bvec[i] = cs_complex_t (0.0, 0.0);
1167 
1168  OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
1169 
1170  for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
1171  {
1172  octave_quit ();
1173 
1174  bvec[j] = cs_complex_t (1.0, 0.0);
1175 
1176  volatile octave_idx_type nm = (nr < nc ? nr : nc);
1177 
1178  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1179  CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec,
1180  reinterpret_cast<cs_complex_t *> (buf),
1181  nr);
1182  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1183 
1184  for (volatile octave_idx_type i = 0; i < nm; i++)
1185  {
1186  octave_quit ();
1187 
1188  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1189  CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i],
1190  reinterpret_cast<cs_complex_t *> (buf));
1191  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1192  }
1193 
1194  for (octave_idx_type i = 0; i < nr; i++)
1195  vec[i+idx] = buf[i];
1196 
1197  bvec[j] = cs_complex_t (0.0, 0.0);
1198  }
1199  }
1200 
1201  return ret.hermitian ();
1202 
1203 #else
1204 
1205  return ComplexMatrix ();
1206 
1207 #endif
1208  }
1209 
1210  template <>
1211  template <>
1212  SparseComplexMatrix
1213  sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix,
1214  SparseComplexMatrix>
1215  (const SparseComplexMatrix& b, octave_idx_type& info) const
1216  {
1217  info = -1;
1218 
1219 #if defined (HAVE_CXSPARSE)
1220 
1221  octave_idx_type nr = nrows;
1222  octave_idx_type nc = ncols;
1223 
1224  octave_idx_type b_nr = b.rows ();
1225  octave_idx_type b_nc = b.cols ();
1226 
1227  SparseComplexMatrix x (nc, b_nc, b.nnz ());
1228  x.xcidx (0) = 0;
1229 
1230  volatile octave_idx_type x_nz = b.nnz ();
1231  volatile octave_idx_type ii = 0;
1232 
1233  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
1234  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
1235  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
1236 
1237  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1238  {
1239  octave_quit ();
1240 
1241  for (octave_idx_type j = 0; j < b_nr; j++)
1242  {
1243  Complex c = b.xelem (j,i);
1244  Xx[j] = c.real ();
1245  Xz[j] = c.imag ();
1246  }
1247 
1248  for (octave_idx_type j = nr; j < S->m2; j++)
1249  buf[j] = 0.;
1250 
1251  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1252  CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr);
1253  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1254 
1255  for (volatile octave_idx_type j = 0; j < nc; j++)
1256  {
1257  octave_quit ();
1258 
1259  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1260  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1261  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1262  }
1263 
1264  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1265  CXSPARSE_DNAME (_usolve) (N->U, buf);
1266  CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc);
1267  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1268 
1269  for (octave_idx_type j = nr; j < S->m2; j++)
1270  buf[j] = 0.;
1271 
1272  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1273  CXSPARSE_DNAME (_ipvec) (S->pinv, Xz, buf, nr);
1274  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1275 
1276  for (volatile octave_idx_type j = 0; j < nc; j++)
1277  {
1278  octave_quit ();
1279 
1280  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1281  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1282  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1283  }
1284 
1285  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1286  CXSPARSE_DNAME (_usolve) (N->U, buf);
1287  CXSPARSE_DNAME (_ipvec) (S->q, buf, Xz, nc);
1288  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1289 
1290  for (octave_idx_type j = 0; j < nc; j++)
1291  {
1292  Complex tmp = Complex (Xx[j], Xz[j]);
1293 
1294  if (tmp != 0.0)
1295  {
1296  if (ii == x_nz)
1297  {
1298  // Resize the sparse matrix
1299  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1300  sz = (sz > 10 ? sz : 10) + x_nz;
1301  x.change_capacity (sz);
1302  x_nz = sz;
1303  }
1304 
1305  x.xdata (ii) = tmp;
1306  x.xridx (ii++) = j;
1307  }
1308  }
1309 
1310  x.xcidx (i+1) = ii;
1311  }
1312 
1313  info = 0;
1314 
1315  return x;
1316 
1317 #else
1318 
1319  octave_unused_parameter (b);
1320 
1321  return SparseComplexMatrix ();
1322 
1323 #endif
1324  }
1325 
1326  template <>
1327  template <>
1328  SparseComplexMatrix
1329  sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix,
1330  SparseComplexMatrix>
1331  (const SparseComplexMatrix& b, octave_idx_type& info) const
1332  {
1333  info = -1;
1334 
1335 #if defined (HAVE_CXSPARSE)
1336 
1337  // These are swapped because the original matrix was transposed in
1338  // sparse_qr<SparseMatrix>::solve.
1339 
1340  octave_idx_type nr = ncols;
1341  octave_idx_type nc = nrows;
1342 
1343  octave_idx_type b_nr = b.rows ();
1344  octave_idx_type b_nc = b.cols ();
1345 
1346  SparseComplexMatrix x (nc, b_nc, b.nnz ());
1347  x.xcidx (0) = 0;
1348 
1349  volatile octave_idx_type x_nz = b.nnz ();
1350  volatile octave_idx_type ii = 0;
1351  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1352 
1353  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
1354  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
1355  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
1356 
1357  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1358  {
1359  octave_quit ();
1360 
1361  for (octave_idx_type j = 0; j < b_nr; j++)
1362  {
1363  Complex c = b.xelem (j,i);
1364  Xx[j] = c.real ();
1365  Xz[j] = c.imag ();
1366  }
1367 
1368  for (octave_idx_type j = nr; j < nbuf; j++)
1369  buf[j] = 0.;
1370 
1371  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1372  CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr);
1373  CXSPARSE_DNAME (_utsolve) (N->U, buf);
1374  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1375 
1376  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
1377  {
1378  octave_quit ();
1379 
1380  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1381  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1382  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1383  }
1384 
1385  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1386  CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc);
1387  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1388 
1389  for (octave_idx_type j = nr; j < nbuf; j++)
1390  buf[j] = 0.;
1391 
1392  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1393  CXSPARSE_DNAME (_pvec) (S->q, Xz, buf, nr);
1394  CXSPARSE_DNAME (_utsolve) (N->U, buf);
1395  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1396 
1397  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
1398  {
1399  octave_quit ();
1400 
1401  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1402  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1403  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1404  }
1405 
1406  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1407  CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xz, nc);
1408  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1409 
1410  for (octave_idx_type j = 0; j < nc; j++)
1411  {
1412  Complex tmp = Complex (Xx[j], Xz[j]);
1413 
1414  if (tmp != 0.0)
1415  {
1416  if (ii == x_nz)
1417  {
1418  // Resize the sparse matrix
1419  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1420  sz = (sz > 10 ? sz : 10) + x_nz;
1421  x.change_capacity (sz);
1422  x_nz = sz;
1423  }
1424 
1425  x.xdata (ii) = tmp;
1426  x.xridx (ii++) = j;
1427  }
1428  }
1429 
1430  x.xcidx (i+1) = ii;
1431  }
1432 
1433  info = 0;
1434 
1435  x.maybe_compress ();
1436 
1437  return x;
1438 
1439 #else
1440 
1441  octave_unused_parameter (b);
1442 
1443  return SparseComplexMatrix ();
1444 
1445 #endif
1446  }
1447 
1448  template <>
1449  template <>
1450  ComplexMatrix
1451  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<double>,
1452  ComplexMatrix>
1453  (const MArray<double>& b, octave_idx_type& info) const
1454  {
1455  info = -1;
1456 
1457 #if defined (HAVE_CXSPARSE)
1458 
1459  octave_idx_type nr = nrows;
1460  octave_idx_type nc = ncols;
1461 
1462  octave_idx_type b_nc = b.cols ();
1463  octave_idx_type b_nr = b.rows ();
1464 
1465  ComplexMatrix x (nc, b_nc);
1466  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
1467 
1468  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
1469  OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
1470 
1471  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1472  {
1473  octave_quit ();
1474 
1475  for (octave_idx_type j = 0; j < b_nr; j++)
1476  Xx[j] = b.xelem (j,i);
1477 
1478  for (octave_idx_type j = nr; j < S->m2; j++)
1479  buf[j] = cs_complex_t (0.0, 0.0);
1480 
1481  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1482  CXSPARSE_ZNAME (_ipvec) (S->pinv,
1483  reinterpret_cast<cs_complex_t *>(Xx),
1484  buf, nr);
1485  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1486 
1487  for (volatile octave_idx_type j = 0; j < nc; j++)
1488  {
1489  octave_quit ();
1490 
1491  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1492  CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
1493  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1494  }
1495 
1496  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1497  CXSPARSE_ZNAME (_usolve) (N->U, buf);
1498  CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc);
1499  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1500  }
1501 
1502  info = 0;
1503 
1504  return x;
1505 
1506 #else
1507 
1508  octave_unused_parameter (b);
1509 
1510  return ComplexMatrix ();
1511 
1512 #endif
1513  }
1514 
1515  template <>
1516  template <>
1517  ComplexMatrix
1518  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<double>,
1519  ComplexMatrix>
1520  (const MArray<double>& b, octave_idx_type& info) const
1521  {
1522  info = -1;
1523 
1524 #if defined (HAVE_CXSPARSE)
1525 
1526  // These are swapped because the original matrix was transposed in
1527  // sparse_qr<SparseComplexMatrix>::solve.
1528 
1529  octave_idx_type nr = ncols;
1530  octave_idx_type nc = nrows;
1531 
1532  octave_idx_type b_nc = b.cols ();
1533  octave_idx_type b_nr = b.rows ();
1534 
1535  ComplexMatrix x (nc, b_nc);
1536  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
1537 
1538  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1539 
1540  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
1541  OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
1542  OCTAVE_LOCAL_BUFFER (double, B, nr);
1543 
1544  for (octave_idx_type i = 0; i < nr; i++)
1545  B[i] = N->B[i];
1546 
1547  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1548  {
1549  octave_quit ();
1550 
1551  for (octave_idx_type j = 0; j < b_nr; j++)
1552  Xx[j] = b.xelem (j,i);
1553 
1554  for (octave_idx_type j = nr; j < nbuf; j++)
1555  buf[j] = cs_complex_t (0.0, 0.0);
1556 
1557  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1558  CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx),
1559  buf, nr);
1560  CXSPARSE_ZNAME (_utsolve) (N->U, buf);
1561  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1562 
1563  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
1564  {
1565  octave_quit ();
1566 
1567  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1568  CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
1569  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1570  }
1571 
1572  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1573  CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc);
1574  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1575  }
1576 
1577  info = 0;
1578 
1579  return x;
1580 
1581 #else
1582 
1583  octave_unused_parameter (b);
1584 
1585  return ComplexMatrix ();
1586 
1587 #endif
1588  }
1589 
1590  template <>
1591  template <>
1592  SparseComplexMatrix
1593  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseMatrix,
1594  SparseComplexMatrix>
1595  (const SparseMatrix& b, octave_idx_type& info) const
1596  {
1597  info = -1;
1598 
1599 #if defined (HAVE_CXSPARSE)
1600 
1601  octave_idx_type nr = nrows;
1602  octave_idx_type nc = ncols;
1603 
1604  octave_idx_type b_nc = b.cols ();
1605  octave_idx_type b_nr = b.rows ();
1606 
1607  SparseComplexMatrix x (nc, b_nc, b.nnz ());
1608  x.xcidx (0) = 0;
1609 
1610  volatile octave_idx_type x_nz = b.nnz ();
1611  volatile octave_idx_type ii = 0;
1612 
1613  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
1614  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
1615 
1616  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1617  {
1618  octave_quit ();
1619 
1620  for (octave_idx_type j = 0; j < b_nr; j++)
1621  Xx[j] = b.xelem (j,i);
1622 
1623  for (octave_idx_type j = nr; j < S->m2; j++)
1624  buf[j] = cs_complex_t (0.0, 0.0);
1625 
1626  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1627  CXSPARSE_ZNAME (_ipvec) (S->pinv,
1628  reinterpret_cast<cs_complex_t *>(Xx),
1629  buf, nr);
1630  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1631 
1632  for (volatile octave_idx_type j = 0; j < nc; j++)
1633  {
1634  octave_quit ();
1635 
1636  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1637  CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
1638  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1639  }
1640 
1641  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1642  CXSPARSE_ZNAME (_usolve) (N->U, buf);
1643  CXSPARSE_ZNAME (_ipvec) (S->q, buf,
1644  reinterpret_cast<cs_complex_t *> (Xx),
1645  nc);
1646  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1647 
1648  for (octave_idx_type j = 0; j < nc; j++)
1649  {
1650  Complex tmp = Xx[j];
1651 
1652  if (tmp != 0.0)
1653  {
1654  if (ii == x_nz)
1655  {
1656  // Resize the sparse matrix
1657  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1658  sz = (sz > 10 ? sz : 10) + x_nz;
1659  x.change_capacity (sz);
1660  x_nz = sz;
1661  }
1662 
1663  x.xdata (ii) = tmp;
1664  x.xridx (ii++) = j;
1665  }
1666  }
1667 
1668  x.xcidx (i+1) = ii;
1669  }
1670 
1671  info = 0;
1672 
1673  x.maybe_compress ();
1674 
1675  return x;
1676 
1677 #else
1678 
1679  octave_unused_parameter (b);
1680 
1681  return SparseComplexMatrix ();
1682 
1683 #endif
1684  }
1685 
1686  template <>
1687  template <>
1688  SparseComplexMatrix
1689  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseMatrix,
1690  SparseComplexMatrix>
1691  (const SparseMatrix& b, octave_idx_type& info) const
1692  {
1693  info = -1;
1694 
1695 #if defined (HAVE_CXSPARSE)
1696 
1697  // These are swapped because the original matrix was transposed in
1698  // sparse_qr<SparseComplexMatrix>::solve.
1699 
1700  octave_idx_type nr = ncols;
1701  octave_idx_type nc = nrows;
1702 
1703  octave_idx_type b_nc = b.cols ();
1704  octave_idx_type b_nr = b.rows ();
1705 
1706  SparseComplexMatrix x (nc, b_nc, b.nnz ());
1707  x.xcidx (0) = 0;
1708 
1709  volatile octave_idx_type x_nz = b.nnz ();
1710  volatile octave_idx_type ii = 0;
1711  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1712 
1713  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
1714  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
1715  OCTAVE_LOCAL_BUFFER (double, B, nr);
1716 
1717  for (octave_idx_type i = 0; i < nr; i++)
1718  B[i] = N->B[i];
1719 
1720  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1721  {
1722  octave_quit ();
1723 
1724  for (octave_idx_type j = 0; j < b_nr; j++)
1725  Xx[j] = b.xelem (j,i);
1726 
1727  for (octave_idx_type j = nr; j < nbuf; j++)
1728  buf[j] = cs_complex_t (0.0, 0.0);
1729 
1730  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1731  CXSPARSE_ZNAME (_pvec) (S->q,
1732  reinterpret_cast<cs_complex_t *>(Xx),
1733  buf, nr);
1734  CXSPARSE_ZNAME (_utsolve) (N->U, buf);
1735  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1736 
1737  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
1738  {
1739  octave_quit ();
1740 
1741  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1742  CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
1743  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1744  }
1745 
1746  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1747  CXSPARSE_ZNAME (_pvec) (S->pinv, buf,
1748  reinterpret_cast<cs_complex_t *> (Xx),
1749  nc);
1750  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1751 
1752  for (octave_idx_type j = 0; j < nc; j++)
1753  {
1754  Complex tmp = Xx[j];
1755 
1756  if (tmp != 0.0)
1757  {
1758  if (ii == x_nz)
1759  {
1760  // Resize the sparse matrix
1761  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1762  sz = (sz > 10 ? sz : 10) + x_nz;
1763  x.change_capacity (sz);
1764  x_nz = sz;
1765  }
1766 
1767  x.xdata (ii) = tmp;
1768  x.xridx (ii++) = j;
1769  }
1770  }
1771 
1772  x.xcidx (i+1) = ii;
1773  }
1774 
1775  info = 0;
1776 
1777  x.maybe_compress ();
1778 
1779  return x;
1780 
1781 #else
1782 
1783  octave_unused_parameter (b);
1784 
1785  return SparseComplexMatrix ();
1786 
1787 #endif
1788  }
1789 
1790  template <>
1791  template <>
1792  ComplexMatrix
1793  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>,
1794  ComplexMatrix>
1795  (const MArray<Complex>& b, octave_idx_type& info) const
1796  {
1797  info = -1;
1798 
1799 #if defined (HAVE_CXSPARSE)
1800 
1801  octave_idx_type nr = nrows;
1802  octave_idx_type nc = ncols;
1803 
1804  octave_idx_type b_nc = b.cols ();
1805  octave_idx_type b_nr = b.rows ();
1806 
1807  const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>
1808  (b.fortran_vec ());
1809 
1810  ComplexMatrix x (nc, b_nc);
1811  cs_complex_t *vec = reinterpret_cast<cs_complex_t *>
1812  (x.fortran_vec ());
1813 
1814  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
1815 
1816  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
1817  i++, idx+=nc, bidx+=b_nr)
1818  {
1819  octave_quit ();
1820 
1821  for (octave_idx_type j = nr; j < S->m2; j++)
1822  buf[j] = cs_complex_t (0.0, 0.0);
1823 
1824  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1825  CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + bidx, buf, nr);
1826  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1827 
1828  for (volatile octave_idx_type j = 0; j < nc; j++)
1829  {
1830  octave_quit ();
1831 
1832  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1833  CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
1834  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1835  }
1836 
1837  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1838  CXSPARSE_ZNAME (_usolve) (N->U, buf);
1839  CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc);
1840  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1841  }
1842 
1843  info = 0;
1844 
1845  return x;
1846 
1847 #else
1848 
1849  octave_unused_parameter (b);
1850 
1851  return ComplexMatrix ();
1852 
1853 #endif
1854  }
1855 
1856  template <>
1857  template <>
1858  ComplexMatrix
1859  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>,
1860  ComplexMatrix>
1861  (const MArray<Complex>& b, octave_idx_type& info) const
1862  {
1863  info = -1;
1864 
1865 #if defined (HAVE_CXSPARSE)
1866 
1867  // These are swapped because the original matrix was transposed in
1868  // sparse_qr<SparseComplexMatrix>::solve.
1869 
1870  octave_idx_type nr = ncols;
1871  octave_idx_type nc = nrows;
1872 
1873  octave_idx_type b_nc = b.cols ();
1874  octave_idx_type b_nr = b.rows ();
1875 
1876  const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>
1877  (b.fortran_vec ());
1878 
1879  ComplexMatrix x (nc, b_nc);
1880  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
1881 
1882  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1883 
1884  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
1885  OCTAVE_LOCAL_BUFFER (double, B, nr);
1886 
1887  for (octave_idx_type i = 0; i < nr; i++)
1888  B[i] = N->B[i];
1889 
1890  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
1891  i++, idx+=nc, bidx+=b_nr)
1892  {
1893  octave_quit ();
1894 
1895  for (octave_idx_type j = nr; j < nbuf; j++)
1896  buf[j] = cs_complex_t (0.0, 0.0);
1897 
1898  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1899  CXSPARSE_ZNAME (_pvec) (S->q, bvec + bidx, buf, nr);
1900  CXSPARSE_ZNAME (_utsolve) (N->U, buf);
1901  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1902 
1903  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
1904  {
1905  octave_quit ();
1906 
1907  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1908  CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
1909  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1910  }
1911 
1912  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1913  CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc);
1914  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1915  }
1916 
1917  info = 0;
1918 
1919  return x;
1920 
1921 #else
1922 
1923  octave_unused_parameter (b);
1924 
1925  return ComplexMatrix ();
1926 
1927 #endif
1928  }
1929 
1930  template <>
1931  template <>
1932  SparseComplexMatrix
1933  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix, SparseComplexMatrix>
1934  (const SparseComplexMatrix& b, octave_idx_type& info) const
1935  {
1936  info = -1;
1937 
1938 #if defined (HAVE_CXSPARSE)
1939 
1940  octave_idx_type nr = nrows;
1941  octave_idx_type nc = ncols;
1942 
1943  octave_idx_type b_nc = b.cols ();
1944  octave_idx_type b_nr = b.rows ();
1945 
1946  SparseComplexMatrix x (nc, b_nc, b.nnz ());
1947  x.xcidx (0) = 0;
1948 
1949  volatile octave_idx_type x_nz = b.nnz ();
1950  volatile octave_idx_type ii = 0;
1951 
1952  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
1953  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
1954 
1955  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1956  {
1957  octave_quit ();
1958 
1959  for (octave_idx_type j = 0; j < b_nr; j++)
1960  Xx[j] = b.xelem (j,i);
1961 
1962  for (octave_idx_type j = nr; j < S->m2; j++)
1963  buf[j] = cs_complex_t (0.0, 0.0);
1964 
1965  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1966  CXSPARSE_ZNAME (_ipvec) (S->pinv,
1967  reinterpret_cast<cs_complex_t *>(Xx),
1968  buf, nr);
1969  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1970 
1971  for (volatile octave_idx_type j = 0; j < nc; j++)
1972  {
1973  octave_quit ();
1974 
1975  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1976  CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
1977  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1978  }
1979 
1980  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1981  CXSPARSE_ZNAME (_usolve) (N->U, buf);
1982  CXSPARSE_ZNAME (_ipvec) (S->q, buf,
1983  reinterpret_cast<cs_complex_t *> (Xx),
1984  nc);
1985  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1986 
1987  for (octave_idx_type j = 0; j < nc; j++)
1988  {
1989  Complex tmp = Xx[j];
1990 
1991  if (tmp != 0.0)
1992  {
1993  if (ii == x_nz)
1994  {
1995  // Resize the sparse matrix
1996  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1997  sz = (sz > 10 ? sz : 10) + x_nz;
1998  x.change_capacity (sz);
1999  x_nz = sz;
2000  }
2001 
2002  x.xdata (ii) = tmp;
2003  x.xridx (ii++) = j;
2004  }
2005  }
2006 
2007  x.xcidx (i+1) = ii;
2008  }
2009 
2010  info = 0;
2011 
2012  x.maybe_compress ();
2013 
2014  return x;
2015 
2016 #else
2017 
2018  octave_unused_parameter (b);
2019 
2020  return SparseComplexMatrix ();
2021 
2022 #endif
2023  }
2024 
2025  template <>
2026  template <>
2027  SparseComplexMatrix
2028  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix, SparseComplexMatrix>
2029  (const SparseComplexMatrix& b, octave_idx_type& info) const
2030  {
2031  info = -1;
2032 
2033 #if defined (HAVE_CXSPARSE)
2034 
2035  // These are swapped because the original matrix was transposed in
2036  // sparse_qr<SparseComplexMatrix>::solve.
2037 
2038  octave_idx_type nr = ncols;
2039  octave_idx_type nc = nrows;
2040 
2041  octave_idx_type b_nc = b.cols ();
2042  octave_idx_type b_nr = b.rows ();
2043 
2044  SparseComplexMatrix x (nc, b_nc, b.nnz ());
2045  x.xcidx (0) = 0;
2046 
2047  volatile octave_idx_type x_nz = b.nnz ();
2048  volatile octave_idx_type ii = 0;
2049  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
2050 
2051  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
2052  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
2053  OCTAVE_LOCAL_BUFFER (double, B, nr);
2054 
2055  for (octave_idx_type i = 0; i < nr; i++)
2056  B[i] = N->B[i];
2057 
2058  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
2059  {
2060  octave_quit ();
2061 
2062  for (octave_idx_type j = 0; j < b_nr; j++)
2063  Xx[j] = b.xelem (j,i);
2064 
2065  for (octave_idx_type j = nr; j < nbuf; j++)
2066  buf[j] = cs_complex_t (0.0, 0.0);
2067 
2068  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2069  CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx),
2070  buf, nr);
2071  CXSPARSE_ZNAME (_utsolve) (N->U, buf);
2072  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2073 
2074  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
2075  {
2076  octave_quit ();
2077 
2078  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2079  CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
2080  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2081  }
2082 
2083  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2084  CXSPARSE_ZNAME (_pvec) (S->pinv, buf,
2085  reinterpret_cast<cs_complex_t *>(Xx), nc);
2086  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2087 
2088  for (octave_idx_type j = 0; j < nc; j++)
2089  {
2090  Complex tmp = Xx[j];
2091 
2092  if (tmp != 0.0)
2093  {
2094  if (ii == x_nz)
2095  {
2096  // Resize the sparse matrix
2097  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
2098  sz = (sz > 10 ? sz : 10) + x_nz;
2099  x.change_capacity (sz);
2100  x_nz = sz;
2101  }
2102 
2103  x.xdata (ii) = tmp;
2104  x.xridx (ii++) = j;
2105  }
2106  }
2107 
2108  x.xcidx (i+1) = ii;
2109  }
2110 
2111  info = 0;
2112 
2113  x.maybe_compress ();
2114 
2115  return x;
2116 
2117 #else
2118 
2119  octave_unused_parameter (b);
2120 
2121  return SparseComplexMatrix ();
2122 
2123 #endif
2124  }
2125 
2126  template <typename SPARSE_T>
2127  sparse_qr<SPARSE_T>::sparse_qr (void)
2128  : rep (new sparse_qr_rep (SPARSE_T (), 0))
2129  { }
2130 
2131  template <typename SPARSE_T>
2132  sparse_qr<SPARSE_T>::sparse_qr (const SPARSE_T& a, int order)
2133  : rep (new sparse_qr_rep (a, order))
2134  { }
2135 
2136  template <typename SPARSE_T>
2137  sparse_qr<SPARSE_T>::sparse_qr (const sparse_qr<SPARSE_T>& a)
2138  : rep (a.rep)
2139  {
2140  rep->count++;
2141  }
2142 
2143  template <typename SPARSE_T>
2144  sparse_qr<SPARSE_T>::~sparse_qr (void)
2145  {
2146  if (--rep->count == 0)
2147  delete rep;
2148  }
2149 
2150  template <typename SPARSE_T>
2151  sparse_qr<SPARSE_T>&
2152  sparse_qr<SPARSE_T>::operator = (const sparse_qr<SPARSE_T>& a)
2153  {
2154  if (this != &a)
2155  {
2156  if (--rep->count == 0)
2157  delete rep;
2158 
2159  rep = a.rep;
2160  rep->count++;
2161  }
2162 
2163  return *this;
2164  }
2165 
2166  template <typename SPARSE_T>
2167  bool
2168  sparse_qr<SPARSE_T>::ok (void) const
2169  {
2170  return rep->ok ();
2171  }
2172 
2173  template <typename SPARSE_T>
2174  SPARSE_T
2175  sparse_qr<SPARSE_T>::V (void) const
2176  {
2177  return rep->V ();
2178  }
2179 
2180  template <typename SPARSE_T>
2181  ColumnVector
2182  sparse_qr<SPARSE_T>::Pinv (void) const
2183  {
2184  return rep->P ();
2185  }
2186 
2187  template <typename SPARSE_T>
2188  ColumnVector
2189  sparse_qr<SPARSE_T>::P (void) const
2190  {
2191  return rep->P ();
2192  }
2193 
2194  template <typename SPARSE_T>
2195  SPARSE_T
2196  sparse_qr<SPARSE_T>::R (bool econ) const
2197  {
2198  return rep->R (econ);
2199  }
2200 
2201  template <typename SPARSE_T>
2202  typename SPARSE_T::dense_matrix_type
2203  sparse_qr<SPARSE_T>::C (const typename SPARSE_T::dense_matrix_type& b) const
2204  {
2205  return rep->C (b);
2206  }
2207 
2208  template <typename SPARSE_T>
2209  typename SPARSE_T::dense_matrix_type
2210  sparse_qr<SPARSE_T>::Q (void) const
2211  {
2212  return rep->Q ();
2213  }
2214 
2215  // FIXME: Why is the "order" of the QR calculation as used in the
2216  // CXSparse function sqr 3 for real matrices and 2 for complex? These
2217  // values seem to be required but there was no explanation in David
2218  // Bateman's original code.
2219 
2220  template <typename SPARSE_T>
2221  class
2222  cxsparse_defaults
2223  {
2224  public:
2225  enum { order = -1 };
2226  };
2227 
2228  template <>
2229  class
2230  cxsparse_defaults<SparseMatrix>
2231  {
2232  public:
2233  enum { order = 3 };
2234  };
2235 
2236  template <>
2237  class
2238  cxsparse_defaults<SparseComplexMatrix>
2239  {
2240  public:
2241  enum { order = 2 };
2242  };
2243 
2244  template <typename SPARSE_T>
2245  template <typename RHS_T, typename RET_T>
2246  RET_T
2247  sparse_qr<SPARSE_T>::solve (const SPARSE_T& a, const RHS_T& b,
2248  octave_idx_type& info)
2249  {
2250  info = -1;
2251 
2252  octave_idx_type nr = a.rows ();
2253  octave_idx_type nc = a.cols ();
2254 
2255  octave_idx_type b_nc = b.cols ();
2256  octave_idx_type b_nr = b.rows ();
2257 
2258  int order = cxsparse_defaults<SPARSE_T>::order;
2259 
2260  if (nr < 0 || nc < 0 || nr != b_nr)
2261  (*current_liboctave_error_handler)
2262  ("matrix dimension mismatch in solution of minimum norm problem");
2263 
2264  if (nr == 0 || nc == 0 || b_nc == 0)
2265  {
2266  info = 0;
2267 
2268  return RET_T (nc, b_nc, 0.0);
2269  }
2270  else if (nr >= nc)
2271  {
2272  sparse_qr<SPARSE_T> q (a, order);
2273 
2274  return q.ok () ? q.tall_solve<RHS_T, RET_T> (b, info) : RET_T ();
2275  }
2276  else
2277  {
2278  sparse_qr<SPARSE_T> q (a.hermitian (), order);
2279 
2280  return q.ok () ? q.wide_solve<RHS_T, RET_T> (b, info) : RET_T ();
2281  }
2282  }
2283 
2284  template <typename SPARSE_T>
2285  template <typename RHS_T, typename RET_T>
2286  RET_T
2287  sparse_qr<SPARSE_T>::tall_solve (const RHS_T& b, octave_idx_type& info) const
2288  {
2289  return rep->template tall_solve<RHS_T, RET_T> (b, info);
2290  }
2291 
2292  template <typename SPARSE_T>
2293  template <typename RHS_T, typename RET_T>
2294  RET_T
2295  sparse_qr<SPARSE_T>::wide_solve (const RHS_T& b, octave_idx_type& info) const
2296  {
2297  return rep->template wide_solve<RHS_T, RET_T> (b, info);
2298  }
2299 
2300  Matrix
2301  qrsolve (const SparseMatrix& a, const MArray<double>& b,
2302  octave_idx_type& info)
2303  {
2304  return sparse_qr<SparseMatrix>::solve<MArray<double>, Matrix> (a, b,
2305  info);
2306  }
2307 
2308  SparseMatrix
2309  qrsolve (const SparseMatrix& a, const SparseMatrix& b,
2310  octave_idx_type& info)
2311  {
2312  return sparse_qr<SparseMatrix>::solve<SparseMatrix, SparseMatrix> (a, b,
2313  info);
2314  }
2315 
2316  ComplexMatrix
2317  qrsolve (const SparseMatrix& a, const MArray<Complex>& b,
2318  octave_idx_type& info)
2319  {
2320  return sparse_qr<SparseMatrix>::solve<MArray<Complex>,
2321  ComplexMatrix> (a, b, info);
2322  }
2323 
2324  SparseComplexMatrix
2325  qrsolve (const SparseMatrix& a, const SparseComplexMatrix& b,
2326  octave_idx_type& info)
2327  {
2328  return sparse_qr<SparseMatrix>::solve<SparseComplexMatrix,
2329  SparseComplexMatrix> (a, b, info);
2330  }
2331 
2332  ComplexMatrix
2333  qrsolve (const SparseComplexMatrix& a, const MArray<double>& b,
2334  octave_idx_type& info)
2335  {
2336  return sparse_qr<SparseComplexMatrix>::solve<MArray<double>,
2337  ComplexMatrix> (a, b, info);
2338  }
2339 
2340  SparseComplexMatrix
2341  qrsolve (const SparseComplexMatrix& a, const SparseMatrix& b,
2342  octave_idx_type& info)
2343  {
2344  return sparse_qr<SparseComplexMatrix>::solve<SparseMatrix,
2345  SparseComplexMatrix>
2346  (a, b, info);
2347  }
2348 
2349  ComplexMatrix
2350  qrsolve (const SparseComplexMatrix& a, const MArray<Complex>& b,
2351  octave_idx_type& info)
2352  {
2353  return sparse_qr<SparseComplexMatrix>::solve<MArray<Complex>,
2354  ComplexMatrix> (a, b, info);
2355  }
2356 
2357  SparseComplexMatrix
2358  qrsolve (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
2359  octave_idx_type& info)
2360  {
2361  return sparse_qr<SparseComplexMatrix>::solve<SparseComplexMatrix,
2362  SparseComplexMatrix>
2363  (a, b, info);
2364  }
2365 
2366  // Instantiations we need.
2367 
2368  template class sparse_qr<SparseMatrix>;
2369 
2370  template class sparse_qr<SparseComplexMatrix>;
2371  }
2372 }
octave::math::sparse_qr::sparse_qr_rep::nrows
octave_idx_type nrows
Definition: sparse-qr.cc:119
octave::math::sparse_qr::sparse_qr_rep::ncols
octave_idx_type ncols
Definition: sparse-qr.cc:120
octave::math::cxsparse_types< SparseMatrix >::numeric_type
cs_din numeric_type
Definition: sparse-qr.cc:58
for
for(octave_idx_type n=0;n< hcv.numel();n++)
Definition: graphics.cc:10831
B
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
Definition: lo-lapack-proto.h:652
octave::math::qrsolve
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
Definition: sparse-qr.cc:2301
octave::math::sparse_qr::Q
SPARSE_T::dense_matrix_type Q(void) const
Definition: sparse-qr.cc:2210
octave::math::sparse_qr::wide_solve
RET_T wide_solve(const RHS_T &b, octave_idx_type &info) const
Definition: sparse-qr.cc:2295
N
F77_RET_T const F77_INT & N
Definition: lo-lapack-proto.h:652
octave::math::cxsparse_types< SparseMatrix >::symbolic_type
cs_dis symbolic_type
Definition: sparse-qr.cc:57
octave::math::cxsparse_types
Definition: sparse-qr.cc:46
octave::math::sparse_qr::Pinv
ColumnVector Pinv(void) const
Definition: sparse-qr.cc:2182
octave::math::sparse_qr::sparse_qr
sparse_qr(void)
Definition: sparse-qr.cc:2127
Array::fortran_vec
const T * fortran_vec(void) const
Definition: Array.h:584
octave::math::sparse_qr::tall_solve
RET_T tall_solve(const RHS_T &b, octave_idx_type &info) const
Definition: sparse-qr.cc:2287
octave::math::sparse_qr::operator=
sparse_qr & operator=(const sparse_qr &a)
Definition: sparse-qr.cc:2152
b_nr
octave_idx_type b_nr
Definition: sylvester.cc:76
octave::math::sparse_qr::sparse_qr_rep::count
refcount< int > count
Definition: sparse-qr.cc:117
c
nd example oindent opens the file binary numeric values will be read assuming they are stored in IEEE format with the least significant bit and then converted to the native representation Opening a file that is already open simply opens it again and returns a separate file id It is not an error to open a file several though writing to the same file through several different file ids may produce unexpected results The possible values of text mode reading and writing automatically converts linefeeds to the appropriate line end character for the you may append a you must also open the file in binary mode The parameter conversions are currently only supported for and permissions will be set to and then everything is written in a single operation This is very efficient and improves performance c
Definition: file-io.cc:587
s
s
Definition: file-io.cc:2729
octave::math::sparse_qr::sparse_qr_rep::P
ColumnVector P(void) const
Definition: sparse-qr.cc:156
octave::math::sparse_qr::sparse_qr_rep::R
SPARSE_T R(bool econ) const
a
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:400
octave::math::cxsparse_types< SparseComplexMatrix >
Definition: sparse-qr.cc:66
octave::math::sparse_qr::sparse_qr_rep
Definition: sparse-qr.cc:80
octave::math::sparse_qr::C
SPARSE_T::dense_matrix_type C(const typename SPARSE_T::dense_matrix_type &b) const
Definition: sparse-qr.cc:2203
octave::math::sparse_qr::sparse_qr_rep::Q
SPARSE_T::dense_matrix_type Q(void) const
octave::math::sparse_qr::~sparse_qr
~sparse_qr(void)
Definition: sparse-qr.cc:2144
Matrix::transpose
Matrix transpose(void) const
Definition: dMatrix.h:132
nm
OCTAVE_EXPORT octave_value_list isdir nd deftypefn *std::string nm
Definition: utils.cc:975
A
F77_RET_T const F77_INT F77_CMPLX * A
Definition: lo-lapack-proto.h:652
octave::suitesparse_integer
int suitesparse_integer
Definition: oct-sparse.h:168
octave::math::sparse_qr::sparse_qr_rep::V
SPARSE_T V(void) const
tmp
double tmp
Definition: data.cc:6252
octave::math::sparse_qr
Definition: sparse-qr.h:46
Sparse::xridx
octave_idx_type * xridx(void)
Definition: Sparse.h:501
octave::math::sparse_qr::sparse_qr_rep::C
SPARSE_T::dense_matrix_type C(const typename SPARSE_T::dense_matrix_type &b) const
octave::math::cxsparse_types< SparseMatrix >
Definition: sparse-qr.cc:52
octave::math::sparse_qr::R
SPARSE_T R(bool econ=false) const
Definition: sparse-qr.cc:2196
Matrix
Definition: dMatrix.h:36
sz
sz
Definition: data.cc:5264
octave::math::cxsparse_types< SparseComplexMatrix >::symbolic_type
cs_cis symbolic_type
Definition: sparse-qr.cc:71
octave::math::cxsparse_defaults
Definition: sparse-qr.cc:2221
Array::xelem
T & xelem(octave_idx_type n)
Definition: Array.h:458
octave::math::sparse_qr::sparse_qr_rep::N
cxsparse_types< SPARSE_T >::numeric_type * N
Definition: sparse-qr.cc:123
CXSPARSE_DNAME
#define CXSPARSE_DNAME(name)
Definition: oct-sparse.h:144
defined
defaults to zero A value of zero computes the digamma a value the trigamma and so on The digamma function is defined
Definition: psi.cc:68
octave::math::sparse_qr::V
SPARSE_T V(void) const
Definition: sparse-qr.cc:2175
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
#define END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
Definition: quit.h:217
octave::math::sparse_qr::solve
static RET_T solve(const SPARSE_T &a, const RHS_T &b, octave_idx_type &info)
Definition: sparse-qr.cc:2247
octave::math::cxsparse_defaults< SparseComplexMatrix >
Definition: sparse-qr.cc:2237
b_nc
octave_idx_type b_nc
Definition: sylvester.cc:77
Sparse::xdata
T * xdata(void)
Definition: Sparse.h:488
octave::math::sparse_qr::sparse_qr_rep::S
cxsparse_types< SPARSE_T >::symbolic_type * S
Definition: sparse-qr.cc:122
octave::math::sparse_qr::sparse_qr_rep::ok
bool ok(void) const
Definition: sparse-qr.cc:94
octave::math::sparse_qr::rep
sparse_qr_rep * rep
Definition: sparse-qr.h:84
b
b
Definition: cellfun.cc:400
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:41
octave::math::sparse_qr::ok
bool ok(void) const
Definition: sparse-qr.cc:2168
octave_idx_type
i
for i
Definition: data.cc:5264
octave::math::sparse_qr::sparse_qr_rep::sparse_qr_rep
sparse_qr_rep(const SPARSE_T &a, int order)
octave::math::cxsparse_defaults< SparseMatrix >
Definition: sparse-qr.cc:2229
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
#define BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
Definition: quit.h:213
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:31
Sparse::xcidx
octave_idx_type * xcidx(void)
Definition: Sparse.h:514
octave::math::sparse_qr::sparse_qr_rep::Pinv
ColumnVector Pinv(void) const
Definition: sparse-qr.cc:136
octave::math::sparse_qr::P
ColumnVector P(void) const
Definition: sparse-qr.cc:2189
octave::math::cxsparse_types< SparseComplexMatrix >::numeric_type
cs_cin numeric_type
Definition: sparse-qr.cc:72
CXSPARSE_ZNAME
#define CXSPARSE_ZNAME(name)
Definition: oct-sparse.h:145
octave::math::sparse_qr::sparse_qr_rep::~sparse_qr_rep
~sparse_qr_rep(void)
x
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 const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x
Definition: lo-slatec-proto.h:55
octave::to_suitesparse_intptr
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
Definition: oct-sparse.cc:48
ComplexMatrix::hermitian
ComplexMatrix hermitian(void) const
Definition: CMatrix.h:166