Mercurial > octave-nkf

comparison liboctave/SparseQR.cc @ 5797:11fcab4c461d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-05-09 06:15:17 by dbateman]
author dbateman
date Tue, 09 May 2006 06:15:18 +0000
parents eb90c83b4f91
children be176b7e110a
comparison
equal deleted inserted replaced
5796:20f4bd627a74 5797:11fcab4c461d
212 } 212 }
213 213
214 Matrix 214 Matrix
215 qrsolve(const SparseMatrix&a, const Matrix &b, octave_idx_type& info) 215 qrsolve(const SparseMatrix&a, const Matrix &b, octave_idx_type& info)
216 { 216 {
217 info = -1;
217 #ifdef HAVE_CXSPARSE 218 #ifdef HAVE_CXSPARSE
218 octave_idx_type nr = a.rows(); 219 octave_idx_type nr = a.rows();
219 octave_idx_type nc = a.cols(); 220 octave_idx_type nc = a.cols();
220 octave_idx_type b_nc = b.cols(); 221 octave_idx_type b_nc = b.cols();
221 octave_idx_type b_nr = b.rows(); 222 octave_idx_type b_nr = b.rows();
222 const double *bvec = b.fortran_vec(); 223 const double *bvec = b.fortran_vec();
223 Matrix x; 224 Matrix x;
224 info = 0;
225 225
226 if (nr < 1 || nc < 1 || nr != b_nr) 226 if (nr < 1 || nc < 1 || nr != b_nr)
227 (*current_liboctave_error_handler) 227 (*current_liboctave_error_handler)
228 ("matrix dimension mismatch in solution of minimum norm problem"); 228 ("matrix dimension mismatch in solution of minimum norm problem");
229 else if (nr >= nc) 229 else if (nr >= nc)
230 { 230 {
231 SparseQR q (a, 3); 231 SparseQR q (a, 3);
232 if (! q.ok ()) 232 if (! q.ok ())
233 { 233 return Matrix();
234 info = -1;
235 return Matrix();
236 }
237 x.resize(nc, b_nc); 234 x.resize(nc, b_nc);
238 double *vec = x.fortran_vec(); 235 double *vec = x.fortran_vec();
239 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 236 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
240 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 237 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
241 i++, idx+=nc, bidx+=b_nr) 238 i++, idx+=nc, bidx+=b_nr)
264 #else 261 #else
265 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx); 262 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx);
266 #endif 263 #endif
267 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 264 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
268 } 265 }
266 info = 0;
269 } 267 }
270 else 268 else
271 { 269 {
272 SparseMatrix at = a.hermitian(); 270 SparseMatrix at = a.hermitian();
273 SparseQR q (at, 3); 271 SparseQR q (at, 3);
274 if (! q.ok ()) 272 if (! q.ok ())
275 { 273 return Matrix();
276 info = -1;
277 return Matrix();
278 }
279 x.resize(nc, b_nc); 274 x.resize(nc, b_nc);
280 double *vec = x.fortran_vec(); 275 double *vec = x.fortran_vec();
281 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 276 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
282 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 277 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
283 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 278 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
307 #else 302 #else
308 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx); 303 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx);
309 #endif 304 #endif
310 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 305 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
311 } 306 }
307 info = 0;
312 } 308 }
313 309
314 return x; 310 return x;
315 #else 311 #else
316 return Matrix (); 312 return Matrix ();
318 } 314 }
319 315
320 SparseMatrix 316 SparseMatrix
321 qrsolve(const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info) 317 qrsolve(const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info)
322 { 318 {
319 info = -1;
323 #ifdef HAVE_CXSPARSE 320 #ifdef HAVE_CXSPARSE
324 octave_idx_type nr = a.rows(); 321 octave_idx_type nr = a.rows();
325 octave_idx_type nc = a.cols(); 322 octave_idx_type nc = a.cols();
326 octave_idx_type b_nr = b.rows(); 323 octave_idx_type b_nr = b.rows();
327 octave_idx_type b_nc = b.cols(); 324 octave_idx_type b_nc = b.cols();
328 SparseMatrix x; 325 SparseMatrix x;
329 volatile octave_idx_type ii, x_nz; 326 volatile octave_idx_type ii, x_nz;
330 info = 0;
331 327
332 if (nr < 1 || nc < 1 || nr != b_nr) 328 if (nr < 1 || nc < 1 || nr != b_nr)
333 (*current_liboctave_error_handler) 329 (*current_liboctave_error_handler)
334 ("matrix dimension mismatch in solution of minimum norm problem"); 330 ("matrix dimension mismatch in solution of minimum norm problem");
335 else if (nr >= nc) 331 else if (nr >= nc)
336 { 332 {
337 SparseQR q (a, 3); 333 SparseQR q (a, 3);
338 if (! q.ok ()) 334 if (! q.ok ())
339 { 335 return SparseMatrix();
340 info = -1;
341 return SparseMatrix();
342 }
343 x = SparseMatrix (nc, b_nc, b.nzmax()); 336 x = SparseMatrix (nc, b_nc, b.nzmax());
344 x.xcidx(0) = 0; 337 x.xcidx(0) = 0;
345 x_nz = b.nzmax(); 338 x_nz = b.nzmax();
346 ii = 0; 339 ii = 0;
347 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 340 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
393 x.xridx(ii++) = j; 386 x.xridx(ii++) = j;
394 } 387 }
395 } 388 }
396 x.xcidx(i+1) = ii; 389 x.xcidx(i+1) = ii;
397 } 390 }
391 info = 0;
398 } 392 }
399 else 393 else
400 { 394 {
401 SparseMatrix at = a.hermitian(); 395 SparseMatrix at = a.hermitian();
402 SparseQR q (at, 3); 396 SparseQR q (at, 3);
403 if (! q.ok ()) 397 if (! q.ok ())
404 { 398 return SparseMatrix();
405 info = -1;
406 return SparseMatrix();
407 }
408 x = SparseMatrix (nc, b_nc, b.nzmax()); 399 x = SparseMatrix (nc, b_nc, b.nzmax());
409 x.xcidx(0) = 0; 400 x.xcidx(0) = 0;
410 x_nz = b.nzmax(); 401 x_nz = b.nzmax();
411 ii = 0; 402 ii = 0;
412 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 403 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
459 x.xridx(ii++) = j; 450 x.xridx(ii++) = j;
460 } 451 }
461 } 452 }
462 x.xcidx(i+1) = ii; 453 x.xcidx(i+1) = ii;
463 } 454 }
455 info = 0;
464 } 456 }
465 457
466 x.maybe_compress (); 458 x.maybe_compress ();
467 return x; 459 return x;
468 #else 460 #else
471 } 463 }
472 464
473 ComplexMatrix 465 ComplexMatrix
474 qrsolve(const SparseMatrix&a, const ComplexMatrix &b, octave_idx_type &info) 466 qrsolve(const SparseMatrix&a, const ComplexMatrix &b, octave_idx_type &info)
475 { 467 {
468 info = -1;
476 #ifdef HAVE_CXSPARSE 469 #ifdef HAVE_CXSPARSE
477 octave_idx_type nr = a.rows(); 470 octave_idx_type nr = a.rows();
478 octave_idx_type nc = a.cols(); 471 octave_idx_type nc = a.cols();
479 octave_idx_type b_nc = b.cols(); 472 octave_idx_type b_nc = b.cols();
480 octave_idx_type b_nr = b.rows(); 473 octave_idx_type b_nr = b.rows();
481 ComplexMatrix x; 474 ComplexMatrix x;
482 info = 0;
483 475
484 if (nr < 1 || nc < 1 || nr != b_nr) 476 if (nr < 1 || nc < 1 || nr != b_nr)
485 (*current_liboctave_error_handler) 477 (*current_liboctave_error_handler)
486 ("matrix dimension mismatch in solution of minimum norm problem"); 478 ("matrix dimension mismatch in solution of minimum norm problem");
487 else if (nr >= nc) 479 else if (nr >= nc)
488 { 480 {
489 SparseQR q (a, 3); 481 SparseQR q (a, 3);
490 if (! q.ok ()) 482 if (! q.ok ())
491 { 483 return ComplexMatrix();
492 info = -1;
493 return ComplexMatrix();
494 }
495 x.resize(nc, b_nc); 484 x.resize(nc, b_nc);
496 Complex *vec = x.fortran_vec(); 485 Complex *vec = x.fortran_vec();
497 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 486 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
498 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 487 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
499 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 488 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
553 #endif 542 #endif
554 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 543 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
555 for (octave_idx_type j = 0; j < nc; j++) 544 for (octave_idx_type j = 0; j < nc; j++)
556 vec[j+idx] = Complex (Xx[j], Xz[j]); 545 vec[j+idx] = Complex (Xx[j], Xz[j]);
557 } 546 }
547 info = 0;
558 } 548 }
559 else 549 else
560 { 550 {
561 SparseMatrix at = a.hermitian(); 551 SparseMatrix at = a.hermitian();
562 SparseQR q (at, 3); 552 SparseQR q (at, 3);
563 if (! q.ok ()) 553 if (! q.ok ())
564 { 554 return ComplexMatrix();
565 info = -1;
566 return ComplexMatrix();
567 }
568 x.resize(nc, b_nc); 555 x.resize(nc, b_nc);
569 Complex *vec = x.fortran_vec(); 556 Complex *vec = x.fortran_vec();
570 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 557 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
571 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 558 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
572 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 559 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
629 #endif 616 #endif
630 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 617 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
631 for (octave_idx_type j = 0; j < nc; j++) 618 for (octave_idx_type j = 0; j < nc; j++)
632 vec[j+idx] = Complex (Xx[j], Xz[j]); 619 vec[j+idx] = Complex (Xx[j], Xz[j]);
633 } 620 }
621 info = 0;
634 } 622 }
635 623
636 return x; 624 return x;
637 #else 625 #else
638 return ComplexMatrix (); 626 return ComplexMatrix ();
640 } 628 }
641 629
642 SparseComplexMatrix 630 SparseComplexMatrix
643 qrsolve(const SparseMatrix&a, const SparseComplexMatrix &b, octave_idx_type &info) 631 qrsolve(const SparseMatrix&a, const SparseComplexMatrix &b, octave_idx_type &info)
644 { 632 {
633 info = -1;
645 #ifdef HAVE_CXSPARSE 634 #ifdef HAVE_CXSPARSE
646 octave_idx_type nr = a.rows(); 635 octave_idx_type nr = a.rows();
647 octave_idx_type nc = a.cols(); 636 octave_idx_type nc = a.cols();
648 octave_idx_type b_nr = b.rows(); 637 octave_idx_type b_nr = b.rows();
649 octave_idx_type b_nc = b.cols(); 638 octave_idx_type b_nc = b.cols();
650 SparseComplexMatrix x; 639 SparseComplexMatrix x;
651 volatile octave_idx_type ii, x_nz; 640 volatile octave_idx_type ii, x_nz;
652 info = 0;
653 641
654 if (nr < 1 || nc < 1 || nr != b_nr) 642 if (nr < 1 || nc < 1 || nr != b_nr)
655 (*current_liboctave_error_handler) 643 (*current_liboctave_error_handler)
656 ("matrix dimension mismatch in solution of minimum norm problem"); 644 ("matrix dimension mismatch in solution of minimum norm problem");
657 else if (nr >= nc) 645 else if (nr >= nc)
658 { 646 {
659 SparseQR q (a, 3); 647 SparseQR q (a, 3);
660 if (! q.ok ()) 648 if (! q.ok ())
661 { 649 return SparseComplexMatrix();
662 info = -1;
663 return SparseComplexMatrix();
664 }
665 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); 650 x = SparseComplexMatrix (nc, b_nc, b.nzmax());
666 x.xcidx(0) = 0; 651 x.xcidx(0) = 0;
667 x_nz = b.nzmax(); 652 x_nz = b.nzmax();
668 ii = 0; 653 ii = 0;
669 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 654 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
744 x.xridx(ii++) = j; 729 x.xridx(ii++) = j;
745 } 730 }
746 } 731 }
747 x.xcidx(i+1) = ii; 732 x.xcidx(i+1) = ii;
748 } 733 }
734 info = 0;
749 } 735 }
750 else 736 else
751 { 737 {
752 SparseMatrix at = a.hermitian(); 738 SparseMatrix at = a.hermitian();
753 SparseQR q (at, 3); 739 SparseQR q (at, 3);
754 if (! q.ok ()) 740 if (! q.ok ())
755 { 741 return SparseComplexMatrix();
756 info = -1;
757 return SparseComplexMatrix();
758 }
759 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); 742 x = SparseComplexMatrix (nc, b_nc, b.nzmax());
760 x.xcidx(0) = 0; 743 x.xcidx(0) = 0;
761 x_nz = b.nzmax(); 744 x_nz = b.nzmax();
762 ii = 0; 745 ii = 0;
763 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 746 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
839 x.xridx(ii++) = j; 822 x.xridx(ii++) = j;
840 } 823 }
841 } 824 }
842 x.xcidx(i+1) = ii; 825 x.xcidx(i+1) = ii;
843 } 826 }
827 info = 0;
844 } 828 }
845 829
846 x.maybe_compress (); 830 x.maybe_compress ();
847 return x; 831 return x;
848 #else 832 #else