Mercurial > hg > octave-nkf

comparison liboctave/SparseQR.cc @ 5610:9761b7d24e9e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-02-09 09:12:02 by dbateman]
author dbateman
date Thu, 09 Feb 2006 09:12:03 +0000
parents
children 69a4f320d95a
comparison
equal deleted inserted replaced
5609:bf96b0f9dbd7 5610:9761b7d24e9e
1 /*
2
3 Copyright (C) 2005 David Bateman
4
5 Octave is free software; you can redistribute it and/or modify it
6 under the terms of the GNU General Public License as published by the
7 Free Software Foundation; either version 2, or (at your option) any
8 later version.
9
10 Octave is distributed in the hope that it will be useful, but WITHOUT
11 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; see the file COPYING. If not, write to the
17 Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
18 Boston, MA 02110-1301, USA.
19
20 */
21
22 #ifdef HAVE_CONFIG_H
23 #include <config.h>
24 #endif
25 #include <vector>
26
27 #include "lo-error.h"
28 #include "SparseQR.h"
29
30 SparseQR::SparseQR_rep::SparseQR_rep (const SparseMatrix& a, int order)
31 {
32 #ifdef HAVE_CXSPARSE
33 CXSPARSE_DNAME (cs) A;
34 A.nzmax = a.nzmax ();
35 A.m = a.rows ();
36 A.n = a.cols ();
37 nrows = A.m;
38 // Cast away const on A, with full knowledge that CSparse won't touch it
39 // Prevents the methods below making a copy of the data.
40 A.p = const_cast<octave_idx_type *>(a.cidx ());
41 A.i = const_cast<octave_idx_type *>(a.ridx ());
42 A.x = const_cast<double *>(a.data ());
43 A.nz = -1;
44 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
45 S = CXSPARSE_DNAME (cs_sqr) (&A, order, 1);
46 N = CXSPARSE_DNAME (cs_qr) (&A, S);
47 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
48 if (!N)
49 (*current_liboctave_error_handler)
50 ("SparseQR: sparse matrix QR factorization filled");
51 count = 1;
52 #else
53 (*current_liboctave_error_handler)
54 ("SparseQR: sparse matrix QR factorization not implemented");
55 #endif
56 }
57
58 SparseQR::SparseQR_rep::~SparseQR_rep (void)
59 {
60 #ifdef HAVE_CXSPARSE
61 CXSPARSE_DNAME (cs_sfree) (S);
62 CXSPARSE_DNAME (cs_nfree) (N);
63 #endif
64 }
65
66 SparseMatrix
67 SparseQR::SparseQR_rep::V (void) const
68 {
69 #ifdef HAVE_CXSPARSE
70 // Drop zeros from V and sort
71 // XXX FIXME XXX Is the double transpose to sort necessary?
72 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
73 CXSPARSE_DNAME (cs_dropzeros) (N->L);
74 CXSPARSE_DNAME (cs) *D = CXSPARSE_DNAME (cs_transpose) (N->L, 1);
75 CXSPARSE_DNAME (cs_spfree) (N->L);
76 N->L = CXSPARSE_DNAME (cs_transpose) (D, 1);
77 CXSPARSE_DNAME (cs_spfree) (D);
78 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
79
80 octave_idx_type nc = N->L->n;
81 octave_idx_type nz = N->L->nzmax;
82 SparseMatrix ret (N->L->m, nc, nz);
83 for (octave_idx_type j = 0; j < nc+1; j++)
84 ret.xcidx (j) = N->L->p[j];
85 for (octave_idx_type j = 0; j < nz; j++)
86 {
87 ret.xridx (j) = N->L->i[j];
88 ret.xdata (j) = N->L->x[j];
89 }
90 return ret;
91 #else
92 return SparseMatrix ();
93 #endif
94 }
95
96 ColumnVector
97 SparseQR::SparseQR_rep::Pinv (void) const
98 {
99 #ifdef HAVE_CXSPARSE
100 ColumnVector ret(N->L->m);
101 for (octave_idx_type i = 0; i < N->L->m; i++)
102 ret.xelem(i) = S->Pinv[i];
103 return ret;
104 #else
105 return ColumnVector ();
106 #endif
107 }
108
109 ColumnVector
110 SparseQR::SparseQR_rep::P (void) const
111 {
112 #ifdef HAVE_CXSPARSE
113 ColumnVector ret(N->L->m);
114 for (octave_idx_type i = 0; i < N->L->m; i++)
115 ret.xelem(S->Pinv[i]) = i;
116 return ret;
117 #else
118 return ColumnVector ();
119 #endif
120 }
121
122 SparseMatrix
123 SparseQR::SparseQR_rep::R (const bool econ) const
124 {
125 #ifdef HAVE_CXSPARSE
126 // Drop zeros from R and sort
127 // XXX FIXME XXX Is the double transpose to sort necessary?
128 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
129 CXSPARSE_DNAME (cs_dropzeros) (N->U);
130 CXSPARSE_DNAME (cs) *D = CXSPARSE_DNAME (cs_transpose) (N->U, 1);
131 CXSPARSE_DNAME (cs_spfree) (N->U);
132 N->U = CXSPARSE_DNAME (cs_transpose) (D, 1);
133 CXSPARSE_DNAME (cs_spfree) (D);
134 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
135
136 octave_idx_type nc = N->U->n;
137 octave_idx_type nz = N->U->nzmax;
138
139 SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
140
141 for (octave_idx_type j = 0; j < nc+1; j++)
142 ret.xcidx (j) = N->U->p[j];
143 for (octave_idx_type j = 0; j < nz; j++)
144 {
145 ret.xridx (j) = N->U->i[j];
146 ret.xdata (j) = N->U->x[j];
147 }
148 return ret;
149 #else
150 return SparseMatrix ();
151 #endif
152 }
153
154 Matrix
155 SparseQR::SparseQR_rep::C (const Matrix &b) const
156 {
157 #ifdef HAVE_CXSPARSE
158 octave_idx_type b_nr = b.rows();
159 octave_idx_type b_nc = b.cols();
160 octave_idx_type nc = N->L->n;
161 octave_idx_type nr = nrows;
162 const double *bvec = b.fortran_vec();
163 Matrix ret(b_nr,b_nc);
164 double *vec = ret.fortran_vec();
165 if (nr < 1 || nc < 1 || nr != b_nr)
166 (*current_liboctave_error_handler) ("matrix dimension mismatch");
167 else
168 {
169 OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
170 for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
171 {
172 OCTAVE_QUIT;
173 volatile octave_idx_type nm = (nr < nc ? nr : nc);
174 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
175 // cast away const on bvec, with full knowledge that CSparse
176 // won't touch it
177 CXSPARSE_DNAME (cs_ipvec) (b_nr, S->Pinv, bvec + idx, buf);
178 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
179
180 for (volatile octave_idx_type i = 0; i < nm; i++)
181 {
182 OCTAVE_QUIT;
183 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
184 CXSPARSE_DNAME (cs_happly) (N->L, i, N->B[i], buf);
185 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
186 }
187 for (octave_idx_type i = 0; i < b_nr; i++)
188 vec[i+idx] = buf[i];
189 }
190 }
191 return ret;
192 #else
193 return Matrix ();
194 #endif
195 }
196
197 Matrix
198 qrsolve(const SparseMatrix&a, const Matrix &b, octave_idx_type& info)
199 {
200 #ifdef HAVE_CXSPARSE
201 octave_idx_type nr = a.rows();
202 octave_idx_type nc = a.cols();
203 octave_idx_type b_nc = b.cols();
204 octave_idx_type b_nr = b.rows();
205 const double *bvec = b.fortran_vec();
206 Matrix x;
207 info = 0;
208
209 if (nr < 1 || nc < 1 || nr != b_nr)
210 (*current_liboctave_error_handler)
211 ("matrix dimension mismatch in solution of minimum norm problem");
212 else if (nr >= nc)
213 {
214 SparseQR q (a, 2);
215 if (! q.ok ())
216 {
217 info = -1;
218 return Matrix();
219 }
220 x.resize(nc, b_nc);
221 double *vec = x.fortran_vec();
222 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
223 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
224 i++, idx+=nc, bidx+=b_nr)
225 {
226 OCTAVE_QUIT;
227 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
228 // cast away const on bvec, with full knowledge that CSparse
229 // won't touch it
230 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf);
231 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
232 for (volatile octave_idx_type j = 0; j < nc; j++)
233 {
234 OCTAVE_QUIT;
235 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
236 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
237 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
238 }
239 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
240 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf);
241 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, vec + idx);
242 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
243 }
244 }
245 else
246 {
247 SparseMatrix at = a.hermitian();
248 SparseQR q (at, 2);
249 if (! q.ok ())
250 {
251 info = -1;
252 return Matrix();
253 }
254 x.resize(nc, b_nc);
255 double *vec = x.fortran_vec();
256 OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2);
257 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
258 i++, idx+=nc, bidx+=b_nr)
259 {
260 OCTAVE_QUIT;
261 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
262 // cast away const on bvec, with full knowledge that CSparse
263 // won't touch it
264 CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, bvec + bidx, buf);
265 CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf);
266 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
267 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
268 {
269 OCTAVE_QUIT;
270 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
271 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
272 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
273 }
274 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
275 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, vec + idx);
276 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
277 }
278 }
279
280 return x;
281 #else
282 return Matrix ();
283 #endif
284 }
285
286 SparseMatrix
287 qrsolve(const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info)
288 {
289 #ifdef HAVE_CXSPARSE
290 octave_idx_type nr = a.rows();
291 octave_idx_type nc = a.cols();
292 octave_idx_type b_nr = b.rows();
293 octave_idx_type b_nc = b.cols();
294 SparseMatrix x;
295 volatile octave_idx_type ii, x_nz;
296 info = 0;
297
298 if (nr < 1 || nc < 1 || nr != b_nr)
299 (*current_liboctave_error_handler)
300 ("matrix dimension mismatch in solution of minimum norm problem");
301 else if (nr >= nc)
302 {
303 SparseQR q (a, 2);
304 if (! q.ok ())
305 {
306 info = -1;
307 return SparseMatrix();
308 }
309 x = SparseMatrix (nc, b_nc, b.nzmax());
310 x.xcidx(0) = 0;
311 x_nz = b.nzmax();
312 ii = 0;
313 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
314 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
315 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
316 {
317 OCTAVE_QUIT;
318 for (octave_idx_type j = 0; j < b_nr; j++)
319 Xx[j] = b.xelem(j,i);
320 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
321 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xx, buf);
322 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
323 for (volatile octave_idx_type j = 0; j < nc; j++)
324 {
325 OCTAVE_QUIT;
326 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
327 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
328 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
329 }
330 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
331 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf);
332 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xx);
333 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
334
335 for (octave_idx_type j = 0; j < nc; j++)
336 {
337 double tmp = Xx[j];
338 if (tmp != 0.0)
339 {
340 if (ii == x_nz)
341 {
342 // Resize the sparse matrix
343 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
344 sz = (sz > 10 ? sz : 10) + x_nz;
345 x.change_capacity (sz);
346 x_nz = sz;
347 }
348 x.xdata(ii) = tmp;
349 x.xridx(ii++) = j;
350 }
351 }
352 x.xcidx(i+1) = ii;
353 }
354 }
355 else
356 {
357 SparseMatrix at = a.hermitian();
358 SparseQR q (at, 2);
359 if (! q.ok ())
360 {
361 info = -1;
362 return SparseMatrix();
363 }
364 x = SparseMatrix (nc, b_nc, b.nzmax());
365 x.xcidx(0) = 0;
366 x_nz = b.nzmax();
367 ii = 0;
368 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
369 OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2);
370 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
371 {
372 OCTAVE_QUIT;
373 for (octave_idx_type j = 0; j < b_nr; j++)
374 Xx[j] = b.xelem(j,i);
375 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
376 CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xx, buf);
377 CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf);
378 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
379 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
380 {
381 OCTAVE_QUIT;
382 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
383 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
384 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
385 }
386 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
387 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xx);
388 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
389
390 for (octave_idx_type j = 0; j < nc; j++)
391 {
392 double tmp = Xx[j];
393 if (tmp != 0.0)
394 {
395 if (ii == x_nz)
396 {
397 // Resize the sparse matrix
398 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
399 sz = (sz > 10 ? sz : 10) + x_nz;
400 x.change_capacity (sz);
401 x_nz = sz;
402 }
403 x.xdata(ii) = tmp;
404 x.xridx(ii++) = j;
405 }
406 }
407 x.xcidx(i+1) = ii;
408 }
409 }
410
411 x.maybe_compress ();
412 return x;
413 #else
414 return SparseMatrix ();
415 #endif
416 }
417
418 ComplexMatrix
419 qrsolve(const SparseMatrix&a, const ComplexMatrix &b, octave_idx_type &info)
420 {
421 #ifdef HAVE_CXSPARSE
422 octave_idx_type nr = a.rows();
423 octave_idx_type nc = a.cols();
424 octave_idx_type b_nc = b.cols();
425 octave_idx_type b_nr = b.rows();
426 ComplexMatrix x;
427 info = 0;
428
429 if (nr < 1 || nc < 1 || nr != b_nr)
430 (*current_liboctave_error_handler)
431 ("matrix dimension mismatch in solution of minimum norm problem");
432 else if (nr >= nc)
433 {
434 SparseQR q (a, 2);
435 if (! q.ok ())
436 {
437 info = -1;
438 return ComplexMatrix();
439 }
440 x.resize(nc, b_nc);
441 Complex *vec = x.fortran_vec();
442 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
443 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
444 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
445 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
446 {
447 OCTAVE_QUIT;
448 for (octave_idx_type j = 0; j < b_nr; j++)
449 {
450 Complex c = b.xelem (j,i);
451 Xx[j] = std::real (c);
452 Xz[j] = std::imag (c);
453 }
454 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
455 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xx, buf);
456 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
457 for (volatile octave_idx_type j = 0; j < nc; j++)
458 {
459 OCTAVE_QUIT;
460 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
461 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
462 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
463 }
464 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
465 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf);
466 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xx);
467
468 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xz, buf);
469 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
470 for (volatile octave_idx_type j = 0; j < nc; j++)
471 {
472 OCTAVE_QUIT;
473 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
474 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
475 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
476 }
477 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
478 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf);
479 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xz);
480 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
481 for (octave_idx_type j = 0; j < nc; j++)
482 vec[j+idx] = Complex (Xx[j], Xz[j]);
483 }
484 }
485 else
486 {
487 SparseMatrix at = a.hermitian();
488 SparseQR q (at, 2);
489 if (! q.ok ())
490 {
491 info = -1;
492 return ComplexMatrix();
493 }
494 x.resize(nc, b_nc);
495 Complex *vec = x.fortran_vec();
496 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
497 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
498 OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2);
499 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
500 {
501 OCTAVE_QUIT;
502 for (octave_idx_type j = 0; j < b_nr; j++)
503 {
504 Complex c = b.xelem (j,i);
505 Xx[j] = std::real (c);
506 Xz[j] = std::imag (c);
507 }
508 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
509 CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xx, buf);
510 CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf);
511 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
512 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
513 {
514 OCTAVE_QUIT;
515 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
516 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
517 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
518 }
519 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
520 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xx);
521 CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xz, buf);
522 CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf);
523 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
524 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
525 {
526 OCTAVE_QUIT;
527 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
528 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
529 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
530 }
531 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
532 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xz);
533 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
534 for (octave_idx_type j = 0; j < nc; j++)
535 vec[j+idx] = Complex (Xx[j], Xz[j]);
536 }
537 }
538
539 return x;
540 #else
541 return ComplexMatrix ();
542 #endif
543 }
544
545 SparseComplexMatrix
546 qrsolve(const SparseMatrix&a, const SparseComplexMatrix &b, octave_idx_type &info)
547 {
548 #ifdef HAVE_CXSPARSE
549 octave_idx_type nr = a.rows();
550 octave_idx_type nc = a.cols();
551 octave_idx_type b_nr = b.rows();
552 octave_idx_type b_nc = b.cols();
553 SparseComplexMatrix x;
554 volatile octave_idx_type ii, x_nz;
555 info = 0;
556
557 if (nr < 1 || nc < 1 || nr != b_nr)
558 (*current_liboctave_error_handler)
559 ("matrix dimension mismatch in solution of minimum norm problem");
560 else if (nr >= nc)
561 {
562 SparseQR q (a, 2);
563 if (! q.ok ())
564 {
565 info = -1;
566 return SparseComplexMatrix();
567 }
568 x = SparseComplexMatrix (nc, b_nc, b.nzmax());
569 x.xcidx(0) = 0;
570 x_nz = b.nzmax();
571 ii = 0;
572 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
573 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
574 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
575 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
576 {
577 OCTAVE_QUIT;
578 for (octave_idx_type j = 0; j < b_nr; j++)
579 {
580 Complex c = b.xelem (j,i);
581 Xx[j] = std::real (c);
582 Xz[j] = std::imag (c);
583 }
584 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
585 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xx, buf);
586 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
587 for (volatile octave_idx_type j = 0; j < nc; j++)
588 {
589 OCTAVE_QUIT;
590 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
591 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
592 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
593 }
594 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
595 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf);
596 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xx);
597 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xz, buf);
598 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
599 for (volatile octave_idx_type j = 0; j < nc; j++)
600 {
601 OCTAVE_QUIT;
602 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
603 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
604 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
605 }
606 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
607 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf);
608 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xz);
609 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
610
611 for (octave_idx_type j = 0; j < nc; j++)
612 {
613 Complex tmp = Complex (Xx[j], Xz[j]);
614 if (tmp != 0.0)
615 {
616 if (ii == x_nz)
617 {
618 // Resize the sparse matrix
619 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
620 sz = (sz > 10 ? sz : 10) + x_nz;
621 x.change_capacity (sz);
622 x_nz = sz;
623 }
624 x.xdata(ii) = tmp;
625 x.xridx(ii++) = j;
626 }
627 }
628 x.xcidx(i+1) = ii;
629 }
630 }
631 else
632 {
633 SparseMatrix at = a.hermitian();
634 SparseQR q (at, 2);
635 if (! q.ok ())
636 {
637 info = -1;
638 return SparseComplexMatrix();
639 }
640 x = SparseComplexMatrix (nc, b_nc, b.nzmax());
641 x.xcidx(0) = 0;
642 x_nz = b.nzmax();
643 ii = 0;
644 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
645 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
646 OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2);
647 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
648 {
649 OCTAVE_QUIT;
650 for (octave_idx_type j = 0; j < b_nr; j++)
651 {
652 Complex c = b.xelem (j,i);
653 Xx[j] = std::real (c);
654 Xz[j] = std::imag (c);
655 }
656 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
657 CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xx, buf);
658 CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf);
659 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
660 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
661 {
662 OCTAVE_QUIT;
663 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
664 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
665 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
666 }
667 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
668 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xx);
669 CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xz, buf);
670 CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf);
671 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
672 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
673 {
674 OCTAVE_QUIT;
675 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
676 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf);
677 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
678 }
679 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
680 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xz);
681 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
682
683 for (octave_idx_type j = 0; j < nc; j++)
684 {
685 Complex tmp = Complex (Xx[j], Xz[j]);
686 if (tmp != 0.0)
687 {
688 if (ii == x_nz)
689 {
690 // Resize the sparse matrix
691 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
692 sz = (sz > 10 ? sz : 10) + x_nz;
693 x.change_capacity (sz);
694 x_nz = sz;
695 }
696 x.xdata(ii) = tmp;
697 x.xridx(ii++) = j;
698 }
699 }
700 x.xcidx(i+1) = ii;
701 }
702 }
703
704 x.maybe_compress ();
705 return x;
706 #else
707 return SparseComplexMatrix ();
708 #endif
709 }
710
711 /*
712 ;;; Local Variables: ***
713 ;;; mode: C++ ***
714 ;;; End: ***
715 */