Mercurial > hg > octave-nkf

comparison liboctave/sparse-dmsolve.cc @ 14888:4315a39da4c9 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Do computations with octave_uint64 to avoid overflow * sparse-dmsolve (MSparse<T>::dmsolve_extract): Perform multiplication and comparison in octave_uint64 to avoid overflow.
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Thu, 19 Jul 2012 18:08:08 -0400
parents 72c96de7a403
children 1316bfc6e260
comparison
equal deleted inserted replaced
14880:9e5df4d6cefa 14888:4315a39da4c9
31 #include "SparseQR.h" 31 #include "SparseQR.h"
32 #include "SparseCmplxQR.h" 32 #include "SparseCmplxQR.h"
33 #include "MatrixType.h" 33 #include "MatrixType.h"
34 #include "oct-sort.h" 34 #include "oct-sort.h"
35 #include "oct-locbuf.h" 35 #include "oct-locbuf.h"
36 #include "oct-inttypes.h"
36 37
37 template <class T> 38 template <class T>
38 static MSparse<T> 39 static MSparse<T>
39 dmsolve_extract (const MSparse<T> &A, const octave_idx_type *Pinv, 40 dmsolve_extract (const MSparse<T> &A, const octave_idx_type *Pinv,
40 const octave_idx_type *Q, octave_idx_type rst, 41 const octave_idx_type *Q, octave_idx_type rst,
41 octave_idx_type rend, octave_idx_type cst, 42 octave_idx_type rend, octave_idx_type cst,
42 octave_idx_type cend, octave_idx_type maxnz = -1, 43 octave_idx_type cend, octave_idx_type maxnz = -1,
43 bool lazy = false) 44 bool lazy = false)
44 { 45 {
45 octave_idx_type nz = (rend - rst) * (cend - cst); 46 octave_idx_type nr = rend - rst, nc = cend - cst;
46 maxnz = (maxnz < 0 ? A.nnz () : maxnz); 47 maxnz = (maxnz < 0 ? A.nnz () : maxnz);
47 MSparse<T> B (rend - rst, cend - cst, (nz < maxnz ? nz : maxnz)); 48 octave_idx_type nz;
49
50 // Cast to uint64 to handle overflow in this multiplication
51 if (octave_uint64 (nr)*octave_uint64 (nc) < octave_uint64 (maxnz))
52 nz = nr*nz;
53 else
54 nz = maxnz;
55
56 MSparse<T> B (nr, nc, (nz < maxnz ? nz : maxnz));
48 // Some sparse functions can support lazy indexing (where elements 57 // Some sparse functions can support lazy indexing (where elements
49 // in the row are in no particular order), even though octave in 58 // in the row are in no particular order), even though octave in
50 // general can't. For those functions that can using it is a big 59 // general can't. For those functions that can using it is a big
51 // win here in terms of speed. 60 // win here in terms of speed.
52 if (lazy) 61 if (lazy)