GNU Octave  6.2.0
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) 2005-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include "CMatrix.h"
31 #include "CSparse.h"
32 #include "MArray.h"
33 #include "dColVector.h"
34 #include "dMatrix.h"
35 #include "dSparse.h"
36 #include "lo-error.h"
37 #include "oct-locbuf.h"
38 #include "oct-refcount.h"
39 #include "oct-sparse.h"
40 #include "quit.h"
41 #include "sparse-qr.h"
42 
43 namespace octave
44 {
45  namespace math
46  {
47  template <typename SPARSE_T>
48  class
49  cxsparse_types
50  {
51  };
52 
53  template <>
54  class
55  cxsparse_types<SparseMatrix>
56  {
57  public:
58 #if defined (HAVE_CXSPARSE)
59  typedef CXSPARSE_DNAME (s) symbolic_type;
60  typedef CXSPARSE_DNAME (n) numeric_type;
61 #else
62  typedef void symbolic_type;
63  typedef void numeric_type;
64 #endif
65  };
66 
67  template <>
68  class
69  cxsparse_types<SparseComplexMatrix>
70  {
71  public:
72 #if defined (HAVE_CXSPARSE)
73  typedef CXSPARSE_ZNAME (s) symbolic_type;
74  typedef CXSPARSE_ZNAME (n) numeric_type;
75 #else
76  typedef void symbolic_type;
77  typedef void numeric_type;
78 #endif
79  };
80 
81  template <typename SPARSE_T>
82  class sparse_qr<SPARSE_T>::sparse_qr_rep
83  {
84  public:
85 
86  sparse_qr_rep (const SPARSE_T& a, int order);
87 
88  // No copying!
89 
90  sparse_qr_rep (const sparse_qr_rep&) = delete;
91 
92  sparse_qr_rep& operator = (const sparse_qr_rep&) = delete;
93 
94  ~sparse_qr_rep (void);
95 
96  bool ok (void) const
97  {
98 #if defined (HAVE_CXSPARSE)
99  return (N && S);
100 #else
101  return false;
102 #endif
103  }
104 
105  SPARSE_T V (void) const;
106 
107  ColumnVector Pinv (void) const;
108 
109  ColumnVector P (void) const;
110 
111  SPARSE_T R (bool econ) const;
112 
113  typename SPARSE_T::dense_matrix_type
114  C (const typename SPARSE_T::dense_matrix_type& b) const;
115 
116  typename SPARSE_T::dense_matrix_type
117  Q (void) const;
118 
119  refcount<octave_idx_type> count;
120 
121  octave_idx_type nrows;
122  octave_idx_type ncols;
123 
124  typename cxsparse_types<SPARSE_T>::symbolic_type *S;
125  typename cxsparse_types<SPARSE_T>::numeric_type *N;
126 
127  template <typename RHS_T, typename RET_T>
128  RET_T
129  tall_solve (const RHS_T& b, octave_idx_type& info) const;
130 
131  template <typename RHS_T, typename RET_T>
132  RET_T
133  wide_solve (const RHS_T& b, octave_idx_type& info) const;
134  };
135 
136  template <typename SPARSE_T>
137  ColumnVector
138  sparse_qr<SPARSE_T>::sparse_qr_rep::Pinv (void) const
139  {
140 #if defined (HAVE_CXSPARSE)
141 
142  ColumnVector ret (N->L->m);
143 
144  for (octave_idx_type i = 0; i < N->L->m; i++)
145  ret.xelem (i) = S->pinv[i];
146 
147  return ret;
148 
149 #else
150 
151  return ColumnVector ();
152 
153 #endif
154  }
155 
156  template <typename SPARSE_T>
157  ColumnVector
158  sparse_qr<SPARSE_T>::sparse_qr_rep::P (void) const
159  {
160 #if defined (HAVE_CXSPARSE)
161 
162  ColumnVector ret (N->L->m);
163 
164  for (octave_idx_type i = 0; i < N->L->m; i++)
165  ret.xelem (S->pinv[i]) = i;
166 
167  return ret;
168 
169 #else
170 
171  return ColumnVector ();
172 
173 #endif
174  }
175 
176  // Specializations.
177 
178  // Real-valued matrices.
179 
180  template <>
181  sparse_qr<SparseMatrix>::sparse_qr_rep::sparse_qr_rep
182  (const SparseMatrix& a, int order)
183  : count (1), nrows (a.rows ()), ncols (a.columns ())
184 #if defined (HAVE_CXSPARSE)
185  , S (nullptr), N (nullptr)
186 #endif
187  {
188 #if defined (HAVE_CXSPARSE)
189 
190  CXSPARSE_DNAME () A;
191 
192  A.nzmax = a.nnz ();
193  A.m = nrows;
194  A.n = ncols;
195  // Cast away const on A, with full knowledge that CSparse won't touch it
196  // Prevents the methods below making a copy of the data.
197  A.p = const_cast<suitesparse_integer *>
198  (to_suitesparse_intptr (a.cidx ()));
199  A.i = const_cast<suitesparse_integer *>
200  (to_suitesparse_intptr (a.ridx ()));
201  A.x = const_cast<double *> (a.data ());
202  A.nz = -1;
203 
204  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
205  S = CXSPARSE_DNAME (_sqr) (order, &A, 1);
206  N = CXSPARSE_DNAME (_qr) (&A, S);
207  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
208 
209  if (! N)
210  (*current_liboctave_error_handler)
211  ("sparse_qr: sparse matrix QR factorization filled");
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 #else
983 
984  octave_unused_parameter (order);
985 
986  (*current_liboctave_error_handler)
987  ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built");
988 
989 #endif
990  }
991 
992  template <>
993  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::~sparse_qr_rep (void)
994  {
995 #if defined (HAVE_CXSPARSE)
996  CXSPARSE_ZNAME (_sfree) (S);
997  CXSPARSE_ZNAME (_nfree) (N);
998 #endif
999  }
1000 
1001  template <>
1002  SparseComplexMatrix
1003  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::V (void) const
1004  {
1005 #if defined (HAVE_CXSPARSE)
1006  // Drop zeros from V and sort
1007  // FIXME: Is the double transpose to sort necessary?
1008 
1009  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1010  CXSPARSE_ZNAME (_dropzeros) (N->L);
1011  CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->L, 1);
1012  CXSPARSE_ZNAME (_spfree) (N->L);
1013  N->L = CXSPARSE_ZNAME (_transpose) (D, 1);
1014  CXSPARSE_ZNAME (_spfree) (D);
1015  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1016 
1017  octave_idx_type nc = N->L->n;
1018  octave_idx_type nz = N->L->nzmax;
1019  SparseComplexMatrix ret (N->L->m, nc, nz);
1020 
1021  for (octave_idx_type j = 0; j < nc+1; j++)
1022  ret.xcidx (j) = N->L->p[j];
1023 
1024  for (octave_idx_type j = 0; j < nz; j++)
1025  {
1026  ret.xridx (j) = N->L->i[j];
1027  ret.xdata (j) = reinterpret_cast<Complex *>(N->L->x)[j];
1028  }
1029 
1030  return ret;
1031 
1032 #else
1033 
1034  return SparseComplexMatrix ();
1035 
1036 #endif
1037  }
1038 
1039  template <>
1040  SparseComplexMatrix
1041  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::R (bool econ) const
1042  {
1043 #if defined (HAVE_CXSPARSE)
1044  // Drop zeros from R and sort
1045  // FIXME: Is the double transpose to sort necessary?
1046 
1047  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1048  CXSPARSE_ZNAME (_dropzeros) (N->U);
1049  CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->U, 1);
1050  CXSPARSE_ZNAME (_spfree) (N->U);
1051  N->U = CXSPARSE_ZNAME (_transpose) (D, 1);
1052  CXSPARSE_ZNAME (_spfree) (D);
1053  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1054 
1055  octave_idx_type nc = N->U->n;
1056  octave_idx_type nz = N->U->nzmax;
1057 
1058  SparseComplexMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows),
1059  nc, nz);
1060 
1061  for (octave_idx_type j = 0; j < nc+1; j++)
1062  ret.xcidx (j) = N->U->p[j];
1063 
1064  for (octave_idx_type j = 0; j < nz; j++)
1065  {
1066  ret.xridx (j) = N->U->i[j];
1067  ret.xdata (j) = reinterpret_cast<Complex *>(N->U->x)[j];
1068  }
1069 
1070  return ret;
1071 
1072 #else
1073 
1074  octave_unused_parameter (econ);
1075 
1076  return SparseComplexMatrix ();
1077 
1078 #endif
1079  }
1080 
1081  template <>
1082  ComplexMatrix
1083  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::C (const ComplexMatrix& b) const
1084  {
1085 #if defined (HAVE_CXSPARSE)
1086  octave_idx_type b_nr = b.rows ();
1087  octave_idx_type b_nc = b.cols ();
1088  octave_idx_type nc = N->L->n;
1089  octave_idx_type nr = nrows;
1090  const cs_complex_t *bvec
1091  = reinterpret_cast<const cs_complex_t *> (b.fortran_vec ());
1092  ComplexMatrix ret (b_nr, b_nc);
1093  Complex *vec = ret.fortran_vec ();
1094 
1095  if (nr < 0 || nc < 0 || nr != b_nr)
1096  (*current_liboctave_error_handler) ("matrix dimension mismatch");
1097 
1098  if (nr == 0 || nc == 0 || b_nc == 0)
1099  ret = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
1100  else
1101  {
1102  OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
1103 
1104  for (volatile octave_idx_type j = 0, idx = 0;
1105  j < b_nc;
1106  j++, idx += b_nr)
1107  {
1108  octave_quit ();
1109 
1110  volatile octave_idx_type nm = (nr < nc ? nr : nc);
1111 
1112  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1113  CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + idx,
1114  reinterpret_cast<cs_complex_t *> (buf),
1115  b_nr);
1116  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1117 
1118  for (volatile octave_idx_type i = 0; i < nm; i++)
1119  {
1120  octave_quit ();
1121 
1122  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1123  CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i],
1124  reinterpret_cast<cs_complex_t *> (buf));
1125  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1126  }
1127 
1128  for (octave_idx_type i = 0; i < b_nr; i++)
1129  vec[i+idx] = buf[i];
1130  }
1131  }
1132 
1133  return ret;
1134 
1135 #else
1136 
1137  octave_unused_parameter (b);
1138 
1139  return ComplexMatrix ();
1140 
1141 #endif
1142  }
1143 
1144  template <>
1145  ComplexMatrix
1146  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::Q (void) const
1147  {
1148 #if defined (HAVE_CXSPARSE)
1149  octave_idx_type nc = N->L->n;
1150  octave_idx_type nr = nrows;
1151  ComplexMatrix ret (nr, nr);
1152  Complex *vec = ret.fortran_vec ();
1153 
1154  if (nr < 0 || nc < 0)
1155  (*current_liboctave_error_handler) ("matrix dimension mismatch");
1156 
1157  if (nr == 0 || nc == 0)
1158  ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0));
1159  else
1160  {
1161  OCTAVE_LOCAL_BUFFER (cs_complex_t, bvec, nr);
1162 
1163  for (octave_idx_type i = 0; i < nr; i++)
1164  bvec[i] = cs_complex_t (0.0, 0.0);
1165 
1166  OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
1167 
1168  for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
1169  {
1170  octave_quit ();
1171 
1172  bvec[j] = cs_complex_t (1.0, 0.0);
1173 
1174  volatile octave_idx_type nm = (nr < nc ? nr : nc);
1175 
1176  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1177  CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec,
1178  reinterpret_cast<cs_complex_t *> (buf),
1179  nr);
1180  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1181 
1182  for (volatile octave_idx_type i = 0; i < nm; i++)
1183  {
1184  octave_quit ();
1185 
1186  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1187  CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i],
1188  reinterpret_cast<cs_complex_t *> (buf));
1189  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1190  }
1191 
1192  for (octave_idx_type i = 0; i < nr; i++)
1193  vec[i+idx] = buf[i];
1194 
1195  bvec[j] = cs_complex_t (0.0, 0.0);
1196  }
1197  }
1198 
1199  return ret.hermitian ();
1200 
1201 #else
1202 
1203  return ComplexMatrix ();
1204 
1205 #endif
1206  }
1207 
1208  template <>
1209  template <>
1210  SparseComplexMatrix
1211  sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix,
1212  SparseComplexMatrix>
1213  (const SparseComplexMatrix& b, octave_idx_type& info) const
1214  {
1215  info = -1;
1216 
1217 #if defined (HAVE_CXSPARSE)
1218 
1219  octave_idx_type nr = nrows;
1220  octave_idx_type nc = ncols;
1221 
1222  octave_idx_type b_nr = b.rows ();
1223  octave_idx_type b_nc = b.cols ();
1224 
1225  SparseComplexMatrix x (nc, b_nc, b.nnz ());
1226  x.xcidx (0) = 0;
1227 
1228  volatile octave_idx_type x_nz = b.nnz ();
1229  volatile octave_idx_type ii = 0;
1230 
1231  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
1232  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
1233  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
1234 
1235  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1236  {
1237  octave_quit ();
1238 
1239  for (octave_idx_type j = 0; j < b_nr; j++)
1240  {
1241  Complex c = b.xelem (j,i);
1242  Xx[j] = c.real ();
1243  Xz[j] = c.imag ();
1244  }
1245 
1246  for (octave_idx_type j = nr; j < S->m2; j++)
1247  buf[j] = 0.;
1248 
1249  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1250  CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr);
1251  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1252 
1253  for (volatile octave_idx_type j = 0; j < nc; j++)
1254  {
1255  octave_quit ();
1256 
1257  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1258  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1259  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1260  }
1261 
1262  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1263  CXSPARSE_DNAME (_usolve) (N->U, buf);
1264  CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc);
1265  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1266 
1267  for (octave_idx_type j = nr; j < S->m2; j++)
1268  buf[j] = 0.;
1269 
1270  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1271  CXSPARSE_DNAME (_ipvec) (S->pinv, Xz, buf, nr);
1272  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1273 
1274  for (volatile octave_idx_type j = 0; j < nc; j++)
1275  {
1276  octave_quit ();
1277 
1278  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1279  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1280  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1281  }
1282 
1283  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1284  CXSPARSE_DNAME (_usolve) (N->U, buf);
1285  CXSPARSE_DNAME (_ipvec) (S->q, buf, Xz, nc);
1286  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1287 
1288  for (octave_idx_type j = 0; j < nc; j++)
1289  {
1290  Complex tmp = Complex (Xx[j], Xz[j]);
1291 
1292  if (tmp != 0.0)
1293  {
1294  if (ii == x_nz)
1295  {
1296  // Resize the sparse matrix
1297  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1298  sz = (sz > 10 ? sz : 10) + x_nz;
1299  x.change_capacity (sz);
1300  x_nz = sz;
1301  }
1302 
1303  x.xdata (ii) = tmp;
1304  x.xridx (ii++) = j;
1305  }
1306  }
1307 
1308  x.xcidx (i+1) = ii;
1309  }
1310 
1311  info = 0;
1312 
1313  return x;
1314 
1315 #else
1316 
1317  octave_unused_parameter (b);
1318 
1319  return SparseComplexMatrix ();
1320 
1321 #endif
1322  }
1323 
1324  template <>
1325  template <>
1326  SparseComplexMatrix
1327  sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix,
1328  SparseComplexMatrix>
1329  (const SparseComplexMatrix& b, octave_idx_type& info) const
1330  {
1331  info = -1;
1332 
1333 #if defined (HAVE_CXSPARSE)
1334 
1335  // These are swapped because the original matrix was transposed in
1336  // sparse_qr<SparseMatrix>::solve.
1337 
1338  octave_idx_type nr = ncols;
1339  octave_idx_type nc = nrows;
1340 
1341  octave_idx_type b_nr = b.rows ();
1342  octave_idx_type b_nc = b.cols ();
1343 
1344  SparseComplexMatrix x (nc, b_nc, b.nnz ());
1345  x.xcidx (0) = 0;
1346 
1347  volatile octave_idx_type x_nz = b.nnz ();
1348  volatile octave_idx_type ii = 0;
1349  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1350 
1351  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
1352  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
1353  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
1354 
1355  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1356  {
1357  octave_quit ();
1358 
1359  for (octave_idx_type j = 0; j < b_nr; j++)
1360  {
1361  Complex c = b.xelem (j,i);
1362  Xx[j] = c.real ();
1363  Xz[j] = c.imag ();
1364  }
1365 
1366  for (octave_idx_type j = nr; j < nbuf; j++)
1367  buf[j] = 0.;
1368 
1369  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1370  CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr);
1371  CXSPARSE_DNAME (_utsolve) (N->U, buf);
1372  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1373 
1374  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
1375  {
1376  octave_quit ();
1377 
1378  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1379  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1380  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1381  }
1382 
1383  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1384  CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc);
1385  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1386 
1387  for (octave_idx_type j = nr; j < nbuf; j++)
1388  buf[j] = 0.;
1389 
1390  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1391  CXSPARSE_DNAME (_pvec) (S->q, Xz, buf, nr);
1392  CXSPARSE_DNAME (_utsolve) (N->U, buf);
1393  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1394 
1395  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
1396  {
1397  octave_quit ();
1398 
1399  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1400  CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1401  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1402  }
1403 
1404  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1405  CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xz, nc);
1406  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1407 
1408  for (octave_idx_type j = 0; j < nc; j++)
1409  {
1410  Complex tmp = Complex (Xx[j], Xz[j]);
1411 
1412  if (tmp != 0.0)
1413  {
1414  if (ii == x_nz)
1415  {
1416  // Resize the sparse matrix
1417  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1418  sz = (sz > 10 ? sz : 10) + x_nz;
1419  x.change_capacity (sz);
1420  x_nz = sz;
1421  }
1422 
1423  x.xdata (ii) = tmp;
1424  x.xridx (ii++) = j;
1425  }
1426  }
1427 
1428  x.xcidx (i+1) = ii;
1429  }
1430 
1431  info = 0;
1432 
1433  x.maybe_compress ();
1434 
1435  return x;
1436 
1437 #else
1438 
1439  octave_unused_parameter (b);
1440 
1441  return SparseComplexMatrix ();
1442 
1443 #endif
1444  }
1445 
1446  template <>
1447  template <>
1448  ComplexMatrix
1449  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<double>,
1450  ComplexMatrix>
1451  (const MArray<double>& b, octave_idx_type& info) const
1452  {
1453  info = -1;
1454 
1455 #if defined (HAVE_CXSPARSE)
1456 
1457  octave_idx_type nr = nrows;
1458  octave_idx_type nc = ncols;
1459 
1460  octave_idx_type b_nc = b.cols ();
1461  octave_idx_type b_nr = b.rows ();
1462 
1463  ComplexMatrix x (nc, b_nc);
1464  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
1465 
1466  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
1467  OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
1468 
1469  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1470  {
1471  octave_quit ();
1472 
1473  for (octave_idx_type j = 0; j < b_nr; j++)
1474  Xx[j] = b.xelem (j,i);
1475 
1476  for (octave_idx_type j = nr; j < S->m2; j++)
1477  buf[j] = cs_complex_t (0.0, 0.0);
1478 
1479  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1480  CXSPARSE_ZNAME (_ipvec) (S->pinv,
1481  reinterpret_cast<cs_complex_t *>(Xx),
1482  buf, nr);
1483  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1484 
1485  for (volatile octave_idx_type j = 0; j < nc; j++)
1486  {
1487  octave_quit ();
1488 
1489  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1490  CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
1491  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1492  }
1493 
1494  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1495  CXSPARSE_ZNAME (_usolve) (N->U, buf);
1496  CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc);
1497  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1498  }
1499 
1500  info = 0;
1501 
1502  return x;
1503 
1504 #else
1505 
1506  octave_unused_parameter (b);
1507 
1508  return ComplexMatrix ();
1509 
1510 #endif
1511  }
1512 
1513  template <>
1514  template <>
1515  ComplexMatrix
1516  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<double>,
1517  ComplexMatrix>
1518  (const MArray<double>& b, octave_idx_type& info) const
1519  {
1520  info = -1;
1521 
1522 #if defined (HAVE_CXSPARSE)
1523 
1524  // These are swapped because the original matrix was transposed in
1525  // sparse_qr<SparseComplexMatrix>::solve.
1526 
1527  octave_idx_type nr = ncols;
1528  octave_idx_type nc = nrows;
1529 
1530  octave_idx_type b_nc = b.cols ();
1531  octave_idx_type b_nr = b.rows ();
1532 
1533  ComplexMatrix x (nc, b_nc);
1534  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
1535 
1536  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1537 
1538  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
1539  OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
1540  OCTAVE_LOCAL_BUFFER (double, B, nr);
1541 
1542  for (octave_idx_type i = 0; i < nr; i++)
1543  B[i] = N->B[i];
1544 
1545  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1546  {
1547  octave_quit ();
1548 
1549  for (octave_idx_type j = 0; j < b_nr; j++)
1550  Xx[j] = b.xelem (j,i);
1551 
1552  for (octave_idx_type j = nr; j < nbuf; j++)
1553  buf[j] = cs_complex_t (0.0, 0.0);
1554 
1555  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1556  CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx),
1557  buf, nr);
1558  CXSPARSE_ZNAME (_utsolve) (N->U, buf);
1559  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1560 
1561  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
1562  {
1563  octave_quit ();
1564 
1565  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1566  CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
1567  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1568  }
1569 
1570  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1571  CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc);
1572  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1573  }
1574 
1575  info = 0;
1576 
1577  return x;
1578 
1579 #else
1580 
1581  octave_unused_parameter (b);
1582 
1583  return ComplexMatrix ();
1584 
1585 #endif
1586  }
1587 
1588  template <>
1589  template <>
1590  SparseComplexMatrix
1591  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseMatrix,
1592  SparseComplexMatrix>
1593  (const SparseMatrix& b, octave_idx_type& info) const
1594  {
1595  info = -1;
1596 
1597 #if defined (HAVE_CXSPARSE)
1598 
1599  octave_idx_type nr = nrows;
1600  octave_idx_type nc = ncols;
1601 
1602  octave_idx_type b_nc = b.cols ();
1603  octave_idx_type b_nr = b.rows ();
1604 
1605  SparseComplexMatrix x (nc, b_nc, b.nnz ());
1606  x.xcidx (0) = 0;
1607 
1608  volatile octave_idx_type x_nz = b.nnz ();
1609  volatile octave_idx_type ii = 0;
1610 
1611  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
1612  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
1613 
1614  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1615  {
1616  octave_quit ();
1617 
1618  for (octave_idx_type j = 0; j < b_nr; j++)
1619  Xx[j] = b.xelem (j,i);
1620 
1621  for (octave_idx_type j = nr; j < S->m2; j++)
1622  buf[j] = cs_complex_t (0.0, 0.0);
1623 
1624  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1625  CXSPARSE_ZNAME (_ipvec) (S->pinv,
1626  reinterpret_cast<cs_complex_t *>(Xx),
1627  buf, nr);
1628  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1629 
1630  for (volatile octave_idx_type j = 0; j < nc; j++)
1631  {
1632  octave_quit ();
1633 
1634  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1635  CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
1636  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1637  }
1638 
1639  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1640  CXSPARSE_ZNAME (_usolve) (N->U, buf);
1641  CXSPARSE_ZNAME (_ipvec) (S->q, buf,
1642  reinterpret_cast<cs_complex_t *> (Xx),
1643  nc);
1644  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1645 
1646  for (octave_idx_type j = 0; j < nc; j++)
1647  {
1648  Complex tmp = Xx[j];
1649 
1650  if (tmp != 0.0)
1651  {
1652  if (ii == x_nz)
1653  {
1654  // Resize the sparse matrix
1655  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1656  sz = (sz > 10 ? sz : 10) + x_nz;
1657  x.change_capacity (sz);
1658  x_nz = sz;
1659  }
1660 
1661  x.xdata (ii) = tmp;
1662  x.xridx (ii++) = j;
1663  }
1664  }
1665 
1666  x.xcidx (i+1) = ii;
1667  }
1668 
1669  info = 0;
1670 
1671  x.maybe_compress ();
1672 
1673  return x;
1674 
1675 #else
1676 
1677  octave_unused_parameter (b);
1678 
1679  return SparseComplexMatrix ();
1680 
1681 #endif
1682  }
1683 
1684  template <>
1685  template <>
1686  SparseComplexMatrix
1687  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseMatrix,
1688  SparseComplexMatrix>
1689  (const SparseMatrix& b, octave_idx_type& info) const
1690  {
1691  info = -1;
1692 
1693 #if defined (HAVE_CXSPARSE)
1694 
1695  // These are swapped because the original matrix was transposed in
1696  // sparse_qr<SparseComplexMatrix>::solve.
1697 
1698  octave_idx_type nr = ncols;
1699  octave_idx_type nc = nrows;
1700 
1701  octave_idx_type b_nc = b.cols ();
1702  octave_idx_type b_nr = b.rows ();
1703 
1704  SparseComplexMatrix x (nc, b_nc, b.nnz ());
1705  x.xcidx (0) = 0;
1706 
1707  volatile octave_idx_type x_nz = b.nnz ();
1708  volatile octave_idx_type ii = 0;
1709  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1710 
1711  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
1712  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
1713  OCTAVE_LOCAL_BUFFER (double, B, nr);
1714 
1715  for (octave_idx_type i = 0; i < nr; i++)
1716  B[i] = N->B[i];
1717 
1718  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1719  {
1720  octave_quit ();
1721 
1722  for (octave_idx_type j = 0; j < b_nr; j++)
1723  Xx[j] = b.xelem (j,i);
1724 
1725  for (octave_idx_type j = nr; j < nbuf; j++)
1726  buf[j] = cs_complex_t (0.0, 0.0);
1727 
1728  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1729  CXSPARSE_ZNAME (_pvec) (S->q,
1730  reinterpret_cast<cs_complex_t *>(Xx),
1731  buf, nr);
1732  CXSPARSE_ZNAME (_utsolve) (N->U, buf);
1733  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1734 
1735  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
1736  {
1737  octave_quit ();
1738 
1739  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1740  CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
1741  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1742  }
1743 
1744  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1745  CXSPARSE_ZNAME (_pvec) (S->pinv, buf,
1746  reinterpret_cast<cs_complex_t *> (Xx),
1747  nc);
1748  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1749 
1750  for (octave_idx_type j = 0; j < nc; j++)
1751  {
1752  Complex tmp = Xx[j];
1753 
1754  if (tmp != 0.0)
1755  {
1756  if (ii == x_nz)
1757  {
1758  // Resize the sparse matrix
1759  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1760  sz = (sz > 10 ? sz : 10) + x_nz;
1761  x.change_capacity (sz);
1762  x_nz = sz;
1763  }
1764 
1765  x.xdata (ii) = tmp;
1766  x.xridx (ii++) = j;
1767  }
1768  }
1769 
1770  x.xcidx (i+1) = ii;
1771  }
1772 
1773  info = 0;
1774 
1775  x.maybe_compress ();
1776 
1777  return x;
1778 
1779 #else
1780 
1781  octave_unused_parameter (b);
1782 
1783  return SparseComplexMatrix ();
1784 
1785 #endif
1786  }
1787 
1788  template <>
1789  template <>
1790  ComplexMatrix
1791  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>,
1792  ComplexMatrix>
1793  (const MArray<Complex>& b, octave_idx_type& info) const
1794  {
1795  info = -1;
1796 
1797 #if defined (HAVE_CXSPARSE)
1798 
1799  octave_idx_type nr = nrows;
1800  octave_idx_type nc = ncols;
1801 
1802  octave_idx_type b_nc = b.cols ();
1803  octave_idx_type b_nr = b.rows ();
1804 
1805  const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>
1806  (b.fortran_vec ());
1807 
1808  ComplexMatrix x (nc, b_nc);
1809  cs_complex_t *vec = reinterpret_cast<cs_complex_t *>
1810  (x.fortran_vec ());
1811 
1812  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
1813 
1814  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
1815  i++, idx+=nc, bidx+=b_nr)
1816  {
1817  octave_quit ();
1818 
1819  for (octave_idx_type j = nr; j < S->m2; j++)
1820  buf[j] = cs_complex_t (0.0, 0.0);
1821 
1822  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1823  CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + bidx, buf, nr);
1824  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1825 
1826  for (volatile octave_idx_type j = 0; j < nc; j++)
1827  {
1828  octave_quit ();
1829 
1830  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1831  CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
1832  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1833  }
1834 
1835  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1836  CXSPARSE_ZNAME (_usolve) (N->U, buf);
1837  CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc);
1838  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1839  }
1840 
1841  info = 0;
1842 
1843  return x;
1844 
1845 #else
1846 
1847  octave_unused_parameter (b);
1848 
1849  return ComplexMatrix ();
1850 
1851 #endif
1852  }
1853 
1854  template <>
1855  template <>
1856  ComplexMatrix
1857  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>,
1858  ComplexMatrix>
1859  (const MArray<Complex>& b, octave_idx_type& info) const
1860  {
1861  info = -1;
1862 
1863 #if defined (HAVE_CXSPARSE)
1864 
1865  // These are swapped because the original matrix was transposed in
1866  // sparse_qr<SparseComplexMatrix>::solve.
1867 
1868  octave_idx_type nr = ncols;
1869  octave_idx_type nc = nrows;
1870 
1871  octave_idx_type b_nc = b.cols ();
1872  octave_idx_type b_nr = b.rows ();
1873 
1874  const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>
1875  (b.fortran_vec ());
1876 
1877  ComplexMatrix x (nc, b_nc);
1878  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
1879 
1880  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1881 
1882  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
1883  OCTAVE_LOCAL_BUFFER (double, B, nr);
1884 
1885  for (octave_idx_type i = 0; i < nr; i++)
1886  B[i] = N->B[i];
1887 
1888  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
1889  i++, idx+=nc, bidx+=b_nr)
1890  {
1891  octave_quit ();
1892 
1893  for (octave_idx_type j = nr; j < nbuf; j++)
1894  buf[j] = cs_complex_t (0.0, 0.0);
1895 
1896  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1897  CXSPARSE_ZNAME (_pvec) (S->q, bvec + bidx, buf, nr);
1898  CXSPARSE_ZNAME (_utsolve) (N->U, buf);
1899  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1900 
1901  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
1902  {
1903  octave_quit ();
1904 
1905  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1906  CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
1907  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1908  }
1909 
1910  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1911  CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc);
1912  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1913  }
1914 
1915  info = 0;
1916 
1917  return x;
1918 
1919 #else
1920 
1921  octave_unused_parameter (b);
1922 
1923  return ComplexMatrix ();
1924 
1925 #endif
1926  }
1927 
1928  template <>
1929  template <>
1930  SparseComplexMatrix
1931  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix, SparseComplexMatrix>
1932  (const SparseComplexMatrix& b, octave_idx_type& info) const
1933  {
1934  info = -1;
1935 
1936 #if defined (HAVE_CXSPARSE)
1937 
1938  octave_idx_type nr = nrows;
1939  octave_idx_type nc = ncols;
1940 
1941  octave_idx_type b_nc = b.cols ();
1942  octave_idx_type b_nr = b.rows ();
1943 
1944  SparseComplexMatrix x (nc, b_nc, b.nnz ());
1945  x.xcidx (0) = 0;
1946 
1947  volatile octave_idx_type x_nz = b.nnz ();
1948  volatile octave_idx_type ii = 0;
1949 
1950  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
1951  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
1952 
1953  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1954  {
1955  octave_quit ();
1956 
1957  for (octave_idx_type j = 0; j < b_nr; j++)
1958  Xx[j] = b.xelem (j,i);
1959 
1960  for (octave_idx_type j = nr; j < S->m2; j++)
1961  buf[j] = cs_complex_t (0.0, 0.0);
1962 
1963  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1964  CXSPARSE_ZNAME (_ipvec) (S->pinv,
1965  reinterpret_cast<cs_complex_t *>(Xx),
1966  buf, nr);
1967  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1968 
1969  for (volatile octave_idx_type j = 0; j < nc; j++)
1970  {
1971  octave_quit ();
1972 
1973  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1974  CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
1975  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1976  }
1977 
1978  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1979  CXSPARSE_ZNAME (_usolve) (N->U, buf);
1980  CXSPARSE_ZNAME (_ipvec) (S->q, buf,
1981  reinterpret_cast<cs_complex_t *> (Xx),
1982  nc);
1983  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
1984 
1985  for (octave_idx_type j = 0; j < nc; j++)
1986  {
1987  Complex tmp = Xx[j];
1988 
1989  if (tmp != 0.0)
1990  {
1991  if (ii == x_nz)
1992  {
1993  // Resize the sparse matrix
1994  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1995  sz = (sz > 10 ? sz : 10) + x_nz;
1996  x.change_capacity (sz);
1997  x_nz = sz;
1998  }
1999 
2000  x.xdata (ii) = tmp;
2001  x.xridx (ii++) = j;
2002  }
2003  }
2004 
2005  x.xcidx (i+1) = ii;
2006  }
2007 
2008  info = 0;
2009 
2010  x.maybe_compress ();
2011 
2012  return x;
2013 
2014 #else
2015 
2016  octave_unused_parameter (b);
2017 
2018  return SparseComplexMatrix ();
2019 
2020 #endif
2021  }
2022 
2023  template <>
2024  template <>
2025  SparseComplexMatrix
2026  sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix, SparseComplexMatrix>
2027  (const SparseComplexMatrix& b, octave_idx_type& info) const
2028  {
2029  info = -1;
2030 
2031 #if defined (HAVE_CXSPARSE)
2032 
2033  // These are swapped because the original matrix was transposed in
2034  // sparse_qr<SparseComplexMatrix>::solve.
2035 
2036  octave_idx_type nr = ncols;
2037  octave_idx_type nc = nrows;
2038 
2039  octave_idx_type b_nc = b.cols ();
2040  octave_idx_type b_nr = b.rows ();
2041 
2042  SparseComplexMatrix x (nc, b_nc, b.nnz ());
2043  x.xcidx (0) = 0;
2044 
2045  volatile octave_idx_type x_nz = b.nnz ();
2046  volatile octave_idx_type ii = 0;
2047  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
2048 
2049  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
2050  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
2051  OCTAVE_LOCAL_BUFFER (double, B, nr);
2052 
2053  for (octave_idx_type i = 0; i < nr; i++)
2054  B[i] = N->B[i];
2055 
2056  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
2057  {
2058  octave_quit ();
2059 
2060  for (octave_idx_type j = 0; j < b_nr; j++)
2061  Xx[j] = b.xelem (j,i);
2062 
2063  for (octave_idx_type j = nr; j < nbuf; j++)
2064  buf[j] = cs_complex_t (0.0, 0.0);
2065 
2066  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2067  CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx),
2068  buf, nr);
2069  CXSPARSE_ZNAME (_utsolve) (N->U, buf);
2070  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2071 
2072  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
2073  {
2074  octave_quit ();
2075 
2076  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2077  CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
2078  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2079  }
2080 
2081  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2082  CXSPARSE_ZNAME (_pvec) (S->pinv, buf,
2083  reinterpret_cast<cs_complex_t *>(Xx), nc);
2084  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
2085 
2086  for (octave_idx_type j = 0; j < nc; j++)
2087  {
2088  Complex tmp = Xx[j];
2089 
2090  if (tmp != 0.0)
2091  {
2092  if (ii == x_nz)
2093  {
2094  // Resize the sparse matrix
2095  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
2096  sz = (sz > 10 ? sz : 10) + x_nz;
2097  x.change_capacity (sz);
2098  x_nz = sz;
2099  }
2100 
2101  x.xdata (ii) = tmp;
2102  x.xridx (ii++) = j;
2103  }
2104  }
2105 
2106  x.xcidx (i+1) = ii;
2107  }
2108 
2109  info = 0;
2110 
2111  x.maybe_compress ();
2112 
2113  return x;
2114 
2115 #else
2116 
2117  octave_unused_parameter (b);
2118 
2119  return SparseComplexMatrix ();
2120 
2121 #endif
2122  }
2123 
2124  template <typename SPARSE_T>
2125  sparse_qr<SPARSE_T>::sparse_qr (void)
2126  : rep (new sparse_qr_rep (SPARSE_T (), 0))
2127  { }
2128 
2129  template <typename SPARSE_T>
2130  sparse_qr<SPARSE_T>::sparse_qr (const SPARSE_T& a, int order)
2131  : rep (new sparse_qr_rep (a, order))
2132  { }
2133 
2134  template <typename SPARSE_T>
2135  sparse_qr<SPARSE_T>::sparse_qr (const sparse_qr<SPARSE_T>& a)
2136  : rep (a.rep)
2137  {
2138  rep->count++;
2139  }
2140 
2141  template <typename SPARSE_T>
2142  sparse_qr<SPARSE_T>::~sparse_qr (void)
2143  {
2144  if (--rep->count == 0)
2145  delete rep;
2146  }
2147 
2148  template <typename SPARSE_T>
2149  sparse_qr<SPARSE_T>&
2150  sparse_qr<SPARSE_T>::operator = (const sparse_qr<SPARSE_T>& a)
2151  {
2152  if (this != &a)
2153  {
2154  if (--rep->count == 0)
2155  delete rep;
2156 
2157  rep = a.rep;
2158  rep->count++;
2159  }
2160 
2161  return *this;
2162  }
2163 
2164  template <typename SPARSE_T>
2165  bool
2166  sparse_qr<SPARSE_T>::ok (void) const
2167  {
2168  return rep->ok ();
2169  }
2170 
2171  template <typename SPARSE_T>
2172  SPARSE_T
2173  sparse_qr<SPARSE_T>::V (void) const
2174  {
2175  return rep->V ();
2176  }
2177 
2178  template <typename SPARSE_T>
2179  ColumnVector
2180  sparse_qr<SPARSE_T>::Pinv (void) const
2181  {
2182  return rep->P ();
2183  }
2184 
2185  template <typename SPARSE_T>
2186  ColumnVector
2187  sparse_qr<SPARSE_T>::P (void) const
2188  {
2189  return rep->P ();
2190  }
2191 
2192  template <typename SPARSE_T>
2193  SPARSE_T
2194  sparse_qr<SPARSE_T>::R (bool econ) const
2195  {
2196  return rep->R (econ);
2197  }
2198 
2199  template <typename SPARSE_T>
2200  typename SPARSE_T::dense_matrix_type
2201  sparse_qr<SPARSE_T>::C (const typename SPARSE_T::dense_matrix_type& b) const
2202  {
2203  return rep->C (b);
2204  }
2205 
2206  template <typename SPARSE_T>
2207  typename SPARSE_T::dense_matrix_type
2208  sparse_qr<SPARSE_T>::Q (void) const
2209  {
2210  return rep->Q ();
2211  }
2212 
2213  // FIXME: Why is the "order" of the QR calculation as used in the
2214  // CXSparse function sqr 3 for real matrices and 2 for complex? These
2215  // values seem to be required but there was no explanation in David
2216  // Bateman's original code.
2217 
2218  template <typename SPARSE_T>
2219  class
2220  cxsparse_defaults
2221  {
2222  public:
2223  enum { order = -1 };
2224  };
2225 
2226  template <>
2227  class
2228  cxsparse_defaults<SparseMatrix>
2229  {
2230  public:
2231  enum { order = 3 };
2232  };
2233 
2234  template <>
2235  class
2236  cxsparse_defaults<SparseComplexMatrix>
2237  {
2238  public:
2239  enum { order = 2 };
2240  };
2241 
2242  template <typename SPARSE_T>
2243  template <typename RHS_T, typename RET_T>
2244  RET_T
2245  sparse_qr<SPARSE_T>::solve (const SPARSE_T& a, const RHS_T& b,
2246  octave_idx_type& info)
2247  {
2248  info = -1;
2249 
2250  octave_idx_type nr = a.rows ();
2251  octave_idx_type nc = a.cols ();
2252 
2253  octave_idx_type b_nc = b.cols ();
2254  octave_idx_type b_nr = b.rows ();
2255 
2256  int order = cxsparse_defaults<SPARSE_T>::order;
2257 
2258  if (nr < 0 || nc < 0 || nr != b_nr)
2259  (*current_liboctave_error_handler)
2260  ("matrix dimension mismatch in solution of minimum norm problem");
2261 
2262  if (nr == 0 || nc == 0 || b_nc == 0)
2263  {
2264  info = 0;
2265 
2266  return RET_T (nc, b_nc, 0.0);
2267  }
2268  else if (nr >= nc)
2269  {
2270  sparse_qr<SPARSE_T> q (a, order);
2271 
2272  return q.ok () ? q.tall_solve<RHS_T, RET_T> (b, info) : RET_T ();
2273  }
2274  else
2275  {
2276  sparse_qr<SPARSE_T> q (a.hermitian (), order);
2277 
2278  return q.ok () ? q.wide_solve<RHS_T, RET_T> (b, info) : RET_T ();
2279  }
2280  }
2281 
2282  template <typename SPARSE_T>
2283  template <typename RHS_T, typename RET_T>
2284  RET_T
2285  sparse_qr<SPARSE_T>::tall_solve (const RHS_T& b, octave_idx_type& info) const
2286  {
2287  return rep->template tall_solve<RHS_T, RET_T> (b, info);
2288  }
2289 
2290  template <typename SPARSE_T>
2291  template <typename RHS_T, typename RET_T>
2292  RET_T
2293  sparse_qr<SPARSE_T>::wide_solve (const RHS_T& b, octave_idx_type& info) const
2294  {
2295  return rep->template wide_solve<RHS_T, RET_T> (b, info);
2296  }
2297 
2298  Matrix
2299  qrsolve (const SparseMatrix& a, const MArray<double>& b,
2300  octave_idx_type& info)
2301  {
2302  return sparse_qr<SparseMatrix>::solve<MArray<double>, Matrix> (a, b,
2303  info);
2304  }
2305 
2306  SparseMatrix
2307  qrsolve (const SparseMatrix& a, const SparseMatrix& b,
2308  octave_idx_type& info)
2309  {
2310  return sparse_qr<SparseMatrix>::solve<SparseMatrix, SparseMatrix> (a, b,
2311  info);
2312  }
2313 
2314  ComplexMatrix
2315  qrsolve (const SparseMatrix& a, const MArray<Complex>& b,
2316  octave_idx_type& info)
2317  {
2318  return sparse_qr<SparseMatrix>::solve<MArray<Complex>,
2319  ComplexMatrix> (a, b, info);
2320  }
2321 
2322  SparseComplexMatrix
2323  qrsolve (const SparseMatrix& a, const SparseComplexMatrix& b,
2324  octave_idx_type& info)
2325  {
2326  return sparse_qr<SparseMatrix>::solve<SparseComplexMatrix,
2327  SparseComplexMatrix> (a, b, info);
2328  }
2329 
2330  ComplexMatrix
2331  qrsolve (const SparseComplexMatrix& a, const MArray<double>& b,
2332  octave_idx_type& info)
2333  {
2334  return sparse_qr<SparseComplexMatrix>::solve<MArray<double>,
2335  ComplexMatrix> (a, b, info);
2336  }
2337 
2338  SparseComplexMatrix
2339  qrsolve (const SparseComplexMatrix& a, const SparseMatrix& b,
2340  octave_idx_type& info)
2341  {
2342  return sparse_qr<SparseComplexMatrix>::solve<SparseMatrix,
2343  SparseComplexMatrix>
2344  (a, b, info);
2345  }
2346 
2347  ComplexMatrix
2348  qrsolve (const SparseComplexMatrix& a, const MArray<Complex>& b,
2349  octave_idx_type& info)
2350  {
2351  return sparse_qr<SparseComplexMatrix>::solve<MArray<Complex>,
2352  ComplexMatrix> (a, b, info);
2353  }
2354 
2355  SparseComplexMatrix
2356  qrsolve (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
2357  octave_idx_type& info)
2358  {
2359  return sparse_qr<SparseComplexMatrix>::solve<SparseComplexMatrix,
2360  SparseComplexMatrix>
2361  (a, b, info);
2362  }
2363 
2364  // Instantiations we need.
2365 
2366  template class sparse_qr<SparseMatrix>;
2367 
2368  template class sparse_qr<SparseComplexMatrix>;
2369  }
2370 }
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:469
Array::cols
octave_idx_type cols(void) const
Definition: Array.h:423
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:415
Array::fortran_vec
const T * fortran_vec(void) const
Size of the specified dimension.
Definition: Array.h:583
ComplexMatrix::hermitian
ComplexMatrix hermitian(void) const
Definition: CMatrix.h:169
Matrix
Definition: dMatrix.h:42
Matrix::transpose
Matrix transpose(void) const
Definition: dMatrix.h:135
Sparse::cols
octave_idx_type cols(void) const
Definition: Sparse.h:251
Sparse::xridx
octave_idx_type * xridx(void)
Definition: Sparse.h:485
Sparse::xelem
T & xelem(octave_idx_type n)
Definition: Sparse.h:293
Sparse::data
T * data(void)
Definition: Sparse.h:470
Sparse::cidx
octave_idx_type * cidx(void)
Definition: Sparse.h:492
Sparse::nnz
octave_idx_type nnz(void) const
Actual number of nonzero terms.
Definition: Sparse.h:238
Sparse::rows
octave_idx_type rows(void) const
Definition: Sparse.h:250
Sparse::xcidx
octave_idx_type * xcidx(void)
Definition: Sparse.h:498
Sparse::xdata
T * xdata(void)
Definition: Sparse.h:472
Sparse::ridx
octave_idx_type * ridx(void)
Definition: Sparse.h:479
octave::math::cxsparse_defaults< SparseComplexMatrix >
Definition: sparse-qr.cc:2237
octave::math::cxsparse_defaults< SparseMatrix >
Definition: sparse-qr.cc:2229
octave::math::cxsparse_defaults
Definition: sparse-qr.cc:2221
octave::math::cxsparse_types< SparseComplexMatrix >
Definition: sparse-qr.cc:70
octave::math::cxsparse_types< SparseComplexMatrix >::symbolic_type
cs_cis symbolic_type
Definition: sparse-qr.cc:73
octave::math::cxsparse_types< SparseComplexMatrix >::numeric_type
cs_cin numeric_type
Definition: sparse-qr.cc:74
octave::math::cxsparse_types< SparseMatrix >
Definition: sparse-qr.cc:56
octave::math::cxsparse_types< SparseMatrix >::symbolic_type
cs_dis symbolic_type
Definition: sparse-qr.cc:59
octave::math::cxsparse_types< SparseMatrix >::numeric_type
cs_din numeric_type
Definition: sparse-qr.cc:60
octave::math::cxsparse_types
Definition: sparse-qr.cc:50
octave::math::sparse_qr::sparse_qr_rep
Definition: sparse-qr.cc:83
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::sparse_qr::sparse_qr_rep::R
SPARSE_T R(bool econ) const
octave::math::sparse_qr::sparse_qr_rep::N
cxsparse_types< SPARSE_T >::numeric_type * N
Definition: sparse-qr.cc:125
octave::math::sparse_qr::sparse_qr_rep::count
refcount< octave_idx_type > count
Definition: sparse-qr.cc:119
octave::math::sparse_qr::sparse_qr_rep::sparse_qr_rep
sparse_qr_rep(const SPARSE_T &a, int order)
octave::math::sparse_qr::sparse_qr_rep::wide_solve
RET_T wide_solve(const RHS_T &b, octave_idx_type &info) const
octave::math::sparse_qr::sparse_qr_rep::P
ColumnVector P(void) const
Definition: sparse-qr.cc:158
octave::math::sparse_qr::sparse_qr_rep::sparse_qr_rep
sparse_qr_rep(const sparse_qr_rep &)=delete
octave::math::sparse_qr::sparse_qr_rep::Q
SPARSE_T::dense_matrix_type Q(void) const
octave::math::sparse_qr::sparse_qr_rep::~sparse_qr_rep
~sparse_qr_rep(void)
octave::math::sparse_qr::sparse_qr_rep::tall_solve
RET_T tall_solve(const RHS_T &b, octave_idx_type &info) const
octave::math::sparse_qr::sparse_qr_rep::Pinv
ColumnVector Pinv(void) const
Definition: sparse-qr.cc:138
octave::math::sparse_qr::sparse_qr_rep::S
cxsparse_types< SPARSE_T >::symbolic_type * S
Definition: sparse-qr.cc:124
octave::math::sparse_qr::sparse_qr_rep::ok
bool ok(void) const
Definition: sparse-qr.cc:96
octave::math::sparse_qr::sparse_qr_rep::nrows
octave_idx_type nrows
Definition: sparse-qr.cc:121
octave::math::sparse_qr::sparse_qr_rep::V
SPARSE_T V(void) const
octave::math::sparse_qr::sparse_qr_rep::ncols
octave_idx_type ncols
Definition: sparse-qr.cc:122
octave::math::sparse_qr
Definition: sparse-qr.h:52
octave::math::sparse_qr::tall_solve
RET_T tall_solve(const RHS_T &b, octave_idx_type &info) const
Definition: sparse-qr.cc:2285
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:2245
octave::math::sparse_qr::R
SPARSE_T R(bool econ=false) const
Definition: sparse-qr.cc:2194
octave::math::sparse_qr::Pinv
ColumnVector Pinv(void) const
Definition: sparse-qr.cc:2180
octave::math::sparse_qr::rep
sparse_qr_rep * rep
Definition: sparse-qr.h:88
octave::math::sparse_qr::wide_solve
RET_T wide_solve(const RHS_T &b, octave_idx_type &info) const
Definition: sparse-qr.cc:2293
octave::math::sparse_qr::~sparse_qr
~sparse_qr(void)
Definition: sparse-qr.cc:2142
octave::math::sparse_qr::sparse_qr
sparse_qr(void)
Definition: sparse-qr.cc:2125
octave::math::sparse_qr::operator=
sparse_qr & operator=(const sparse_qr &a)
Definition: sparse-qr.cc:2150
octave::math::sparse_qr::V
SPARSE_T V(void) const
Definition: sparse-qr.cc:2173
octave::math::sparse_qr::ok
bool ok(void) const
Definition: sparse-qr.cc:2166
octave::math::sparse_qr::P
ColumnVector P(void) const
Definition: sparse-qr.cc:2187
octave::math::sparse_qr::C
SPARSE_T::dense_matrix_type C(const typename SPARSE_T::dense_matrix_type &b) const
Definition: sparse-qr.cc:2201
octave::math::sparse_qr::Q
SPARSE_T::dense_matrix_type Q(void) const
Definition: sparse-qr.cc:2208
octave_idx_type
B
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
Definition: lo-lapack-proto.h:657
N
F77_RET_T const F77_INT & N
Definition: lo-lapack-proto.h:655
A
F77_RET_T const F77_INT F77_CMPLX * A
Definition: lo-lapack-proto.h:656
x
F77_RET_T const F77_DBLE * x
Definition: lo-slatec-proto.h:38
for
for(octave_idx_type i=0;i< n;i++) ac+
Definition: mx-inlines.cc:756
n
octave_idx_type n
Definition: mx-inlines.cc:753
octave::math::qrsolve
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
Definition: sparse-qr.cc:2299
octave::suitesparse_integer
int suitesparse_integer
Definition: oct-sparse.h:171
octave::to_suitesparse_intptr
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
Definition: oct-sparse.cc:51
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:33
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:44
CXSPARSE_DNAME
#define CXSPARSE_DNAME(name)
Definition: oct-sparse.h:147
CXSPARSE_ZNAME
#define CXSPARSE_ZNAME(name)
Definition: oct-sparse.h:148
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
#define BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
Definition: quit.h:280
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
#define END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
Definition: quit.h:284