Mercurial > jwe > octave

comparison liboctave/numeric/oct-norm.cc @ 21139:538b57866b90

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
consistently use "typename" intead of "class" in template declarations * Object.h, QtHandlesUtils.cc, QtHandlesUtils.h, ToolBarButton.cc, ToolBarButton.h, Cell.h, __lin_interpn__.cc, bitfcns.cc, bsxfun.cc, cellfun.cc, data.cc, filter.cc, gcd.cc, graphics.cc, help.cc, kron.cc, lookup.cc, ls-mat5.cc, ls-oct-text.h, lu.cc, max.cc, mgorth.cc, oct-map.cc, oct-map.h, oct-stream.cc, oct-stream.h, octave-link.h, pr-output.cc, profiler.h, schur.cc, sparse-xdiv.cc, sparse-xpow.cc, sqrtm.cc, symtab.h, tril.cc, typecast.cc, variables.cc, xdiv.cc, zfstream.h, __init_fltk__.cc, __magick_read__.cc, chol.cc, qr.cc, ov-base-diag.cc, ov-base-diag.h, ov-base-int.cc, ov-base-int.h, ov-base-mat.cc, ov-base-mat.h, ov-base-scalar.cc, ov-base-scalar.h, ov-base-sparse.cc, ov-base-sparse.h, ov-base.h, ov-classdef.cc, ov-int-traits.h, ov-java.h, ov-usr-fcn.h, ov.cc, ov.h, op-dms-template.cc, oct-parse.in.yy, parse.h, pt-mat.cc, Array-b.cc, Array.cc, Array.h, CDiagMatrix.h, CMatrix.h, CNDArray.h, DiagArray2.cc, DiagArray2.h, MArray.cc, MArray.h, MDiagArray2.cc, MDiagArray2.h, MSparse.cc, MSparse.h, MatrixType.cc, Sparse.cc, Sparse.h, dDiagMatrix.h, dMatrix.h, dNDArray.h, fCDiagMatrix.h, fCMatrix.h, fCNDArray.h, fDiagMatrix.h, fMatrix.h, fNDArray.h, idx-vector.cc, idx-vector.h, intNDArray.cc, intNDArray.h, DET.h, base-aepbal.h, base-lu.cc, base-lu.h, base-qr.cc, base-qr.h, bsxfun-defs.cc, eigs-base.cc, lo-mappers.h, lo-specfun.cc, lo-specfun.h, oct-convn.cc, oct-fftw.cc, oct-norm.cc, sparse-base-chol.cc, sparse-base-chol.h, sparse-base-lu.cc, sparse-base-lu.h, sparse-dmsolve.cc, mx-inlines.cc, action-container.h, base-list.h, lo-traits.h, lo-utils.h, oct-base64.h, oct-binmap.h, oct-cmplx.h, oct-inttypes.cc, oct-inttypes.h, oct-locbuf.h, oct-refcount.h, oct-sort.cc, oct-sort.h: Use "typename" instead of "class" in template declarations.
author John W. Eaton <jwe@octave.org>
date Sun, 24 Jan 2016 13:50:04 -0500
parents 7b6d8c19dab0
children f7121e111991
comparison
equal deleted inserted replaced
21138:e2fca7d79169 21139:538b57866b90
65 // to handle both real and complex element, and a cast operator 65 // to handle both real and complex element, and a cast operator
66 // returning the intermediate norm. Reference: Higham, N. "Estimating 66 // returning the intermediate norm. Reference: Higham, N. "Estimating
67 // the Matrix p-Norm." Numer. Math. 62, 539-555, 1992. 67 // the Matrix p-Norm." Numer. Math. 62, 539-555, 1992.
68 68
69 // norm accumulator for the p-norm 69 // norm accumulator for the p-norm
70 template <class R> 70 template <typename R>
71 class norm_accumulator_p 71 class norm_accumulator_p
72 { 72 {
73 R p,scl,sum; 73 R p,scl,sum;
74 public: 74 public:
75 norm_accumulator_p () {} // we need this one for Array 75 norm_accumulator_p () {} // we need this one for Array
76 norm_accumulator_p (R pp) : p(pp), scl(0), sum(1) {} 76 norm_accumulator_p (R pp) : p(pp), scl(0), sum(1) {}
77 77
78 template<class U> 78 template <typename U>
79 void accum (U val) 79 void accum (U val)
80 { 80 {
81 octave_quit (); 81 octave_quit ();
82 R t = std::abs (val); 82 R t = std::abs (val);
83 if (scl == t) // we need this to handle Infs properly 83 if (scl == t) // we need this to handle Infs properly
93 } 93 }
94 operator R () { return scl * std::pow (sum, 1/p); } 94 operator R () { return scl * std::pow (sum, 1/p); }
95 }; 95 };
96 96
97 // norm accumulator for the minus p-pseudonorm 97 // norm accumulator for the minus p-pseudonorm
98 template <class R> 98 template <typename R>
99 class norm_accumulator_mp 99 class norm_accumulator_mp
100 { 100 {
101 R p,scl,sum; 101 R p,scl,sum;
102 public: 102 public:
103 norm_accumulator_mp () {} // we need this one for Array 103 norm_accumulator_mp () {} // we need this one for Array
104 norm_accumulator_mp (R pp) : p(pp), scl(0), sum(1) {} 104 norm_accumulator_mp (R pp) : p(pp), scl(0), sum(1) {}
105 105
106 template<class U> 106 template <typename U>
107 void accum (U val) 107 void accum (U val)
108 { 108 {
109 octave_quit (); 109 octave_quit ();
110 R t = 1 / std::abs (val); 110 R t = 1 / std::abs (val);
111 if (scl == t) 111 if (scl == t)
121 } 121 }
122 operator R () { return scl * std::pow (sum, -1/p); } 122 operator R () { return scl * std::pow (sum, -1/p); }
123 }; 123 };
124 124
125 // norm accumulator for the 2-norm (euclidean) 125 // norm accumulator for the 2-norm (euclidean)
126 template <class R> 126 template <typename R>
127 class norm_accumulator_2 127 class norm_accumulator_2
128 { 128 {
129 R scl,sum; 129 R scl,sum;
130 static R pow2 (R x) { return x*x; } 130 static R pow2 (R x) { return x*x; }
131 public: 131 public:
154 154
155 operator R () { return scl * std::sqrt (sum); } 155 operator R () { return scl * std::sqrt (sum); }
156 }; 156 };
157 157
158 // norm accumulator for the 1-norm (city metric) 158 // norm accumulator for the 1-norm (city metric)
159 template <class R> 159 template <typename R>
160 class norm_accumulator_1 160 class norm_accumulator_1
161 { 161 {
162 R sum; 162 R sum;
163 public: 163 public:
164 norm_accumulator_1 () : sum (0) {} 164 norm_accumulator_1 () : sum (0) {}
165 template<class U> 165 template <typename U>
166 void accum (U val) 166 void accum (U val)
167 { 167 {
168 sum += std::abs (val); 168 sum += std::abs (val);
169 } 169 }
170 operator R () { return sum; } 170 operator R () { return sum; }
171 }; 171 };
172 172
173 // norm accumulator for the inf-norm (max metric) 173 // norm accumulator for the inf-norm (max metric)
174 template <class R> 174 template <typename R>
175 class norm_accumulator_inf 175 class norm_accumulator_inf
176 { 176 {
177 R max; 177 R max;
178 public: 178 public:
179 norm_accumulator_inf () : max (0) {} 179 norm_accumulator_inf () : max (0) {}
180 template<class U> 180 template <typename U>
181 void accum (U val) 181 void accum (U val)
182 { 182 {
183 max = std::max (max, std::abs (val)); 183 max = std::max (max, std::abs (val));
184 } 184 }
185 operator R () { return max; } 185 operator R () { return max; }
186 }; 186 };
187 187
188 // norm accumulator for the -inf pseudonorm (min abs value) 188 // norm accumulator for the -inf pseudonorm (min abs value)
189 template <class R> 189 template <typename R>
190 class norm_accumulator_minf 190 class norm_accumulator_minf
191 { 191 {
192 R min; 192 R min;
193 public: 193 public:
194 norm_accumulator_minf () : min (octave_Inf) {} 194 norm_accumulator_minf () : min (octave_Inf) {}
195 template<class U> 195 template <typename U>
196 void accum (U val) 196 void accum (U val)
197 { 197 {
198 min = std::min (min, std::abs (val)); 198 min = std::min (min, std::abs (val));
199 } 199 }
200 operator R () { return min; } 200 operator R () { return min; }
201 }; 201 };
202 202
203 // norm accumulator for the 0-pseudonorm (hamming distance) 203 // norm accumulator for the 0-pseudonorm (hamming distance)
204 template <class R> 204 template <typename R>
205 class norm_accumulator_0 205 class norm_accumulator_0
206 { 206 {
207 unsigned int num; 207 unsigned int num;
208 public: 208 public:
209 norm_accumulator_0 () : num (0) {} 209 norm_accumulator_0 () : num (0) {}
210 template<class U> 210 template <typename U>
211 void accum (U val) 211 void accum (U val)
212 { 212 {
213 if (val != static_cast<U> (0)) ++num; 213 if (val != static_cast<U> (0)) ++num;
214 } 214 }
215 operator R () { return num; } 215 operator R () { return num; }
216 }; 216 };
217 217
218 218
219 // OK, we're armed :) Now let's go for the fun 219 // OK, we're armed :) Now let's go for the fun
220 220
221 template <class T, class R, class ACC> 221 template <typename T, typename R, typename ACC>
222 inline void vector_norm (const Array<T>& v, R& res, ACC acc) 222 inline void vector_norm (const Array<T>& v, R& res, ACC acc)
223 { 223 {
224 for (octave_idx_type i = 0; i < v.numel (); i++) 224 for (octave_idx_type i = 0; i < v.numel (); i++)
225 acc.accum (v(i)); 225 acc.accum (v(i));
226 226
227 res = acc; 227 res = acc;
228 } 228 }
229 229
230 // dense versions 230 // dense versions
231 template <class T, class R, class ACC> 231 template <typename T, typename R, typename ACC>
232 void column_norms (const MArray<T>& m, MArray<R>& res, ACC acc) 232 void column_norms (const MArray<T>& m, MArray<R>& res, ACC acc)
233 { 233 {
234 res = MArray<R> (dim_vector (1, m.columns ())); 234 res = MArray<R> (dim_vector (1, m.columns ()));
235 for (octave_idx_type j = 0; j < m.columns (); j++) 235 for (octave_idx_type j = 0; j < m.columns (); j++)
236 { 236 {
240 240
241 res.xelem (j) = accj; 241 res.xelem (j) = accj;
242 } 242 }
243 } 243 }
244 244
245 template <class T, class R, class ACC> 245 template <typename T, typename R, typename ACC>
246 void row_norms (const MArray<T>& m, MArray<R>& res, ACC acc) 246 void row_norms (const MArray<T>& m, MArray<R>& res, ACC acc)
247 { 247 {
248 res = MArray<R> (dim_vector (m.rows (), 1)); 248 res = MArray<R> (dim_vector (m.rows (), 1));
249 std::vector<ACC> acci (m.rows (), acc); 249 std::vector<ACC> acci (m.rows (), acc);
250 for (octave_idx_type j = 0; j < m.columns (); j++) 250 for (octave_idx_type j = 0; j < m.columns (); j++)
256 for (octave_idx_type i = 0; i < m.rows (); i++) 256 for (octave_idx_type i = 0; i < m.rows (); i++)
257 res.xelem (i) = acci[i]; 257 res.xelem (i) = acci[i];
258 } 258 }
259 259
260 // sparse versions 260 // sparse versions
261 template <class T, class R, class ACC> 261 template <typename T, typename R, typename ACC>
262 void column_norms (const MSparse<T>& m, MArray<R>& res, ACC acc) 262 void column_norms (const MSparse<T>& m, MArray<R>& res, ACC acc)
263 { 263 {
264 res = MArray<R> (dim_vector (1, m.columns ())); 264 res = MArray<R> (dim_vector (1, m.columns ()));
265 for (octave_idx_type j = 0; j < m.columns (); j++) 265 for (octave_idx_type j = 0; j < m.columns (); j++)
266 { 266 {
270 270
271 res.xelem (j) = accj; 271 res.xelem (j) = accj;
272 } 272 }
273 } 273 }
274 274
275 template <class T, class R, class ACC> 275 template <typename T, typename R, typename ACC>
276 void row_norms (const MSparse<T>& m, MArray<R>& res, ACC acc) 276 void row_norms (const MSparse<T>& m, MArray<R>& res, ACC acc)
277 { 277 {
278 res = MArray<R> (dim_vector (m.rows (), 1)); 278 res = MArray<R> (dim_vector (m.rows (), 1));
279 std::vector<ACC> acci (m.rows (), acc); 279 std::vector<ACC> acci (m.rows (), acc);
280 for (octave_idx_type j = 0; j < m.columns (); j++) 280 for (octave_idx_type j = 0; j < m.columns (); j++)
287 res.xelem (i) = acci[i]; 287 res.xelem (i) = acci[i];
288 } 288 }
289 289
290 // now the dispatchers 290 // now the dispatchers
291 #define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE) \ 291 #define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE) \
292 template <class T, class R> \ 292 template <typename T, typename R> \
293 RES_TYPE FUNC_NAME (const ARG_TYPE& v, R p) \ 293 RES_TYPE FUNC_NAME (const ARG_TYPE& v, R p) \
294 { \ 294 { \
295 RES_TYPE res; \ 295 RES_TYPE res; \
296 if (p == 2) \ 296 if (p == 2) \
297 FUNC_NAME (v, res, norm_accumulator_2<R> ()); \ 297 FUNC_NAME (v, res, norm_accumulator_2<R> ()); \
320 DEFINE_DISPATCHER (row_norms, MSparse<T>, MArray<R>) 320 DEFINE_DISPATCHER (row_norms, MSparse<T>, MArray<R>)
321 321
322 // The approximate subproblem in Higham's method. Find lambda and mu such that 322 // The approximate subproblem in Higham's method. Find lambda and mu such that
323 // norm ([lambda, mu], p) == 1 and norm (y*lambda + col*mu, p) is maximized. 323 // norm ([lambda, mu], p) == 1 and norm (y*lambda + col*mu, p) is maximized.
324 // Real version. As in Higham's paper. 324 // Real version. As in Higham's paper.
325 template <class ColVectorT, class R> 325 template <typename ColVectorT, typename R>
326 static void 326 static void
327 higham_subp (const ColVectorT& y, const ColVectorT& col, 327 higham_subp (const ColVectorT& y, const ColVectorT& col,
328 octave_idx_type nsamp, R p, R& lambda, R& mu) 328 octave_idx_type nsamp, R p, R& lambda, R& mu)
329 { 329 {
330 R nrm = 0; 330 R nrm = 0;
348 } 348 }
349 349
350 // Complex version. Higham's paper does not deal with complex case, so we use a 350 // Complex version. Higham's paper does not deal with complex case, so we use a
351 // simple extension. First, guess the magnitudes as in real version, then try 351 // simple extension. First, guess the magnitudes as in real version, then try
352 // to rotate lambda to improve further. 352 // to rotate lambda to improve further.
353 template <class ColVectorT, class R> 353 template <typename ColVectorT, typename R>
354 static void 354 static void
355 higham_subp (const ColVectorT& y, const ColVectorT& col, 355 higham_subp (const ColVectorT& y, const ColVectorT& col,
356 octave_idx_type nsamp, R p, 356 octave_idx_type nsamp, R p,
357 std::complex<R>& lambda, std::complex<R>& mu) 357 std::complex<R>& lambda, std::complex<R>& mu)
358 { 358 {
393 } 393 }
394 } 394 }
395 } 395 }
396 396
397 // the p-dual element (should work for both real and complex) 397 // the p-dual element (should work for both real and complex)
398 template <class T, class R> 398 template <typename T, typename R>
399 inline T elem_dual_p (T x, R p) 399 inline T elem_dual_p (T x, R p)
400 { 400 {
401 return signum (x) * std::pow (std::abs (x), p-1); 401 return signum (x) * std::pow (std::abs (x), p-1);
402 } 402 }
403 403
404 // the VectorT is used for vectors, but actually it has to be 404 // the VectorT is used for vectors, but actually it has to be
405 // a Matrix type to allow all the operations. For instance SparseMatrix 405 // a Matrix type to allow all the operations. For instance SparseMatrix
406 // does not support multiplication with column/row vectors. 406 // does not support multiplication with column/row vectors.
407 // the dual vector 407 // the dual vector
408 template <class VectorT, class R> 408 template <typename VectorT, typename R>
409 VectorT dual_p (const VectorT& x, R p, R q) 409 VectorT dual_p (const VectorT& x, R p, R q)
410 { 410 {
411 VectorT res (x.dims ()); 411 VectorT res (x.dims ());
412 for (octave_idx_type i = 0; i < x.numel (); i++) 412 for (octave_idx_type i = 0; i < x.numel (); i++)
413 res.xelem (i) = elem_dual_p (x(i), p); 413 res.xelem (i) = elem_dual_p (x(i), p);
414 return res / vector_norm (res, q); 414 return res / vector_norm (res, q);
415 } 415 }
416 416
417 // Higham's hybrid method 417 // Higham's hybrid method
418 template <class MatrixT, class VectorT, class R> 418 template <typename MatrixT, typename VectorT, typename R>
419 R higham (const MatrixT& m, R p, R tol, int maxiter, 419 R higham (const MatrixT& m, R p, R tol, int maxiter,
420 VectorT& x) 420 VectorT& x)
421 { 421 {
422 x.resize (m.columns (), 1); 422 x.resize (m.columns (), 1);
423 // the OSE part 423 // the OSE part
473 // be a good value. Eventually, we can provide a means to change this 473 // be a good value. Eventually, we can provide a means to change this
474 // constant from Octave. 474 // constant from Octave.
475 static int max_norm_iter = 100; 475 static int max_norm_iter = 100;
476 476
477 // version with SVD for dense matrices 477 // version with SVD for dense matrices
478 template <class MatrixT, class VectorT, class SVDT, class R> 478 template <typename MatrixT, typename VectorT, typename SVDT, typename R>
479 R matrix_norm (const MatrixT& m, R p, VectorT, SVDT) 479 R matrix_norm (const MatrixT& m, R p, VectorT, SVDT)
480 { 480 {
481 R res = 0; 481 R res = 0;
482 if (p == 2) 482 if (p == 2)
483 { 483 {
499 499
500 return res; 500 return res;
501 } 501 }
502 502
503 // SVD-free version for sparse matrices 503 // SVD-free version for sparse matrices
504 template <class MatrixT, class VectorT, class R> 504 template <typename MatrixT, typename VectorT, typename R>
505 R matrix_norm (const MatrixT& m, R p, VectorT) 505 R matrix_norm (const MatrixT& m, R p, VectorT)
506 { 506 {
507 R res = 0; 507 R res = 0;
508 if (p == 1) 508 if (p == 1)
509 res = xcolnorms (m, 1).max (); 509 res = xcolnorms (m, 1).max ();
537 DEFINE_XNORM_FUNCS(Complex, double) 537 DEFINE_XNORM_FUNCS(Complex, double)
538 DEFINE_XNORM_FUNCS(Float, float) 538 DEFINE_XNORM_FUNCS(Float, float)
539 DEFINE_XNORM_FUNCS(FloatComplex, float) 539 DEFINE_XNORM_FUNCS(FloatComplex, float)
540 540
541 // this is needed to avoid copying the sparse matrix for xfrobnorm 541 // this is needed to avoid copying the sparse matrix for xfrobnorm
542 template <class T, class R> 542 template <typename T, typename R>
543 inline void array_norm_2 (const T* v, octave_idx_type n, R& res) 543 inline void array_norm_2 (const T* v, octave_idx_type n, R& res)
544 { 544 {
545 norm_accumulator_2<R> acc; 545 norm_accumulator_2<R> acc;
546 for (octave_idx_type i = 0; i < n; i++) 546 for (octave_idx_type i = 0; i < n; i++)
547 acc.accum (v[i]); 547 acc.accum (v[i]);