Mercurial > octave-nkf

comparison liboctave/SparseCmplxQR.cc @ 5792:eb90c83b4f91

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-05-04 20:14:49 by dbateman]
author dbateman
date Thu, 04 May 2006 20:14:50 +0000
parents ace8d8d26933
children 11fcab4c461d
comparison
equal deleted inserted replaced
5791:70215aff5ccf 5792:eb90c83b4f91
51 A.i = const_cast<octave_idx_type *>(a.ridx ()); 51 A.i = const_cast<octave_idx_type *>(a.ridx ());
52 A.x = const_cast<double _Complex *>(reinterpret_cast<const double _Complex *> 52 A.x = const_cast<double _Complex *>(reinterpret_cast<const double _Complex *>
53 (a.data ())); 53 (a.data ()));
54 A.nz = -1; 54 A.nz = -1;
55 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 55 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
56 S = CXSPARSE_ZNAME (_sqr) (&A, order, 1); 56 #if defined(CS_VER) && (CS_VER >= 2)
57 S = CXSPARSE_ZNAME (_sqr) (order, &A, 1);
58 #else
59 S = CXSPARSE_ZNAME (_sqr) (&A, order - 1, 1);
60 #endif
57 N = CXSPARSE_ZNAME (_qr) (&A, S); 61 N = CXSPARSE_ZNAME (_qr) (&A, S);
58 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 62 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
59 if (!N) 63 if (!N)
60 (*current_liboctave_error_handler) 64 (*current_liboctave_error_handler)
61 ("SparseComplexQR: sparse matrix QR factorization filled"); 65 ("SparseComplexQR: sparse matrix QR factorization filled");
108 SparseComplexQR::SparseComplexQR_rep::Pinv (void) const 112 SparseComplexQR::SparseComplexQR_rep::Pinv (void) const
109 { 113 {
110 #ifdef HAVE_CXSPARSE 114 #ifdef HAVE_CXSPARSE
111 ColumnVector ret(N->L->m); 115 ColumnVector ret(N->L->m);
112 for (octave_idx_type i = 0; i < N->L->m; i++) 116 for (octave_idx_type i = 0; i < N->L->m; i++)
117 #if defined(CS_VER) && (CS_VER >= 2)
118 ret.xelem(i) = S->pinv[i];
119 #else
113 ret.xelem(i) = S->Pinv[i]; 120 ret.xelem(i) = S->Pinv[i];
121 #endif
114 return ret; 122 return ret;
115 #else 123 #else
116 return ColumnVector (); 124 return ColumnVector ();
117 #endif 125 #endif
118 } 126 }
121 SparseComplexQR::SparseComplexQR_rep::P (void) const 129 SparseComplexQR::SparseComplexQR_rep::P (void) const
122 { 130 {
123 #ifdef HAVE_CXSPARSE 131 #ifdef HAVE_CXSPARSE
124 ColumnVector ret(N->L->m); 132 ColumnVector ret(N->L->m);
125 for (octave_idx_type i = 0; i < N->L->m; i++) 133 for (octave_idx_type i = 0; i < N->L->m; i++)
134 #if defined(CS_VER) && (CS_VER >= 2)
135 ret.xelem(S->pinv[i]) = i;
136 #else
126 ret.xelem(S->Pinv[i]) = i; 137 ret.xelem(S->Pinv[i]) = i;
138 #endif
127 return ret; 139 return ret;
128 #else 140 #else
129 return ColumnVector (); 141 return ColumnVector ();
130 #endif 142 #endif
131 } 143 }
180 for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr) 192 for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
181 { 193 {
182 OCTAVE_QUIT; 194 OCTAVE_QUIT;
183 volatile octave_idx_type nm = (nr < nc ? nr : nc); 195 volatile octave_idx_type nm = (nr < nc ? nr : nc);
184 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 196 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
185 CXSPARSE_ZNAME (_ipvec) (b_nr, S->Pinv, bvec + idx, 197 #if defined(CS_VER) && (CS_VER >= 2)
186 reinterpret_cast<double _Complex *>(buf)); 198 CXSPARSE_ZNAME (_ipvec)
199 (S->pinv, bvec + idx, reinterpret_cast<double _Complex *>(buf), b_nr);
200 #else
201 CXSPARSE_ZNAME (_ipvec)
202 (b_nr, S->Pinv, bvec + idx, reinterpret_cast<double _Complex *>(buf));
203 #endif
187 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 204 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
188 for (volatile octave_idx_type i = 0; i < nm; i++) 205 for (volatile octave_idx_type i = 0; i < nm; i++)
189 { 206 {
190 OCTAVE_QUIT; 207 OCTAVE_QUIT;
191 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 208 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
236 for (octave_idx_type j = 0; j < b_nr; j++) 253 for (octave_idx_type j = 0; j < b_nr; j++)
237 Xx[j] = b.xelem(j,i); 254 Xx[j] = b.xelem(j,i);
238 for (octave_idx_type j = nr; j < q.S()->m2; j++) 255 for (octave_idx_type j = nr; j < q.S()->m2; j++)
239 buf[j] = OCTAVE_C99_ZERO; 256 buf[j] = OCTAVE_C99_ZERO;
240 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 257 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
258 #if defined(CS_VER) && (CS_VER >= 2)
259 CXSPARSE_ZNAME (_ipvec)
260 (q.S()->pinv, reinterpret_cast<double _Complex *>(Xx), buf, nr);
261 #else
241 CXSPARSE_ZNAME (_ipvec) 262 CXSPARSE_ZNAME (_ipvec)
242 (nr, q.S()->Pinv, reinterpret_cast<double _Complex *>(Xx), buf); 263 (nr, q.S()->Pinv, reinterpret_cast<double _Complex *>(Xx), buf);
264 #endif
243 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 265 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
244 for (volatile octave_idx_type j = 0; j < nc; j++) 266 for (volatile octave_idx_type j = 0; j < nc; j++)
245 { 267 {
246 OCTAVE_QUIT; 268 OCTAVE_QUIT;
247 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 269 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
248 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 270 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
249 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 271 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
250 } 272 }
251 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 273 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
252 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf); 274 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf);
275 #if defined(CS_VER) && (CS_VER >= 2)
276 CXSPARSE_ZNAME (_ipvec) (q.S()->q, buf, vec + idx, nc);
277 #else
253 CXSPARSE_ZNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx); 278 CXSPARSE_ZNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx);
279 #endif
254 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 280 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
255 } 281 }
256 } 282 }
257 else 283 else
258 { 284 {
278 for (octave_idx_type j = 0; j < b_nr; j++) 304 for (octave_idx_type j = 0; j < b_nr; j++)
279 Xx[j] = b.xelem(j,i); 305 Xx[j] = b.xelem(j,i);
280 for (octave_idx_type j = nr; j < nbuf; j++) 306 for (octave_idx_type j = nr; j < nbuf; j++)
281 buf[j] = OCTAVE_C99_ZERO; 307 buf[j] = OCTAVE_C99_ZERO;
282 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 308 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
309 #if defined(CS_VER) && (CS_VER >= 2)
310 CXSPARSE_ZNAME (_pvec)
311 (q.S()->q, reinterpret_cast<double _Complex *>(Xx), buf, nr);
312 #else
283 CXSPARSE_ZNAME (_pvec) 313 CXSPARSE_ZNAME (_pvec)
284 (nr, q.S()->Q, reinterpret_cast<double _Complex *>(Xx), buf); 314 (nr, q.S()->Q, reinterpret_cast<double _Complex *>(Xx), buf);
315 #endif
285 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf); 316 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf);
286 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 317 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
287 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 318 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
288 { 319 {
289 OCTAVE_QUIT; 320 OCTAVE_QUIT;
292 CXSPARSE_ZNAME (_happly) 323 CXSPARSE_ZNAME (_happly)
293 (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf); 324 (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf);
294 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 325 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
295 } 326 }
296 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 327 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
328 #if defined(CS_VER) && (CS_VER >= 2)
329 CXSPARSE_ZNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc);
330 #else
297 CXSPARSE_ZNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx); 331 CXSPARSE_ZNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx);
332 #endif
298 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 333 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
299 } 334 }
300 } 335 }
301 336
302 return x; 337 return x;
340 for (octave_idx_type j = 0; j < b_nr; j++) 375 for (octave_idx_type j = 0; j < b_nr; j++)
341 Xx[j] = b.xelem(j,i); 376 Xx[j] = b.xelem(j,i);
342 for (octave_idx_type j = nr; j < q.S()->m2; j++) 377 for (octave_idx_type j = nr; j < q.S()->m2; j++)
343 buf[j] = OCTAVE_C99_ZERO; 378 buf[j] = OCTAVE_C99_ZERO;
344 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 379 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
380 #if defined(CS_VER) && (CS_VER >= 2)
381 CXSPARSE_ZNAME (_ipvec)
382 (q.S()->pinv, reinterpret_cast<double _Complex *>(Xx), buf, nr);
383 #else
345 CXSPARSE_ZNAME (_ipvec) 384 CXSPARSE_ZNAME (_ipvec)
346 (nr, q.S()->Pinv, reinterpret_cast<double _Complex *>(Xx), buf); 385 (nr, q.S()->Pinv, reinterpret_cast<double _Complex *>(Xx), buf);
386 #endif
347 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 387 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
348 for (volatile octave_idx_type j = 0; j < nc; j++) 388 for (volatile octave_idx_type j = 0; j < nc; j++)
349 { 389 {
350 OCTAVE_QUIT; 390 OCTAVE_QUIT;
351 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 391 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
352 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 392 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
353 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 393 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
354 } 394 }
355 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 395 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
356 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf); 396 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf);
357 CXSPARSE_ZNAME (_ipvec) (nc, q.S()->Q, buf, 397 #if defined(CS_VER) && (CS_VER >= 2)
358 reinterpret_cast<double _Complex *>(Xx)); 398 CXSPARSE_ZNAME (_ipvec)
399 (q.S()->q, buf, reinterpret_cast<double _Complex *>(Xx), nc);
400 #else
401 CXSPARSE_ZNAME (_ipvec)
402 (nc, q.S()->Q, buf, reinterpret_cast<double _Complex *>(Xx));
403 #endif
359 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 404 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
360 405
361 for (octave_idx_type j = 0; j < nc; j++) 406 for (octave_idx_type j = 0; j < nc; j++)
362 { 407 {
363 Complex tmp = Xx[j]; 408 Complex tmp = Xx[j];
403 for (octave_idx_type j = 0; j < b_nr; j++) 448 for (octave_idx_type j = 0; j < b_nr; j++)
404 Xx[j] = b.xelem(j,i); 449 Xx[j] = b.xelem(j,i);
405 for (octave_idx_type j = nr; j < nbuf; j++) 450 for (octave_idx_type j = nr; j < nbuf; j++)
406 buf[j] = OCTAVE_C99_ZERO; 451 buf[j] = OCTAVE_C99_ZERO;
407 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 452 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
453 #if defined(CS_VER) && (CS_VER >= 2)
454 CXSPARSE_ZNAME (_pvec)
455 (q.S()->q, reinterpret_cast<double _Complex *>(Xx), buf, nr);
456 #else
408 CXSPARSE_ZNAME (_pvec) 457 CXSPARSE_ZNAME (_pvec)
409 (nr, q.S()->Q, reinterpret_cast<double _Complex *>(Xx), buf); 458 (nr, q.S()->Q, reinterpret_cast<double _Complex *>(Xx), buf);
459 #endif
410 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf); 460 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf);
411 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 461 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
412 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 462 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
413 { 463 {
414 OCTAVE_QUIT; 464 OCTAVE_QUIT;
416 CXSPARSE_ZNAME (_happly) 466 CXSPARSE_ZNAME (_happly)
417 (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf); 467 (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf);
418 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 468 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
419 } 469 }
420 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 470 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
421 CXSPARSE_ZNAME (_pvec) (nc, q.S()->Pinv, buf, 471 #if defined(CS_VER) && (CS_VER >= 2)
422 reinterpret_cast<double _Complex *>(Xx)); 472 CXSPARSE_ZNAME (_pvec)
473 (q.S()->pinv, buf, reinterpret_cast<double _Complex *>(Xx), nc);
474 #else
475 CXSPARSE_ZNAME (_pvec)
476 (nc, q.S()->Pinv, buf, reinterpret_cast<double _Complex *>(Xx));
477 #endif
423 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 478 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
424 479
425 for (octave_idx_type j = 0; j < nc; j++) 480 for (octave_idx_type j = 0; j < nc; j++)
426 { 481 {
427 Complex tmp = Xx[j]; 482 Complex tmp = Xx[j];
483 { 538 {
484 OCTAVE_QUIT; 539 OCTAVE_QUIT;
485 for (octave_idx_type j = nr; j < q.S()->m2; j++) 540 for (octave_idx_type j = nr; j < q.S()->m2; j++)
486 buf[j] = OCTAVE_C99_ZERO; 541 buf[j] = OCTAVE_C99_ZERO;
487 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 542 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
543 #if defined(CS_VER) && (CS_VER >= 2)
544 CXSPARSE_ZNAME (_ipvec) (q.S()->pinv, bvec + bidx, buf, nr);
545 #else
488 CXSPARSE_ZNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf); 546 CXSPARSE_ZNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf);
547 #endif
489 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 548 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
490 for (volatile octave_idx_type j = 0; j < nc; j++) 549 for (volatile octave_idx_type j = 0; j < nc; j++)
491 { 550 {
492 OCTAVE_QUIT; 551 OCTAVE_QUIT;
493 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 552 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
494 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 553 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
495 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 554 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
496 } 555 }
497 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 556 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
498 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf); 557 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf);
558 #if defined(CS_VER) && (CS_VER >= 2)
559 CXSPARSE_ZNAME (_ipvec) (q.S()->q, buf, vec + idx, nc);
560 #else
499 CXSPARSE_ZNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx); 561 CXSPARSE_ZNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx);
562 #endif
500 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 563 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
501 } 564 }
502 } 565 }
503 else 566 else
504 { 567 {
522 { 585 {
523 OCTAVE_QUIT; 586 OCTAVE_QUIT;
524 for (octave_idx_type j = nr; j < nbuf; j++) 587 for (octave_idx_type j = nr; j < nbuf; j++)
525 buf[j] = OCTAVE_C99_ZERO; 588 buf[j] = OCTAVE_C99_ZERO;
526 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 589 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
590 #if defined(CS_VER) && (CS_VER >= 2)
591 CXSPARSE_ZNAME (_pvec) (q.S()->q, bvec + bidx, buf, nr);
592 #else
527 CXSPARSE_ZNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf); 593 CXSPARSE_ZNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf);
594 #endif
528 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf); 595 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf);
529 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 596 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
530 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 597 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
531 { 598 {
532 OCTAVE_QUIT; 599 OCTAVE_QUIT;
534 CXSPARSE_ZNAME (_happly) 601 CXSPARSE_ZNAME (_happly)
535 (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf); 602 (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf);
536 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 603 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
537 } 604 }
538 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 605 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
606 #if defined(CS_VER) && (CS_VER >= 2)
607 CXSPARSE_ZNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc);
608 #else
539 CXSPARSE_ZNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx); 609 CXSPARSE_ZNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx);
610 #endif
540 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 611 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
541 } 612 }
542 } 613 }
543 614
544 return x; 615 return x;
582 for (octave_idx_type j = 0; j < b_nr; j++) 653 for (octave_idx_type j = 0; j < b_nr; j++)
583 Xx[j] = b.xelem(j,i); 654 Xx[j] = b.xelem(j,i);
584 for (octave_idx_type j = nr; j < q.S()->m2; j++) 655 for (octave_idx_type j = nr; j < q.S()->m2; j++)
585 buf[j] = OCTAVE_C99_ZERO; 656 buf[j] = OCTAVE_C99_ZERO;
586 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 657 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
658 #if defined(CS_VER) && (CS_VER >= 2)
659 CXSPARSE_ZNAME (_ipvec)
660 (q.S()->pinv, reinterpret_cast<double _Complex *>(Xx), buf, nr);
661 #else
587 CXSPARSE_ZNAME (_ipvec) 662 CXSPARSE_ZNAME (_ipvec)
588 (nr, q.S()->Pinv, reinterpret_cast<double _Complex *>(Xx), buf); 663 (nr, q.S()->Pinv, reinterpret_cast<double _Complex *>(Xx), buf);
664 #endif
589 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 665 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
590 for (volatile octave_idx_type j = 0; j < nc; j++) 666 for (volatile octave_idx_type j = 0; j < nc; j++)
591 { 667 {
592 OCTAVE_QUIT; 668 OCTAVE_QUIT;
593 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 669 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
594 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 670 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
595 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 671 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
596 } 672 }
597 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 673 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
598 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf); 674 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf);
599 CXSPARSE_ZNAME (_ipvec) (nc, q.S()->Q, buf, 675 #if defined(CS_VER) && (CS_VER >= 2)
600 reinterpret_cast<double _Complex *>(Xx)); 676 CXSPARSE_ZNAME (_ipvec)
677 (q.S()->q, buf, reinterpret_cast<double _Complex *>(Xx), nc);
678 #else
679 CXSPARSE_ZNAME (_ipvec)
680 (nc, q.S()->Q, buf, reinterpret_cast<double _Complex *>(Xx));
681 #endif
601 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 682 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
602 683
603 for (octave_idx_type j = 0; j < nc; j++) 684 for (octave_idx_type j = 0; j < nc; j++)
604 { 685 {
605 Complex tmp = Xx[j]; 686 Complex tmp = Xx[j];
645 for (octave_idx_type j = 0; j < b_nr; j++) 726 for (octave_idx_type j = 0; j < b_nr; j++)
646 Xx[j] = b.xelem(j,i); 727 Xx[j] = b.xelem(j,i);
647 for (octave_idx_type j = nr; j < nbuf; j++) 728 for (octave_idx_type j = nr; j < nbuf; j++)
648 buf[j] = OCTAVE_C99_ZERO; 729 buf[j] = OCTAVE_C99_ZERO;
649 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 730 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
731 #if defined(CS_VER) && (CS_VER >= 2)
732 CXSPARSE_ZNAME (_pvec)
733 (q.S()->q, reinterpret_cast<double _Complex *>(Xx), buf, nr);
734 #else
650 CXSPARSE_ZNAME (_pvec) 735 CXSPARSE_ZNAME (_pvec)
651 (nr, q.S()->Q, reinterpret_cast<double _Complex *>(Xx), buf); 736 (nr, q.S()->Q, reinterpret_cast<double _Complex *>(Xx), buf);
737 #endif
652 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf); 738 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf);
653 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 739 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
654 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 740 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
655 { 741 {
656 OCTAVE_QUIT; 742 OCTAVE_QUIT;
658 CXSPARSE_ZNAME (_happly) 744 CXSPARSE_ZNAME (_happly)
659 (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf); 745 (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf);
660 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 746 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
661 } 747 }
662 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 748 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
663 CXSPARSE_ZNAME (_pvec) (nc, q.S()->Pinv, buf, 749 #if defined(CS_VER) && (CS_VER >= 2)
664 reinterpret_cast<double _Complex *>(Xx)); 750 CXSPARSE_ZNAME (_pvec)
751 (q.S()->pinv, buf, reinterpret_cast<double _Complex *>(Xx), nc);
752 #else
753 CXSPARSE_ZNAME (_pvec)
754 (nc, q.S()->Pinv, buf, reinterpret_cast<double _Complex *>(Xx));
755 #endif
665 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 756 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
666 757
667 for (octave_idx_type j = 0; j < nc; j++) 758 for (octave_idx_type j = 0; j < nc; j++)
668 { 759 {
669 Complex tmp = Xx[j]; 760 Complex tmp = Xx[j];