octave-nkf: liboctave/numeric/oct-norm.cc annotate

annotate liboctave/numeric/oct-norm.cc @ 17769:49a5a4be04a1

maint: Use GNU style coding conventions for code in liboctave/ * liboctave/array/Array-C.cc, liboctave/array/Array-b.cc, liboctave/array/Array-ch.cc, liboctave/array/Array-d.cc, liboctave/array/Array-f.cc, liboctave/array/Array-fC.cc, liboctave/array/Array-util.cc, liboctave/array/Array-util.h, liboctave/array/Array.cc, liboctave/array/Array.h, liboctave/array/Array3.h, liboctave/array/CColVector.cc, liboctave/array/CColVector.h, liboctave/array/CDiagMatrix.cc, liboctave/array/CDiagMatrix.h, liboctave/array/CMatrix.cc, liboctave/array/CMatrix.h, liboctave/array/CNDArray.cc, liboctave/array/CNDArray.h, liboctave/array/CRowVector.cc, liboctave/array/CRowVector.h, liboctave/array/CSparse.cc, liboctave/array/CSparse.h, liboctave/array/DiagArray2.h, liboctave/array/MArray.cc, liboctave/array/MArray.h, liboctave/array/MDiagArray2.cc, liboctave/array/MDiagArray2.h, liboctave/array/MSparse.cc, liboctave/array/MSparse.h, liboctave/array/MatrixType.cc, liboctave/array/MatrixType.h, liboctave/array/PermMatrix.h, liboctave/array/Range.cc, liboctave/array/Range.h, liboctave/array/Sparse.cc, liboctave/array/Sparse.h, liboctave/array/boolMatrix.cc, liboctave/array/boolMatrix.h, liboctave/array/boolNDArray.cc, liboctave/array/boolNDArray.h, liboctave/array/boolSparse.cc, liboctave/array/boolSparse.h, liboctave/array/chMatrix.cc, liboctave/array/chMatrix.h, liboctave/array/chNDArray.cc, liboctave/array/chNDArray.h, liboctave/array/dColVector.h, liboctave/array/dDiagMatrix.cc, liboctave/array/dDiagMatrix.h, liboctave/array/dMatrix.cc, liboctave/array/dMatrix.h, liboctave/array/dNDArray.cc, liboctave/array/dNDArray.h, liboctave/array/dRowVector.h, liboctave/array/dSparse.cc, liboctave/array/dSparse.h, liboctave/array/dim-vector.cc, liboctave/array/dim-vector.h, liboctave/array/fCColVector.cc, liboctave/array/fCColVector.h, liboctave/array/fCDiagMatrix.cc, liboctave/array/fCDiagMatrix.h, liboctave/array/fCMatrix.cc, liboctave/array/fCMatrix.h, liboctave/array/fCNDArray.cc, liboctave/array/fCNDArray.h, liboctave/array/fCRowVector.cc, liboctave/array/fCRowVector.h, liboctave/array/fColVector.h, liboctave/array/fDiagMatrix.cc, liboctave/array/fDiagMatrix.h, liboctave/array/fMatrix.cc, liboctave/array/fMatrix.h, liboctave/array/fNDArray.cc, liboctave/array/fNDArray.h, liboctave/array/fRowVector.h, liboctave/array/idx-vector.cc, liboctave/array/idx-vector.h, liboctave/array/intNDArray.cc, liboctave/array/intNDArray.h, liboctave/cruft/misc/blaswrap.c, liboctave/cruft/misc/quit.cc, liboctave/numeric/CmplxCHOL.cc, liboctave/numeric/CmplxCHOL.h, liboctave/numeric/CmplxGEPBAL.cc, liboctave/numeric/CmplxGEPBAL.h, liboctave/numeric/CmplxHESS.h, liboctave/numeric/CmplxLU.cc, liboctave/numeric/CmplxLU.h, liboctave/numeric/CmplxQR.cc, liboctave/numeric/CmplxQRP.cc, liboctave/numeric/CmplxQRP.h, liboctave/numeric/CmplxSCHUR.h, liboctave/numeric/CmplxSVD.cc, liboctave/numeric/CmplxSVD.h, liboctave/numeric/CollocWt.h, liboctave/numeric/DAE.h, liboctave/numeric/DAEFunc.h, liboctave/numeric/DAERT.h, liboctave/numeric/DAERTFunc.h, liboctave/numeric/DASPK.cc, liboctave/numeric/DASRT.cc, liboctave/numeric/DASRT.h, liboctave/numeric/DASSL.cc, liboctave/numeric/DET.h, liboctave/numeric/EIG.cc, liboctave/numeric/EIG.h, liboctave/numeric/LSODE.cc, liboctave/numeric/ODE.h, liboctave/numeric/ODEFunc.h, liboctave/numeric/ODES.h, liboctave/numeric/ODESFunc.h, liboctave/numeric/Quad.cc, liboctave/numeric/Quad.h, liboctave/numeric/SparseCmplxCHOL.h, liboctave/numeric/SparseCmplxLU.cc, liboctave/numeric/SparseCmplxLU.h, liboctave/numeric/SparseCmplxQR.cc, liboctave/numeric/SparseCmplxQR.h, liboctave/numeric/SparseQR.cc, liboctave/numeric/SparseQR.h, liboctave/numeric/SparsedbleCHOL.h, liboctave/numeric/SparsedbleLU.cc, liboctave/numeric/SparsedbleLU.h, liboctave/numeric/base-aepbal.h, liboctave/numeric/base-dae.h, liboctave/numeric/base-de.h, liboctave/numeric/base-lu.cc, liboctave/numeric/base-lu.h, liboctave/numeric/base-min.h, liboctave/numeric/base-qr.h, liboctave/numeric/bsxfun.h, liboctave/numeric/dbleCHOL.cc, liboctave/numeric/dbleCHOL.h, liboctave/numeric/dbleGEPBAL.h, liboctave/numeric/dbleHESS.h, liboctave/numeric/dbleLU.cc, liboctave/numeric/dbleLU.h, liboctave/numeric/dbleQR.cc, liboctave/numeric/dbleQRP.cc, liboctave/numeric/dbleQRP.h, liboctave/numeric/dbleSCHUR.cc, liboctave/numeric/dbleSCHUR.h, liboctave/numeric/dbleSVD.cc, liboctave/numeric/dbleSVD.h, liboctave/numeric/eigs-base.cc, liboctave/numeric/fCmplxAEPBAL.cc, liboctave/numeric/fCmplxAEPBAL.h, liboctave/numeric/fCmplxCHOL.cc, liboctave/numeric/fCmplxCHOL.h, liboctave/numeric/fCmplxGEPBAL.cc, liboctave/numeric/fCmplxGEPBAL.h, liboctave/numeric/fCmplxHESS.h, liboctave/numeric/fCmplxLU.cc, liboctave/numeric/fCmplxLU.h, liboctave/numeric/fCmplxQR.cc, liboctave/numeric/fCmplxQR.h, liboctave/numeric/fCmplxQRP.cc, liboctave/numeric/fCmplxQRP.h, liboctave/numeric/fCmplxSCHUR.cc, liboctave/numeric/fCmplxSCHUR.h, liboctave/numeric/fCmplxSVD.h, liboctave/numeric/fEIG.cc, liboctave/numeric/fEIG.h, liboctave/numeric/floatCHOL.cc, liboctave/numeric/floatCHOL.h, liboctave/numeric/floatGEPBAL.cc, liboctave/numeric/floatGEPBAL.h, liboctave/numeric/floatHESS.h, liboctave/numeric/floatLU.cc, liboctave/numeric/floatLU.h, liboctave/numeric/floatQR.cc, liboctave/numeric/floatQRP.cc, liboctave/numeric/floatQRP.h, liboctave/numeric/floatSCHUR.cc, liboctave/numeric/floatSCHUR.h, liboctave/numeric/floatSVD.cc, liboctave/numeric/floatSVD.h, liboctave/numeric/lo-mappers.cc, liboctave/numeric/lo-mappers.h, liboctave/numeric/lo-specfun.cc, liboctave/numeric/lo-specfun.h, liboctave/numeric/oct-convn.cc, liboctave/numeric/oct-fftw.cc, liboctave/numeric/oct-fftw.h, liboctave/numeric/oct-norm.cc, liboctave/numeric/oct-rand.cc, liboctave/numeric/oct-rand.h, liboctave/numeric/randgamma.c, liboctave/numeric/randgamma.h, liboctave/numeric/randmtzig.c, liboctave/numeric/randpoisson.c, liboctave/numeric/randpoisson.h, liboctave/numeric/sparse-base-chol.h, liboctave/numeric/sparse-base-lu.h, liboctave/numeric/sparse-dmsolve.cc, liboctave/operators/Sparse-diag-op-defs.h, liboctave/operators/Sparse-op-defs.h, liboctave/operators/mx-inlines.cc, liboctave/system/dir-ops.h, liboctave/system/file-ops.cc, liboctave/system/file-stat.cc, liboctave/system/file-stat.h, liboctave/system/lo-sysdep.cc, liboctave/system/lo-sysdep.h, liboctave/system/mach-info.cc, liboctave/system/mach-info.h, liboctave/system/oct-env.cc, liboctave/system/oct-group.cc, liboctave/system/oct-syscalls.cc, liboctave/system/oct-syscalls.h, liboctave/system/oct-time.h, liboctave/system/tempname.c, liboctave/util/action-container.h, liboctave/util/base-list.h, liboctave/util/cmd-edit.cc, liboctave/util/cmd-edit.h, liboctave/util/cmd-hist.cc, liboctave/util/cmd-hist.h, liboctave/util/data-conv.cc, liboctave/util/data-conv.h, liboctave/util/kpse.cc, liboctave/util/lo-array-gripes.cc, liboctave/util/lo-cieee.c, liboctave/util/lo-regexp.cc, liboctave/util/lo-utils.cc, liboctave/util/oct-alloc.cc, liboctave/util/oct-base64.cc, liboctave/util/oct-binmap.h, liboctave/util/oct-cmplx.h, liboctave/util/oct-glob.cc, liboctave/util/oct-inttypes.cc, liboctave/util/oct-inttypes.h, liboctave/util/oct-locbuf.cc, liboctave/util/oct-locbuf.h, liboctave/util/oct-mem.h, liboctave/util/oct-mutex.cc, liboctave/util/oct-refcount.h, liboctave/util/oct-shlib.cc, liboctave/util/oct-shlib.h, liboctave/util/oct-sort.cc, liboctave/util/oct-sort.h, liboctave/util/pathsearch.cc, liboctave/util/pathsearch.h, liboctave/util/sparse-util.cc, liboctave/util/str-vec.cc, liboctave/util/str-vec.h, liboctave/util/unwind-prot.h, liboctave/util/url-transfer.cc, liboctave/util/url-transfer.h: Use GNU style coding conventions.

author	Rik <rik@octave.org>
date	Sat, 26 Oct 2013 18:57:05 -0700
parents	d63878346099
children	8e056300994b

rev	line source
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	1 /*
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	2
17744 d63878346099 maint: Update copyright notices for release. John W. Eaton <jwe@octave.org> parents: 15271 diff changeset	3 Copyright (C) 2008-2013 VZLU Prague, a.s.
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	4
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	5 This file is part of Octave.
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	6
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	7 Octave is free software; you can redistribute it and/or modify it
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	8 under the terms of the GNU General Public License as published by the
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	9 Free Software Foundation; either version 3 of the License, or (at your
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	10 option) any later version.
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	11
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	12 Octave is distributed in the hope that it will be useful, but WITHOUT
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	15 for more details.
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	16
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	17 You should have received a copy of the GNU General Public License
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	18 along with Octave; see the file COPYING. If not, see
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	19 <http://www.gnu.org/licenses/>.
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	20
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	21 */
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	22
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	23 // author: Jaroslav Hajek <highegg@gmail.com>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	24
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	25 #ifdef HAVE_CONFIG_H
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	26 #include <config.h>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	27 #endif
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	28
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	29 #include <cassert>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	30 #include <cfloat>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	31 #include <cmath>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	32
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	33 #include <iostream>
8308 5fe0f4dfdbec use std::vector as a simple linear container in oct-norm.cc to avoid problems with instantiating Array<T> Jaroslav Hajek <highegg@gmail.com> parents: 8303 diff changeset	34 #include <vector>
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	35
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	36 #include "oct-cmplx.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	37 #include "lo-error.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	38 #include "lo-ieee.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	39 #include "Array.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	40 #include "Array-util.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	41 #include "CMatrix.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	42 #include "dMatrix.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	43 #include "fCMatrix.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	44 #include "fMatrix.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	45 #include "CColVector.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	46 #include "dColVector.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	47 #include "CRowVector.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	48 #include "dRowVector.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	49 #include "fCColVector.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	50 #include "fColVector.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	51 #include "fCRowVector.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	52 #include "fRowVector.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	53 #include "CSparse.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	54 #include "dSparse.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	55 #include "dbleSVD.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	56 #include "CmplxSVD.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	57 #include "floatSVD.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	58 #include "fCmplxSVD.h"
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	59
13934 7a756a7e145b Add a citation to Higham's matrix norm paper Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 11586 diff changeset	60 // Theory: norm accumulator is an object that has an accum method able
7a756a7e145b Add a citation to Higham's matrix norm paper Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 11586 diff changeset	61 // to handle both real and complex element, and a cast operator
7a756a7e145b Add a citation to Higham's matrix norm paper Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 11586 diff changeset	62 // returning the intermediate norm. Reference: Higham, N. "Estimating
7a756a7e145b Add a citation to Higham's matrix norm paper Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 11586 diff changeset	63 // the Matrix p-Norm." Numer. Math. 62, 539-555, 1992.
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	64
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	65 // norm accumulator for the p-norm
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	66 template <class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	67 class norm_accumulator_p
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	68 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	69 R p,scl,sum;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	70 public:
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	71 norm_accumulator_p () {} // we need this one for Array
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	72 norm_accumulator_p (R pp) : p(pp), scl(0), sum(1) {}
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	73
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	74 template<class U>
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	75 void accum (U val)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	76 {
10145 fa01c1670b3e make p-norms breakable Jaroslav Hajek <highegg@gmail.com> parents: 10142 diff changeset	77 octave_quit ();
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	78 R t = std::abs (val);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	79 if (scl == t) // we need this to handle Infs properly
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	80 sum += 1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	81 else if (scl < t)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	82 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	83 sum *= std::pow (scl/t, p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	84 sum += 1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	85 scl = t;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	86 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	87 else if (t != 0)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	88 sum += std::pow (t/scl, p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	89 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	90 operator R () { return scl * std::pow (sum, 1/p); }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	91 };
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	92
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	93 // norm accumulator for the minus p-pseudonorm
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	94 template <class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	95 class norm_accumulator_mp
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	96 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	97 R p,scl,sum;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	98 public:
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	99 norm_accumulator_mp () {} // we need this one for Array
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	100 norm_accumulator_mp (R pp) : p(pp), scl(0), sum(1) {}
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	101
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	102 template<class U>
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	103 void accum (U val)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	104 {
10145 fa01c1670b3e make p-norms breakable Jaroslav Hajek <highegg@gmail.com> parents: 10142 diff changeset	105 octave_quit ();
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	106 R t = 1 / std::abs (val);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	107 if (scl == t)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	108 sum += 1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	109 else if (scl < t)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	110 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	111 sum *= std::pow (scl/t, p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	112 sum += 1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	113 scl = t;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	114 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	115 else if (t != 0)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	116 sum += std::pow (t/scl, p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	117 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	118 operator R () { return scl * std::pow (sum, -1/p); }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	119 };
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	120
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	121 // norm accumulator for the 2-norm (euclidean)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	122 template <class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	123 class norm_accumulator_2
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	124 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	125 R scl,sum;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	126 static R pow2 (R x) { return x*x; }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	127 public:
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	128 norm_accumulator_2 () : scl(0), sum(1) {}
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	129
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	130 void accum (R val)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	131 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	132 R t = std::abs (val);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	133 if (scl == t)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	134 sum += 1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	135 else if (scl < t)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	136 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	137 sum *= pow2 (scl/t);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	138 sum += 1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	139 scl = t;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	140 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	141 else if (t != 0)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	142 sum += pow2 (t/scl);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	143 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	144
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	145 void accum (std::complex<R> val)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	146 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	147 accum (val.real ());
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	148 accum (val.imag ());
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	149 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	150
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	151 operator R () { return scl * std::sqrt (sum); }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	152 };
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	153
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	154 // norm accumulator for the 1-norm (city metric)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	155 template <class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	156 class norm_accumulator_1
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	157 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	158 R sum;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	159 public:
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	160 norm_accumulator_1 () : sum (0) {}
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	161 template<class U>
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	162 void accum (U val)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	163 {
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	164 sum += std::abs (val);
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	165 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	166 operator R () { return sum; }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	167 };
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	168
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	169 // norm accumulator for the inf-norm (max metric)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	170 template <class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	171 class norm_accumulator_inf
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	172 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	173 R max;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	174 public:
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	175 norm_accumulator_inf () : max (0) {}
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	176 template<class U>
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	177 void accum (U val)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	178 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	179 max = std::max (max, std::abs (val));
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	180 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	181 operator R () { return max; }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	182 };
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	183
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	184 // norm accumulator for the -inf pseudonorm (min abs value)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	185 template <class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	186 class norm_accumulator_minf
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	187 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	188 R min;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	189 public:
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	190 norm_accumulator_minf () : min (octave_Inf) {}
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	191 template<class U>
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	192 void accum (U val)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	193 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	194 min = std::min (min, std::abs (val));
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	195 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	196 operator R () { return min; }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	197 };
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	198
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	199 // norm accumulator for the 0-pseudonorm (hamming distance)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	200 template <class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	201 class norm_accumulator_0
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	202 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	203 unsigned int num;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	204 public:
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	205 norm_accumulator_0 () : num (0) {}
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	206 template<class U>
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	207 void accum (U val)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	208 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	209 if (val != static_cast<U> (0)) ++num;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	210 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	211 operator R () { return num; }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	212 };
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	213
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	214
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	215 // OK, we're armed :) Now let's go for the fun
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	216
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	217 template <class T, class R, class ACC>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	218 inline void vector_norm (const Array<T>& v, R& res, ACC acc)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	219 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	220 for (octave_idx_type i = 0; i < v.numel (); i++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	221 acc.accum (v(i));
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	222
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	223 res = acc;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	224 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	225
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	226 // dense versions
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	227 template <class T, class R, class ACC>
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10145 diff changeset	228 void column_norms (const MArray<T>& m, MArray<R>& res, ACC acc)
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	229 {
11570 57632dea2446 attempt better backward compatibility for Array constructors John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	230 res = MArray<R> (dim_vector (1, m.columns ()));
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	231 for (octave_idx_type j = 0; j < m.columns (); j++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	232 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	233 ACC accj = acc;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	234 for (octave_idx_type i = 0; i < m.rows (); i++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	235 accj.accum (m(i, j));
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	236
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	237 res.xelem (j) = accj;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	238 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	239 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	240
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	241 template <class T, class R, class ACC>
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10145 diff changeset	242 void row_norms (const MArray<T>& m, MArray<R>& res, ACC acc)
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	243 {
11570 57632dea2446 attempt better backward compatibility for Array constructors John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	244 res = MArray<R> (dim_vector (m.rows (), 1));
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	245 std::vector<ACC> acci (m.rows (), acc);
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	246 for (octave_idx_type j = 0; j < m.columns (); j++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	247 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	248 for (octave_idx_type i = 0; i < m.rows (); i++)
8308 5fe0f4dfdbec use std::vector as a simple linear container in oct-norm.cc to avoid problems with instantiating Array<T> Jaroslav Hajek <highegg@gmail.com> parents: 8303 diff changeset	249 acci[i].accum (m(i, j));
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	250 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	251
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	252 for (octave_idx_type i = 0; i < m.rows (); i++)
8308 5fe0f4dfdbec use std::vector as a simple linear container in oct-norm.cc to avoid problems with instantiating Array<T> Jaroslav Hajek <highegg@gmail.com> parents: 8303 diff changeset	253 res.xelem (i) = acci[i];
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	254 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	255
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	256 // sparse versions
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	257 template <class T, class R, class ACC>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	258 void column_norms (const MSparse<T>& m, MArray<R>& res, ACC acc)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	259 {
11570 57632dea2446 attempt better backward compatibility for Array constructors John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	260 res = MArray<R> (dim_vector (1, m.columns ()));
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	261 for (octave_idx_type j = 0; j < m.columns (); j++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	262 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	263 ACC accj = acc;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	264 for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	265 accj.accum (m.data (k));
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	266
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	267 res.xelem (j) = accj;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	268 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	269 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	270
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	271 template <class T, class R, class ACC>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	272 void row_norms (const MSparse<T>& m, MArray<R>& res, ACC acc)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	273 {
11570 57632dea2446 attempt better backward compatibility for Array constructors John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	274 res = MArray<R> (dim_vector (m.rows (), 1));
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	275 std::vector<ACC> acci (m.rows (), acc);
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	276 for (octave_idx_type j = 0; j < m.columns (); j++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	277 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	278 for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++)
8308 5fe0f4dfdbec use std::vector as a simple linear container in oct-norm.cc to avoid problems with instantiating Array<T> Jaroslav Hajek <highegg@gmail.com> parents: 8303 diff changeset	279 acci[m.ridx (k)].accum (m.data (k));
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	280 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	281
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	282 for (octave_idx_type i = 0; i < m.rows (); i++)
8309 05b7a8ffec65 omission from last patch Jaroslav Hajek <highegg@gmail.com> parents: 8308 diff changeset	283 res.xelem (i) = acci[i];
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	284 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	285
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	286 // now the dispatchers
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	287 #define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE) \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	288 template <class T, class R> \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	289 RES_TYPE FUNC_NAME (const ARG_TYPE& v, R p) \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	290 { \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	291 RES_TYPE res; \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	292 if (p == 2) \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	293 FUNC_NAME (v, res, norm_accumulator_2<R> ()); \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	294 else if (p == 1) \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	295 FUNC_NAME (v, res, norm_accumulator_1<R> ()); \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	296 else if (lo_ieee_isinf (p)) \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	297 { \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	298 if (p > 0) \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	299 FUNC_NAME (v, res, norm_accumulator_inf<R> ()); \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	300 else \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	301 FUNC_NAME (v, res, norm_accumulator_minf<R> ()); \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	302 } \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	303 else if (p == 0) \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	304 FUNC_NAME (v, res, norm_accumulator_0<R> ()); \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	305 else if (p > 0) \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	306 FUNC_NAME (v, res, norm_accumulator_p<R> (p)); \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	307 else \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	308 FUNC_NAME (v, res, norm_accumulator_mp<R> (p)); \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	309 return res; \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	310 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	311
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	312 DEFINE_DISPATCHER (vector_norm, MArray<T>, R)
10350 12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10145 diff changeset	313 DEFINE_DISPATCHER (column_norms, MArray<T>, MArray<R>)
12884915a8e4 merge MArray classes & improve Array interface Jaroslav Hajek <highegg@gmail.com> parents: 10145 diff changeset	314 DEFINE_DISPATCHER (row_norms, MArray<T>, MArray<R>)
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	315 DEFINE_DISPATCHER (column_norms, MSparse<T>, MArray<R>)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	316 DEFINE_DISPATCHER (row_norms, MSparse<T>, MArray<R>)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	317
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	318 // The approximate subproblem in Higham's method. Find lambda and mu such that
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	319 // norm ([lambda, mu], p) == 1 and norm (ylambda + colmu, p) is maximized.
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	320 // Real version. As in Higham's paper.
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	321 template <class ColVectorT, class R>
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	322 static void
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	323 higham_subp (const ColVectorT& y, const ColVectorT& col,
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	324 octave_idx_type nsamp, R p, R& lambda, R& mu)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	325 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	326 R nrm = 0;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	327 for (octave_idx_type i = 0; i < nsamp; i++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	328 {
10145 fa01c1670b3e make p-norms breakable Jaroslav Hajek <highegg@gmail.com> parents: 10142 diff changeset	329 octave_quit ();
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	330 R fi = i*M_PI/nsamp, lambda1 = cos (fi), mu1 = sin (fi);
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	331 R lmnr = std::pow (std::pow (std::abs (lambda1), p) +
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	332 std::pow (std::abs (mu1), p), 1/p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	333 lambda1 /= lmnr; mu1 /= lmnr;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	334 R nrm1 = vector_norm (lambda1 * y + mu1 * col, p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	335 if (nrm1 > nrm)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	336 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	337 lambda = lambda1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	338 mu = mu1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	339 nrm = nrm1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	340 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	341 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	342 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	343
17769 49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	344 // Complex version. Higham's paper does not deal with complex case, so we use a
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	345 // simple extension. First, guess the magnitudes as in real version, then try
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	346 // to rotate lambda to improve further.
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	347 template <class ColVectorT, class R>
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	348 static void
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	349 higham_subp (const ColVectorT& y, const ColVectorT& col,
12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	350 octave_idx_type nsamp, R p,
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	351 std::complex<R>& lambda, std::complex<R>& mu)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	352 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	353 typedef std::complex<R> CR;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	354 R nrm = 0;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	355 lambda = 1.0;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	356 CR lamcu = lambda / std::abs (lambda);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	357 // Probe magnitudes
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	358 for (octave_idx_type i = 0; i < nsamp; i++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	359 {
10145 fa01c1670b3e make p-norms breakable Jaroslav Hajek <highegg@gmail.com> parents: 10142 diff changeset	360 octave_quit ();
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	361 R fi = i*M_PI/nsamp, lambda1 = cos (fi), mu1 = sin (fi);
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	362 R lmnr = std::pow (std::pow (std::abs (lambda1), p) +
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	363 std::pow (std::abs (mu1), p), 1/p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	364 lambda1 /= lmnr; mu1 /= lmnr;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	365 R nrm1 = vector_norm (lambda1 * lamcu * y + mu1 * col, p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	366 if (nrm1 > nrm)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	367 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	368 lambda = lambda1 * lamcu;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	369 mu = mu1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	370 nrm = nrm1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	371 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	372 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	373 R lama = std::abs (lambda);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	374 // Probe orientation
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	375 for (octave_idx_type i = 0; i < nsamp; i++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	376 {
10145 fa01c1670b3e make p-norms breakable Jaroslav Hajek <highegg@gmail.com> parents: 10142 diff changeset	377 octave_quit ();
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	378 R fi = i*M_PI/nsamp;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	379 lamcu = CR (cos (fi), sin (fi));
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	380 R nrm1 = vector_norm (lama * lamcu * y + mu * col, p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	381 if (nrm1 > nrm)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	382 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	383 lambda = lama * lamcu;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	384 nrm = nrm1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	385 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	386 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	387 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	388
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	389 // the p-dual element (should work for both real and complex)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	390 template <class T, class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	391 inline T elem_dual_p (T x, R p)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	392 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	393 return signum (x) * std::pow (std::abs (x), p-1);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	394 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	395
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	396 // the VectorT is used for vectors, but actually it has to be
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	397 // a Matrix type to allow all the operations. For instance SparseMatrix
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	398 // does not support multiplication with column/row vectors.
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	399 // the dual vector
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	400 template <class VectorT, class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	401 VectorT dual_p (const VectorT& x, R p, R q)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	402 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	403 VectorT res (x.dims ());
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	404 for (octave_idx_type i = 0; i < x.numel (); i++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	405 res.xelem (i) = elem_dual_p (x(i), p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	406 return res / vector_norm (res, q);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	407 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	408
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	409 // Higham's hybrid method
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	410 template <class MatrixT, class VectorT, class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	411 R higham (const MatrixT& m, R p, R tol, int maxiter,
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	412 VectorT& x)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	413 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	414 x.resize (m.columns (), 1);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	415 // the OSE part
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	416 VectorT y(m.rows (), 1, 0), z(m.rows (), 1);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	417 typedef typename VectorT::element_type RR;
14520 c4ed0fdf2b62 Fix initialization in Higham's method for norm estimate (bug #36031). Marco Caliari <marco.caliari@univr.it> parents: 14138 diff changeset	418 RR lambda = 0, mu = 1;
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	419 for (octave_idx_type k = 0; k < m.columns (); k++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	420 {
10142 829e69ec3110 make OCTAVE_QUIT a function Jaroslav Hajek <highegg@gmail.com> parents: 9245 diff changeset	421 octave_quit ();
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	422 VectorT col (m.column (k));
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	423 if (k > 0)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	424 higham_subp (y, col, 4*k, p, lambda, mu);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	425 for (octave_idx_type i = 0; i < k; i++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	426 x(i) *= lambda;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	427 x(k) = mu;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	428 y = lambda * y + mu * col;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	429 }
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	430
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	431 // the PM part
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	432 x = x / vector_norm (x, p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	433 R q = p/(p-1);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	434
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	435 R gamma = 0, gamma1;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	436 int iter = 0;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	437 while (iter < maxiter)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	438 {
10142 829e69ec3110 make OCTAVE_QUIT a function Jaroslav Hajek <highegg@gmail.com> parents: 9245 diff changeset	439 octave_quit ();
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	440 y = m*x;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	441 gamma1 = gamma;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	442 gamma = vector_norm (y, p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	443 z = dual_p (y, p, q);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	444 z = z.hermitian ();
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	445 z = z * m;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	446
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	447 if (iter > 0 && (vector_norm (z, q) <= gamma
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	448 \|\| (gamma - gamma1) <= tol*gamma))
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	449 break;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	450
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	451 z = z.hermitian ();
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	452 x = dual_p (z, q, p);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	453 iter ++;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	454 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	455
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	456 return gamma;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	457 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	458
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	459 // derive column vector and SVD types
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	460
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	461 static const char *p_less1_gripe = "xnorm: p must be at least 1";
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	462
17769 49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	463 // Static constant to control the maximum number of iterations. 100 seems to
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	464 // be a good value. Eventually, we can provide a means to change this
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/ Rik <rik@octave.org> parents: 17744 diff changeset	465 // constant from Octave.
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	466 static int max_norm_iter = 100;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	467
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	468 // version with SVD for dense matrices
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	469 template <class MatrixT, class VectorT, class SVDT, class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	470 R matrix_norm (const MatrixT& m, R p, VectorT, SVDT)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	471 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	472 R res = 0;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	473 if (p == 2)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	474 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	475 SVDT svd (m, SVD::sigma_only);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	476 res = svd.singular_values () (0,0);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	477 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	478 else if (p == 1)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	479 res = xcolnorms (m, 1).max ();
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	480 else if (lo_ieee_isinf (p))
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	481 res = xrownorms (m, 1).max ();
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	482 else if (p > 1)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	483 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	484 VectorT x;
8995 1b097d86a61a remove a TODO in oct-norm.cc Jaroslav Hajek <highegg@gmail.com> parents: 8319 diff changeset	485 const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
1b097d86a61a remove a TODO in oct-norm.cc Jaroslav Hajek <highegg@gmail.com> parents: 8319 diff changeset	486 res = higham (m, p, sqrteps, max_norm_iter, x);
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	487 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	488 else
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	489 (*current_liboctave_error_handler) (p_less1_gripe);
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	490
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	491 return res;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	492 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	493
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	494 // SVD-free version for sparse matrices
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	495 template <class MatrixT, class VectorT, class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	496 R matrix_norm (const MatrixT& m, R p, VectorT)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	497 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	498 R res = 0;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	499 if (p == 1)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	500 res = xcolnorms (m, 1).max ();
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	501 else if (lo_ieee_isinf (p))
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	502 res = xrownorms (m, 1).max ();
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	503 else if (p > 1)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	504 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	505 VectorT x;
8995 1b097d86a61a remove a TODO in oct-norm.cc Jaroslav Hajek <highegg@gmail.com> parents: 8319 diff changeset	506 const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
1b097d86a61a remove a TODO in oct-norm.cc Jaroslav Hajek <highegg@gmail.com> parents: 8319 diff changeset	507 res = higham (m, p, sqrteps, max_norm_iter, x);
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	508 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	509 else
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11570 diff changeset	510 (*current_liboctave_error_handler) (p_less1_gripe);
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	511
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	512 return res;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	513 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	514
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	515 // and finally, here's what we've promised in the header file
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	516
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	517 #define DEFINE_XNORM_FUNCS(PREFIX, RTYPE) \
8319 c374691576f6 Fix for MSVC compilation Michael Goffioul <michael.goffioul@gmail.com> parents: 8309 diff changeset	518 OCTAVE_API RTYPE xnorm (const PREFIX##ColumnVector& x, RTYPE p) \
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	519 { return vector_norm (x, p); } \
8319 c374691576f6 Fix for MSVC compilation Michael Goffioul <michael.goffioul@gmail.com> parents: 8309 diff changeset	520 OCTAVE_API RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	521 { return vector_norm (x, p); } \
8319 c374691576f6 Fix for MSVC compilation Michael Goffioul <michael.goffioul@gmail.com> parents: 8309 diff changeset	522 OCTAVE_API RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	523 { return matrix_norm (x, p, PREFIX##Matrix (), PREFIX##SVD ()); } \
8319 c374691576f6 Fix for MSVC compilation Michael Goffioul <michael.goffioul@gmail.com> parents: 8309 diff changeset	524 OCTAVE_API RTYPE xfrobnorm (const PREFIX##Matrix& x) \
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	525 { return vector_norm (x, static_cast<RTYPE> (2)); }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	526
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	527 DEFINE_XNORM_FUNCS(, double)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	528 DEFINE_XNORM_FUNCS(Complex, double)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	529 DEFINE_XNORM_FUNCS(Float, float)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	530 DEFINE_XNORM_FUNCS(FloatComplex, float)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	531
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	532 // this is needed to avoid copying the sparse matrix for xfrobnorm
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	533 template <class T, class R>
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	534 inline void array_norm_2 (const T* v, octave_idx_type n, R& res)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	535 {
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	536 norm_accumulator_2<R> acc;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	537 for (octave_idx_type i = 0; i < n; i++)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	538 acc.accum (v[i]);
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	539
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	540 res = acc;
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	541 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	542
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	543 #define DEFINE_XNORM_SPARSE_FUNCS(PREFIX, RTYPE) \
8319 c374691576f6 Fix for MSVC compilation Michael Goffioul <michael.goffioul@gmail.com> parents: 8309 diff changeset	544 OCTAVE_API RTYPE xnorm (const Sparse##PREFIX##Matrix& x, RTYPE p) \
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	545 { return matrix_norm (x, p, PREFIX##Matrix ()); } \
8319 c374691576f6 Fix for MSVC compilation Michael Goffioul <michael.goffioul@gmail.com> parents: 8309 diff changeset	546 OCTAVE_API RTYPE xfrobnorm (const Sparse##PREFIX##Matrix& x) \
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	547 { \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	548 RTYPE res; \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	549 array_norm_2 (x.data (), x.nnz (), res); \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	550 return res; \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	551 }
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	552
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	553 DEFINE_XNORM_SPARSE_FUNCS(, double)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	554 DEFINE_XNORM_SPARSE_FUNCS(Complex, double)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	555
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	556 #define DEFINE_COLROW_NORM_FUNCS(PREFIX, RPREFIX, RTYPE) \
8319 c374691576f6 Fix for MSVC compilation Michael Goffioul <michael.goffioul@gmail.com> parents: 8309 diff changeset	557 extern OCTAVE_API RPREFIX##RowVector xcolnorms (const PREFIX##Matrix& m, RTYPE p) \
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	558 { return column_norms (m, p); } \
8319 c374691576f6 Fix for MSVC compilation Michael Goffioul <michael.goffioul@gmail.com> parents: 8309 diff changeset	559 extern OCTAVE_API RPREFIX##ColumnVector xrownorms (const PREFIX##Matrix& m, RTYPE p) \
8303 b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	560 { return row_norms (m, p); } \
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	561
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	562 DEFINE_COLROW_NORM_FUNCS(, , double)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	563 DEFINE_COLROW_NORM_FUNCS(Complex, , double)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	564 DEFINE_COLROW_NORM_FUNCS(Float, Float, float)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	565 DEFINE_COLROW_NORM_FUNCS(FloatComplex, Float, float)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	566
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	567 DEFINE_COLROW_NORM_FUNCS(Sparse, , double)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	568 DEFINE_COLROW_NORM_FUNCS(SparseComplex, , double)
b11c31849b44 improve norm computation capabilities Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	569

Mercurial > hg > octave-nkf

annotate liboctave/numeric/oct-norm.cc @ 17769:49a5a4be04a1