Mercurial > hg > octave-nkf

comparison liboctave/SparseQR.cc @ 5648:69a4f320d95a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-03-08 20:17:37 by dbateman]
author dbateman
date Wed, 08 Mar 2006 20:17:38 +0000
parents 9761b7d24e9e
children 233d98d95659
comparison
equal deleted inserted replaced
5647:9e3a2d1e5e72 5648:69a4f320d95a
28 #include "SparseQR.h" 28 #include "SparseQR.h"
29 29
30 SparseQR::SparseQR_rep::SparseQR_rep (const SparseMatrix& a, int order) 30 SparseQR::SparseQR_rep::SparseQR_rep (const SparseMatrix& a, int order)
31 { 31 {
32 #ifdef HAVE_CXSPARSE 32 #ifdef HAVE_CXSPARSE
33 CXSPARSE_DNAME (cs) A; 33 CXSPARSE_DNAME () A;
34 A.nzmax = a.nzmax (); 34 A.nzmax = a.nzmax ();
35 A.m = a.rows (); 35 A.m = a.rows ();
36 A.n = a.cols (); 36 A.n = a.cols ();
37 nrows = A.m; 37 nrows = A.m;
38 // Cast away const on A, with full knowledge that CSparse won't touch it 38 // Cast away const on A, with full knowledge that CSparse won't touch it
40 A.p = const_cast<octave_idx_type *>(a.cidx ()); 40 A.p = const_cast<octave_idx_type *>(a.cidx ());
41 A.i = const_cast<octave_idx_type *>(a.ridx ()); 41 A.i = const_cast<octave_idx_type *>(a.ridx ());
42 A.x = const_cast<double *>(a.data ()); 42 A.x = const_cast<double *>(a.data ());
43 A.nz = -1; 43 A.nz = -1;
44 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 44 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
45 S = CXSPARSE_DNAME (cs_sqr) (&A, order, 1); 45 S = CXSPARSE_DNAME (_sqr) (&A, order, 1);
46 N = CXSPARSE_DNAME (cs_qr) (&A, S); 46 N = CXSPARSE_DNAME (_qr) (&A, S);
47 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 47 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
48 if (!N) 48 if (!N)
49 (*current_liboctave_error_handler) 49 (*current_liboctave_error_handler)
50 ("SparseQR: sparse matrix QR factorization filled"); 50 ("SparseQR: sparse matrix QR factorization filled");
51 count = 1; 51 count = 1;
56 } 56 }
57 57
58 SparseQR::SparseQR_rep::~SparseQR_rep (void) 58 SparseQR::SparseQR_rep::~SparseQR_rep (void)
59 { 59 {
60 #ifdef HAVE_CXSPARSE 60 #ifdef HAVE_CXSPARSE
61 CXSPARSE_DNAME (cs_sfree) (S); 61 CXSPARSE_DNAME (_sfree) (S);
62 CXSPARSE_DNAME (cs_nfree) (N); 62 CXSPARSE_DNAME (_nfree) (N);
63 #endif 63 #endif
64 } 64 }
65 65
66 SparseMatrix 66 SparseMatrix
67 SparseQR::SparseQR_rep::V (void) const 67 SparseQR::SparseQR_rep::V (void) const
68 { 68 {
69 #ifdef HAVE_CXSPARSE 69 #ifdef HAVE_CXSPARSE
70 // Drop zeros from V and sort 70 // Drop zeros from V and sort
71 // XXX FIXME XXX Is the double transpose to sort necessary? 71 // XXX FIXME XXX Is the double transpose to sort necessary?
72 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 72 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
73 CXSPARSE_DNAME (cs_dropzeros) (N->L); 73 CXSPARSE_DNAME (_dropzeros) (N->L);
74 CXSPARSE_DNAME (cs) *D = CXSPARSE_DNAME (cs_transpose) (N->L, 1); 74 CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->L, 1);
75 CXSPARSE_DNAME (cs_spfree) (N->L); 75 CXSPARSE_DNAME (_spfree) (N->L);
76 N->L = CXSPARSE_DNAME (cs_transpose) (D, 1); 76 N->L = CXSPARSE_DNAME (_transpose) (D, 1);
77 CXSPARSE_DNAME (cs_spfree) (D); 77 CXSPARSE_DNAME (_spfree) (D);
78 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 78 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
79 79
80 octave_idx_type nc = N->L->n; 80 octave_idx_type nc = N->L->n;
81 octave_idx_type nz = N->L->nzmax; 81 octave_idx_type nz = N->L->nzmax;
82 SparseMatrix ret (N->L->m, nc, nz); 82 SparseMatrix ret (N->L->m, nc, nz);
124 { 124 {
125 #ifdef HAVE_CXSPARSE 125 #ifdef HAVE_CXSPARSE
126 // Drop zeros from R and sort 126 // Drop zeros from R and sort
127 // XXX FIXME XXX Is the double transpose to sort necessary? 127 // XXX FIXME XXX Is the double transpose to sort necessary?
128 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 128 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
129 CXSPARSE_DNAME (cs_dropzeros) (N->U); 129 CXSPARSE_DNAME (_dropzeros) (N->U);
130 CXSPARSE_DNAME (cs) *D = CXSPARSE_DNAME (cs_transpose) (N->U, 1); 130 CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->U, 1);
131 CXSPARSE_DNAME (cs_spfree) (N->U); 131 CXSPARSE_DNAME (_spfree) (N->U);
132 N->U = CXSPARSE_DNAME (cs_transpose) (D, 1); 132 N->U = CXSPARSE_DNAME (_transpose) (D, 1);
133 CXSPARSE_DNAME (cs_spfree) (D); 133 CXSPARSE_DNAME (_spfree) (D);
134 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 134 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
135 135
136 octave_idx_type nc = N->U->n; 136 octave_idx_type nc = N->U->n;
137 octave_idx_type nz = N->U->nzmax; 137 octave_idx_type nz = N->U->nzmax;
138 138
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 volatile octave_idx_type nm = (nr < nc ? nr : nc); 173 volatile octave_idx_type nm = (nr < nc ? nr : nc);
174 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 174 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
175 // cast away const on bvec, with full knowledge that CSparse 175 CXSPARSE_DNAME (_ipvec) (b_nr, S->Pinv, bvec + idx, buf);
176 // won't touch it
177 CXSPARSE_DNAME (cs_ipvec) (b_nr, S->Pinv, bvec + idx, buf);
178 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 176 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
179 177
180 for (volatile octave_idx_type i = 0; i < nm; i++) 178 for (volatile octave_idx_type i = 0; i < nm; i++)
181 { 179 {
182 OCTAVE_QUIT; 180 OCTAVE_QUIT;
183 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 181 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
184 CXSPARSE_DNAME (cs_happly) (N->L, i, N->B[i], buf); 182 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
185 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 183 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
186 } 184 }
187 for (octave_idx_type i = 0; i < b_nr; i++) 185 for (octave_idx_type i = 0; i < b_nr; i++)
188 vec[i+idx] = buf[i]; 186 vec[i+idx] = buf[i];
189 } 187 }
223 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 221 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
224 i++, idx+=nc, bidx+=b_nr) 222 i++, idx+=nc, bidx+=b_nr)
225 { 223 {
226 OCTAVE_QUIT; 224 OCTAVE_QUIT;
227 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 225 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
228 // cast away const on bvec, with full knowledge that CSparse 226 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf);
229 // won't touch it
230 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf);
231 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 227 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
232 for (volatile octave_idx_type j = 0; j < nc; j++) 228 for (volatile octave_idx_type j = 0; j < nc; j++)
233 { 229 {
234 OCTAVE_QUIT; 230 OCTAVE_QUIT;
235 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 231 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
236 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 232 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
237 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 233 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
238 } 234 }
239 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 235 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
240 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); 236 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
241 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, vec + idx); 237 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx);
242 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 238 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
243 } 239 }
244 } 240 }
245 else 241 else
246 { 242 {
257 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 253 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
258 i++, idx+=nc, bidx+=b_nr) 254 i++, idx+=nc, bidx+=b_nr)
259 { 255 {
260 OCTAVE_QUIT; 256 OCTAVE_QUIT;
261 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 257 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
262 // cast away const on bvec, with full knowledge that CSparse 258 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf);
263 // won't touch it 259 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
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; 260 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
267 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 261 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
268 { 262 {
269 OCTAVE_QUIT; 263 OCTAVE_QUIT;
270 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 264 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
271 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 265 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
272 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 266 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
273 } 267 }
274 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 268 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
275 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, vec + idx); 269 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx);
276 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 270 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
277 } 271 }
278 } 272 }
279 273
280 return x; 274 return x;
316 { 310 {
317 OCTAVE_QUIT; 311 OCTAVE_QUIT;
318 for (octave_idx_type j = 0; j < b_nr; j++) 312 for (octave_idx_type j = 0; j < b_nr; j++)
319 Xx[j] = b.xelem(j,i); 313 Xx[j] = b.xelem(j,i);
320 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 314 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
321 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xx, buf); 315 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
322 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 316 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
323 for (volatile octave_idx_type j = 0; j < nc; j++) 317 for (volatile octave_idx_type j = 0; j < nc; j++)
324 { 318 {
325 OCTAVE_QUIT; 319 OCTAVE_QUIT;
326 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 320 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
327 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 321 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
328 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 322 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
329 } 323 }
330 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 324 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
331 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); 325 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
332 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xx); 326 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx);
333 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 327 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
334 328
335 for (octave_idx_type j = 0; j < nc; j++) 329 for (octave_idx_type j = 0; j < nc; j++)
336 { 330 {
337 double tmp = Xx[j]; 331 double tmp = Xx[j];
371 { 365 {
372 OCTAVE_QUIT; 366 OCTAVE_QUIT;
373 for (octave_idx_type j = 0; j < b_nr; j++) 367 for (octave_idx_type j = 0; j < b_nr; j++)
374 Xx[j] = b.xelem(j,i); 368 Xx[j] = b.xelem(j,i);
375 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 369 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
376 CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xx, buf); 370 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
377 CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf); 371 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
378 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 372 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
379 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 373 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
380 { 374 {
381 OCTAVE_QUIT; 375 OCTAVE_QUIT;
382 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 376 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
383 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 377 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
384 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 378 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
385 } 379 }
386 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 380 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
387 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xx); 381 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx);
388 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 382 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
389 383
390 for (octave_idx_type j = 0; j < nc; j++) 384 for (octave_idx_type j = 0; j < nc; j++)
391 { 385 {
392 double tmp = Xx[j]; 386 double tmp = Xx[j];
450 Complex c = b.xelem (j,i); 444 Complex c = b.xelem (j,i);
451 Xx[j] = std::real (c); 445 Xx[j] = std::real (c);
452 Xz[j] = std::imag (c); 446 Xz[j] = std::imag (c);
453 } 447 }
454 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 448 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
455 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xx, buf); 449 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
456 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 450 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
457 for (volatile octave_idx_type j = 0; j < nc; j++) 451 for (volatile octave_idx_type j = 0; j < nc; j++)
458 { 452 {
459 OCTAVE_QUIT; 453 OCTAVE_QUIT;
460 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 454 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
461 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 455 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
462 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 456 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
463 } 457 }
464 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 458 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
465 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); 459 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
466 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xx); 460 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx);
467 461
468 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xz, buf); 462 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf);
469 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 463 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
470 for (volatile octave_idx_type j = 0; j < nc; j++) 464 for (volatile octave_idx_type j = 0; j < nc; j++)
471 { 465 {
472 OCTAVE_QUIT; 466 OCTAVE_QUIT;
473 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 467 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
474 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 468 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
475 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 469 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
476 } 470 }
477 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 471 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
478 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); 472 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
479 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xz); 473 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz);
480 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 474 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
481 for (octave_idx_type j = 0; j < nc; j++) 475 for (octave_idx_type j = 0; j < nc; j++)
482 vec[j+idx] = Complex (Xx[j], Xz[j]); 476 vec[j+idx] = Complex (Xx[j], Xz[j]);
483 } 477 }
484 } 478 }
504 Complex c = b.xelem (j,i); 498 Complex c = b.xelem (j,i);
505 Xx[j] = std::real (c); 499 Xx[j] = std::real (c);
506 Xz[j] = std::imag (c); 500 Xz[j] = std::imag (c);
507 } 501 }
508 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 502 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
509 CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xx, buf); 503 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
510 CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf); 504 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
511 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 505 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
512 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 506 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
513 { 507 {
514 OCTAVE_QUIT; 508 OCTAVE_QUIT;
515 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 509 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
516 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 510 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
517 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 511 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
518 } 512 }
519 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 513 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
520 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xx); 514 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx);
521 CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xz, buf); 515 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf);
522 CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf); 516 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
523 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 517 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
524 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 518 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
525 { 519 {
526 OCTAVE_QUIT; 520 OCTAVE_QUIT;
527 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 521 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
528 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 522 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
529 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 523 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
530 } 524 }
531 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 525 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
532 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xz); 526 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz);
533 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 527 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
534 for (octave_idx_type j = 0; j < nc; j++) 528 for (octave_idx_type j = 0; j < nc; j++)
535 vec[j+idx] = Complex (Xx[j], Xz[j]); 529 vec[j+idx] = Complex (Xx[j], Xz[j]);
536 } 530 }
537 } 531 }
580 Complex c = b.xelem (j,i); 574 Complex c = b.xelem (j,i);
581 Xx[j] = std::real (c); 575 Xx[j] = std::real (c);
582 Xz[j] = std::imag (c); 576 Xz[j] = std::imag (c);
583 } 577 }
584 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 578 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
585 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xx, buf); 579 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf);
586 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 580 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
587 for (volatile octave_idx_type j = 0; j < nc; j++) 581 for (volatile octave_idx_type j = 0; j < nc; j++)
588 { 582 {
589 OCTAVE_QUIT; 583 OCTAVE_QUIT;
590 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 584 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
591 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 585 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
592 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 586 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
593 } 587 }
594 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 588 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
595 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); 589 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
596 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xx); 590 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx);
597 CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xz, buf); 591 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf);
598 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 592 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
599 for (volatile octave_idx_type j = 0; j < nc; j++) 593 for (volatile octave_idx_type j = 0; j < nc; j++)
600 { 594 {
601 OCTAVE_QUIT; 595 OCTAVE_QUIT;
602 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 596 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
603 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 597 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
604 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 598 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
605 } 599 }
606 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 600 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
607 CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); 601 CXSPARSE_DNAME (_usolve) (q.N()->U, buf);
608 CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xz); 602 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz);
609 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 603 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
610 604
611 for (octave_idx_type j = 0; j < nc; j++) 605 for (octave_idx_type j = 0; j < nc; j++)
612 { 606 {
613 Complex tmp = Complex (Xx[j], Xz[j]); 607 Complex tmp = Complex (Xx[j], Xz[j]);
652 Complex c = b.xelem (j,i); 646 Complex c = b.xelem (j,i);
653 Xx[j] = std::real (c); 647 Xx[j] = std::real (c);
654 Xz[j] = std::imag (c); 648 Xz[j] = std::imag (c);
655 } 649 }
656 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 650 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
657 CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xx, buf); 651 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf);
658 CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf); 652 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
659 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 653 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
660 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 654 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
661 { 655 {
662 OCTAVE_QUIT; 656 OCTAVE_QUIT;
663 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 657 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
664 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 658 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
665 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 659 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
666 } 660 }
667 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 661 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
668 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xx); 662 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx);
669 CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xz, buf); 663 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf);
670 CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf); 664 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf);
671 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 665 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
672 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 666 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
673 { 667 {
674 OCTAVE_QUIT; 668 OCTAVE_QUIT;
675 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 669 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
676 CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); 670 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
677 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 671 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
678 } 672 }
679 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 673 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
680 CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xz); 674 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz);
681 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 675 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
682 676
683 for (octave_idx_type j = 0; j < nc; j++) 677 for (octave_idx_type j = 0; j < nc; j++)
684 { 678 {
685 Complex tmp = Complex (Xx[j], Xz[j]); 679 Complex tmp = Complex (Xx[j], Xz[j]);