Mercurial > hg > octave-lyh

comparison liboctave/sparse-dmsolve.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
71 } 71 }
72 else 72 else
73 { 73 {
74 OCTAVE_LOCAL_BUFFER (T, X, rend - rst); 74 OCTAVE_LOCAL_BUFFER (T, X, rend - rst);
75 octave_sort<octave_idx_type> sort; 75 octave_sort<octave_idx_type> sort;
76 octave_idx_type *ri = B.xridx(); 76 octave_idx_type *ri = B.xridx ();
77 nz = 0; 77 nz = 0;
78 for (octave_idx_type j = cst ; j < cend ; j++) 78 for (octave_idx_type j = cst ; j < cend ; j++)
79 { 79 {
80 octave_idx_type qq = (Q ? Q [j] : j); 80 octave_idx_type qq = (Q ? Q [j] : j);
81 B.xcidx (j - cst) = nz; 81 B.xcidx (j - cst) = nz;
156 template <class T> 156 template <class T>
157 static void 157 static void
158 dmsolve_insert (MArray<T> &a, const MArray<T> &b, const octave_idx_type *Q, 158 dmsolve_insert (MArray<T> &a, const MArray<T> &b, const octave_idx_type *Q,
159 octave_idx_type r, octave_idx_type c) 159 octave_idx_type r, octave_idx_type c)
160 { 160 {
161 T *ax = a.fortran_vec(); 161 T *ax = a.fortran_vec ();
162 const T *bx = b.fortran_vec(); 162 const T *bx = b.fortran_vec ();
163 octave_idx_type anr = a.rows(); 163 octave_idx_type anr = a.rows ();
164 octave_idx_type nr = b.rows(); 164 octave_idx_type nr = b.rows ();
165 octave_idx_type nc = b.cols(); 165 octave_idx_type nc = b.cols ();
166 for (octave_idx_type j = 0; j < nc; j++) 166 for (octave_idx_type j = 0; j < nc; j++)
167 { 167 {
168 octave_idx_type aoff = (c + j) * anr; 168 octave_idx_type aoff = (c + j) * anr;
169 octave_idx_type boff = j * nr; 169 octave_idx_type boff = j * nr;
170 for (octave_idx_type i = 0; i < nr; i++) 170 for (octave_idx_type i = 0; i < nr; i++)
212 212
213 OCTAVE_LOCAL_BUFFER (T, X, nr); 213 OCTAVE_LOCAL_BUFFER (T, X, nr);
214 octave_sort<octave_idx_type> sort; 214 octave_sort<octave_idx_type> sort;
215 MSparse<T> tmp (a); 215 MSparse<T> tmp (a);
216 a = MSparse<T> (nr, nc, nel); 216 a = MSparse<T> (nr, nc, nel);
217 octave_idx_type *ri = a.xridx(); 217 octave_idx_type *ri = a.xridx ();
218 218
219 for (octave_idx_type i = 0; i < tmp.cidx(c); i++) 219 for (octave_idx_type i = 0; i < tmp.cidx(c); i++)
220 { 220 {
221 a.xdata(i) = tmp.xdata(i); 221 a.xdata(i) = tmp.xdata(i);
222 a.xridx(i) = tmp.xridx(i); 222 a.xridx(i) = tmp.xridx(i);
276 static void 276 static void
277 dmsolve_permute (MArray<RT> &a, const MArray<T>& b, const octave_idx_type *p) 277 dmsolve_permute (MArray<RT> &a, const MArray<T>& b, const octave_idx_type *p)
278 { 278 {
279 octave_idx_type b_nr = b.rows (); 279 octave_idx_type b_nr = b.rows ();
280 octave_idx_type b_nc = b.cols (); 280 octave_idx_type b_nc = b.cols ();
281 const T *Bx = b.fortran_vec(); 281 const T *Bx = b.fortran_vec ();
282 a.resize (dim_vector (b_nr, b_nc)); 282 a.resize (dim_vector (b_nr, b_nc));
283 RT *Btx = a.fortran_vec(); 283 RT *Btx = a.fortran_vec ();
284 for (octave_idx_type j = 0; j < b_nc; j++) 284 for (octave_idx_type j = 0; j < b_nc; j++)
285 { 285 {
286 octave_idx_type off = j * b_nr; 286 octave_idx_type off = j * b_nr;
287 for (octave_idx_type i = 0; i < b_nr; i++) 287 for (octave_idx_type i = 0; i < b_nr; i++)
288 { 288 {
314 octave_idx_type b_nc = b.cols (); 314 octave_idx_type b_nc = b.cols ();
315 octave_idx_type b_nz = b.nnz (); 315 octave_idx_type b_nz = b.nnz ();
316 octave_idx_type nz = 0; 316 octave_idx_type nz = 0;
317 a = MSparse<RT> (b_nr, b_nc, b_nz); 317 a = MSparse<RT> (b_nr, b_nc, b_nz);
318 octave_sort<octave_idx_type> sort; 318 octave_sort<octave_idx_type> sort;
319 octave_idx_type *ri = a.xridx(); 319 octave_idx_type *ri = a.xridx ();
320 OCTAVE_LOCAL_BUFFER (RT, X, b_nr); 320 OCTAVE_LOCAL_BUFFER (RT, X, b_nr);
321 a.xcidx(0) = 0; 321 a.xcidx(0) = 0;
322 for (octave_idx_type j = 0; j < b_nc; j++) 322 for (octave_idx_type j = 0; j < b_nc; j++)
323 { 323 {
324 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 324 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
409 // Leading over-determined block 409 // Leading over-determined block
410 if (dm->rr [2] < nr && dm->cc [3] < nc) 410 if (dm->rr [2] < nr && dm->cc [3] < nc)
411 { 411 {
412 ST m = dmsolve_extract (a, pinv, q, dm->rr [2], nr, dm->cc [3], nc, 412 ST m = dmsolve_extract (a, pinv, q, dm->rr [2], nr, dm->cc [3], nc,
413 nnz_remaining, true); 413 nnz_remaining, true);
414 nnz_remaining -= m.nnz(); 414 nnz_remaining -= m.nnz ();
415 RT mtmp = 415 RT mtmp =
416 qrsolve (m, dmsolve_extract (btmp, 0, 0, dm->rr[2], b_nr, 0, 416 qrsolve (m, dmsolve_extract (btmp, 0, 0, dm->rr[2], b_nr, 0,
417 b_nc), info); 417 b_nc), info);
418 dmsolve_insert (retval, mtmp, q, dm->cc [3], 0); 418 dmsolve_insert (retval, mtmp, q, dm->cc [3], 0);
419 if (dm->rr [2] > 0 && !info) 419 if (dm->rr [2] > 0 && !info)
420 { 420 {
421 m = dmsolve_extract (a, pinv, q, 0, dm->rr [2], 421 m = dmsolve_extract (a, pinv, q, 0, dm->rr [2],
422 dm->cc [3], nc, nnz_remaining, true); 422 dm->cc [3], nc, nnz_remaining, true);
423 nnz_remaining -= m.nnz(); 423 nnz_remaining -= m.nnz ();
424 RT ctmp = dmsolve_extract (btmp, 0, 0, 0, 424 RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
425 dm->rr[2], 0, b_nc); 425 dm->rr[2], 0, b_nc);
426 btmp.insert (ctmp - m * mtmp, 0, 0); 426 btmp.insert (ctmp - m * mtmp, 0, 0);
427 } 427 }
428 } 428 }
431 // FIXME Should use fine Dulmange-Mendelsohn decomposition here. 431 // FIXME Should use fine Dulmange-Mendelsohn decomposition here.
432 if (dm->rr [1] < dm->rr [2] && dm->cc [2] < dm->cc [3] && !info) 432 if (dm->rr [1] < dm->rr [2] && dm->cc [2] < dm->cc [3] && !info)
433 { 433 {
434 ST m = dmsolve_extract (a, pinv, q, dm->rr [1], dm->rr [2], 434 ST m = dmsolve_extract (a, pinv, q, dm->rr [1], dm->rr [2],
435 dm->cc [2], dm->cc [3], nnz_remaining, false); 435 dm->cc [2], dm->cc [3], nnz_remaining, false);
436 nnz_remaining -= m.nnz(); 436 nnz_remaining -= m.nnz ();
437 RT btmp2 = dmsolve_extract (btmp, 0, 0, dm->rr [1], dm->rr [2], 437 RT btmp2 = dmsolve_extract (btmp, 0, 0, dm->rr [1], dm->rr [2],
438 0, b_nc); 438 0, b_nc);
439 double rcond = 0.0; 439 double rcond = 0.0;
440 MatrixType mtyp (MatrixType::Full); 440 MatrixType mtyp (MatrixType::Full);
441 RT mtmp = m.solve (mtyp, btmp2, info, rcond, 441 RT mtmp = m.solve (mtyp, btmp2, info, rcond,
449 dmsolve_insert (retval, mtmp, q, dm->cc [2], 0); 449 dmsolve_insert (retval, mtmp, q, dm->cc [2], 0);
450 if (dm->rr [1] > 0 && !info) 450 if (dm->rr [1] > 0 && !info)
451 { 451 {
452 m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], dm->cc [2], 452 m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], dm->cc [2],
453 dm->cc [3], nnz_remaining, true); 453 dm->cc [3], nnz_remaining, true);
454 nnz_remaining -= m.nnz(); 454 nnz_remaining -= m.nnz ();
455 RT ctmp = dmsolve_extract (btmp, 0, 0, 0, 455 RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
456 dm->rr[1], 0, b_nc); 456 dm->rr[1], 0, b_nc);
457 btmp.insert (ctmp - m * mtmp, 0, 0); 457 btmp.insert (ctmp - m * mtmp, 0, 0);
458 } 458 }
459 } 459 }