Mercurial > hg > octave-nkf

comparison liboctave/SparseCmplxQR.cc @ 10314:07ebe522dac2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
untabify liboctave C++ sources
author John W. Eaton <jwe@octave.org>
date Thu, 11 Feb 2010 12:23:32 -0500
parents 4c0cdbe0acca
children 12884915a8e4
comparison
equal deleted inserted replaced
10313:f3b65e1ae355 10314:07ebe522dac2
60 // Cast away const on A, with full knowledge that CSparse won't touch it 60 // Cast away const on A, with full knowledge that CSparse won't touch it
61 // Prevents the methods below making a copy of the data. 61 // Prevents the methods below making a copy of the data.
62 A.p = const_cast<octave_idx_type *>(a.cidx ()); 62 A.p = const_cast<octave_idx_type *>(a.cidx ());
63 A.i = const_cast<octave_idx_type *>(a.ridx ()); 63 A.i = const_cast<octave_idx_type *>(a.ridx ());
64 A.x = const_cast<cs_complex_t *>(reinterpret_cast<const cs_complex_t *> 64 A.x = const_cast<cs_complex_t *>(reinterpret_cast<const cs_complex_t *>
65 (a.data ())); 65 (a.data ()));
66 A.nz = -1; 66 A.nz = -1;
67 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 67 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
68 #if defined(CS_VER) && (CS_VER >= 2) 68 #if defined(CS_VER) && (CS_VER >= 2)
69 S = CXSPARSE_ZNAME (_sqr) (order, &A, 1); 69 S = CXSPARSE_ZNAME (_sqr) (order, &A, 1);
70 #else 70 #else
202 ret = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); 202 ret = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
203 else 203 else
204 { 204 {
205 OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2); 205 OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
206 for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr) 206 for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
207 { 207 {
208 octave_quit (); 208 octave_quit ();
209 volatile octave_idx_type nm = (nr < nc ? nr : nc); 209 volatile octave_idx_type nm = (nr < nc ? nr : nc);
210 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 210 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
211 #if defined(CS_VER) && (CS_VER >= 2) 211 #if defined(CS_VER) && (CS_VER >= 2)
212 CXSPARSE_ZNAME (_ipvec) 212 CXSPARSE_ZNAME (_ipvec)
213 (S->pinv, bvec + idx, reinterpret_cast<cs_complex_t *>(buf), b_nr); 213 (S->pinv, bvec + idx, reinterpret_cast<cs_complex_t *>(buf), b_nr);
214 #else 214 #else
215 CXSPARSE_ZNAME (_ipvec) 215 CXSPARSE_ZNAME (_ipvec)
216 (b_nr, S->Pinv, bvec + idx, reinterpret_cast<cs_complex_t *>(buf)); 216 (b_nr, S->Pinv, bvec + idx, reinterpret_cast<cs_complex_t *>(buf));
217 #endif 217 #endif
218 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 218 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
219 for (volatile octave_idx_type i = 0; i < nm; i++) 219 for (volatile octave_idx_type i = 0; i < nm; i++)
220 { 220 {
221 octave_quit (); 221 octave_quit ();
222 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 222 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
223 CXSPARSE_ZNAME (_happly) 223 CXSPARSE_ZNAME (_happly)
224 (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf)); 224 (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf));
225 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 225 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
226 } 226 }
227 for (octave_idx_type i = 0; i < b_nr; i++) 227 for (octave_idx_type i = 0; i < b_nr; i++)
228 vec[i+idx] = buf[i]; 228 vec[i+idx] = buf[i];
229 } 229 }
230 } 230 }
231 return ret; 231 return ret;
232 #else 232 #else
233 return ComplexMatrix (); 233 return ComplexMatrix ();
234 #endif 234 #endif
248 ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0)); 248 ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0));
249 else 249 else
250 { 250 {
251 OCTAVE_C99_COMPLEX (bvec, nr); 251 OCTAVE_C99_COMPLEX (bvec, nr);
252 for (octave_idx_type i = 0; i < nr; i++) 252 for (octave_idx_type i = 0; i < nr; i++)
253 bvec[i] = OCTAVE_C99_ZERO; 253 bvec[i] = OCTAVE_C99_ZERO;
254 OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2); 254 OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
255 for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr) 255 for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
256 { 256 {
257 octave_quit (); 257 octave_quit ();
258 bvec[j] = OCTAVE_C99_ONE; 258 bvec[j] = OCTAVE_C99_ONE;
259 volatile octave_idx_type nm = (nr < nc ? nr : nc); 259 volatile octave_idx_type nm = (nr < nc ? nr : nc);
260 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 260 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
261 #if defined(CS_VER) && (CS_VER >= 2) 261 #if defined(CS_VER) && (CS_VER >= 2)
262 CXSPARSE_ZNAME (_ipvec) 262 CXSPARSE_ZNAME (_ipvec)
263 (S->pinv, bvec, reinterpret_cast<cs_complex_t *>(buf), nr); 263 (S->pinv, bvec, reinterpret_cast<cs_complex_t *>(buf), nr);
264 #else 264 #else
265 CXSPARSE_ZNAME (_ipvec) 265 CXSPARSE_ZNAME (_ipvec)
266 (nr, S->Pinv, bvec, reinterpret_cast<cs_complex_t *>(buf)); 266 (nr, S->Pinv, bvec, reinterpret_cast<cs_complex_t *>(buf));
267 #endif 267 #endif
268 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 268 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
269 for (volatile octave_idx_type i = 0; i < nm; i++) 269 for (volatile octave_idx_type i = 0; i < nm; i++)
270 { 270 {
271 octave_quit (); 271 octave_quit ();
272 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 272 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
273 CXSPARSE_ZNAME (_happly) 273 CXSPARSE_ZNAME (_happly)
274 (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf)); 274 (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf));
275 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 275 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
276 } 276 }
277 for (octave_idx_type i = 0; i < nr; i++) 277 for (octave_idx_type i = 0; i < nr; i++)
278 vec[i+idx] = buf[i]; 278 vec[i+idx] = buf[i];
279 bvec[j] = OCTAVE_C99_ZERO; 279 bvec[j] = OCTAVE_C99_ZERO;
280 } 280 }
281 } 281 }
282 return ret.hermitian (); 282 return ret.hermitian ();
283 #else 283 #else
284 return ComplexMatrix (); 284 return ComplexMatrix ();
285 #endif 285 #endif
303 x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); 303 x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
304 else if (nr >= nc) 304 else if (nr >= nc)
305 { 305 {
306 SparseComplexQR q (a, 2); 306 SparseComplexQR q (a, 2);
307 if (! q.ok ()) 307 if (! q.ok ())
308 return ComplexMatrix(); 308 return ComplexMatrix();
309 x.resize(nc, b_nc); 309 x.resize(nc, b_nc);
310 cs_complex_t *vec = reinterpret_cast<cs_complex_t *> 310 cs_complex_t *vec = reinterpret_cast<cs_complex_t *>
311 (x.fortran_vec()); 311 (x.fortran_vec());
312 OCTAVE_C99_COMPLEX (buf, q.S()->m2); 312 OCTAVE_C99_COMPLEX (buf, q.S()->m2);
313 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); 313 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
314 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 314 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
315 { 315 {
316 octave_quit (); 316 octave_quit ();
317 for (octave_idx_type j = 0; j < b_nr; j++) 317 for (octave_idx_type j = 0; j < b_nr; j++)
318 Xx[j] = b.xelem(j,i); 318 Xx[j] = b.xelem(j,i);
319 for (octave_idx_type j = nr; j < q.S()->m2; j++) 319 for (octave_idx_type j = nr; j < q.S()->m2; j++)
320 buf[j] = OCTAVE_C99_ZERO; 320 buf[j] = OCTAVE_C99_ZERO;
321 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 321 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
322 #if defined(CS_VER) && (CS_VER >= 2) 322 #if defined(CS_VER) && (CS_VER >= 2)
323 CXSPARSE_ZNAME (_ipvec) 323 CXSPARSE_ZNAME (_ipvec)
324 (q.S()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); 324 (q.S()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
325 #else 325 #else
326 CXSPARSE_ZNAME (_ipvec) 326 CXSPARSE_ZNAME (_ipvec)
327 (nr, q.S()->Pinv, reinterpret_cast<cs_complex_t *>(Xx), buf); 327 (nr, q.S()->Pinv, reinterpret_cast<cs_complex_t *>(Xx), buf);
328 #endif 328 #endif
329 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 329 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
330 for (volatile octave_idx_type j = 0; j < nc; j++) 330 for (volatile octave_idx_type j = 0; j < nc; j++)
331 { 331 {
332 octave_quit (); 332 octave_quit ();
333 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 333 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
334 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 334 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
335 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 335 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
336 } 336 }
337 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 337 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
338 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf); 338 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf);
339 #if defined(CS_VER) && (CS_VER >= 2) 339 #if defined(CS_VER) && (CS_VER >= 2)
340 CXSPARSE_ZNAME (_ipvec) (q.S()->q, buf, vec + idx, nc); 340 CXSPARSE_ZNAME (_ipvec) (q.S()->q, buf, vec + idx, nc);
341 #else 341 #else
342 CXSPARSE_ZNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx); 342 CXSPARSE_ZNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx);
343 #endif 343 #endif
344 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 344 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
345 } 345 }
346 info = 0; 346 info = 0;
347 } 347 }
348 else 348 else
349 { 349 {
350 SparseComplexMatrix at = a.hermitian(); 350 SparseComplexMatrix at = a.hermitian();
351 SparseComplexQR q (at, 2); 351 SparseComplexQR q (at, 2);
352 if (! q.ok ()) 352 if (! q.ok ())
353 return ComplexMatrix(); 353 return ComplexMatrix();
354 x.resize(nc, b_nc); 354 x.resize(nc, b_nc);
355 cs_complex_t *vec = reinterpret_cast<cs_complex_t *> 355 cs_complex_t *vec = reinterpret_cast<cs_complex_t *>
356 (x.fortran_vec()); 356 (x.fortran_vec());
357 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 357 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
358 OCTAVE_C99_COMPLEX (buf, nbuf); 358 OCTAVE_C99_COMPLEX (buf, nbuf);
359 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); 359 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
360 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2)) 360 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
361 OCTAVE_LOCAL_BUFFER (double, B, nr); 361 OCTAVE_LOCAL_BUFFER (double, B, nr);
362 for (octave_idx_type i = 0; i < nr; i++) 362 for (octave_idx_type i = 0; i < nr; i++)
363 B[i] = q.N()->B [i]; 363 B[i] = q.N()->B [i];
364 #else 364 #else
365 OCTAVE_LOCAL_BUFFER (Complex, B, nr); 365 OCTAVE_LOCAL_BUFFER (Complex, B, nr);
366 for (octave_idx_type i = 0; i < nr; i++) 366 for (octave_idx_type i = 0; i < nr; i++)
367 B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]); 367 B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]);
368 #endif 368 #endif
369 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 369 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
370 { 370 {
371 octave_quit (); 371 octave_quit ();
372 for (octave_idx_type j = 0; j < b_nr; j++) 372 for (octave_idx_type j = 0; j < b_nr; j++)
373 Xx[j] = b.xelem(j,i); 373 Xx[j] = b.xelem(j,i);
374 for (octave_idx_type j = nr; j < nbuf; j++) 374 for (octave_idx_type j = nr; j < nbuf; j++)
375 buf[j] = OCTAVE_C99_ZERO; 375 buf[j] = OCTAVE_C99_ZERO;
376 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 376 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
377 #if defined(CS_VER) && (CS_VER >= 2) 377 #if defined(CS_VER) && (CS_VER >= 2)
378 CXSPARSE_ZNAME (_pvec) 378 CXSPARSE_ZNAME (_pvec)
379 (q.S()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); 379 (q.S()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
380 #else 380 #else
381 CXSPARSE_ZNAME (_pvec) 381 CXSPARSE_ZNAME (_pvec)
382 (nr, q.S()->Q, reinterpret_cast<cs_complex_t *>(Xx), buf); 382 (nr, q.S()->Q, reinterpret_cast<cs_complex_t *>(Xx), buf);
383 #endif 383 #endif
384 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf); 384 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf);
385 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 385 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
386 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 386 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
387 { 387 {
388 octave_quit (); 388 octave_quit ();
389 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 389 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
390 390
391 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2)) 391 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
392 CXSPARSE_ZNAME (_happly) (q.N()->L, j, B[j], buf); 392 CXSPARSE_ZNAME (_happly) (q.N()->L, j, B[j], buf);
393 #else 393 #else
394 CXSPARSE_ZNAME (_happly) 394 CXSPARSE_ZNAME (_happly)
395 (q.N()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf); 395 (q.N()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf);
396 #endif 396 #endif
397 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 397 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
398 } 398 }
399 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 399 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
400 #if defined(CS_VER) && (CS_VER >= 2) 400 #if defined(CS_VER) && (CS_VER >= 2)
401 CXSPARSE_ZNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc); 401 CXSPARSE_ZNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc);
402 #else 402 #else
403 CXSPARSE_ZNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx); 403 CXSPARSE_ZNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx);
404 #endif 404 #endif
405 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 405 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
406 } 406 }
407 info = 0; 407 info = 0;
408 } 408 }
409 409
410 return x; 410 return x;
411 #else 411 #else
432 x = SparseComplexMatrix (nc, b_nc); 432 x = SparseComplexMatrix (nc, b_nc);
433 else if (nr >= nc) 433 else if (nr >= nc)
434 { 434 {
435 SparseComplexQR q (a, 2); 435 SparseComplexQR q (a, 2);
436 if (! q.ok ()) 436 if (! q.ok ())
437 return SparseComplexMatrix(); 437 return SparseComplexMatrix();
438 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); 438 x = SparseComplexMatrix (nc, b_nc, b.nzmax());
439 x.xcidx(0) = 0; 439 x.xcidx(0) = 0;
440 x_nz = b.nzmax(); 440 x_nz = b.nzmax();
441 ii = 0; 441 ii = 0;
442 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); 442 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
443 OCTAVE_C99_COMPLEX (buf, q.S()->m2); 443 OCTAVE_C99_COMPLEX (buf, q.S()->m2);
444 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 444 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
445 { 445 {
446 octave_quit (); 446 octave_quit ();
447 for (octave_idx_type j = 0; j < b_nr; j++) 447 for (octave_idx_type j = 0; j < b_nr; j++)
448 Xx[j] = b.xelem(j,i); 448 Xx[j] = b.xelem(j,i);
449 for (octave_idx_type j = nr; j < q.S()->m2; j++) 449 for (octave_idx_type j = nr; j < q.S()->m2; j++)
450 buf[j] = OCTAVE_C99_ZERO; 450 buf[j] = OCTAVE_C99_ZERO;
451 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 451 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
452 #if defined(CS_VER) && (CS_VER >= 2) 452 #if defined(CS_VER) && (CS_VER >= 2)
453 CXSPARSE_ZNAME (_ipvec) 453 CXSPARSE_ZNAME (_ipvec)
454 (q.S()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); 454 (q.S()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
455 #else 455 #else
456 CXSPARSE_ZNAME (_ipvec) 456 CXSPARSE_ZNAME (_ipvec)
457 (nr, q.S()->Pinv, reinterpret_cast<cs_complex_t *>(Xx), buf); 457 (nr, q.S()->Pinv, reinterpret_cast<cs_complex_t *>(Xx), buf);
458 #endif 458 #endif
459 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 459 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
460 for (volatile octave_idx_type j = 0; j < nc; j++) 460 for (volatile octave_idx_type j = 0; j < nc; j++)
461 { 461 {
462 octave_quit (); 462 octave_quit ();
463 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 463 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
464 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 464 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
465 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 465 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
466 } 466 }
467 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 467 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
468 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf); 468 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf);
469 #if defined(CS_VER) && (CS_VER >= 2) 469 #if defined(CS_VER) && (CS_VER >= 2)
470 CXSPARSE_ZNAME (_ipvec) 470 CXSPARSE_ZNAME (_ipvec)
471 (q.S()->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); 471 (q.S()->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
472 #else 472 #else
473 CXSPARSE_ZNAME (_ipvec) 473 CXSPARSE_ZNAME (_ipvec)
474 (nc, q.S()->Q, buf, reinterpret_cast<cs_complex_t *>(Xx)); 474 (nc, q.S()->Q, buf, reinterpret_cast<cs_complex_t *>(Xx));
475 #endif 475 #endif
476 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 476 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
477 477
478 for (octave_idx_type j = 0; j < nc; j++) 478 for (octave_idx_type j = 0; j < nc; j++)
479 { 479 {
480 Complex tmp = Xx[j]; 480 Complex tmp = Xx[j];
481 if (tmp != 0.0) 481 if (tmp != 0.0)
482 { 482 {
483 if (ii == x_nz) 483 if (ii == x_nz)
484 { 484 {
485 // Resize the sparse matrix 485 // Resize the sparse matrix
486 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 486 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
487 sz = (sz > 10 ? sz : 10) + x_nz; 487 sz = (sz > 10 ? sz : 10) + x_nz;
488 x.change_capacity (sz); 488 x.change_capacity (sz);
489 x_nz = sz; 489 x_nz = sz;
490 } 490 }
491 x.xdata(ii) = tmp; 491 x.xdata(ii) = tmp;
492 x.xridx(ii++) = j; 492 x.xridx(ii++) = j;
493 } 493 }
494 } 494 }
495 x.xcidx(i+1) = ii; 495 x.xcidx(i+1) = ii;
496 } 496 }
497 info = 0; 497 info = 0;
498 } 498 }
499 else 499 else
500 { 500 {
501 SparseComplexMatrix at = a.hermitian(); 501 SparseComplexMatrix at = a.hermitian();
502 SparseComplexQR q (at, 2); 502 SparseComplexQR q (at, 2);
503 if (! q.ok ()) 503 if (! q.ok ())
504 return SparseComplexMatrix(); 504 return SparseComplexMatrix();
505 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); 505 x = SparseComplexMatrix (nc, b_nc, b.nzmax());
506 x.xcidx(0) = 0; 506 x.xcidx(0) = 0;
507 x_nz = b.nzmax(); 507 x_nz = b.nzmax();
508 ii = 0; 508 ii = 0;
509 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 509 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
511 OCTAVE_C99_COMPLEX (buf, nbuf); 511 OCTAVE_C99_COMPLEX (buf, nbuf);
512 512
513 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2)) 513 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
514 OCTAVE_LOCAL_BUFFER (double, B, nr); 514 OCTAVE_LOCAL_BUFFER (double, B, nr);
515 for (octave_idx_type i = 0; i < nr; i++) 515 for (octave_idx_type i = 0; i < nr; i++)
516 B[i] = q.N()->B [i]; 516 B[i] = q.N()->B [i];
517 #else 517 #else
518 OCTAVE_LOCAL_BUFFER (Complex, B, nr); 518 OCTAVE_LOCAL_BUFFER (Complex, B, nr);
519 for (octave_idx_type i = 0; i < nr; i++) 519 for (octave_idx_type i = 0; i < nr; i++)
520 B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]); 520 B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]);
521 #endif 521 #endif
522 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 522 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
523 { 523 {
524 octave_quit (); 524 octave_quit ();
525 for (octave_idx_type j = 0; j < b_nr; j++) 525 for (octave_idx_type j = 0; j < b_nr; j++)
526 Xx[j] = b.xelem(j,i); 526 Xx[j] = b.xelem(j,i);
527 for (octave_idx_type j = nr; j < nbuf; j++) 527 for (octave_idx_type j = nr; j < nbuf; j++)
528 buf[j] = OCTAVE_C99_ZERO; 528 buf[j] = OCTAVE_C99_ZERO;
529 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 529 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
530 #if defined(CS_VER) && (CS_VER >= 2) 530 #if defined(CS_VER) && (CS_VER >= 2)
531 CXSPARSE_ZNAME (_pvec) 531 CXSPARSE_ZNAME (_pvec)
532 (q.S()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); 532 (q.S()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
533 #else 533 #else
534 CXSPARSE_ZNAME (_pvec) 534 CXSPARSE_ZNAME (_pvec)
535 (nr, q.S()->Q, reinterpret_cast<cs_complex_t *>(Xx), buf); 535 (nr, q.S()->Q, reinterpret_cast<cs_complex_t *>(Xx), buf);
536 #endif 536 #endif
537 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf); 537 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf);
538 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 538 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
539 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 539 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
540 { 540 {
541 octave_quit (); 541 octave_quit ();
542 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 542 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
543 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2)) 543 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
544 CXSPARSE_ZNAME (_happly) (q.N()->L, j, B[j], buf); 544 CXSPARSE_ZNAME (_happly) (q.N()->L, j, B[j], buf);
545 #else 545 #else
546 CXSPARSE_ZNAME (_happly) 546 CXSPARSE_ZNAME (_happly)
547 (q.N()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf); 547 (q.N()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf);
548 #endif 548 #endif
549 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 549 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
550 } 550 }
551 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 551 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
552 #if defined(CS_VER) && (CS_VER >= 2) 552 #if defined(CS_VER) && (CS_VER >= 2)
553 CXSPARSE_ZNAME (_pvec) 553 CXSPARSE_ZNAME (_pvec)
554 (q.S()->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); 554 (q.S()->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
555 #else 555 #else
556 CXSPARSE_ZNAME (_pvec) 556 CXSPARSE_ZNAME (_pvec)
557 (nc, q.S()->Pinv, buf, reinterpret_cast<cs_complex_t *>(Xx)); 557 (nc, q.S()->Pinv, buf, reinterpret_cast<cs_complex_t *>(Xx));
558 #endif 558 #endif
559 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 559 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
560 560
561 for (octave_idx_type j = 0; j < nc; j++) 561 for (octave_idx_type j = 0; j < nc; j++)
562 { 562 {
563 Complex tmp = Xx[j]; 563 Complex tmp = Xx[j];
564 if (tmp != 0.0) 564 if (tmp != 0.0)
565 { 565 {
566 if (ii == x_nz) 566 if (ii == x_nz)
567 { 567 {
568 // Resize the sparse matrix 568 // Resize the sparse matrix
569 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 569 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
570 sz = (sz > 10 ? sz : 10) + x_nz; 570 sz = (sz > 10 ? sz : 10) + x_nz;
571 x.change_capacity (sz); 571 x.change_capacity (sz);
572 x_nz = sz; 572 x_nz = sz;
573 } 573 }
574 x.xdata(ii) = tmp; 574 x.xdata(ii) = tmp;
575 x.xridx(ii++) = j; 575 x.xridx(ii++) = j;
576 } 576 }
577 } 577 }
578 x.xcidx(i+1) = ii; 578 x.xcidx(i+1) = ii;
579 } 579 }
580 info = 0; 580 info = 0;
581 } 581 }
582 582
583 x.maybe_compress (); 583 x.maybe_compress ();
584 return x; 584 return x;
607 x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); 607 x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
608 else if (nr >= nc) 608 else if (nr >= nc)
609 { 609 {
610 SparseComplexQR q (a, 2); 610 SparseComplexQR q (a, 2);
611 if (! q.ok ()) 611 if (! q.ok ())
612 return ComplexMatrix(); 612 return ComplexMatrix();
613 x.resize(nc, b_nc); 613 x.resize(nc, b_nc);
614 cs_complex_t *vec = reinterpret_cast<cs_complex_t *> 614 cs_complex_t *vec = reinterpret_cast<cs_complex_t *>
615 (x.fortran_vec()); 615 (x.fortran_vec());
616 OCTAVE_C99_COMPLEX (buf, q.S()->m2); 616 OCTAVE_C99_COMPLEX (buf, q.S()->m2);
617 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 617 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
618 i++, idx+=nc, bidx+=b_nr) 618 i++, idx+=nc, bidx+=b_nr)
619 { 619 {
620 octave_quit (); 620 octave_quit ();
621 for (octave_idx_type j = nr; j < q.S()->m2; j++) 621 for (octave_idx_type j = nr; j < q.S()->m2; j++)
622 buf[j] = OCTAVE_C99_ZERO; 622 buf[j] = OCTAVE_C99_ZERO;
623 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 623 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
624 #if defined(CS_VER) && (CS_VER >= 2) 624 #if defined(CS_VER) && (CS_VER >= 2)
625 CXSPARSE_ZNAME (_ipvec) (q.S()->pinv, bvec + bidx, buf, nr); 625 CXSPARSE_ZNAME (_ipvec) (q.S()->pinv, bvec + bidx, buf, nr);
626 #else 626 #else
627 CXSPARSE_ZNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf); 627 CXSPARSE_ZNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf);
628 #endif 628 #endif
629 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 629 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
630 for (volatile octave_idx_type j = 0; j < nc; j++) 630 for (volatile octave_idx_type j = 0; j < nc; j++)
631 { 631 {
632 octave_quit (); 632 octave_quit ();
633 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 633 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
634 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 634 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
635 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 635 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
636 } 636 }
637 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 637 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
638 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf); 638 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf);
639 #if defined(CS_VER) && (CS_VER >= 2) 639 #if defined(CS_VER) && (CS_VER >= 2)
640 CXSPARSE_ZNAME (_ipvec) (q.S()->q, buf, vec + idx, nc); 640 CXSPARSE_ZNAME (_ipvec) (q.S()->q, buf, vec + idx, nc);
641 #else 641 #else
642 CXSPARSE_ZNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx); 642 CXSPARSE_ZNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx);
643 #endif 643 #endif
644 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 644 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
645 } 645 }
646 info = 0; 646 info = 0;
647 } 647 }
648 else 648 else
649 { 649 {
650 SparseComplexMatrix at = a.hermitian(); 650 SparseComplexMatrix at = a.hermitian();
651 SparseComplexQR q (at, 2); 651 SparseComplexQR q (at, 2);
652 if (! q.ok ()) 652 if (! q.ok ())
653 return ComplexMatrix(); 653 return ComplexMatrix();
654 x.resize(nc, b_nc); 654 x.resize(nc, b_nc);
655 cs_complex_t *vec = reinterpret_cast<cs_complex_t *> 655 cs_complex_t *vec = reinterpret_cast<cs_complex_t *>
656 (x.fortran_vec()); 656 (x.fortran_vec());
657 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 657 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
658 OCTAVE_C99_COMPLEX (buf, nbuf); 658 OCTAVE_C99_COMPLEX (buf, nbuf);
659 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2)) 659 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
660 OCTAVE_LOCAL_BUFFER (double, B, nr); 660 OCTAVE_LOCAL_BUFFER (double, B, nr);
661 for (octave_idx_type i = 0; i < nr; i++) 661 for (octave_idx_type i = 0; i < nr; i++)
662 B[i] = q.N()->B [i]; 662 B[i] = q.N()->B [i];
663 #else 663 #else
664 OCTAVE_LOCAL_BUFFER (Complex, B, nr); 664 OCTAVE_LOCAL_BUFFER (Complex, B, nr);
665 for (octave_idx_type i = 0; i < nr; i++) 665 for (octave_idx_type i = 0; i < nr; i++)
666 B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]); 666 B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]);
667 #endif 667 #endif
668 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 668 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
669 i++, idx+=nc, bidx+=b_nr) 669 i++, idx+=nc, bidx+=b_nr)
670 { 670 {
671 octave_quit (); 671 octave_quit ();
672 for (octave_idx_type j = nr; j < nbuf; j++) 672 for (octave_idx_type j = nr; j < nbuf; j++)
673 buf[j] = OCTAVE_C99_ZERO; 673 buf[j] = OCTAVE_C99_ZERO;
674 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 674 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
675 #if defined(CS_VER) && (CS_VER >= 2) 675 #if defined(CS_VER) && (CS_VER >= 2)
676 CXSPARSE_ZNAME (_pvec) (q.S()->q, bvec + bidx, buf, nr); 676 CXSPARSE_ZNAME (_pvec) (q.S()->q, bvec + bidx, buf, nr);
677 #else 677 #else
678 CXSPARSE_ZNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf); 678 CXSPARSE_ZNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf);
679 #endif 679 #endif
680 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf); 680 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf);
681 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 681 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
682 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 682 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
683 { 683 {
684 octave_quit (); 684 octave_quit ();
685 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 685 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
686 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2)) 686 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
687 CXSPARSE_ZNAME (_happly) (q.N()->L, j, B[j], buf); 687 CXSPARSE_ZNAME (_happly) (q.N()->L, j, B[j], buf);
688 #else 688 #else
689 CXSPARSE_ZNAME (_happly) 689 CXSPARSE_ZNAME (_happly)
690 (q.N()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf); 690 (q.N()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf);
691 #endif 691 #endif
692 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 692 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
693 } 693 }
694 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 694 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
695 #if defined(CS_VER) && (CS_VER >= 2) 695 #if defined(CS_VER) && (CS_VER >= 2)
696 CXSPARSE_ZNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc); 696 CXSPARSE_ZNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc);
697 #else 697 #else
698 CXSPARSE_ZNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx); 698 CXSPARSE_ZNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx);
699 #endif 699 #endif
700 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 700 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
701 } 701 }
702 info = 0; 702 info = 0;
703 } 703 }
704 704
705 return x; 705 return x;
706 #else 706 #else
727 x = SparseComplexMatrix (nc, b_nc); 727 x = SparseComplexMatrix (nc, b_nc);
728 else if (nr >= nc) 728 else if (nr >= nc)
729 { 729 {
730 SparseComplexQR q (a, 2); 730 SparseComplexQR q (a, 2);
731 if (! q.ok ()) 731 if (! q.ok ())
732 return SparseComplexMatrix(); 732 return SparseComplexMatrix();
733 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); 733 x = SparseComplexMatrix (nc, b_nc, b.nzmax());
734 x.xcidx(0) = 0; 734 x.xcidx(0) = 0;
735 x_nz = b.nzmax(); 735 x_nz = b.nzmax();
736 ii = 0; 736 ii = 0;
737 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); 737 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
738 OCTAVE_C99_COMPLEX (buf, q.S()->m2); 738 OCTAVE_C99_COMPLEX (buf, q.S()->m2);
739 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 739 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
740 { 740 {
741 octave_quit (); 741 octave_quit ();
742 for (octave_idx_type j = 0; j < b_nr; j++) 742 for (octave_idx_type j = 0; j < b_nr; j++)
743 Xx[j] = b.xelem(j,i); 743 Xx[j] = b.xelem(j,i);
744 for (octave_idx_type j = nr; j < q.S()->m2; j++) 744 for (octave_idx_type j = nr; j < q.S()->m2; j++)
745 buf[j] = OCTAVE_C99_ZERO; 745 buf[j] = OCTAVE_C99_ZERO;
746 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 746 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
747 #if defined(CS_VER) && (CS_VER >= 2) 747 #if defined(CS_VER) && (CS_VER >= 2)
748 CXSPARSE_ZNAME (_ipvec) 748 CXSPARSE_ZNAME (_ipvec)
749 (q.S()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); 749 (q.S()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
750 #else 750 #else
751 CXSPARSE_ZNAME (_ipvec) 751 CXSPARSE_ZNAME (_ipvec)
752 (nr, q.S()->Pinv, reinterpret_cast<cs_complex_t *>(Xx), buf); 752 (nr, q.S()->Pinv, reinterpret_cast<cs_complex_t *>(Xx), buf);
753 #endif 753 #endif
754 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 754 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
755 for (volatile octave_idx_type j = 0; j < nc; j++) 755 for (volatile octave_idx_type j = 0; j < nc; j++)
756 { 756 {
757 octave_quit (); 757 octave_quit ();
758 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 758 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
759 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 759 CXSPARSE_ZNAME (_happly) (q.N()->L, j, q.N()->B[j], buf);
760 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 760 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
761 } 761 }
762 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 762 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
763 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf); 763 CXSPARSE_ZNAME (_usolve) (q.N()->U, buf);
764 #if defined(CS_VER) && (CS_VER >= 2) 764 #if defined(CS_VER) && (CS_VER >= 2)
765 CXSPARSE_ZNAME (_ipvec) 765 CXSPARSE_ZNAME (_ipvec)
766 (q.S()->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); 766 (q.S()->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
767 #else 767 #else
768 CXSPARSE_ZNAME (_ipvec) 768 CXSPARSE_ZNAME (_ipvec)
769 (nc, q.S()->Q, buf, reinterpret_cast<cs_complex_t *>(Xx)); 769 (nc, q.S()->Q, buf, reinterpret_cast<cs_complex_t *>(Xx));
770 #endif 770 #endif
771 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 771 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
772 772
773 for (octave_idx_type j = 0; j < nc; j++) 773 for (octave_idx_type j = 0; j < nc; j++)
774 { 774 {
775 Complex tmp = Xx[j]; 775 Complex tmp = Xx[j];
776 if (tmp != 0.0) 776 if (tmp != 0.0)
777 { 777 {
778 if (ii == x_nz) 778 if (ii == x_nz)
779 { 779 {
780 // Resize the sparse matrix 780 // Resize the sparse matrix
781 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 781 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
782 sz = (sz > 10 ? sz : 10) + x_nz; 782 sz = (sz > 10 ? sz : 10) + x_nz;
783 x.change_capacity (sz); 783 x.change_capacity (sz);
784 x_nz = sz; 784 x_nz = sz;
785 } 785 }
786 x.xdata(ii) = tmp; 786 x.xdata(ii) = tmp;
787 x.xridx(ii++) = j; 787 x.xridx(ii++) = j;
788 } 788 }
789 } 789 }
790 x.xcidx(i+1) = ii; 790 x.xcidx(i+1) = ii;
791 } 791 }
792 info = 0; 792 info = 0;
793 } 793 }
794 else 794 else
795 { 795 {
796 SparseComplexMatrix at = a.hermitian(); 796 SparseComplexMatrix at = a.hermitian();
797 SparseComplexQR q (at, 2); 797 SparseComplexQR q (at, 2);
798 if (! q.ok ()) 798 if (! q.ok ())
799 return SparseComplexMatrix(); 799 return SparseComplexMatrix();
800 x = SparseComplexMatrix (nc, b_nc, b.nzmax()); 800 x = SparseComplexMatrix (nc, b_nc, b.nzmax());
801 x.xcidx(0) = 0; 801 x.xcidx(0) = 0;
802 x_nz = b.nzmax(); 802 x_nz = b.nzmax();
803 ii = 0; 803 ii = 0;
804 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 804 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
805 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); 805 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
806 OCTAVE_C99_COMPLEX (buf, nbuf); 806 OCTAVE_C99_COMPLEX (buf, nbuf);
807 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2)) 807 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
808 OCTAVE_LOCAL_BUFFER (double, B, nr); 808 OCTAVE_LOCAL_BUFFER (double, B, nr);
809 for (octave_idx_type i = 0; i < nr; i++) 809 for (octave_idx_type i = 0; i < nr; i++)
810 B[i] = q.N()->B [i]; 810 B[i] = q.N()->B [i];
811 #else 811 #else
812 OCTAVE_LOCAL_BUFFER (Complex, B, nr); 812 OCTAVE_LOCAL_BUFFER (Complex, B, nr);
813 for (octave_idx_type i = 0; i < nr; i++) 813 for (octave_idx_type i = 0; i < nr; i++)
814 B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]); 814 B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]);
815 #endif 815 #endif
816 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 816 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
817 { 817 {
818 octave_quit (); 818 octave_quit ();
819 for (octave_idx_type j = 0; j < b_nr; j++) 819 for (octave_idx_type j = 0; j < b_nr; j++)
820 Xx[j] = b.xelem(j,i); 820 Xx[j] = b.xelem(j,i);
821 for (octave_idx_type j = nr; j < nbuf; j++) 821 for (octave_idx_type j = nr; j < nbuf; j++)
822 buf[j] = OCTAVE_C99_ZERO; 822 buf[j] = OCTAVE_C99_ZERO;
823 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 823 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
824 #if defined(CS_VER) && (CS_VER >= 2) 824 #if defined(CS_VER) && (CS_VER >= 2)
825 CXSPARSE_ZNAME (_pvec) 825 CXSPARSE_ZNAME (_pvec)
826 (q.S()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); 826 (q.S()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
827 #else 827 #else
828 CXSPARSE_ZNAME (_pvec) 828 CXSPARSE_ZNAME (_pvec)
829 (nr, q.S()->Q, reinterpret_cast<cs_complex_t *>(Xx), buf); 829 (nr, q.S()->Q, reinterpret_cast<cs_complex_t *>(Xx), buf);
830 #endif 830 #endif
831 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf); 831 CXSPARSE_ZNAME (_utsolve) (q.N()->U, buf);
832 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 832 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
833 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 833 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
834 { 834 {
835 octave_quit (); 835 octave_quit ();
836 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 836 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
837 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2)) 837 #if defined(CS_VER) && (((CS_VER == 2) && (CS_SUBVER >= 2)) || (CS_VER > 2))
838 CXSPARSE_ZNAME (_happly) (q.N()->L, j, B[j], buf); 838 CXSPARSE_ZNAME (_happly) (q.N()->L, j, B[j], buf);
839 #else 839 #else
840 CXSPARSE_ZNAME (_happly) 840 CXSPARSE_ZNAME (_happly)
841 (q.N()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf); 841 (q.N()->L, j, reinterpret_cast<cs_complex_t *>(B)[j], buf);
842 #endif 842 #endif
843 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 843 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
844 } 844 }
845 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 845 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
846 #if defined(CS_VER) && (CS_VER >= 2) 846 #if defined(CS_VER) && (CS_VER >= 2)
847 CXSPARSE_ZNAME (_pvec) 847 CXSPARSE_ZNAME (_pvec)
848 (q.S()->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); 848 (q.S()->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
849 #else 849 #else
850 CXSPARSE_ZNAME (_pvec) 850 CXSPARSE_ZNAME (_pvec)
851 (nc, q.S()->Pinv, buf, reinterpret_cast<cs_complex_t *>(Xx)); 851 (nc, q.S()->Pinv, buf, reinterpret_cast<cs_complex_t *>(Xx));
852 #endif 852 #endif
853 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 853 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
854 854
855 for (octave_idx_type j = 0; j < nc; j++) 855 for (octave_idx_type j = 0; j < nc; j++)
856 { 856 {
857 Complex tmp = Xx[j]; 857 Complex tmp = Xx[j];
858 if (tmp != 0.0) 858 if (tmp != 0.0)
859 { 859 {
860 if (ii == x_nz) 860 if (ii == x_nz)
861 { 861 {
862 // Resize the sparse matrix 862 // Resize the sparse matrix
863 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 863 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
864 sz = (sz > 10 ? sz : 10) + x_nz; 864 sz = (sz > 10 ? sz : 10) + x_nz;
865 x.change_capacity (sz); 865 x.change_capacity (sz);
866 x_nz = sz; 866 x_nz = sz;
867 } 867 }
868 x.xdata(ii) = tmp; 868 x.xdata(ii) = tmp;
869 x.xridx(ii++) = j; 869 x.xridx(ii++) = j;
870 } 870 }
871 } 871 }
872 x.xcidx(i+1) = ii; 872 x.xcidx(i+1) = ii;
873 } 873 }
874 info = 0; 874 info = 0;
875 } 875 }
876 876
877 x.maybe_compress (); 877 x.maybe_compress ();
878 return x; 878 return x;
881 #endif 881 #endif
882 } 882 }
883 883
884 ComplexMatrix 884 ComplexMatrix
885 qrsolve (const SparseComplexMatrix &a, const MArray2<double> &b, 885 qrsolve (const SparseComplexMatrix &a, const MArray2<double> &b,
886 octave_idx_type &info) 886 octave_idx_type &info)
887 { 887 {
888 return qrsolve (a, Matrix (b), info); 888 return qrsolve (a, Matrix (b), info);
889 } 889 }
890 890
891 ComplexMatrix 891 ComplexMatrix
892 qrsolve (const SparseComplexMatrix &a, const MArray2<Complex> &b, 892 qrsolve (const SparseComplexMatrix &a, const MArray2<Complex> &b,
893 octave_idx_type &info) 893 octave_idx_type &info)
894 { 894 {
895 return qrsolve (a, ComplexMatrix (b), info); 895 return qrsolve (a, ComplexMatrix (b), info);
896 } 896 }