GNU Octave  4.0.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
SparseCmplxQR.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2005-2015 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 "SparseCmplxQR.h"
30 #include "oct-locbuf.h"
31 
32 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER < 2)) || (CS_VER < 2))
33 typedef double _Complex cs_complex_t;
34 
35 // Why did g++ 4.x stl_vector.h make
36 // OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, n)
37 // an error ?
38 #define OCTAVE_C99_COMPLEX(buf, n) \
39  OCTAVE_LOCAL_BUFFER (double, buf ## tmp, (2 * (n))); \
40  cs_complex_t *buf = reinterpret_cast<cs_complex_t *> (buf ## tmp);
41 
42 #define OCTAVE_C99_ZERO (0. + 0.iF)
43 #define OCTAVE_C99_ONE (1. + 0.iF)
44 #else
45 #define OCTAVE_C99_COMPLEX(buf, n) \
46  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, (n));
47 #define OCTAVE_C99_ZERO cs_complex_t(0., 0.);
48 #define OCTAVE_C99_ONE cs_complex_t(1., 0.);
49 #endif
50 
52  (GCC_ATTR_UNUSED const SparseComplexMatrix& a, GCC_ATTR_UNUSED int order)
53  : count (1), nrows (0)
54 #ifdef HAVE_CXSPARSE
55  , S (0), N (0)
56 #endif
57 {
58 #ifdef HAVE_CXSPARSE
59  CXSPARSE_ZNAME () A;
60  A.nzmax = a.nnz ();
61  A.m = a.rows ();
62  A.n = a.cols ();
63  nrows = A.m;
64  // Cast away const on A, with full knowledge that CSparse won't touch it
65  // Prevents the methods below making a copy of the data.
66  A.p = const_cast<octave_idx_type *>(a.cidx ());
67  A.i = const_cast<octave_idx_type *>(a.ridx ());
68  A.x = const_cast<cs_complex_t *>(reinterpret_cast<const cs_complex_t *>
69  (a.data ()));
70  A.nz = -1;
72 #if defined (CS_VER) && (CS_VER >= 2)
73  S = CXSPARSE_ZNAME (_sqr) (order, &A, 1);
74 #else
75  S = CXSPARSE_ZNAME (_sqr) (&A, order - 1, 1);
76 #endif
77  N = CXSPARSE_ZNAME (_qr) (&A, S);
79  if (!N)
80  (*current_liboctave_error_handler)
81  ("SparseComplexQR: sparse matrix QR factorization filled");
82  count = 1;
83 #else
84  (*current_liboctave_error_handler)
85  ("SparseComplexQR: sparse matrix QR factorization not implemented");
86 #endif
87 }
88 
90 {
91 #ifdef HAVE_CXSPARSE
92  CXSPARSE_ZNAME (_sfree) (S);
93  CXSPARSE_ZNAME (_nfree) (N);
94 #endif
95 }
96 
99 {
100 #ifdef HAVE_CXSPARSE
101  // Drop zeros from V and sort
102  // FIXME: Is the double transpose to sort necessary?
104  CXSPARSE_ZNAME (_dropzeros) (N->L);
105  CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->L, 1);
106  CXSPARSE_ZNAME (_spfree) (N->L);
107  N->L = CXSPARSE_ZNAME (_transpose) (D, 1);
108  CXSPARSE_ZNAME (_spfree) (D);
110 
111  octave_idx_type nc = N->L->n;
112  octave_idx_type nz = N->L->nzmax;
113  SparseComplexMatrix ret (N->L->m, nc, nz);
114  for (octave_idx_type j = 0; j < nc+1; j++)
115  ret.xcidx (j) = N->L->p[j];
116  for (octave_idx_type j = 0; j < nz; j++)
117  {
118  ret.xridx (j) = N->L->i[j];
119  ret.xdata (j) = reinterpret_cast<Complex *>(N->L->x)[j];
120  }
121  return ret;
122 #else
123  return SparseComplexMatrix ();
124 #endif
125 }
126 
129 {
130 #ifdef HAVE_CXSPARSE
131  ColumnVector ret(N->L->m);
132  for (octave_idx_type i = 0; i < N->L->m; i++)
133 #if defined (CS_VER) && (CS_VER >= 2)
134  ret.xelem (i) = S->pinv[i];
135 #else
136  ret.xelem (i) = S->Pinv[i];
137 #endif
138  return ret;
139 #else
140  return ColumnVector ();
141 #endif
142 }
143 
146 {
147 #ifdef HAVE_CXSPARSE
148  ColumnVector ret(N->L->m);
149  for (octave_idx_type i = 0; i < N->L->m; i++)
150 #if defined (CS_VER) && (CS_VER >= 2)
151  ret.xelem (S->pinv[i]) = i;
152 #else
153  ret.xelem (S->Pinv[i]) = i;
154 #endif
155  return ret;
156 #else
157  return ColumnVector ();
158 #endif
159 }
160 
163 {
164 #ifdef HAVE_CXSPARSE
165  // Drop zeros from R and sort
166  // FIXME: Is the double transpose to sort necessary?
168  CXSPARSE_ZNAME (_dropzeros) (N->U);
169  CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->U, 1);
170  CXSPARSE_ZNAME (_spfree) (N->U);
171  N->U = CXSPARSE_ZNAME (_transpose) (D, 1);
172  CXSPARSE_ZNAME (_spfree) (D);
174 
175  octave_idx_type nc = N->U->n;
176  octave_idx_type nz = N->U->nzmax;
177  SparseComplexMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
178  for (octave_idx_type j = 0; j < nc+1; j++)
179  ret.xcidx (j) = N->U->p[j];
180  for (octave_idx_type j = 0; j < nz; j++)
181  {
182  ret.xridx (j) = N->U->i[j];
183  ret.xdata (j) = reinterpret_cast<Complex *>(N->U->x)[j];
184  }
185  return ret;
186 #else
187  return SparseComplexMatrix ();
188 #endif
189 }
190 
193 {
194 #ifdef HAVE_CXSPARSE
195  octave_idx_type b_nr = b.rows ();
196  octave_idx_type b_nc = b.cols ();
197  octave_idx_type nc = N->L->n;
198  octave_idx_type nr = nrows;
199  const cs_complex_t *bvec =
200  reinterpret_cast<const cs_complex_t *>(b.fortran_vec ());
201  ComplexMatrix ret(b_nr, b_nc);
202  Complex *vec = ret.fortran_vec ();
203  if (nr < 0 || nc < 0 || nr != b_nr)
204  (*current_liboctave_error_handler) ("matrix dimension mismatch");
205  else if (nr == 0 || nc == 0 || b_nc == 0)
206  ret = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
207  else
208  {
209  OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
210  for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
211  {
212  octave_quit ();
213  volatile octave_idx_type nm = (nr < nc ? nr : nc);
215 #if defined (CS_VER) && (CS_VER >= 2)
216  CXSPARSE_ZNAME (_ipvec)
217  (S->pinv, bvec + idx, reinterpret_cast<cs_complex_t *>(buf), b_nr);
218 #else
219  CXSPARSE_ZNAME (_ipvec)
220  (b_nr, S->Pinv, bvec + idx, reinterpret_cast<cs_complex_t *>(buf));
221 #endif
223  for (volatile octave_idx_type i = 0; i < nm; i++)
224  {
225  octave_quit ();
227  CXSPARSE_ZNAME (_happly)
228  (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf));
230  }
231  for (octave_idx_type i = 0; i < b_nr; i++)
232  vec[i+idx] = buf[i];
233  }
234  }
235  return ret;
236 #else
237  return ComplexMatrix ();
238 #endif
239 }
240 
243 {
244 #ifdef HAVE_CXSPARSE
245  octave_idx_type nc = N->L->n;
246  octave_idx_type nr = nrows;
247  ComplexMatrix ret(nr, nr);
248  Complex *vec = ret.fortran_vec ();
249  if (nr < 0 || nc < 0)
250  (*current_liboctave_error_handler) ("matrix dimension mismatch");
251  else if (nr == 0 || nc == 0)
252  ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0));
253  else
254  {
255  OCTAVE_C99_COMPLEX (bvec, nr);
256  for (octave_idx_type i = 0; i < nr; i++)
257  bvec[i] = OCTAVE_C99_ZERO;
258  OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
259  for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
260  {
261  octave_quit ();
262  bvec[j] = OCTAVE_C99_ONE;
263  volatile octave_idx_type nm = (nr < nc ? nr : nc);
265 #if defined (CS_VER) && (CS_VER >= 2)
266  CXSPARSE_ZNAME (_ipvec)
267  (S->pinv, bvec, reinterpret_cast<cs_complex_t *>(buf), nr);
268 #else
269  CXSPARSE_ZNAME (_ipvec)
270  (nr, S->Pinv, bvec, reinterpret_cast<cs_complex_t *>(buf));
271 #endif
273  for (volatile octave_idx_type i = 0; i < nm; i++)
274  {
275  octave_quit ();
277  CXSPARSE_ZNAME (_happly)
278  (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf));
280  }
281  for (octave_idx_type i = 0; i < nr; i++)
282  vec[i+idx] = buf[i];
283  bvec[j] = OCTAVE_C99_ZERO;
284  }
285  }
286  return ret.hermitian ();
287 #else
288  return ComplexMatrix ();
289 #endif
290 }
291 
294 {
295  info = -1;
296 #ifdef HAVE_CXSPARSE
297  octave_idx_type nr = a.rows ();
298  octave_idx_type nc = a.cols ();
299  octave_idx_type b_nc = b.cols ();
300  octave_idx_type b_nr = b.rows ();
302 
303  if (nr < 0 || nc < 0 || nr != b_nr)
304  (*current_liboctave_error_handler)
305  ("matrix dimension mismatch in solution of minimum norm problem");
306  else if (nr == 0 || nc == 0 || b_nc == 0)
307  x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
308  else if (nr >= nc)
309  {
310  SparseComplexQR q (a, 2);
311  if (! q.ok ())
312  return ComplexMatrix ();
313  x.resize (nc, b_nc);
314  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
315  OCTAVE_C99_COMPLEX (buf, q.S ()->m2);
316  OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
317  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
318  {
319  octave_quit ();
320  for (octave_idx_type j = 0; j < b_nr; j++)
321  Xx[j] = b.xelem (j,i);
322  for (octave_idx_type j = nr; j < q.S ()->m2; j++)
323  buf[j] = OCTAVE_C99_ZERO;
325 #if defined (CS_VER) && (CS_VER >= 2)
326  CXSPARSE_ZNAME (_ipvec)
327  (q.S ()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
328 #else
329  CXSPARSE_ZNAME (_ipvec)
330  (nr, q.S ()->Pinv, reinterpret_cast<cs_complex_t *>(Xx), buf);
331 #endif
333  for (volatile octave_idx_type j = 0; j < nc; j++)
334  {
335  octave_quit ();
337  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
339  }
341  CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf);
342 #if defined (CS_VER) && (CS_VER >= 2)
343  CXSPARSE_ZNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc);
344 #else
345  CXSPARSE_ZNAME (_ipvec) (nc, q.S ()->Q, buf, vec + idx);
346 #endif
348  }
349  info = 0;
350  }
351  else
352  {
353  SparseComplexMatrix at = a.hermitian ();
354  SparseComplexQR q (at, 2);
355  if (! q.ok ())
356  return ComplexMatrix ();
357  x.resize (nc, b_nc);
358  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
359  volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
360  OCTAVE_C99_COMPLEX (buf, nbuf);
361  OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
362 #if defined (CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
363  OCTAVE_LOCAL_BUFFER (double, B, nr);
364  for (octave_idx_type i = 0; i < nr; i++)
365  B[i] = q.N ()->B[i];
366 #else
368  for (octave_idx_type i = 0; i < nr; i++)
369  B[i] = conj (reinterpret_cast<Complex *>(q.N ()->B)[i]);
370 #endif
371  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
372  {
373  octave_quit ();
374  for (octave_idx_type j = 0; j < b_nr; j++)
375  Xx[j] = b.xelem (j,i);
376  for (octave_idx_type j = nr; j < nbuf; j++)
377  buf[j] = OCTAVE_C99_ZERO;
379 #if defined (CS_VER) && (CS_VER >= 2)
380  CXSPARSE_ZNAME (_pvec)
381  (q.S ()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
382 #else
383  CXSPARSE_ZNAME (_pvec)
384  (nr, q.S ()->Q, reinterpret_cast<cs_complex_t *>(Xx), buf);
385 #endif
386  CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf);
388  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
389  {
390  octave_quit ();
392 
393 #if defined (CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
394  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf);
395 #else
396  CXSPARSE_ZNAME (_happly)
397  (q.N ()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf);
398 #endif
400  }
402 #if defined (CS_VER) && (CS_VER >= 2)
403  CXSPARSE_ZNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc);
404 #else
405  CXSPARSE_ZNAME (_pvec) (nc, q.S ()->Pinv, buf, vec + idx);
406 #endif
408  }
409  info = 0;
410  }
411 
412  return x;
413 #else
414  return ComplexMatrix ();
415 #endif
416 }
417 
420  octave_idx_type &info)
421 {
422  info = -1;
423 #ifdef HAVE_CXSPARSE
424  octave_idx_type nr = a.rows ();
425  octave_idx_type nc = a.cols ();
426  octave_idx_type b_nc = b.cols ();
427  octave_idx_type b_nr = b.rows ();
429  volatile octave_idx_type ii, x_nz;
430 
431  if (nr < 0 || nc < 0 || nr != b_nr)
432  (*current_liboctave_error_handler)
433  ("matrix dimension mismatch in solution of minimum norm problem");
434  else if (nr == 0 || nc == 0 || b_nc == 0)
435  x = SparseComplexMatrix (nc, b_nc);
436  else if (nr >= nc)
437  {
438  SparseComplexQR q (a, 2);
439  if (! q.ok ())
440  return SparseComplexMatrix ();
441  x = SparseComplexMatrix (nc, b_nc, b.nnz ());
442  x.xcidx (0) = 0;
443  x_nz = b.nnz ();
444  ii = 0;
445  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
446  OCTAVE_C99_COMPLEX (buf, q.S ()->m2);
447  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
448  {
449  octave_quit ();
450  for (octave_idx_type j = 0; j < b_nr; j++)
451  Xx[j] = b.xelem (j,i);
452  for (octave_idx_type j = nr; j < q.S ()->m2; j++)
453  buf[j] = OCTAVE_C99_ZERO;
455 #if defined (CS_VER) && (CS_VER >= 2)
456  CXSPARSE_ZNAME (_ipvec)
457  (q.S ()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
458 #else
459  CXSPARSE_ZNAME (_ipvec)
460  (nr, q.S ()->Pinv, reinterpret_cast<cs_complex_t *>(Xx), buf);
461 #endif
463  for (volatile octave_idx_type j = 0; j < nc; j++)
464  {
465  octave_quit ();
467  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
469  }
471  CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf);
472 #if defined (CS_VER) && (CS_VER >= 2)
473  CXSPARSE_ZNAME (_ipvec)
474  (q.S ()->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
475 #else
476  CXSPARSE_ZNAME (_ipvec)
477  (nc, q.S ()->Q, buf, reinterpret_cast<cs_complex_t *>(Xx));
478 #endif
480 
481  for (octave_idx_type j = 0; j < nc; j++)
482  {
483  Complex tmp = Xx[j];
484  if (tmp != 0.0)
485  {
486  if (ii == x_nz)
487  {
488  // Resize the sparse matrix
489  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
490  sz = (sz > 10 ? sz : 10) + x_nz;
491  x.change_capacity (sz);
492  x_nz = sz;
493  }
494  x.xdata (ii) = tmp;
495  x.xridx (ii++) = j;
496  }
497  }
498  x.xcidx (i+1) = ii;
499  }
500  info = 0;
501  }
502  else
503  {
504  SparseComplexMatrix at = a.hermitian ();
505  SparseComplexQR q (at, 2);
506  if (! q.ok ())
507  return SparseComplexMatrix ();
508  x = SparseComplexMatrix (nc, b_nc, b.nnz ());
509  x.xcidx (0) = 0;
510  x_nz = b.nnz ();
511  ii = 0;
512  volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
513  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
514  OCTAVE_C99_COMPLEX (buf, nbuf);
515 
516 #if defined (CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
517  OCTAVE_LOCAL_BUFFER (double, B, nr);
518  for (octave_idx_type i = 0; i < nr; i++)
519  B[i] = q.N ()->B[i];
520 #else
522  for (octave_idx_type i = 0; i < nr; i++)
523  B[i] = conj (reinterpret_cast<Complex *>(q.N ()->B)[i]);
524 #endif
525  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
526  {
527  octave_quit ();
528  for (octave_idx_type j = 0; j < b_nr; j++)
529  Xx[j] = b.xelem (j,i);
530  for (octave_idx_type j = nr; j < nbuf; j++)
531  buf[j] = OCTAVE_C99_ZERO;
533 #if defined (CS_VER) && (CS_VER >= 2)
534  CXSPARSE_ZNAME (_pvec)
535  (q.S ()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
536 #else
537  CXSPARSE_ZNAME (_pvec)
538  (nr, q.S ()->Q, reinterpret_cast<cs_complex_t *>(Xx), buf);
539 #endif
540  CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf);
542  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
543  {
544  octave_quit ();
546 #if defined (CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
547  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf);
548 #else
549  CXSPARSE_ZNAME (_happly)
550  (q.N ()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf);
551 #endif
553  }
555 #if defined (CS_VER) && (CS_VER >= 2)
556  CXSPARSE_ZNAME (_pvec)
557  (q.S ()->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
558 #else
559  CXSPARSE_ZNAME (_pvec)
560  (nc, q.S ()->Pinv, buf, reinterpret_cast<cs_complex_t *>(Xx));
561 #endif
563 
564  for (octave_idx_type j = 0; j < nc; j++)
565  {
566  Complex tmp = Xx[j];
567  if (tmp != 0.0)
568  {
569  if (ii == x_nz)
570  {
571  // Resize the sparse matrix
572  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
573  sz = (sz > 10 ? sz : 10) + x_nz;
574  x.change_capacity (sz);
575  x_nz = sz;
576  }
577  x.xdata (ii) = tmp;
578  x.xridx (ii++) = j;
579  }
580  }
581  x.xcidx (i+1) = ii;
582  }
583  info = 0;
584  }
585 
586  x.maybe_compress ();
587  return x;
588 #else
589  return SparseComplexMatrix ();
590 #endif
591 }
592 
595  octave_idx_type &info)
596 {
597  info = -1;
598 #ifdef HAVE_CXSPARSE
599  octave_idx_type nr = a.rows ();
600  octave_idx_type nc = a.cols ();
601  octave_idx_type b_nc = b.cols ();
602  octave_idx_type b_nr = b.rows ();
603  const cs_complex_t *bvec =
604  reinterpret_cast<const cs_complex_t *>(b.fortran_vec ());
606 
607  if (nr < 0 || nc < 0 || nr != b_nr)
608  (*current_liboctave_error_handler)
609  ("matrix dimension mismatch in solution of minimum norm problem");
610  else if (nr == 0 || nc == 0 || b_nc == 0)
611  x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
612  else if (nr >= nc)
613  {
614  SparseComplexQR q (a, 2);
615  if (! q.ok ())
616  return ComplexMatrix ();
617  x.resize (nc, b_nc);
618  cs_complex_t *vec = reinterpret_cast<cs_complex_t *>
619  (x.fortran_vec ());
620  OCTAVE_C99_COMPLEX (buf, q.S ()->m2);
621  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
622  i++, idx+=nc, bidx+=b_nr)
623  {
624  octave_quit ();
625  for (octave_idx_type j = nr; j < q.S ()->m2; j++)
626  buf[j] = OCTAVE_C99_ZERO;
628 #if defined (CS_VER) && (CS_VER >= 2)
629  CXSPARSE_ZNAME (_ipvec) (q.S ()->pinv, bvec + bidx, buf, nr);
630 #else
631  CXSPARSE_ZNAME (_ipvec) (nr, q.S ()->Pinv, bvec + bidx, buf);
632 #endif
634  for (volatile octave_idx_type j = 0; j < nc; j++)
635  {
636  octave_quit ();
638  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
640  }
642  CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf);
643 #if defined (CS_VER) && (CS_VER >= 2)
644  CXSPARSE_ZNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc);
645 #else
646  CXSPARSE_ZNAME (_ipvec) (nc, q.S ()->Q, buf, vec + idx);
647 #endif
649  }
650  info = 0;
651  }
652  else
653  {
654  SparseComplexMatrix at = a.hermitian ();
655  SparseComplexQR q (at, 2);
656  if (! q.ok ())
657  return ComplexMatrix ();
658  x.resize (nc, b_nc);
659  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
660  volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
661  OCTAVE_C99_COMPLEX (buf, nbuf);
662 #if defined (CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
663  OCTAVE_LOCAL_BUFFER (double, B, nr);
664  for (octave_idx_type i = 0; i < nr; i++)
665  B[i] = q.N ()->B[i];
666 #else
668  for (octave_idx_type i = 0; i < nr; i++)
669  B[i] = conj (reinterpret_cast<Complex *>(q.N ()->B)[i]);
670 #endif
671  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
672  i++, idx+=nc, bidx+=b_nr)
673  {
674  octave_quit ();
675  for (octave_idx_type j = nr; j < nbuf; j++)
676  buf[j] = OCTAVE_C99_ZERO;
678 #if defined (CS_VER) && (CS_VER >= 2)
679  CXSPARSE_ZNAME (_pvec) (q.S ()->q, bvec + bidx, buf, nr);
680 #else
681  CXSPARSE_ZNAME (_pvec) (nr, q.S ()->Q, bvec + bidx, buf);
682 #endif
683  CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf);
685  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
686  {
687  octave_quit ();
689 #if defined (CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
690  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf);
691 #else
692  CXSPARSE_ZNAME (_happly)
693  (q.N ()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf);
694 #endif
696  }
698 #if defined (CS_VER) && (CS_VER >= 2)
699  CXSPARSE_ZNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc);
700 #else
701  CXSPARSE_ZNAME (_pvec) (nc, q.S ()->Pinv, buf, vec + idx);
702 #endif
704  }
705  info = 0;
706  }
707 
708  return x;
709 #else
710  return ComplexMatrix ();
711 #endif
712 }
713 
716  octave_idx_type &info)
717 {
718  info = -1;
719 #ifdef HAVE_CXSPARSE
720  octave_idx_type nr = a.rows ();
721  octave_idx_type nc = a.cols ();
722  octave_idx_type b_nc = b.cols ();
723  octave_idx_type b_nr = b.rows ();
725  volatile octave_idx_type ii, x_nz;
726 
727  if (nr < 0 || nc < 0 || nr != b_nr)
728  (*current_liboctave_error_handler)
729  ("matrix dimension mismatch in solution of minimum norm problem");
730  else if (nr == 0 || nc == 0 || b_nc == 0)
731  x = SparseComplexMatrix (nc, b_nc);
732  else if (nr >= nc)
733  {
734  SparseComplexQR q (a, 2);
735  if (! q.ok ())
736  return SparseComplexMatrix ();
737  x = SparseComplexMatrix (nc, b_nc, b.nnz ());
738  x.xcidx (0) = 0;
739  x_nz = b.nnz ();
740  ii = 0;
741  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
742  OCTAVE_C99_COMPLEX (buf, q.S ()->m2);
743  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
744  {
745  octave_quit ();
746  for (octave_idx_type j = 0; j < b_nr; j++)
747  Xx[j] = b.xelem (j,i);
748  for (octave_idx_type j = nr; j < q.S ()->m2; j++)
749  buf[j] = OCTAVE_C99_ZERO;
751 #if defined (CS_VER) && (CS_VER >= 2)
752  CXSPARSE_ZNAME (_ipvec)
753  (q.S ()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
754 #else
755  CXSPARSE_ZNAME (_ipvec)
756  (nr, q.S ()->Pinv, reinterpret_cast<cs_complex_t *>(Xx), buf);
757 #endif
759  for (volatile octave_idx_type j = 0; j < nc; j++)
760  {
761  octave_quit ();
763  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
765  }
767  CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf);
768 #if defined (CS_VER) && (CS_VER >= 2)
769  CXSPARSE_ZNAME (_ipvec)
770  (q.S ()->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
771 #else
772  CXSPARSE_ZNAME (_ipvec)
773  (nc, q.S ()->Q, buf, reinterpret_cast<cs_complex_t *>(Xx));
774 #endif
776 
777  for (octave_idx_type j = 0; j < nc; j++)
778  {
779  Complex tmp = Xx[j];
780  if (tmp != 0.0)
781  {
782  if (ii == x_nz)
783  {
784  // Resize the sparse matrix
785  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
786  sz = (sz > 10 ? sz : 10) + x_nz;
787  x.change_capacity (sz);
788  x_nz = sz;
789  }
790  x.xdata (ii) = tmp;
791  x.xridx (ii++) = j;
792  }
793  }
794  x.xcidx (i+1) = ii;
795  }
796  info = 0;
797  }
798  else
799  {
800  SparseComplexMatrix at = a.hermitian ();
801  SparseComplexQR q (at, 2);
802  if (! q.ok ())
803  return SparseComplexMatrix ();
804  x = SparseComplexMatrix (nc, b_nc, b.nnz ());
805  x.xcidx (0) = 0;
806  x_nz = b.nnz ();
807  ii = 0;
808  volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
809  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
810  OCTAVE_C99_COMPLEX (buf, nbuf);
811 #if defined (CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
812  OCTAVE_LOCAL_BUFFER (double, B, nr);
813  for (octave_idx_type i = 0; i < nr; i++)
814  B[i] = q.N ()->B[i];
815 #else
817  for (octave_idx_type i = 0; i < nr; i++)
818  B[i] = conj (reinterpret_cast<Complex *>(q.N ()->B)[i]);
819 #endif
820  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
821  {
822  octave_quit ();
823  for (octave_idx_type j = 0; j < b_nr; j++)
824  Xx[j] = b.xelem (j,i);
825  for (octave_idx_type j = nr; j < nbuf; j++)
826  buf[j] = OCTAVE_C99_ZERO;
828 #if defined (CS_VER) && (CS_VER >= 2)
829  CXSPARSE_ZNAME (_pvec)
830  (q.S ()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
831 #else
832  CXSPARSE_ZNAME (_pvec)
833  (nr, q.S ()->Q, reinterpret_cast<cs_complex_t *>(Xx), buf);
834 #endif
835  CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf);
837  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
838  {
839  octave_quit ();
841 #if defined (CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
842  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf);
843 #else
844  CXSPARSE_ZNAME (_happly)
845  (q.N ()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf);
846 #endif
848  }
850 #if defined (CS_VER) && (CS_VER >= 2)
851  CXSPARSE_ZNAME (_pvec)
852  (q.S ()->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
853 #else
854  CXSPARSE_ZNAME (_pvec)
855  (nc, q.S ()->Pinv, buf, reinterpret_cast<cs_complex_t *>(Xx));
856 #endif
858 
859  for (octave_idx_type j = 0; j < nc; j++)
860  {
861  Complex tmp = Xx[j];
862  if (tmp != 0.0)
863  {
864  if (ii == x_nz)
865  {
866  // Resize the sparse matrix
867  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
868  sz = (sz > 10 ? sz : 10) + x_nz;
869  x.change_capacity (sz);
870  x_nz = sz;
871  }
872  x.xdata (ii) = tmp;
873  x.xridx (ii++) = j;
874  }
875  }
876  x.xcidx (i+1) = ii;
877  }
878  info = 0;
879  }
880 
881  x.maybe_compress ();
882  return x;
883 #else
884  return SparseComplexMatrix ();
885 #endif
886 }
887 
890  octave_idx_type &info)
891 {
892  return qrsolve (a, Matrix (b), info);
893 }
894 
897  octave_idx_type &info)
898 {
899  return qrsolve (a, ComplexMatrix (b), info);
900 }
SparseComplexMatrix V(void) const
octave_idx_type * xridx(void)
Definition: Sparse.h:524
octave_idx_type cols(void) const
Definition: Sparse.h:264
octave_idx_type rows(void) const
Definition: Sparse.h:263
T & xelem(octave_idx_type n)
Definition: Sparse.h:300
SparseComplexQR_rep(const SparseComplexMatrix &a, int order)
#define CXSPARSE_ZNAME(name)
Definition: SparseCmplxQR.h:37
F77_RET_T const octave_idx_type Complex * A
Definition: CmplxGEPBAL.cc:39
#define OCTAVE_C99_COMPLEX(buf, n)
octave_idx_type * xcidx(void)
Definition: Sparse.h:537
#define OCTAVE_C99_ZERO
SparseComplexMatrix R(const bool econ) const
octave_idx_type rows(void) const
Definition: Array.h:313
octave_idx_type nnz(void) const
Definition: Sparse.h:248
#define OCTAVE_C99_ONE
ComplexColumnVector conj(const ComplexColumnVector &a)
Definition: CColVector.cc:244
ComplexMatrix hermitian(void) const
Definition: CMatrix.h:149
cs_cis * S(void)
void change_capacity(octave_idx_type nz)
Definition: Sparse.h:493
Sparse< T > maybe_compress(bool remove_zeros=false)
Definition: Sparse.h:469
ComplexMatrix C(const ComplexMatrix &b) const
SparseComplexMatrix hermitian(void) const
Definition: CSparse.cc:679
Definition: dMatrix.h:35
cs_cin * N(void)
T & xelem(octave_idx_type n)
Definition: Array.h:353
F77_RET_T const octave_idx_type Complex const octave_idx_type Complex * B
Definition: CmplxGEPBAL.cc:39
F77_RET_T const octave_idx_type & N
Definition: CmplxGEPBAL.cc:39
#define END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
Definition: quit.h:180
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Definition: CMatrix.h:170
T * xdata(void)
Definition: Sparse.h:511
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
ComplexMatrix Q(void) const
friend ComplexMatrix qrsolve(const SparseComplexMatrix &a, const Matrix &b, octave_idx_type &info)
#define BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE
Definition: quit.h:158
std::complex< double > Complex
Definition: oct-cmplx.h:29
const T * fortran_vec(void) const
Definition: Array.h:481
ColumnVector Pinv(void) const
octave_idx_type cols(void) const
Definition: Array.h:321
bool ok(void) const
F77_RET_T const double * x