Mercurial > hg > octave-nkf

comparison liboctave/SparseQR.cc @ 14846:460a3c6d8bf1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists. Example: func() => func () * dynamic.txi, func.txi, oop.txi, var.txi, embedded.cc, fortdemo.cc, funcdemo.cc, paramdemo.cc, stringdemo.cc, unwinddemo.cc, Array.cc, Array.h, CColVector.cc, CDiagMatrix.h, CMatrix.cc, CNDArray.cc, CRowVector.cc, CSparse.cc, CmplxGEPBAL.cc, EIG.cc, MSparse.cc, MatrixType.cc, Sparse-op-defs.h, Sparse-perm-op-defs.h, Sparse.cc, Sparse.h, SparseCmplxCHOL.cc, SparseCmplxCHOL.h, SparseCmplxLU.cc, SparseCmplxQR.cc, SparseCmplxQR.h, SparseQR.cc, SparseQR.h, SparsedbleCHOL.cc, SparsedbleCHOL.h, SparsedbleLU.cc, SparsedbleLU.h, base-lu.cc, cmd-hist.cc, dColVector.cc, dDiagMatrix.h, dMatrix.cc, dNDArray.cc, dRowVector.cc, dSparse.cc, dbleCHOL.cc, dbleGEPBAL.cc, dim-vector.cc, eigs-base.cc, f2c-main.c, fCColVector.cc, fCDiagMatrix.h, fCMatrix.cc, fCNDArray.cc, fCRowVector.cc, fCmplxGEPBAL.cc, fColVector.cc, fDiagMatrix.h, fEIG.cc, fMatrix.cc, fNDArray.cc, fRowVector.cc, file-ops.cc, file-stat.cc, floatCHOL.cc, floatGEPBAL.cc, idx-vector.h, lo-specfun.cc, lo-sysdep.cc, mx-inlines.cc, oct-binmap.h, oct-convn.cc, oct-md5.cc, oct-mem.h, oct-rand.cc, oct-syscalls.cc, randgamma.c, randmtzig.c, sparse-base-chol.cc, sparse-base-chol.h, sparse-base-lu.cc, sparse-dmsolve.cc, tempname.c, curl.m, divergence.m, randi.m, dlmwrite.m, edit.m, getappdata.m, what.m, getarchdir.m, install.m, installed_packages.m, repackage.m, unload_packages.m, colorbar.m, figure.m, isosurface.m, legend.m, loglog.m, plot.m, plot3.m, plotyy.m, polar.m, __errplot__.m, __ghostscript__.m, __marching_cube__.m, __plt__.m, __scatter__.m, semilogx.m, semilogy.m, trimesh.m, trisurf.m, demo.m, test.m, datetick.m, __delaunayn__.cc, __dsearchn__.cc, __fltk_uigetfile__.cc, __glpk__.cc, __init_fltk__.cc, __lin_interpn__.cc, __magick_read__.cc, __pchip_deriv__.cc, balance.cc, bsxfun.cc, ccolamd.cc, cellfun.cc, chol.cc, daspk.cc, dasrt.cc, dassl.cc, dmperm.cc, eig.cc, eigs.cc, fftw.cc, filter.cc, find.cc, kron.cc, lookup.cc, lsode.cc, matrix_type.cc, md5sum.cc, mgorth.cc, qr.cc, quad.cc, rand.cc, regexp.cc, symbfact.cc, tril.cc, urlwrite.cc, op-bm-bm.cc, op-cdm-cdm.cc, op-cell.cc, op-chm.cc, op-cm-cm.cc, op-cm-scm.cc, op-cm-sm.cc, op-cs-scm.cc, op-cs-sm.cc, op-dm-dm.cc, op-dm-scm.cc, op-dm-sm.cc, op-fcdm-fcdm.cc, op-fcm-fcm.cc, op-fdm-fdm.cc, op-fm-fm.cc, op-int.h, op-m-m.cc, op-m-scm.cc, op-m-sm.cc, op-pm-pm.cc, op-pm-scm.cc, op-pm-sm.cc, op-range.cc, op-s-scm.cc, op-s-sm.cc, op-sbm-sbm.cc, op-scm-cm.cc, op-scm-cs.cc, op-scm-m.cc, op-scm-s.cc, op-scm-scm.cc, op-scm-sm.cc, op-sm-cm.cc, op-sm-cs.cc, op-sm-m.cc, op-sm-s.cc, op-sm-scm.cc, op-sm-sm.cc, op-str-str.cc, op-struct.cc, bitfcns.cc, data.cc, debug.cc, dynamic-ld.cc, error.cc, gl-render.cc, graphics.cc, graphics.in.h, load-path.cc, ls-hdf5.cc, ls-mat5.cc, ls-mat5.h, ls-oct-ascii.cc, ls-oct-ascii.h, mex.cc, mk-errno-list, oct-map.cc, oct-obj.h, oct-parse.yy, octave-config.in.cc, ov-base-int.cc, ov-base-mat.cc, ov-base.cc, ov-bool-mat.cc, ov-bool-sparse.cc, ov-bool.cc, ov-cell.cc, ov-class.cc, ov-class.h, ov-cx-mat.cc, ov-cx-sparse.cc, ov-fcn-handle.cc, ov-flt-cx-mat.cc, ov-flt-re-mat.cc, ov-intx.h, ov-range.h, ov-re-mat.cc, ov-re-sparse.cc, ov-str-mat.cc, ov-struct.cc, ov-usr-fcn.h, ov.h, pr-output.cc, pt-id.cc, pt-id.h, pt-mat.cc, pt-select.cc, sparse.cc, symtab.cc, symtab.h, syscalls.cc, toplev.cc, txt-eng-ft.cc, variables.cc, zfstream.cc, zfstream.h, Dork.m, getStash.m, myStash.m, Gork.m, Pork.m, myStash.m, getStash.m, myStash.m, getStash.m, myStash.m, fntests.m: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
author Rik <octave@nomad.inbox5.com>
date Sun, 08 Jul 2012 11:28:50 -0700
parents 72c96de7a403
children 3d8ace26c5b4
comparison
equal deleted inserted replaced
14844:5bc9b9cb4362 14846:460a3c6d8bf1
172 172
173 Matrix 173 Matrix
174 SparseQR::SparseQR_rep::C (const Matrix &b) const 174 SparseQR::SparseQR_rep::C (const Matrix &b) const
175 { 175 {
176 #ifdef HAVE_CXSPARSE 176 #ifdef HAVE_CXSPARSE
177 octave_idx_type b_nr = b.rows(); 177 octave_idx_type b_nr = b.rows ();
178 octave_idx_type b_nc = b.cols(); 178 octave_idx_type b_nc = b.cols ();
179 octave_idx_type nc = N->L->n; 179 octave_idx_type nc = N->L->n;
180 octave_idx_type nr = nrows; 180 octave_idx_type nr = nrows;
181 const double *bvec = b.fortran_vec(); 181 const double *bvec = b.fortran_vec ();
182 Matrix ret (b_nr, b_nc); 182 Matrix ret (b_nr, b_nc);
183 double *vec = ret.fortran_vec(); 183 double *vec = ret.fortran_vec ();
184 if (nr < 0 || nc < 0 || nr != b_nr) 184 if (nr < 0 || nc < 0 || nr != b_nr)
185 (*current_liboctave_error_handler) ("matrix dimension mismatch"); 185 (*current_liboctave_error_handler) ("matrix dimension mismatch");
186 else if (nr == 0 || nc == 0 || b_nc == 0) 186 else if (nr == 0 || nc == 0 || b_nc == 0)
187 ret = Matrix (nc, b_nc, 0.0); 187 ret = Matrix (nc, b_nc, 0.0);
188 else 188 else
224 { 224 {
225 #ifdef HAVE_CXSPARSE 225 #ifdef HAVE_CXSPARSE
226 octave_idx_type nc = N->L->n; 226 octave_idx_type nc = N->L->n;
227 octave_idx_type nr = nrows; 227 octave_idx_type nr = nrows;
228 Matrix ret (nr, nr); 228 Matrix ret (nr, nr);
229 double *vec = ret.fortran_vec(); 229 double *vec = ret.fortran_vec ();
230 if (nr < 0 || nc < 0) 230 if (nr < 0 || nc < 0)
231 (*current_liboctave_error_handler) ("matrix dimension mismatch"); 231 (*current_liboctave_error_handler) ("matrix dimension mismatch");
232 else if (nr == 0 || nc == 0) 232 else if (nr == 0 || nc == 0)
233 ret = Matrix (nc, nr, 0.0); 233 ret = Matrix (nc, nr, 0.0);
234 else 234 else
273 Matrix 273 Matrix
274 qrsolve(const SparseMatrix&a, const Matrix &b, octave_idx_type& info) 274 qrsolve(const SparseMatrix&a, const Matrix &b, octave_idx_type& info)
275 { 275 {
276 info = -1; 276 info = -1;
277 #ifdef HAVE_CXSPARSE 277 #ifdef HAVE_CXSPARSE
278 octave_idx_type nr = a.rows(); 278 octave_idx_type nr = a.rows ();
279 octave_idx_type nc = a.cols(); 279 octave_idx_type nc = a.cols ();
280 octave_idx_type b_nc = b.cols(); 280 octave_idx_type b_nc = b.cols ();
281 octave_idx_type b_nr = b.rows(); 281 octave_idx_type b_nr = b.rows ();
282 const double *bvec = b.fortran_vec(); 282 const double *bvec = b.fortran_vec ();
283 Matrix x; 283 Matrix x;
284 284
285 if (nr < 0 || nc < 0 || nr != b_nr) 285 if (nr < 0 || nc < 0 || nr != b_nr)
286 (*current_liboctave_error_handler) 286 (*current_liboctave_error_handler)
287 ("matrix dimension mismatch in solution of minimum norm problem"); 287 ("matrix dimension mismatch in solution of minimum norm problem");
289 x = Matrix (nc, b_nc, 0.0); 289 x = Matrix (nc, b_nc, 0.0);
290 else if (nr >= nc) 290 else if (nr >= nc)
291 { 291 {
292 SparseQR q (a, 3); 292 SparseQR q (a, 3);
293 if (! q.ok ()) 293 if (! q.ok ())
294 return Matrix(); 294 return Matrix ();
295 x.resize(nc, b_nc); 295 x.resize(nc, b_nc);
296 double *vec = x.fortran_vec(); 296 double *vec = x.fortran_vec ();
297 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 297 OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
298 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 298 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
299 i++, idx+=nc, bidx+=b_nr) 299 i++, idx+=nc, bidx+=b_nr)
300 { 300 {
301 octave_quit (); 301 octave_quit ();
302 for (octave_idx_type j = nr; j < q.S()->m2; j++) 302 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
303 buf[j] = 0.; 303 buf[j] = 0.;
304 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 304 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
305 #if defined(CS_VER) && (CS_VER >= 2) 305 #if defined(CS_VER) && (CS_VER >= 2)
306 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, bvec + bidx, buf, nr); 306 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, bvec + bidx, buf, nr);
307 #else 307 #else
308 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf); 308 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, bvec + bidx, buf);
309 #endif 309 #endif
310 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 310 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
311 for (volatile octave_idx_type j = 0; j < nc; j++) 311 for (volatile octave_idx_type j = 0; j < nc; j++)
312 { 312 {
313 octave_quit (); 313 octave_quit ();
314 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 314 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
315 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 315 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
316 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 316 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
317 } 317 }
318 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 318 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
319 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 319 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
320 #if defined(CS_VER) && (CS_VER >= 2) 320 #if defined(CS_VER) && (CS_VER >= 2)
321 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, vec + idx, nc); 321 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc);
322 #else 322 #else
323 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, vec + idx); 323 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, vec + idx);
324 #endif 324 #endif
325 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 325 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
326 } 326 }
327 info = 0; 327 info = 0;
328 } 328 }
329 else 329 else
330 { 330 {
331 SparseMatrix at = a.hermitian(); 331 SparseMatrix at = a.hermitian ();
332 SparseQR q (at, 3); 332 SparseQR q (at, 3);
333 if (! q.ok ()) 333 if (! q.ok ())
334 return Matrix(); 334 return Matrix ();
335 x.resize(nc, b_nc); 335 x.resize(nc, b_nc);
336 double *vec = x.fortran_vec(); 336 double *vec = x.fortran_vec ();
337 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 337 volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
338 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 338 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
339 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 339 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
340 i++, idx+=nc, bidx+=b_nr) 340 i++, idx+=nc, bidx+=b_nr)
341 { 341 {
342 octave_quit (); 342 octave_quit ();
343 for (octave_idx_type j = nr; j < nbuf; j++) 343 for (octave_idx_type j = nr; j < nbuf; j++)
344 buf[j] = 0.; 344 buf[j] = 0.;
345 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 345 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
346 #if defined(CS_VER) && (CS_VER >= 2) 346 #if defined(CS_VER) && (CS_VER >= 2)
347 CXSPARSE_DNAME (_pvec) (q.S()->q, bvec + bidx, buf, nr); 347 CXSPARSE_DNAME (_pvec) (q.S ()->q, bvec + bidx, buf, nr);
348 #else 348 #else
349 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, bvec + bidx, buf); 349 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, bvec + bidx, buf);
350 #endif 350 #endif
351 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 351 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
352 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 352 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
353 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 353 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
354 { 354 {
355 octave_quit (); 355 octave_quit ();
356 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 356 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
357 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 357 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
358 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 358 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
359 } 359 }
360 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 360 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
361 #if defined(CS_VER) && (CS_VER >= 2) 361 #if defined(CS_VER) && (CS_VER >= 2)
362 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, vec + idx, nc); 362 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc);
363 #else 363 #else
364 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, vec + idx); 364 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, vec + idx);
365 #endif 365 #endif
366 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 366 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
367 } 367 }
368 info = 0; 368 info = 0;
369 } 369 }
377 SparseMatrix 377 SparseMatrix
378 qrsolve(const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info) 378 qrsolve(const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info)
379 { 379 {
380 info = -1; 380 info = -1;
381 #ifdef HAVE_CXSPARSE 381 #ifdef HAVE_CXSPARSE
382 octave_idx_type nr = a.rows(); 382 octave_idx_type nr = a.rows ();
383 octave_idx_type nc = a.cols(); 383 octave_idx_type nc = a.cols ();
384 octave_idx_type b_nr = b.rows(); 384 octave_idx_type b_nr = b.rows ();
385 octave_idx_type b_nc = b.cols(); 385 octave_idx_type b_nc = b.cols ();
386 SparseMatrix x; 386 SparseMatrix x;
387 volatile octave_idx_type ii, x_nz; 387 volatile octave_idx_type ii, x_nz;
388 388
389 if (nr < 0 || nc < 0 || nr != b_nr) 389 if (nr < 0 || nc < 0 || nr != b_nr)
390 (*current_liboctave_error_handler) 390 (*current_liboctave_error_handler)
393 x = SparseMatrix (nc, b_nc); 393 x = SparseMatrix (nc, b_nc);
394 else if (nr >= nc) 394 else if (nr >= nc)
395 { 395 {
396 SparseQR q (a, 3); 396 SparseQR q (a, 3);
397 if (! q.ok ()) 397 if (! q.ok ())
398 return SparseMatrix(); 398 return SparseMatrix ();
399 x = SparseMatrix (nc, b_nc, b.nnz()); 399 x = SparseMatrix (nc, b_nc, b.nnz ());
400 x.xcidx(0) = 0; 400 x.xcidx(0) = 0;
401 x_nz = b.nnz(); 401 x_nz = b.nnz ();
402 ii = 0; 402 ii = 0;
403 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 403 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
404 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 404 OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
405 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 405 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
406 { 406 {
407 octave_quit (); 407 octave_quit ();
408 for (octave_idx_type j = 0; j < b_nr; j++) 408 for (octave_idx_type j = 0; j < b_nr; j++)
409 Xx[j] = b.xelem(j,i); 409 Xx[j] = b.xelem(j,i);
410 for (octave_idx_type j = nr; j < q.S()->m2; j++) 410 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
411 buf[j] = 0.; 411 buf[j] = 0.;
412 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 412 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
413 #if defined(CS_VER) && (CS_VER >= 2) 413 #if defined(CS_VER) && (CS_VER >= 2)
414 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); 414 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
415 #else 415 #else
416 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); 416 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xx, buf);
417 #endif 417 #endif
418 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 418 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
419 for (volatile octave_idx_type j = 0; j < nc; j++) 419 for (volatile octave_idx_type j = 0; j < nc; j++)
420 { 420 {
421 octave_quit (); 421 octave_quit ();
422 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 422 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
423 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 423 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
424 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 424 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
425 } 425 }
426 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 426 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
427 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 427 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
428 #if defined(CS_VER) && (CS_VER >= 2) 428 #if defined(CS_VER) && (CS_VER >= 2)
429 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); 429 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
430 #else 430 #else
431 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); 431 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xx);
432 #endif 432 #endif
433 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 433 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
434 434
435 for (octave_idx_type j = 0; j < nc; j++) 435 for (octave_idx_type j = 0; j < nc; j++)
436 { 436 {
453 } 453 }
454 info = 0; 454 info = 0;
455 } 455 }
456 else 456 else
457 { 457 {
458 SparseMatrix at = a.hermitian(); 458 SparseMatrix at = a.hermitian ();
459 SparseQR q (at, 3); 459 SparseQR q (at, 3);
460 if (! q.ok ()) 460 if (! q.ok ())
461 return SparseMatrix(); 461 return SparseMatrix ();
462 x = SparseMatrix (nc, b_nc, b.nnz()); 462 x = SparseMatrix (nc, b_nc, b.nnz ());
463 x.xcidx(0) = 0; 463 x.xcidx(0) = 0;
464 x_nz = b.nnz(); 464 x_nz = b.nnz ();
465 ii = 0; 465 ii = 0;
466 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 466 volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
467 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 467 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
468 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 468 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
469 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 469 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
470 { 470 {
471 octave_quit (); 471 octave_quit ();
473 Xx[j] = b.xelem(j,i); 473 Xx[j] = b.xelem(j,i);
474 for (octave_idx_type j = nr; j < nbuf; j++) 474 for (octave_idx_type j = nr; j < nbuf; j++)
475 buf[j] = 0.; 475 buf[j] = 0.;
476 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 476 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
477 #if defined(CS_VER) && (CS_VER >= 2) 477 #if defined(CS_VER) && (CS_VER >= 2)
478 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); 478 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
479 #else 479 #else
480 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); 480 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xx, buf);
481 #endif 481 #endif
482 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 482 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
483 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 483 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
484 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 484 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
485 { 485 {
486 octave_quit (); 486 octave_quit ();
487 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 487 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
488 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 488 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
489 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 489 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
490 } 490 }
491 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 491 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
492 #if defined(CS_VER) && (CS_VER >= 2) 492 #if defined(CS_VER) && (CS_VER >= 2)
493 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); 493 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
494 #else 494 #else
495 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); 495 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xx);
496 #endif 496 #endif
497 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 497 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
498 498
499 for (octave_idx_type j = 0; j < nc; j++) 499 for (octave_idx_type j = 0; j < nc; j++)
500 { 500 {
528 ComplexMatrix 528 ComplexMatrix
529 qrsolve(const SparseMatrix&a, const ComplexMatrix &b, octave_idx_type &info) 529 qrsolve(const SparseMatrix&a, const ComplexMatrix &b, octave_idx_type &info)
530 { 530 {
531 info = -1; 531 info = -1;
532 #ifdef HAVE_CXSPARSE 532 #ifdef HAVE_CXSPARSE
533 octave_idx_type nr = a.rows(); 533 octave_idx_type nr = a.rows ();
534 octave_idx_type nc = a.cols(); 534 octave_idx_type nc = a.cols ();
535 octave_idx_type b_nc = b.cols(); 535 octave_idx_type b_nc = b.cols ();
536 octave_idx_type b_nr = b.rows(); 536 octave_idx_type b_nr = b.rows ();
537 ComplexMatrix x; 537 ComplexMatrix x;
538 538
539 if (nr < 0 || nc < 0 || nr != b_nr) 539 if (nr < 0 || nc < 0 || nr != b_nr)
540 (*current_liboctave_error_handler) 540 (*current_liboctave_error_handler)
541 ("matrix dimension mismatch in solution of minimum norm problem"); 541 ("matrix dimension mismatch in solution of minimum norm problem");
543 x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); 543 x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
544 else if (nr >= nc) 544 else if (nr >= nc)
545 { 545 {
546 SparseQR q (a, 3); 546 SparseQR q (a, 3);
547 if (! q.ok ()) 547 if (! q.ok ())
548 return ComplexMatrix(); 548 return ComplexMatrix ();
549 x.resize(nc, b_nc); 549 x.resize(nc, b_nc);
550 Complex *vec = x.fortran_vec(); 550 Complex *vec = x.fortran_vec ();
551 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 551 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
552 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 552 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
553 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 553 OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
554 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 554 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
555 { 555 {
556 octave_quit (); 556 octave_quit ();
557 for (octave_idx_type j = 0; j < b_nr; j++) 557 for (octave_idx_type j = 0; j < b_nr; j++)
558 { 558 {
559 Complex c = b.xelem (j,i); 559 Complex c = b.xelem (j,i);
560 Xx[j] = std::real (c); 560 Xx[j] = std::real (c);
561 Xz[j] = std::imag (c); 561 Xz[j] = std::imag (c);
562 } 562 }
563 for (octave_idx_type j = nr; j < q.S()->m2; j++) 563 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
564 buf[j] = 0.; 564 buf[j] = 0.;
565 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 565 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
566 #if defined(CS_VER) && (CS_VER >= 2) 566 #if defined(CS_VER) && (CS_VER >= 2)
567 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); 567 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
568 #else 568 #else
569 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); 569 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xx, buf);
570 #endif 570 #endif
571 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 571 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
572 for (volatile octave_idx_type j = 0; j < nc; j++) 572 for (volatile octave_idx_type j = 0; j < nc; j++)
573 { 573 {
574 octave_quit (); 574 octave_quit ();
575 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 575 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
576 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 576 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
577 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 577 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
578 } 578 }
579 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 579 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
580 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 580 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
581 #if defined(CS_VER) && (CS_VER >= 2) 581 #if defined(CS_VER) && (CS_VER >= 2)
582 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); 582 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
583 #else 583 #else
584 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); 584 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xx);
585 #endif 585 #endif
586 for (octave_idx_type j = nr; j < q.S()->m2; j++) 586 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
587 buf[j] = 0.; 587 buf[j] = 0.;
588 #if defined(CS_VER) && (CS_VER >= 2) 588 #if defined(CS_VER) && (CS_VER >= 2)
589 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr); 589 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr);
590 #else 590 #else
591 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf); 591 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xz, buf);
592 #endif 592 #endif
593 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 593 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
594 for (volatile octave_idx_type j = 0; j < nc; j++) 594 for (volatile octave_idx_type j = 0; j < nc; j++)
595 { 595 {
596 octave_quit (); 596 octave_quit ();
597 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 597 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
598 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 598 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
599 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 599 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
600 } 600 }
601 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 601 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
602 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 602 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
603 #if defined(CS_VER) && (CS_VER >= 2) 603 #if defined(CS_VER) && (CS_VER >= 2)
604 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc); 604 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc);
605 #else 605 #else
606 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz); 606 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xz);
607 #endif 607 #endif
608 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 608 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
609 for (octave_idx_type j = 0; j < nc; j++) 609 for (octave_idx_type j = 0; j < nc; j++)
610 vec[j+idx] = Complex (Xx[j], Xz[j]); 610 vec[j+idx] = Complex (Xx[j], Xz[j]);
611 } 611 }
612 info = 0; 612 info = 0;
613 } 613 }
614 else 614 else
615 { 615 {
616 SparseMatrix at = a.hermitian(); 616 SparseMatrix at = a.hermitian ();
617 SparseQR q (at, 3); 617 SparseQR q (at, 3);
618 if (! q.ok ()) 618 if (! q.ok ())
619 return ComplexMatrix(); 619 return ComplexMatrix ();
620 x.resize(nc, b_nc); 620 x.resize(nc, b_nc);
621 Complex *vec = x.fortran_vec(); 621 Complex *vec = x.fortran_vec ();
622 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 622 volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
623 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 623 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
624 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 624 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
625 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 625 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
626 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 626 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
627 { 627 {
634 } 634 }
635 for (octave_idx_type j = nr; j < nbuf; j++) 635 for (octave_idx_type j = nr; j < nbuf; j++)
636 buf[j] = 0.; 636 buf[j] = 0.;
637 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 637 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
638 #if defined(CS_VER) && (CS_VER >= 2) 638 #if defined(CS_VER) && (CS_VER >= 2)
639 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); 639 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
640 #else 640 #else
641 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); 641 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xx, buf);
642 #endif 642 #endif
643 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 643 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
644 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 644 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
645 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 645 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
646 { 646 {
647 octave_quit (); 647 octave_quit ();
648 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 648 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
649 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 649 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
650 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 650 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
651 } 651 }
652 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 652 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
653 #if defined(CS_VER) && (CS_VER >= 2) 653 #if defined(CS_VER) && (CS_VER >= 2)
654 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); 654 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
655 #else 655 #else
656 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); 656 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xx);
657 #endif 657 #endif
658 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 658 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
659 for (octave_idx_type j = nr; j < nbuf; j++) 659 for (octave_idx_type j = nr; j < nbuf; j++)
660 buf[j] = 0.; 660 buf[j] = 0.;
661 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 661 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
662 #if defined(CS_VER) && (CS_VER >= 2) 662 #if defined(CS_VER) && (CS_VER >= 2)
663 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr); 663 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr);
664 #else 664 #else
665 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf); 665 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xz, buf);
666 #endif 666 #endif
667 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 667 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
668 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 668 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
669 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 669 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
670 { 670 {
671 octave_quit (); 671 octave_quit ();
672 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 672 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
673 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 673 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
674 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 674 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
675 } 675 }
676 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 676 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
677 #if defined(CS_VER) && (CS_VER >= 2) 677 #if defined(CS_VER) && (CS_VER >= 2)
678 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc); 678 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc);
679 #else 679 #else
680 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz); 680 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xz);
681 #endif 681 #endif
682 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 682 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
683 for (octave_idx_type j = 0; j < nc; j++) 683 for (octave_idx_type j = 0; j < nc; j++)
684 vec[j+idx] = Complex (Xx[j], Xz[j]); 684 vec[j+idx] = Complex (Xx[j], Xz[j]);
685 } 685 }
695 SparseComplexMatrix 695 SparseComplexMatrix
696 qrsolve(const SparseMatrix&a, const SparseComplexMatrix &b, octave_idx_type &info) 696 qrsolve(const SparseMatrix&a, const SparseComplexMatrix &b, octave_idx_type &info)
697 { 697 {
698 info = -1; 698 info = -1;
699 #ifdef HAVE_CXSPARSE 699 #ifdef HAVE_CXSPARSE
700 octave_idx_type nr = a.rows(); 700 octave_idx_type nr = a.rows ();
701 octave_idx_type nc = a.cols(); 701 octave_idx_type nc = a.cols ();
702 octave_idx_type b_nr = b.rows(); 702 octave_idx_type b_nr = b.rows ();
703 octave_idx_type b_nc = b.cols(); 703 octave_idx_type b_nc = b.cols ();
704 SparseComplexMatrix x; 704 SparseComplexMatrix x;
705 volatile octave_idx_type ii, x_nz; 705 volatile octave_idx_type ii, x_nz;
706 706
707 if (nr < 0 || nc < 0 || nr != b_nr) 707 if (nr < 0 || nc < 0 || nr != b_nr)
708 (*current_liboctave_error_handler) 708 (*current_liboctave_error_handler)
711 x = SparseComplexMatrix (nc, b_nc); 711 x = SparseComplexMatrix (nc, b_nc);
712 else if (nr >= nc) 712 else if (nr >= nc)
713 { 713 {
714 SparseQR q (a, 3); 714 SparseQR q (a, 3);
715 if (! q.ok ()) 715 if (! q.ok ())
716 return SparseComplexMatrix(); 716 return SparseComplexMatrix ();
717 x = SparseComplexMatrix (nc, b_nc, b.nnz()); 717 x = SparseComplexMatrix (nc, b_nc, b.nnz ());
718 x.xcidx(0) = 0; 718 x.xcidx(0) = 0;
719 x_nz = b.nnz(); 719 x_nz = b.nnz ();
720 ii = 0; 720 ii = 0;
721 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 721 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
722 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 722 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
723 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 723 OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
724 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 724 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
725 { 725 {
726 octave_quit (); 726 octave_quit ();
727 for (octave_idx_type j = 0; j < b_nr; j++) 727 for (octave_idx_type j = 0; j < b_nr; j++)
728 { 728 {
729 Complex c = b.xelem (j,i); 729 Complex c = b.xelem (j,i);
730 Xx[j] = std::real (c); 730 Xx[j] = std::real (c);
731 Xz[j] = std::imag (c); 731 Xz[j] = std::imag (c);
732 } 732 }
733 for (octave_idx_type j = nr; j < q.S()->m2; j++) 733 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
734 buf[j] = 0.; 734 buf[j] = 0.;
735 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 735 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
736 #if defined(CS_VER) && (CS_VER >= 2) 736 #if defined(CS_VER) && (CS_VER >= 2)
737 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xx, buf, nr); 737 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
738 #else 738 #else
739 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xx, buf); 739 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xx, buf);
740 #endif 740 #endif
741 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 741 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
742 for (volatile octave_idx_type j = 0; j < nc; j++) 742 for (volatile octave_idx_type j = 0; j < nc; j++)
743 { 743 {
744 octave_quit (); 744 octave_quit ();
745 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 745 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
746 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 746 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
747 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 747 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
748 } 748 }
749 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 749 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
750 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 750 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
751 #if defined(CS_VER) && (CS_VER >= 2) 751 #if defined(CS_VER) && (CS_VER >= 2)
752 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xx, nc); 752 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
753 #else 753 #else
754 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xx); 754 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xx);
755 #endif 755 #endif
756 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 756 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
757 for (octave_idx_type j = nr; j < q.S()->m2; j++) 757 for (octave_idx_type j = nr; j < q.S ()->m2; j++)
758 buf[j] = 0.; 758 buf[j] = 0.;
759 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 759 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
760 #if defined(CS_VER) && (CS_VER >= 2) 760 #if defined(CS_VER) && (CS_VER >= 2)
761 CXSPARSE_DNAME (_ipvec) (q.S()->pinv, Xz, buf, nr); 761 CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr);
762 #else 762 #else
763 CXSPARSE_DNAME (_ipvec) (nr, q.S()->Pinv, Xz, buf); 763 CXSPARSE_DNAME (_ipvec) (nr, q.S ()->Pinv, Xz, buf);
764 #endif 764 #endif
765 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 765 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
766 for (volatile octave_idx_type j = 0; j < nc; j++) 766 for (volatile octave_idx_type j = 0; j < nc; j++)
767 { 767 {
768 octave_quit (); 768 octave_quit ();
769 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 769 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
770 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 770 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
771 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 771 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
772 } 772 }
773 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 773 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
774 CXSPARSE_DNAME (_usolve) (q.N()->U, buf); 774 CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
775 #if defined(CS_VER) && (CS_VER >= 2) 775 #if defined(CS_VER) && (CS_VER >= 2)
776 CXSPARSE_DNAME (_ipvec) (q.S()->q, buf, Xz, nc); 776 CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc);
777 #else 777 #else
778 CXSPARSE_DNAME (_ipvec) (nc, q.S()->Q, buf, Xz); 778 CXSPARSE_DNAME (_ipvec) (nc, q.S ()->Q, buf, Xz);
779 #endif 779 #endif
780 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 780 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
781 781
782 for (octave_idx_type j = 0; j < nc; j++) 782 for (octave_idx_type j = 0; j < nc; j++)
783 { 783 {
800 } 800 }
801 info = 0; 801 info = 0;
802 } 802 }
803 else 803 else
804 { 804 {
805 SparseMatrix at = a.hermitian(); 805 SparseMatrix at = a.hermitian ();
806 SparseQR q (at, 3); 806 SparseQR q (at, 3);
807 if (! q.ok ()) 807 if (! q.ok ())
808 return SparseComplexMatrix(); 808 return SparseComplexMatrix ();
809 x = SparseComplexMatrix (nc, b_nc, b.nnz()); 809 x = SparseComplexMatrix (nc, b_nc, b.nnz ());
810 x.xcidx(0) = 0; 810 x.xcidx(0) = 0;
811 x_nz = b.nnz(); 811 x_nz = b.nnz ();
812 ii = 0; 812 ii = 0;
813 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 813 volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
814 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 814 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
815 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 815 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
816 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 816 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
817 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 817 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
818 { 818 {
825 } 825 }
826 for (octave_idx_type j = nr; j < nbuf; j++) 826 for (octave_idx_type j = nr; j < nbuf; j++)
827 buf[j] = 0.; 827 buf[j] = 0.;
828 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 828 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
829 #if defined(CS_VER) && (CS_VER >= 2) 829 #if defined(CS_VER) && (CS_VER >= 2)
830 CXSPARSE_DNAME (_pvec) (q.S()->q, Xx, buf, nr); 830 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
831 #else 831 #else
832 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xx, buf); 832 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xx, buf);
833 #endif 833 #endif
834 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 834 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
835 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 835 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
836 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 836 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
837 { 837 {
838 octave_quit (); 838 octave_quit ();
839 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 839 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
840 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 840 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
841 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 841 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
842 } 842 }
843 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 843 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
844 #if defined(CS_VER) && (CS_VER >= 2) 844 #if defined(CS_VER) && (CS_VER >= 2)
845 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xx, nc); 845 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
846 #else 846 #else
847 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xx); 847 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xx);
848 #endif 848 #endif
849 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 849 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
850 for (octave_idx_type j = nr; j < nbuf; j++) 850 for (octave_idx_type j = nr; j < nbuf; j++)
851 buf[j] = 0.; 851 buf[j] = 0.;
852 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 852 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
853 #if defined(CS_VER) && (CS_VER >= 2) 853 #if defined(CS_VER) && (CS_VER >= 2)
854 CXSPARSE_DNAME (_pvec) (q.S()->q, Xz, buf, nr); 854 CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr);
855 #else 855 #else
856 CXSPARSE_DNAME (_pvec) (nr, q.S()->Q, Xz, buf); 856 CXSPARSE_DNAME (_pvec) (nr, q.S ()->Q, Xz, buf);
857 #endif 857 #endif
858 CXSPARSE_DNAME (_utsolve) (q.N()->U, buf); 858 CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
859 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 859 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
860 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 860 for (volatile octave_idx_type j = nr-1; j >= 0; j--)
861 { 861 {
862 octave_quit (); 862 octave_quit ();
863 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 863 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
864 CXSPARSE_DNAME (_happly) (q.N()->L, j, q.N()->B[j], buf); 864 CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
865 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 865 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
866 } 866 }
867 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 867 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
868 #if defined(CS_VER) && (CS_VER >= 2) 868 #if defined(CS_VER) && (CS_VER >= 2)
869 CXSPARSE_DNAME (_pvec) (q.S()->pinv, buf, Xz, nc); 869 CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc);
870 #else 870 #else
871 CXSPARSE_DNAME (_pvec) (nc, q.S()->Pinv, buf, Xz); 871 CXSPARSE_DNAME (_pvec) (nc, q.S ()->Pinv, buf, Xz);
872 #endif 872 #endif
873 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 873 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
874 874
875 for (octave_idx_type j = 0; j < nc; j++) 875 for (octave_idx_type j = 0; j < nc; j++)
876 { 876 {