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
SparseCmplxQR.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 "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 
51 SparseComplexQR::SparseComplexQR_rep::SparseComplexQR_rep
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;
71  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
78  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
89 SparseComplexQR::SparseComplexQR_rep::~SparseComplexQR_rep (void)
90 {
91 #ifdef HAVE_CXSPARSE
92  CXSPARSE_ZNAME (_sfree) (S);
93  CXSPARSE_ZNAME (_nfree) (N);
94 #endif
95 }
96 
97 SparseComplexMatrix
98 SparseComplexQR::SparseComplexQR_rep::V (void) const
99 {
100 #ifdef HAVE_CXSPARSE
101  // Drop zeros from V and sort
102  // FIXME: Is the double transpose to sort necessary?
103  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
109  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
127 ColumnVector
128 SparseComplexQR::SparseComplexQR_rep::Pinv (void) const
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 
144 ColumnVector
145 SparseComplexQR::SparseComplexQR_rep::P (void) const
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 
161 SparseComplexMatrix
162 SparseComplexQR::SparseComplexQR_rep::R (const bool econ) const
163 {
164 #ifdef HAVE_CXSPARSE
165  // Drop zeros from R and sort
166  // FIXME: Is the double transpose to sort necessary?
167  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
173  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
191 ComplexMatrix
192 SparseComplexQR::SparseComplexQR_rep::C (const ComplexMatrix &b) const
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);
214  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
222  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
223  for (volatile octave_idx_type i = 0; i < nm; i++)
224  {
225  octave_quit ();
226  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
227  CXSPARSE_ZNAME (_happly)
228  (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf));
229  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
241 ComplexMatrix
242 SparseComplexQR::SparseComplexQR_rep::Q (void) const
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);
264  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
272  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
273  for (volatile octave_idx_type i = 0; i < nm; i++)
274  {
275  octave_quit ();
276  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
277  CXSPARSE_ZNAME (_happly)
278  (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf));
279  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
292 ComplexMatrix
293 qrsolve (const SparseComplexMatrix&a, const Matrix &b, octave_idx_type &info)
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 ();
301  ComplexMatrix x;
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;
324  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
332  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
333  for (volatile octave_idx_type j = 0; j < nc; j++)
334  {
335  octave_quit ();
336  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
337  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
338  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
339  }
340  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
347  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
367  OCTAVE_LOCAL_BUFFER (Complex, B, nr);
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;
378  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
387  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
388  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
389  {
390  octave_quit ();
391  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
399  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
400  }
401  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
407  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
408  }
409  info = 0;
410  }
411 
412  return x;
413 #else
414  return ComplexMatrix ();
415 #endif
416 }
417 
418 SparseComplexMatrix
419 qrsolve (const SparseComplexMatrix&a, const SparseMatrix &b,
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 ();
428  SparseComplexMatrix x;
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;
454  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
462  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
463  for (volatile octave_idx_type j = 0; j < nc; j++)
464  {
465  octave_quit ();
466  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
467  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
468  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
469  }
470  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
479  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
521  OCTAVE_LOCAL_BUFFER (Complex, B, nr);
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;
532  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
541  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
542  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
543  {
544  octave_quit ();
545  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
552  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
553  }
554  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
562  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
593 ComplexMatrix
594 qrsolve (const SparseComplexMatrix&a, const ComplexMatrix &b,
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 ());
605  ComplexMatrix x;
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;
627  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
633  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
634  for (volatile octave_idx_type j = 0; j < nc; j++)
635  {
636  octave_quit ();
637  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
638  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
639  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
640  }
641  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
648  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
667  OCTAVE_LOCAL_BUFFER (Complex, B, nr);
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;
677  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
684  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
685  for (volatile octave_idx_type j = nr-1; j >= 0; j--)
686  {
687  octave_quit ();
688  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
695  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
696  }
697  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
703  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
704  }
705  info = 0;
706  }
707 
708  return x;
709 #else
710  return ComplexMatrix ();
711 #endif
712 }
713 
714 SparseComplexMatrix
715 qrsolve (const SparseComplexMatrix&a, const SparseComplexMatrix &b,
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 ();
724  SparseComplexMatrix x;
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;
750  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
758  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
759  for (volatile octave_idx_type j = 0; j < nc; j++)
760  {
761  octave_quit ();
762  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
763  CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
764  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
765  }
766  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
775  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
816  OCTAVE_LOCAL_BUFFER (Complex, B, nr);
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;
827  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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);
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 #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
847  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
848  }
849  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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
857  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
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 
888 ComplexMatrix
889 qrsolve (const SparseComplexMatrix &a, const MArray<double> &b,
890  octave_idx_type &info)
891 {
892  return qrsolve (a, Matrix (b), info);
893 }
894 
895 ComplexMatrix
896 qrsolve (const SparseComplexMatrix &a, const MArray<Complex> &b,
897  octave_idx_type &info)
898 {
899  return qrsolve (a, ComplexMatrix (b), info);
900 }
Generated on Mon Dec 30 2013 03:04:50 for GNU Octave by   doxygen 1.8.1.2