Mercurial > hg > octave-thorsten

comparison liboctave/SparseQR.cc @ 5681:233d98d95659

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-03-16 17:48:55 by dbateman]
author dbateman
date Thu, 16 Mar 2006 17:48:56 +0000
parents 69a4f320d95a
children ace8d8d26933
comparison
equal deleted inserted replaced
5680:cc6a965ae4ca 5681:233d98d95659
168 { 168 {
169 OCTAVE_LOCAL_BUFFER (double, buf, S->m2); 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) 170 for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
171 { 171 {
172 OCTAVE_QUIT; 172 OCTAVE_QUIT;
173 for (octave_idx_type i = nr; i < S->m2; i++)
174 buf[i] = 0.;
173 volatile octave_idx_type nm = (nr < nc ? nr : nc); 175 volatile octave_idx_type nm = (nr < nc ? nr : nc);
174 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 176 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
175 CXSPARSE_DNAME (_ipvec) (b_nr, S->Pinv, bvec + idx, buf); 177 CXSPARSE_DNAME (_ipvec) (b_nr, S->Pinv, bvec + idx, buf);
176 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 178 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
177 179
220 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 222 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
221 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 223 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
222 i++, idx+=nc, bidx+=b_nr) 224 i++, idx+=nc, bidx+=b_nr)
223 { 225 {
224 OCTAVE_QUIT; 226 OCTAVE_QUIT;
227 for (octave_idx_type j = nr; j < q.S()->m2; j++)
228 buf[j] = 0.;
225 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 229 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
226 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf); 230 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf);
227 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 231 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
228 for (volatile octave_idx_type j = 0; j < nc; j++) 232 for (volatile octave_idx_type j = 0; j < nc; j++)
229 { 233 {
247 info = -1; 251 info = -1;
248 return Matrix(); 252 return Matrix();
249 } 253 }
250 x.resize(nc, b_nc); 254 x.resize(nc, b_nc);
251 double *vec = x.fortran_vec(); 255 double *vec = x.fortran_vec();
252 OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2); 256 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
257 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
253 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 258 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
254 i++, idx+=nc, bidx+=b_nr) 259 i++, idx+=nc, bidx+=b_nr)
255 { 260 {
256 OCTAVE_QUIT; 261 OCTAVE_QUIT;
262 for (octave_idx_type j = nr; j < nbuf; j++)
263 buf[j] = 0.;
257 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 264 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
258 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf); 265 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf);
259 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 266 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
260 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 267 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
261 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 268 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
309 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 316 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
310 { 317 {
311 OCTAVE_QUIT; 318 OCTAVE_QUIT;
312 for (octave_idx_type j = 0; j < b_nr; j++) 319 for (octave_idx_type j = 0; j < b_nr; j++)
313 Xx[j] = b.xelem(j,i); 320 Xx[j] = b.xelem(j,i);
321 for (octave_idx_type j = nr; j < q.S()->m2; j++)
322 buf[j] = 0.;
314 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 323 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
315 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); 324 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
316 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 325 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
317 for (volatile octave_idx_type j = 0; j < nc; j++) 326 for (volatile octave_idx_type j = 0; j < nc; j++)
318 { 327 {
357 } 366 }
358 x = SparseMatrix (nc, b_nc, b.nzmax()); 367 x = SparseMatrix (nc, b_nc, b.nzmax());
359 x.xcidx(0) = 0; 368 x.xcidx(0) = 0;
360 x_nz = b.nzmax(); 369 x_nz = b.nzmax();
361 ii = 0; 370 ii = 0;
371 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
362 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 372 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
363 OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2); 373 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
364 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 374 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
365 { 375 {
366 OCTAVE_QUIT; 376 OCTAVE_QUIT;
367 for (octave_idx_type j = 0; j < b_nr; j++) 377 for (octave_idx_type j = 0; j < b_nr; j++)
368 Xx[j] = b.xelem(j,i); 378 Xx[j] = b.xelem(j,i);
379 for (octave_idx_type j = nr; j < nbuf; j++)
380 buf[j] = 0.;
369 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 381 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
370 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); 382 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
371 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 383 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
372 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 384 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
373 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 385 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
443 { 455 {
444 Complex c = b.xelem (j,i); 456 Complex c = b.xelem (j,i);
445 Xx[j] = std::real (c); 457 Xx[j] = std::real (c);
446 Xz[j] = std::imag (c); 458 Xz[j] = std::imag (c);
447 } 459 }
460 for (octave_idx_type j = nr; j < q.S()->m2; j++)
461 buf[j] = 0.;
448 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 462 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
449 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); 463 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
450 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 464 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
451 for (volatile octave_idx_type j = 0; j < nc; j++) 465 for (volatile octave_idx_type j = 0; j < nc; j++)
452 { 466 {
457 } 471 }
458 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 472 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
459 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 473 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
460 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); 474 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx);
461 475
476 for (octave_idx_type j = nr; j < q.S()->m2; j++)
477 buf[j] = 0.;
462 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf); 478 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf);
463 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 479 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
464 for (volatile octave_idx_type j = 0; j < nc; j++) 480 for (volatile octave_idx_type j = 0; j < nc; j++)
465 { 481 {
466 OCTAVE_QUIT; 482 OCTAVE_QUIT;
485 info = -1; 501 info = -1;
486 return ComplexMatrix(); 502 return ComplexMatrix();
487 } 503 }
488 x.resize(nc, b_nc); 504 x.resize(nc, b_nc);
489 Complex *vec = x.fortran_vec(); 505 Complex *vec = x.fortran_vec();
506 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
490 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 507 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
491 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 508 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
492 OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2); 509 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
493 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 510 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
494 { 511 {
495 OCTAVE_QUIT; 512 OCTAVE_QUIT;
496 for (octave_idx_type j = 0; j < b_nr; j++) 513 for (octave_idx_type j = 0; j < b_nr; j++)
497 { 514 {
498 Complex c = b.xelem (j,i); 515 Complex c = b.xelem (j,i);
499 Xx[j] = std::real (c); 516 Xx[j] = std::real (c);
500 Xz[j] = std::imag (c); 517 Xz[j] = std::imag (c);
501 } 518 }
519 for (octave_idx_type j = nr; j < nbuf; j++)
520 buf[j] = 0.;
502 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 521 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
503 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); 522 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
504 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 523 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
505 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 524 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
506 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 525 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
510 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 529 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
511 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 530 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
512 } 531 }
513 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 532 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
514 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); 533 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx);
534 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
535 for (octave_idx_type j = nr; j < nbuf; j++)
536 buf[j] = 0.;
537 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
515 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf); 538 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf);
516 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 539 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
517 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 540 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
518 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 541 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
519 { 542 {
573 { 596 {
574 Complex c = b.xelem (j,i); 597 Complex c = b.xelem (j,i);
575 Xx[j] = std::real (c); 598 Xx[j] = std::real (c);
576 Xz[j] = std::imag (c); 599 Xz[j] = std::imag (c);
577 } 600 }
601 for (octave_idx_type j = nr; j < q.S()->m2; j++)
602 buf[j] = 0.;
578 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 603 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
579 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); 604 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
580 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 605 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
581 for (volatile octave_idx_type j = 0; j < nc; j++) 606 for (volatile octave_idx_type j = 0; j < nc; j++)
582 { 607 {
586 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 611 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
587 } 612 }
588 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 613 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
589 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 614 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
590 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); 615 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx);
616 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
617 for (octave_idx_type j = nr; j < q.S()->m2; j++)
618 buf[j] = 0.;
619 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
591 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf); 620 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf);
592 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 621 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
593 for (volatile octave_idx_type j = 0; j < nc; j++) 622 for (volatile octave_idx_type j = 0; j < nc; j++)
594 { 623 {
595 OCTAVE_QUIT; 624 OCTAVE_QUIT;
633 } 662 }
634 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); 663 x = SparseComplexMatrix (nc, b_nc, b.nzmax());
635 x.xcidx(0) = 0; 664 x.xcidx(0) = 0;
636 x_nz = b.nzmax(); 665 x_nz = b.nzmax();
637 ii = 0; 666 ii = 0;
667 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
638 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 668 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
639 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 669 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
640 OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2); 670 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
641 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 671 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
642 { 672 {
643 OCTAVE_QUIT; 673 OCTAVE_QUIT;
644 for (octave_idx_type j = 0; j < b_nr; j++) 674 for (octave_idx_type j = 0; j < b_nr; j++)
645 { 675 {
646 Complex c = b.xelem (j,i); 676 Complex c = b.xelem (j,i);
647 Xx[j] = std::real (c); 677 Xx[j] = std::real (c);
648 Xz[j] = std::imag (c); 678 Xz[j] = std::imag (c);
649 } 679 }
680 for (octave_idx_type j = nr; j < nbuf; j++)
681 buf[j] = 0.;
650 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 682 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
651 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); 683 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
652 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 684 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
653 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 685 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
654 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 686 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
658 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 690 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
659 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 691 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
660 } 692 }
661 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 693 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
662 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); 694 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx);
695 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
696 for (octave_idx_type j = nr; j < nbuf; j++)
697 buf[j] = 0.;
698 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
663 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf); 699 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf);
664 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 700 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
665 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 701 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
666 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 702 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
667 { 703 {