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 
32 SparseQR::SparseQR_rep::SparseQR_rep (const SparseMatrix& a, int order)
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;
50  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
58  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
69 SparseQR::SparseQR_rep::~SparseQR_rep (void)
70 {
71 #ifdef HAVE_CXSPARSE
72  CXSPARSE_DNAME (_sfree) (S);
73  CXSPARSE_DNAME (_nfree) (N);
74 #endif
75 }
76 
77 SparseMatrix
78 SparseQR::SparseQR_rep::V (void) const
79 {
80 #ifdef HAVE_CXSPARSE
81  // Drop zeros from V and sort
82  // FIXME: Is the double transpose to sort necessary?
83  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
89  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
107 ColumnVector
108 SparseQR::SparseQR_rep::Pinv (void) const
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 
124 ColumnVector
125 SparseQR::SparseQR_rep::P (void) const
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 
141 SparseMatrix
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?
147  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
153  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
174 SparseQR::SparseQR_rep::C (const Matrix &b) const
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);
197  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
203  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
204 
205  for (volatile octave_idx_type i = 0; i < nm; i++)
206  {
207  octave_quit ();
208  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
209  CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
210  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
223 SparseQR::SparseQR_rep::Q (void) const
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);
247  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
253  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
254 
255  for (volatile octave_idx_type i = 0; i < nm; i++)
256  {
257  octave_quit ();
258  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
259  CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
260  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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.;
304  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
310  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
311  for (volatile octave_idx_type j = 0; j < nc; j++)
312  {
313  octave_quit ();
314  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
315  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
316  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
317  }
318  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
325  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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.;
345  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
352  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
353  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
354  {
355  octave_quit ();
356  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
357  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
358  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
359  }
360  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
366  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
367  }
368  info = 0;
369  }
370 
371  return x;
372 #else
373  return Matrix ();
374 #endif
375 }
376 
377 SparseMatrix
378 qrsolve (const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info)
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.;
412  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
418  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
419  for (volatile octave_idx_type j = 0; j < nc; j++)
420  {
421  octave_quit ();
422  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
423  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
424  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
425  }
426  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
433  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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.;
476  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
483  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
484  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
485  {
486  octave_quit ();
487  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
488  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
489  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
490  }
491  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
497  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
528 ComplexMatrix
529 qrsolve (const SparseMatrix&a, const ComplexMatrix &b, octave_idx_type &info)
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 ();
537  ComplexMatrix x;
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.;
565  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
571  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
572  for (volatile octave_idx_type j = 0; j < nc; j++)
573  {
574  octave_quit ();
575  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
576  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
577  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
578  }
579  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
593  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
594  for (volatile octave_idx_type j = 0; j < nc; j++)
595  {
596  octave_quit ();
597  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
598  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
599  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
600  }
601  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
608  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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.;
637  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
644  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
645  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
646  {
647  octave_quit ();
648  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
649  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
650  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
651  }
652  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
658  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
659  for (octave_idx_type j = nr; j < nbuf; j++)
660  buf[j] = 0.;
661  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
668  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
669  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
670  {
671  octave_quit ();
672  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
673  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
674  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
675  }
676  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
682  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
695 SparseComplexMatrix
696 qrsolve (const SparseMatrix&a, const SparseComplexMatrix &b,
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 ();
705  SparseComplexMatrix x;
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.;
736  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
742  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
743  for (volatile octave_idx_type j = 0; j < nc; j++)
744  {
745  octave_quit ();
746  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
747  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
748  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
749  }
750  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
757  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
758  for (octave_idx_type j = nr; j < q.S ()->m2; j++)
759  buf[j] = 0.;
760  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
766  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
767  for (volatile octave_idx_type j = 0; j < nc; j++)
768  {
769  octave_quit ();
770  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
771  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
772  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
773  }
774  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
781  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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.;
829  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
836  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
837  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
838  {
839  octave_quit ();
840  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
841  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
842  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
843  }
844  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
850  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
851  for (octave_idx_type j = nr; j < nbuf; j++)
852  buf[j] = 0.;
853  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
860  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
861  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
862  {
863  octave_quit ();
864  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
865  CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
866  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
867  }
868  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
874  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
912 ComplexMatrix
913 qrsolve (const SparseMatrix &a, const MArray<Complex> &b,
914  octave_idx_type &info)
915 {
916  return qrsolve (a, ComplexMatrix (b), info);
917 }
Generated on Mon Dec 30 2013 03:04:50 for GNU Octave by   doxygen 1.8.1.2