Mercurial > octave

comparison liboctave/numeric/SparseQR.cc @ 21152:8ad3907b8fad

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
require CXSparse 2.2 or later * oct-sparse.h: Require CXSparse 2.2 or later. * dmperm.cc, SparseCmplxQR.cc, SparseQR.cc, sparse-dmsolve.cc: Assume CXSparse 2.2 or later.
author John W. Eaton <jwe@octave.org>
date Fri, 29 Jan 2016 14:55:20 -0500
parents 7cac4e7458f2
children 9c61ab1f7588
comparison
equal deleted inserted replaced
21151:bfd5e48c41a1 21152:8ad3907b8fad
46 A.p = const_cast<octave_idx_type *>(a.cidx ()); 46 A.p = const_cast<octave_idx_type *>(a.cidx ());
47 A.i = const_cast<octave_idx_type *>(a.ridx ()); 47 A.i = const_cast<octave_idx_type *>(a.ridx ());
48 A.x = const_cast<double *>(a.data ()); 48 A.x = const_cast<double *>(a.data ());
49 A.nz = -1; 49 A.nz = -1;
50 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 50 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
51 #if defined (CS_VER) && (CS_VER >= 2)
52 S = CXSPARSE_DNAME (_sqr) (order, &A, 1); 51 S = CXSPARSE_DNAME (_sqr) (order, &A, 1);
53 #else
54 S = CXSPARSE_DNAME (_sqr) (&A, order - 1, 1);
55 #endif
56 52
57 N = CXSPARSE_DNAME (_qr) (&A, S); 53 N = CXSPARSE_DNAME (_qr) (&A, S);
58 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 54 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
59 if (! N) 55 if (! N)
60 (*current_liboctave_error_handler) 56 (*current_liboctave_error_handler)
110 SparseQR::SparseQR_rep::Pinv (void) const 106 SparseQR::SparseQR_rep::Pinv (void) const
111 { 107 {
112 #ifdef HAVE_CXSPARSE 108 #ifdef HAVE_CXSPARSE
113 ColumnVector ret(N->L->m); 109 ColumnVector ret(N->L->m);
114 for (octave_idx_type i = 0; i < N->L->m; i++) 110 for (octave_idx_type i = 0; i < N->L->m; i++)
115 #if defined (CS_VER) && (CS_VER >= 2)
116 ret.xelem (i) = S->pinv[i]; 111 ret.xelem (i) = S->pinv[i];
117 #else
118 ret.xelem (i) = S->Pinv[i];
119 #endif
120 return ret; 112 return ret;
121 #else 113 #else
122 return ColumnVector (); 114 return ColumnVector ();
123 #endif 115 #endif
124 } 116 }
127 SparseQR::SparseQR_rep::P (void) const 119 SparseQR::SparseQR_rep::P (void) const
128 { 120 {
129 #ifdef HAVE_CXSPARSE 121 #ifdef HAVE_CXSPARSE
130 ColumnVector ret(N->L->m); 122 ColumnVector ret(N->L->m);
131 for (octave_idx_type i = 0; i < N->L->m; i++) 123 for (octave_idx_type i = 0; i < N->L->m; i++)
132 #if defined (CS_VER) && (CS_VER >= 2)
133 ret.xelem (S->pinv[i]) = i; 124 ret.xelem (S->pinv[i]) = i;
134 #else
135 ret.xelem (S->Pinv[i]) = i;
136 #endif
137 return ret; 125 return ret;
138 #else 126 #else
139 return ColumnVector (); 127 return ColumnVector ();
140 #endif 128 #endif
141 } 129 }
196 octave_quit (); 184 octave_quit ();
197 for (octave_idx_type i = nr; i < S->m2; i++) 185 for (octave_idx_type i = nr; i < S->m2; i++)
198 buf[i] = 0.; 186 buf[i] = 0.;
199 volatile octave_idx_type nm = (nr < nc ? nr : nc); 187 volatile octave_idx_type nm = (nr < nc ? nr : nc);
200 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 188 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
201 #if defined (CS_VER) && (CS_VER >= 2)
202 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr); 189 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr);
203 #else
204 CXSPARSE_DNAME (_ipvec) (b_nr, S->Pinv, bvec + idx, buf);
205 #endif
206 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 190 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
207 191
208 for (volatile octave_idx_type i = 0; i < nm; i++) 192 for (volatile octave_idx_type i = 0; i < nm; i++)
209 { 193 {
210 octave_quit (); 194 octave_quit ();
247 bvec[j] = 1.0; 231 bvec[j] = 1.0;
248 for (octave_idx_type i = nr; i < S->m2; i++) 232 for (octave_idx_type i = nr; i < S->m2; i++)
249 buf[i] = 0.; 233 buf[i] = 0.;
250 volatile octave_idx_type nm = (nr < nc ? nr : nc); 234 volatile octave_idx_type nm = (nr < nc ? nr : nc);
251 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 235 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
252 #if defined (CS_VER) && (CS_VER >= 2)
253 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr); 236 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr);
254 #else
255 CXSPARSE_DNAME (_ipvec) (nr, S->Pinv, bvec, buf);
256 #endif
257 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 237 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
258 238
259 for (volatile octave_idx_type i = 0; i < nm; i++) 239 for (volatile octave_idx_type i = 0; i < nm; i++)
260 { 240 {
261 octave_quit (); 241 octave_quit ();
305 { 285 {
306 octave_quit (); 286 octave_quit ();
307 for (octave_idx_type j = nr; j < q.S ()->m2; j++) 287 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
308 buf[j] = 0.; 288 buf[j] = 0.;
309 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 289 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
310 #if defined (CS_VER) && (CS_VER >= 2)
311 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, bvec + bidx, buf, nr); 290 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, bvec + bidx, buf, nr);
312 #else
313 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, bvec + bidx, buf);
314 #endif
315 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 291 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
316 for (volatile octave_idx_type j = 0; j < nc; j++) 292 for (volatile octave_idx_type j = 0; j < nc; j++)
317 { 293 {
318 octave_quit (); 294 octave_quit ();
319 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 295 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
320 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 296 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
321 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 297 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
322 } 298 }
323 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 299 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
324 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); 300 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
325 #if defined (CS_VER) && (CS_VER >= 2)
326 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc); 301 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc);
327 #else
328 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, vec + idx);
329 #endif
330 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 302 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
331 } 303 }
332 info = 0; 304 info = 0;
333 } 305 }
334 else 306 else
346 { 318 {
347 octave_quit (); 319 octave_quit ();
348 for (octave_idx_type j = nr; j < nbuf; j++) 320 for (octave_idx_type j = nr; j < nbuf; j++)
349 buf[j] = 0.; 321 buf[j] = 0.;
350 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 322 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
351 #if defined (CS_VER) && (CS_VER >= 2)
352 CXSPARSE_DNAME (_pvec) (q.S ()->q, bvec + bidx, buf, nr); 323 CXSPARSE_DNAME (_pvec) (q.S ()->q, bvec + bidx, buf, nr);
353 #else
354 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, bvec + bidx, buf);
355 #endif
356 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); 324 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
357 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 325 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
358 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 326 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
359 { 327 {
360 octave_quit (); 328 octave_quit ();
361 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 329 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
362 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 330 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
363 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 331 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
364 } 332 }
365 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 333 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
366 #if defined (CS_VER) && (CS_VER >= 2)
367 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc); 334 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc);
368 #else
369 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, vec + idx);
370 #endif
371 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 335 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
372 } 336 }
373 info = 0; 337 info = 0;
374 } 338 }
375 339
414 for (octave_idx_type j = 0; j < b_nr; j++) 378 for (octave_idx_type j = 0; j < b_nr; j++)
415 Xx[j] = b.xelem (j,i); 379 Xx[j] = b.xelem (j,i);
416 for (octave_idx_type j = nr; j < q.S ()->m2; j++) 380 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
417 buf[j] = 0.; 381 buf[j] = 0.;
418 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 382 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
419 #if defined (CS_VER) && (CS_VER >= 2)
420 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr); 383 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
421 #else
422 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xx, buf);
423 #endif
424 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 384 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
425 for (volatile octave_idx_type j = 0; j < nc; j++) 385 for (volatile octave_idx_type j = 0; j < nc; j++)
426 { 386 {
427 octave_quit (); 387 octave_quit ();
428 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 388 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
429 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 389 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
430 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 390 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
431 } 391 }
432 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 392 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
433 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); 393 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
434 #if defined (CS_VER) && (CS_VER >= 2)
435 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc); 394 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
436 #else
437 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xx);
438 #endif
439 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 395 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
440 396
441 for (octave_idx_type j = 0; j < nc; j++) 397 for (octave_idx_type j = 0; j < nc; j++)
442 { 398 {
443 double tmp = Xx[j]; 399 double tmp = Xx[j];
478 for (octave_idx_type j = 0; j < b_nr; j++) 434 for (octave_idx_type j = 0; j < b_nr; j++)
479 Xx[j] = b.xelem (j,i); 435 Xx[j] = b.xelem (j,i);
480 for (octave_idx_type j = nr; j < nbuf; j++) 436 for (octave_idx_type j = nr; j < nbuf; j++)
481 buf[j] = 0.; 437 buf[j] = 0.;
482 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 438 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
483 #if defined (CS_VER) && (CS_VER >= 2)
484 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr); 439 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
485 #else
486 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xx, buf);
487 #endif
488 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); 440 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
489 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 441 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
490 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 442 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
491 { 443 {
492 octave_quit (); 444 octave_quit ();
493 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 445 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
494 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 446 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
495 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 447 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
496 } 448 }
497 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 449 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
498 #if defined (CS_VER) && (CS_VER >= 2)
499 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc); 450 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
500 #else
501 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xx);
502 #endif
503 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 451 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
504 452
505 for (octave_idx_type j = 0; j < nc; j++) 453 for (octave_idx_type j = 0; j < nc; j++)
506 { 454 {
507 double tmp = Xx[j]; 455 double tmp = Xx[j];
568 Xz[j] = std::imag (c); 516 Xz[j] = std::imag (c);
569 } 517 }
570 for (octave_idx_type j = nr; j < q.S ()->m2; j++) 518 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
571 buf[j] = 0.; 519 buf[j] = 0.;
572 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 520 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
573 #if defined (CS_VER) && (CS_VER >= 2)
574 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr); 521 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
575 #else
576 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xx, buf);
577 #endif
578 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 522 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
579 for (volatile octave_idx_type j = 0; j < nc; j++) 523 for (volatile octave_idx_type j = 0; j < nc; j++)
580 { 524 {
581 octave_quit (); 525 octave_quit ();
582 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 526 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
583 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 527 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
584 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 528 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
585 } 529 }
586 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 530 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
587 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); 531 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
588 #if defined (CS_VER) && (CS_VER >= 2)
589 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc); 532 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
590 #else
591 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xx);
592 #endif
593 for (octave_idx_type j = nr; j < q.S ()->m2; j++) 533 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
594 buf[j] = 0.; 534 buf[j] = 0.;
595 #if defined (CS_VER) && (CS_VER >= 2)
596 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr); 535 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr);
597 #else
598 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xz, buf);
599 #endif
600 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 536 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
601 for (volatile octave_idx_type j = 0; j < nc; j++) 537 for (volatile octave_idx_type j = 0; j < nc; j++)
602 { 538 {
603 octave_quit (); 539 octave_quit ();
604 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 540 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
605 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 541 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
606 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 542 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
607 } 543 }
608 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 544 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
609 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); 545 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
610 #if defined (CS_VER) && (CS_VER >= 2)
611 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc); 546 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc);
612 #else
613 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xz);
614 #endif
615 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 547 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
616 for (octave_idx_type j = 0; j < nc; j++) 548 for (octave_idx_type j = 0; j < nc; j++)
617 vec[j+idx] = Complex (Xx[j], Xz[j]); 549 vec[j+idx] = Complex (Xx[j], Xz[j]);
618 } 550 }
619 info = 0; 551 info = 0;
640 Xz[j] = std::imag (c); 572 Xz[j] = std::imag (c);
641 } 573 }
642 for (octave_idx_type j = nr; j < nbuf; j++) 574 for (octave_idx_type j = nr; j < nbuf; j++)
643 buf[j] = 0.; 575 buf[j] = 0.;
644 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 576 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
645 #if defined (CS_VER) && (CS_VER >= 2)
646 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr); 577 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
647 #else
648 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xx, buf);
649 #endif
650 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); 578 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
651 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 579 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
652 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 580 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
653 { 581 {
654 octave_quit (); 582 octave_quit ();
655 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 583 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
656 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 584 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
657 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 585 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
658 } 586 }
659 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 587 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
660 #if defined (CS_VER) && (CS_VER >= 2)
661 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc); 588 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
662 #else
663 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xx);
664 #endif
665 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 589 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
666 for (octave_idx_type j = nr; j < nbuf; j++) 590 for (octave_idx_type j = nr; j < nbuf; j++)
667 buf[j] = 0.; 591 buf[j] = 0.;
668 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 592 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
669 #if defined (CS_VER) && (CS_VER >= 2)
670 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr); 593 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr);
671 #else
672 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xz, buf);
673 #endif
674 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); 594 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
675 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 595 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
676 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 596 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
677 { 597 {
678 octave_quit (); 598 octave_quit ();
679 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 599 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
680 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 600 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
681 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 601 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
682 } 602 }
683 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 603 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
684 #if defined (CS_VER) && (CS_VER >= 2)
685 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc); 604 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc);
686 #else
687 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xz);
688 #endif
689 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 605 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
690 for (octave_idx_type j = 0; j < nc; j++) 606 for (octave_idx_type j = 0; j < nc; j++)
691 vec[j+idx] = Complex (Xx[j], Xz[j]); 607 vec[j+idx] = Complex (Xx[j], Xz[j]);
692 } 608 }
693 info = 0; 609 info = 0;
740 Xz[j] = std::imag (c); 656 Xz[j] = std::imag (c);
741 } 657 }
742 for (octave_idx_type j = nr; j < q.S ()->m2; j++) 658 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
743 buf[j] = 0.; 659 buf[j] = 0.;
744 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 660 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
745 #if defined (CS_VER) && (CS_VER >= 2)
746 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr); 661 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
747 #else
748 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xx, buf);
749 #endif
750 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 662 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
751 for (volatile octave_idx_type j = 0; j < nc; j++) 663 for (volatile octave_idx_type j = 0; j < nc; j++)
752 { 664 {
753 octave_quit (); 665 octave_quit ();
754 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 666 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
755 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 667 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
756 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 668 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
757 } 669 }
758 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 670 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
759 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); 671 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
760 #if defined (CS_VER) && (CS_VER >= 2)
761 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc); 672 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
762 #else
763 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xx);
764 #endif
765 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 673 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
766 for (octave_idx_type j = nr; j < q.S ()->m2; j++) 674 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
767 buf[j] = 0.; 675 buf[j] = 0.;
768 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 676 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
769 #if defined (CS_VER) && (CS_VER >= 2)
770 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr); 677 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr);
771 #else
772 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xz, buf);
773 #endif
774 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 678 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
775 for (volatile octave_idx_type j = 0; j < nc; j++) 679 for (volatile octave_idx_type j = 0; j < nc; j++)
776 { 680 {
777 octave_quit (); 681 octave_quit ();
778 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 682 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
779 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 683 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
780 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 684 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
781 } 685 }
782 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 686 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
783 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); 687 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
784 #if defined (CS_VER) && (CS_VER >= 2)
785 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc); 688 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc);
786 #else
787 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xz);
788 #endif
789 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 689 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
790 690
791 for (octave_idx_type j = 0; j < nc; j++) 691 for (octave_idx_type j = 0; j < nc; j++)
792 { 692 {
793 Complex tmp = Complex (Xx[j], Xz[j]); 693 Complex tmp = Complex (Xx[j], Xz[j]);
833 Xz[j] = std::imag (c); 733 Xz[j] = std::imag (c);
834 } 734 }
835 for (octave_idx_type j = nr; j < nbuf; j++) 735 for (octave_idx_type j = nr; j < nbuf; j++)
836 buf[j] = 0.; 736 buf[j] = 0.;
837 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 737 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
838 #if defined (CS_VER) && (CS_VER >= 2)
839 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr); 738 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
840 #else
841 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xx, buf);
842 #endif
843 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); 739 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
844 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 740 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
845 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 741 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
846 { 742 {
847 octave_quit (); 743 octave_quit ();
848 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 744 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
849 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 745 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
850 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 746 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
851 } 747 }
852 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 748 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
853 #if defined (CS_VER) && (CS_VER >= 2)
854 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc); 749 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
855 #else
856 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xx);
857 #endif
858 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 750 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
859 for (octave_idx_type j = nr; j < nbuf; j++) 751 for (octave_idx_type j = nr; j < nbuf; j++)
860 buf[j] = 0.; 752 buf[j] = 0.;
861 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 753 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
862 #if defined (CS_VER) && (CS_VER >= 2)
863 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr); 754 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr);
864 #else
865 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xz, buf);
866 #endif
867 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); 755 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
868 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 756 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
869 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 757 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
870 { 758 {
871 octave_quit (); 759 octave_quit ();
872 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 760 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
873 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); 761 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
874 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 762 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
875 } 763 }
876 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 764 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
877 #if defined (CS_VER) && (CS_VER >= 2)
878 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc); 765 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc);
879 #else
880 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xz);
881 #endif
882 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 766 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
883 767
884 for (octave_idx_type j = 0; j < nc; j++) 768 for (octave_idx_type j = 0; j < nc; j++)
885 { 769 {
886 Complex tmp = Complex (Xx[j], Xz[j]); 770 Complex tmp = Complex (Xx[j], Xz[j]);