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