GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SparseQR.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2005-2013 David Bateman
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 #include <vector>
27 
28 #include "lo-error.h"
29 #include "SparseQR.h"
30 #include "oct-locbuf.h"
31 
33  : count (1), nrows (0)
34 #ifdef HAVE_CXSPARSE
35  , S (0), N (0)
36 #endif
37 {
38 #ifdef HAVE_CXSPARSE
39  CXSPARSE_DNAME () A;
40  A.nzmax = a.nnz ();
41  A.m = a.rows ();
42  A.n = a.cols ();
43  nrows = A.m;
44  // Cast away const on A, with full knowledge that CSparse won't touch it
45  // Prevents the methods below making a copy of the data.
46  A.p = const_cast<octave_idx_type *>(a.cidx ());
47  A.i = const_cast<octave_idx_type *>(a.ridx ());
48  A.x = const_cast<double *>(a.data ());
49  A.nz = -1;
51 #if defined (CS_VER) && (CS_VER >= 2)
52  S = CXSPARSE_DNAME (_sqr) (order, &A, 1);
53 #else
54  S = CXSPARSE_DNAME (_sqr) (&A, order - 1, 1);
55 #endif
56 
57  N = CXSPARSE_DNAME (_qr) (&A, S);
59  if (!N)
60  (*current_liboctave_error_handler)
61  ("SparseQR: sparse matrix QR factorization filled");
62  count = 1;
63 #else
64  (*current_liboctave_error_handler)
65  ("SparseQR: sparse matrix QR factorization not implemented");
66 #endif
67 }
68 
70 {
71 #ifdef HAVE_CXSPARSE
72  CXSPARSE_DNAME (_sfree) (S);
73  CXSPARSE_DNAME (_nfree) (N);
74 #endif
75 }
76 
79 {
80 #ifdef HAVE_CXSPARSE
81  // Drop zeros from V and sort
82  // FIXME: Is the double transpose to sort necessary?
84  CXSPARSE_DNAME (_dropzeros) (N->L);
85  CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->L, 1);
86  CXSPARSE_DNAME (_spfree) (N->L);
87  N->L = CXSPARSE_DNAME (_transpose) (D, 1);
88  CXSPARSE_DNAME (_spfree) (D);
90 
91  octave_idx_type nc = N->L->n;
92  octave_idx_type nz = N->L->nzmax;
93  SparseMatrix ret (N->L->m, nc, nz);
94  for (octave_idx_type j = 0; j < nc+1; j++)
95  ret.xcidx (j) = N->L->p[j];
96  for (octave_idx_type j = 0; j < nz; j++)
97  {
98  ret.xridx (j) = N->L->i[j];
99  ret.xdata (j) = N->L->x[j];
100  }
101  return ret;
102 #else
103  return SparseMatrix ();
104 #endif
105 }
106 
109 {
110 #ifdef HAVE_CXSPARSE
111  ColumnVector ret(N->L->m);
112  for (octave_idx_type i = 0; i < N->L->m; i++)
113 #if defined (CS_VER) && (CS_VER >= 2)
114  ret.xelem (i) = S->pinv[i];
115 #else
116  ret.xelem (i) = S->Pinv[i];
117 #endif
118  return ret;
119 #else
120  return ColumnVector ();
121 #endif
122 }
123 
126 {
127 #ifdef HAVE_CXSPARSE
128  ColumnVector ret(N->L->m);
129  for (octave_idx_type i = 0; i < N->L->m; i++)
130 #if defined (CS_VER) && (CS_VER >= 2)
131  ret.xelem (S->pinv[i]) = i;
132 #else
133  ret.xelem (S->Pinv[i]) = i;
134 #endif
135  return ret;
136 #else
137  return ColumnVector ();
138 #endif
139 }
140 
142 SparseQR::SparseQR_rep::R (const bool econ) const
143 {
144 #ifdef HAVE_CXSPARSE
145  // Drop zeros from R and sort
146  // FIXME: Is the double transpose to sort necessary?
148  CXSPARSE_DNAME (_dropzeros) (N->U);
149  CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->U, 1);
150  CXSPARSE_DNAME (_spfree) (N->U);
151  N->U = CXSPARSE_DNAME (_transpose) (D, 1);
152  CXSPARSE_DNAME (_spfree) (D);
154 
155  octave_idx_type nc = N->U->n;
156  octave_idx_type nz = N->U->nzmax;
157 
158  SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
159 
160  for (octave_idx_type j = 0; j < nc+1; j++)
161  ret.xcidx (j) = N->U->p[j];
162  for (octave_idx_type j = 0; j < nz; j++)
163  {
164  ret.xridx (j) = N->U->i[j];
165  ret.xdata (j) = N->U->x[j];
166  }
167  return ret;
168 #else
169  return SparseMatrix ();
170 #endif
171 }
172 
173 Matrix
175 {
176 #ifdef HAVE_CXSPARSE
177  octave_idx_type b_nr = b.rows ();
178  octave_idx_type b_nc = b.cols ();
179  octave_idx_type nc = N->L->n;
180  octave_idx_type nr = nrows;
181  const double *bvec = b.fortran_vec ();
182  Matrix ret (b_nr, b_nc);
183  double *vec = ret.fortran_vec ();
184  if (nr < 0 || nc < 0 || nr != b_nr)
185  (*current_liboctave_error_handler) ("matrix dimension mismatch");
186  else if (nr == 0 || nc == 0 || b_nc == 0)
187  ret = Matrix (nc, b_nc, 0.0);
188  else
189  {
190  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
191  for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
192  {
193  octave_quit ();
194  for (octave_idx_type i = nr; i < S->m2; i++)
195  buf[i] = 0.;
196  volatile octave_idx_type nm = (nr < nc ? nr : nc);
198 #if defined (CS_VER) && (CS_VER >= 2)
199  CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr);
200 #else
201  CXSPARSE_DNAME (_ipvec) (b_nr, S->Pinv, bvec + idx, buf);
202 #endif
204 
205  for (volatile octave_idx_type i = 0; i < nm; i++)
206  {
207  octave_quit ();
209  CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
211  }
212  for (octave_idx_type i = 0; i < b_nr; i++)
213  vec[i+idx] = buf[i];
214  }
215  }
216  return ret;
217 #else
218  return Matrix ();
219 #endif
220 }
221 
222 Matrix
224 {
225 #ifdef HAVE_CXSPARSE
226  octave_idx_type nc = N->L->n;
227  octave_idx_type nr = nrows;
228  Matrix ret (nr, nr);
229  double *vec = ret.fortran_vec ();
230  if (nr < 0 || nc < 0)
231  (*current_liboctave_error_handler) ("matrix dimension mismatch");
232  else if (nr == 0 || nc == 0)
233  ret = Matrix (nc, nr, 0.0);
234  else
235  {
236  OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1);
237  for (octave_idx_type i = 0; i < nr; i++)
238  bvec[i] = 0.;
239  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
240  for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
241  {
242  octave_quit ();
243  bvec[j] = 1.0;
244  for (octave_idx_type i = nr; i < S->m2; i++)
245  buf[i] = 0.;
246  volatile octave_idx_type nm = (nr < nc ? nr : nc);
248 #if defined (CS_VER) && (CS_VER >= 2)
249  CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr);
250 #else
251  CXSPARSE_DNAME (_ipvec) (nr, S->Pinv, bvec, buf);
252 #endif
254 
255  for (volatile octave_idx_type i = 0; i < nm; i++)
256  {
257  octave_quit ();
259  CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
261  }
262  for (octave_idx_type i = 0; i < nr; i++)
263  vec[i+idx] = buf[i];
264  bvec[j] = 0.0;
265  }
266  }
267  return ret.transpose ();
268 #else
269  return Matrix ();
270 #endif
271 }
272 
273 Matrix
274 qrsolve (const SparseMatrix&a, const Matrix &b, octave_idx_type& info)
275 {
276  info = -1;
277 #ifdef HAVE_CXSPARSE
278  octave_idx_type nr = a.rows ();
279  octave_idx_type nc = a.cols ();
280  octave_idx_type b_nc = b.cols ();
281  octave_idx_type b_nr = b.rows ();
282  const double *bvec = b.fortran_vec ();
283  Matrix x;
284 
285  if (nr < 0 || nc < 0 || nr != b_nr)
286  (*current_liboctave_error_handler)
287  ("matrix dimension mismatch in solution of minimum norm problem");
288  else if (nr == 0 || nc == 0 || b_nc == 0)
289  x = Matrix (nc, b_nc, 0.0);
290  else if (nr >= nc)
291  {
292  SparseQR q (a, 3);
293  if (! q.ok ())
294  return Matrix ();
295  x.resize (nc, b_nc);
296  double *vec = x.fortran_vec ();
297  OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
298  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
299  i++, idx+=nc, bidx+=b_nr)
300  {
301  octave_quit ();
302  for (octave_idx_type j = nr; j < q.S ()->m2; j++)
303  buf[j] = 0.;
305 #if defined (CS_VER) && (CS_VER >= 2)
306  CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, bvec + bidx, buf, nr);
307 #else
308  CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, bvec + bidx, buf);
309 #endif
311  for (volatile octave_idx_type j = 0; j < nc; j++)
312  {
313  octave_quit ();
315  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
317  }
319  CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
320 #if defined (CS_VER) && (CS_VER >= 2)
321  CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc);
322 #else
323  CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, vec + idx);
324 #endif
326  }
327  info = 0;
328  }
329  else
330  {
331  SparseMatrix at = a.hermitian ();
332  SparseQR q (at, 3);
333  if (! q.ok ())
334  return Matrix ();
335  x.resize (nc, b_nc);
336  double *vec = x.fortran_vec ();
337  volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
338  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
339  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
340  i++, idx+=nc, bidx+=b_nr)
341  {
342  octave_quit ();
343  for (octave_idx_type j = nr; j < nbuf; j++)
344  buf[j] = 0.;
346 #if defined (CS_VER) && (CS_VER >= 2)
347  CXSPARSE_DNAME (_pvec) (q.S ()->q, bvec + bidx, buf, nr);
348 #else
349  CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, bvec + bidx, buf);
350 #endif
351  CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
353  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
354  {
355  octave_quit ();
357  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
359  }
361 #if defined (CS_VER) && (CS_VER >= 2)
362  CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc);
363 #else
364  CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, vec + idx);
365 #endif
367  }
368  info = 0;
369  }
370 
371  return x;
372 #else
373  return Matrix ();
374 #endif
375 }
376 
379 {
380  info = -1;
381 #ifdef HAVE_CXSPARSE
382  octave_idx_type nr = a.rows ();
383  octave_idx_type nc = a.cols ();
384  octave_idx_type b_nr = b.rows ();
385  octave_idx_type b_nc = b.cols ();
386  SparseMatrix x;
387  volatile octave_idx_type ii, x_nz;
388 
389  if (nr < 0 || nc < 0 || nr != b_nr)
390  (*current_liboctave_error_handler)
391  ("matrix dimension mismatch in solution of minimum norm problem");
392  else if (nr == 0 || nc == 0 || b_nc == 0)
393  x = SparseMatrix (nc, b_nc);
394  else if (nr >= nc)
395  {
396  SparseQR q (a, 3);
397  if (! q.ok ())
398  return SparseMatrix ();
399  x = SparseMatrix (nc, b_nc, b.nnz ());
400  x.xcidx (0) = 0;
401  x_nz = b.nnz ();
402  ii = 0;
403  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
404  OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
405  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
406  {
407  octave_quit ();
408  for (octave_idx_type j = 0; j < b_nr; j++)
409  Xx[j] = b.xelem (j,i);
410  for (octave_idx_type j = nr; j < q.S ()->m2; j++)
411  buf[j] = 0.;
413 #if defined (CS_VER) && (CS_VER >= 2)
414  CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
415 #else
416  CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xx, buf);
417 #endif
419  for (volatile octave_idx_type j = 0; j < nc; j++)
420  {
421  octave_quit ();
423  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
425  }
427  CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
428 #if defined (CS_VER) && (CS_VER >= 2)
429  CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
430 #else
431  CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xx);
432 #endif
434 
435  for (octave_idx_type j = 0; j < nc; j++)
436  {
437  double tmp = Xx[j];
438  if (tmp != 0.0)
439  {
440  if (ii == x_nz)
441  {
442  // Resize the sparse matrix
443  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
444  sz = (sz > 10 ? sz : 10) + x_nz;
445  x.change_capacity (sz);
446  x_nz = sz;
447  }
448  x.xdata (ii) = tmp;
449  x.xridx (ii++) = j;
450  }
451  }
452  x.xcidx (i+1) = ii;
453  }
454  info = 0;
455  }
456  else
457  {
458  SparseMatrix at = a.hermitian ();
459  SparseQR q (at, 3);
460  if (! q.ok ())
461  return SparseMatrix ();
462  x = SparseMatrix (nc, b_nc, b.nnz ());
463  x.xcidx (0) = 0;
464  x_nz = b.nnz ();
465  ii = 0;
466  volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
467  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
468  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
469  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
470  {
471  octave_quit ();
472  for (octave_idx_type j = 0; j < b_nr; j++)
473  Xx[j] = b.xelem (j,i);
474  for (octave_idx_type j = nr; j < nbuf; j++)
475  buf[j] = 0.;
477 #if defined (CS_VER) && (CS_VER >= 2)
478  CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
479 #else
480  CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xx, buf);
481 #endif
482  CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
484  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
485  {
486  octave_quit ();
488  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
490  }
492 #if defined (CS_VER) && (CS_VER >= 2)
493  CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
494 #else
495  CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xx);
496 #endif
498 
499  for (octave_idx_type j = 0; j < nc; j++)
500  {
501  double tmp = Xx[j];
502  if (tmp != 0.0)
503  {
504  if (ii == x_nz)
505  {
506  // Resize the sparse matrix
507  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
508  sz = (sz > 10 ? sz : 10) + x_nz;
509  x.change_capacity (sz);
510  x_nz = sz;
511  }
512  x.xdata (ii) = tmp;
513  x.xridx (ii++) = j;
514  }
515  }
516  x.xcidx (i+1) = ii;
517  }
518  info = 0;
519  }
520 
521  x.maybe_compress ();
522  return x;
523 #else
524  return SparseMatrix ();
525 #endif
526 }
527 
530 {
531  info = -1;
532 #ifdef HAVE_CXSPARSE
533  octave_idx_type nr = a.rows ();
534  octave_idx_type nc = a.cols ();
535  octave_idx_type b_nc = b.cols ();
536  octave_idx_type b_nr = b.rows ();
538 
539  if (nr < 0 || nc < 0 || nr != b_nr)
540  (*current_liboctave_error_handler)
541  ("matrix dimension mismatch in solution of minimum norm problem");
542  else if (nr == 0 || nc == 0 || b_nc == 0)
543  x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
544  else if (nr >= nc)
545  {
546  SparseQR q (a, 3);
547  if (! q.ok ())
548  return ComplexMatrix ();
549  x.resize (nc, b_nc);
550  Complex *vec = x.fortran_vec ();
551  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
552  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
553  OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
554  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
555  {
556  octave_quit ();
557  for (octave_idx_type j = 0; j < b_nr; j++)
558  {
559  Complex c = b.xelem (j,i);
560  Xx[j] = std::real (c);
561  Xz[j] = std::imag (c);
562  }
563  for (octave_idx_type j = nr; j < q.S ()->m2; j++)
564  buf[j] = 0.;
566 #if defined (CS_VER) && (CS_VER >= 2)
567  CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
568 #else
569  CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xx, buf);
570 #endif
572  for (volatile octave_idx_type j = 0; j < nc; j++)
573  {
574  octave_quit ();
576  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
578  }
580  CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
581 #if defined (CS_VER) && (CS_VER >= 2)
582  CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
583 #else
584  CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xx);
585 #endif
586  for (octave_idx_type j = nr; j < q.S ()->m2; j++)
587  buf[j] = 0.;
588 #if defined (CS_VER) && (CS_VER >= 2)
589  CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr);
590 #else
591  CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xz, buf);
592 #endif
594  for (volatile octave_idx_type j = 0; j < nc; j++)
595  {
596  octave_quit ();
598  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
600  }
602  CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
603 #if defined (CS_VER) && (CS_VER >= 2)
604  CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc);
605 #else
606  CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xz);
607 #endif
609  for (octave_idx_type j = 0; j < nc; j++)
610  vec[j+idx] = Complex (Xx[j], Xz[j]);
611  }
612  info = 0;
613  }
614  else
615  {
616  SparseMatrix at = a.hermitian ();
617  SparseQR q (at, 3);
618  if (! q.ok ())
619  return ComplexMatrix ();
620  x.resize (nc, b_nc);
621  Complex *vec = x.fortran_vec ();
622  volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
623  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
624  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
625  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
626  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
627  {
628  octave_quit ();
629  for (octave_idx_type j = 0; j < b_nr; j++)
630  {
631  Complex c = b.xelem (j,i);
632  Xx[j] = std::real (c);
633  Xz[j] = std::imag (c);
634  }
635  for (octave_idx_type j = nr; j < nbuf; j++)
636  buf[j] = 0.;
638 #if defined (CS_VER) && (CS_VER >= 2)
639  CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
640 #else
641  CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xx, buf);
642 #endif
643  CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
645  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
646  {
647  octave_quit ();
649  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
651  }
653 #if defined (CS_VER) && (CS_VER >= 2)
654  CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
655 #else
656  CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xx);
657 #endif
659  for (octave_idx_type j = nr; j < nbuf; j++)
660  buf[j] = 0.;
662 #if defined (CS_VER) && (CS_VER >= 2)
663  CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr);
664 #else
665  CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xz, buf);
666 #endif
667  CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
669  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
670  {
671  octave_quit ();
673  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
675  }
677 #if defined (CS_VER) && (CS_VER >= 2)
678  CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc);
679 #else
680  CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xz);
681 #endif
683  for (octave_idx_type j = 0; j < nc; j++)
684  vec[j+idx] = Complex (Xx[j], Xz[j]);
685  }
686  info = 0;
687  }
688 
689  return x;
690 #else
691  return ComplexMatrix ();
692 #endif
693 }
694 
697  octave_idx_type &info)
698 {
699  info = -1;
700 #ifdef HAVE_CXSPARSE
701  octave_idx_type nr = a.rows ();
702  octave_idx_type nc = a.cols ();
703  octave_idx_type b_nr = b.rows ();
704  octave_idx_type b_nc = b.cols ();
706  volatile octave_idx_type ii, x_nz;
707 
708  if (nr < 0 || nc < 0 || nr != b_nr)
709  (*current_liboctave_error_handler)
710  ("matrix dimension mismatch in solution of minimum norm problem");
711  else if (nr == 0 || nc == 0 || b_nc == 0)
712  x = SparseComplexMatrix (nc, b_nc);
713  else if (nr >= nc)
714  {
715  SparseQR q (a, 3);
716  if (! q.ok ())
717  return SparseComplexMatrix ();
718  x = SparseComplexMatrix (nc, b_nc, b.nnz ());
719  x.xcidx (0) = 0;
720  x_nz = b.nnz ();
721  ii = 0;
722  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
723  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
724  OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
725  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
726  {
727  octave_quit ();
728  for (octave_idx_type j = 0; j < b_nr; j++)
729  {
730  Complex c = b.xelem (j,i);
731  Xx[j] = std::real (c);
732  Xz[j] = std::imag (c);
733  }
734  for (octave_idx_type j = nr; j < q.S ()->m2; j++)
735  buf[j] = 0.;
737 #if defined (CS_VER) && (CS_VER >= 2)
738  CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
739 #else
740  CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xx, buf);
741 #endif
743  for (volatile octave_idx_type j = 0; j < nc; j++)
744  {
745  octave_quit ();
747  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
749  }
751  CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
752 #if defined (CS_VER) && (CS_VER >= 2)
753  CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
754 #else
755  CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xx);
756 #endif
758  for (octave_idx_type j = nr; j < q.S ()->m2; j++)
759  buf[j] = 0.;
761 #if defined (CS_VER) && (CS_VER >= 2)
762  CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr);
763 #else
764  CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xz, buf);
765 #endif
767  for (volatile octave_idx_type j = 0; j < nc; j++)
768  {
769  octave_quit ();
771  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
773  }
775  CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
776 #if defined (CS_VER) && (CS_VER >= 2)
777  CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc);
778 #else
779  CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xz);
780 #endif
782 
783  for (octave_idx_type j = 0; j < nc; j++)
784  {
785  Complex tmp = Complex (Xx[j], Xz[j]);
786  if (tmp != 0.0)
787  {
788  if (ii == x_nz)
789  {
790  // Resize the sparse matrix
791  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
792  sz = (sz > 10 ? sz : 10) + x_nz;
793  x.change_capacity (sz);
794  x_nz = sz;
795  }
796  x.xdata (ii) = tmp;
797  x.xridx (ii++) = j;
798  }
799  }
800  x.xcidx (i+1) = ii;
801  }
802  info = 0;
803  }
804  else
805  {
806  SparseMatrix at = a.hermitian ();
807  SparseQR q (at, 3);
808  if (! q.ok ())
809  return SparseComplexMatrix ();
810  x = SparseComplexMatrix (nc, b_nc, b.nnz ());
811  x.xcidx (0) = 0;
812  x_nz = b.nnz ();
813  ii = 0;
814  volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
815  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
816  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
817  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
818  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
819  {
820  octave_quit ();
821  for (octave_idx_type j = 0; j < b_nr; j++)
822  {
823  Complex c = b.xelem (j,i);
824  Xx[j] = std::real (c);
825  Xz[j] = std::imag (c);
826  }
827  for (octave_idx_type j = nr; j < nbuf; j++)
828  buf[j] = 0.;
830 #if defined (CS_VER) && (CS_VER >= 2)
831  CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
832 #else
833  CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xx, buf);
834 #endif
835  CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
837  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
838  {
839  octave_quit ();
841  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
843  }
845 #if defined (CS_VER) && (CS_VER >= 2)
846  CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
847 #else
848  CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xx);
849 #endif
851  for (octave_idx_type j = nr; j < nbuf; j++)
852  buf[j] = 0.;
854 #if defined (CS_VER) && (CS_VER >= 2)
855  CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr);
856 #else
857  CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xz, buf);
858 #endif
859  CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
861  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
862  {
863  octave_quit ();
865  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
867  }
869 #if defined (CS_VER) && (CS_VER >= 2)
870  CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc);
871 #else
872  CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xz);
873 #endif
875 
876  for (octave_idx_type j = 0; j < nc; j++)
877  {
878  Complex tmp = Complex (Xx[j], Xz[j]);
879  if (tmp != 0.0)
880  {
881  if (ii == x_nz)
882  {
883  // Resize the sparse matrix
884  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
885  sz = (sz > 10 ? sz : 10) + x_nz;
886  x.change_capacity (sz);
887  x_nz = sz;
888  }
889  x.xdata (ii) = tmp;
890  x.xridx (ii++) = j;
891  }
892  }
893  x.xcidx (i+1) = ii;
894  }
895  info = 0;
896  }
897 
898  x.maybe_compress ();
899  return x;
900 #else
901  return SparseComplexMatrix ();
902 #endif
903 }
904 
905 Matrix
906 qrsolve (const SparseMatrix &a, const MArray<double> &b,
907  octave_idx_type &info)
908 {
909  return qrsolve (a, Matrix (b), info);
910 }
911 
914  octave_idx_type &info)
915 {
916  return qrsolve (a, ComplexMatrix (b), info);
917 }