annotate liboctave/numeric/oct-norm.cc @ 22204:469c817eb256

svd: reduce code duplication with more use of template and macro. * liboctave/numeric/svd.cc, liboctave/numeric/svd.h: remove unused constructor with reference for int (info). This allows to move all of the constructor into a single template, so remove init(). Two new methods, gesvd and gesdd, are fully specialized but the main hunck of code are the long list of arguments. Scope type and drive enums to the svd class for clarity, and rename member names. Add a new member for the drive used. * libinterp/corefcn/svd.cc: fix typenames for the svd enums which are now scoped. * CMatrix.cc, dMatrix.cc, fCMatrix.cc, fMatrix.cc: fix typenames for the svd enums which are now scoped.
author Carnë Draug <carandraug@octave.org>
date Thu, 04 Aug 2016 20:20:27 +0100
parents e43d83253e28
children 6ca3acf5fad8
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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
19697
4197fc428c7d maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents: 19269
diff changeset
3 Copyright (C) 2008-2015 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
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21723
diff changeset
25 #if defined (HAVE_CONFIG_H)
21301
40de9f8f23a6 Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents: 21275
diff changeset
26 # include "config.h"
8303
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
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
36 #include "Array-util.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
37 #include "Array.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
38 #include "CColVector.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
39 #include "CMatrix.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
40 #include "CRowVector.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
41 #include "CSparse.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
42 #include "dColVector.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
43 #include "dDiagMatrix.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
44 #include "dMatrix.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
45 #include "dRowVector.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
46 #include "dSparse.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
47 #include "fCColVector.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
48 #include "fCMatrix.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
49 #include "fCRowVector.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
50 #include "fColVector.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
51 #include "fDiagMatrix.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
52 #include "fMatrix.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
53 #include "fRowVector.h"
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
54 #include "lo-error.h"
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
55 #include "lo-ieee.h"
19269
65554f5847ac don't include oct-locbuf.h in header files unnecessarily
John W. Eaton <jwe@octave.org>
parents: 18084
diff changeset
56 #include "mx-cm-s.h"
65554f5847ac don't include oct-locbuf.h in header files unnecessarily
John W. Eaton <jwe@octave.org>
parents: 18084
diff changeset
57 #include "mx-fcm-fs.h"
65554f5847ac don't include oct-locbuf.h in header files unnecessarily
John W. Eaton <jwe@octave.org>
parents: 18084
diff changeset
58 #include "mx-fs-fcm.h"
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
59 #include "mx-s-cm.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
60 #include "oct-cmplx.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
61 #include "svd.h"
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
62
13934
7a756a7e145b Add a citation to Higham's matrix norm paper
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 11586
diff changeset
63 // 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
64 // to handle both real and complex element, and a cast operator
21751
b571fc85953f maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents: 21724
diff changeset
65 // returning the intermediate norm. Reference: Higham, N. "Estimating
13934
7a756a7e145b Add a citation to Higham's matrix norm paper
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 11586
diff changeset
66 // the Matrix p-Norm." Numer. Math. 62, 539-555, 1992.
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
67
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
68 // norm accumulator for the p-norm
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
69 template <typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
70 class norm_accumulator_p
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
71 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
72 R p,scl,sum;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
73 public:
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
74 norm_accumulator_p () {} // we need this one for Array
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
75 norm_accumulator_p (R pp) : p(pp), scl(0), sum(1) {}
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
76
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
77 template <typename U>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
78 void accum (U val)
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
79 {
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
80 octave_quit ();
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
81 R t = std::abs (val);
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
82 if (scl == t) // we need this to handle Infs properly
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
83 sum += 1;
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
84 else if (scl < t)
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
85 {
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
86 sum *= std::pow (scl/t, p);
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
87 sum += 1;
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
88 scl = t;
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
89 }
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
90 else if (t != 0)
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
91 sum += std::pow (t/scl, p);
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
92 }
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
93 operator R () { return scl * std::pow (sum, 1/p); }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
94 };
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
95
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
96 // norm accumulator for the minus p-pseudonorm
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
97 template <typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
98 class norm_accumulator_mp
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
99 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
100 R p,scl,sum;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
101 public:
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
102 norm_accumulator_mp () {} // we need this one for Array
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
103 norm_accumulator_mp (R pp) : p(pp), scl(0), sum(1) {}
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
104
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
105 template <typename U>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
106 void accum (U val)
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
107 {
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
108 octave_quit ();
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
109 R t = 1 / std::abs (val);
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
110 if (scl == t)
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
111 sum += 1;
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
112 else if (scl < t)
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
113 {
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
114 sum *= std::pow (scl/t, p);
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
115 sum += 1;
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
116 scl = t;
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
117 }
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
118 else if (t != 0)
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
119 sum += std::pow (t/scl, p);
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
120 }
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
121 operator R () { return scl * std::pow (sum, -1/p); }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
122 };
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
123
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
124 // norm accumulator for the 2-norm (euclidean)
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
125 template <typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
126 class norm_accumulator_2
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
127 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
128 R scl,sum;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
129 static R pow2 (R x) { return x*x; }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
130 public:
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
131 norm_accumulator_2 () : scl(0), sum(1) {}
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
132
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
133 void accum (R val)
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
134 {
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
135 R t = std::abs (val);
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
136 if (scl == t)
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
137 sum += 1;
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
138 else if (scl < t)
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
139 {
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
140 sum *= pow2 (scl/t);
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
141 sum += 1;
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
142 scl = t;
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
143 }
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
144 else if (t != 0)
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
145 sum += pow2 (t/scl);
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
146 }
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
147
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
148 void accum (std::complex<R> val)
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
149 {
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
150 accum (val.real ());
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
151 accum (val.imag ());
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
152 }
8303
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 operator R () { return scl * std::sqrt (sum); }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
155 };
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
156
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
157 // norm accumulator for the 1-norm (city metric)
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
158 template <typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
159 class norm_accumulator_1
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
160 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
161 R sum;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
162 public:
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
163 norm_accumulator_1 () : sum (0) {}
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
164 template <typename U>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
165 void accum (U val)
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
166 {
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
167 sum += std::abs (val);
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
168 }
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
169 operator R () { return sum; }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
170 };
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
171
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
172 // norm accumulator for the inf-norm (max metric)
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
173 template <typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
174 class norm_accumulator_inf
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
175 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
176 R max;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
177 public:
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
178 norm_accumulator_inf () : max (0) {}
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
179 template <typename U>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
180 void accum (U val)
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
181 {
21782
2aef506f3fec use namespace for lo-mappers.h functions
John W. Eaton <jwe@octave.org>
parents: 21751
diff changeset
182 if (octave::math::isnan (val))
21723
bae585228161 use namespace for numeric_limits
John W. Eaton <jwe@octave.org>
parents: 21721
diff changeset
183 max = octave::numeric_limits<R>::NaN ();
21275
85d8280c64f4 Return NaN for norm (..., +/-Inf) if input contains NaN (bug #32855).
Rik <rik@octave.org>
parents: 21273
diff changeset
184 else
85d8280c64f4 Return NaN for norm (..., +/-Inf) if input contains NaN (bug #32855).
Rik <rik@octave.org>
parents: 21273
diff changeset
185 max = std::max (max, std::abs (val));
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
186 }
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
187 operator R () { return max; }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
188 };
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
189
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
190 // norm accumulator for the -inf pseudonorm (min abs value)
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
191 template <typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
192 class norm_accumulator_minf
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 R min;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
195 public:
21723
bae585228161 use namespace for numeric_limits
John W. Eaton <jwe@octave.org>
parents: 21721
diff changeset
196 norm_accumulator_minf () : min (octave::numeric_limits<R>::Inf ()) {}
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
197 template <typename U>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
198 void accum (U val)
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
199 {
21782
2aef506f3fec use namespace for lo-mappers.h functions
John W. Eaton <jwe@octave.org>
parents: 21751
diff changeset
200 if (octave::math::isnan (val))
21723
bae585228161 use namespace for numeric_limits
John W. Eaton <jwe@octave.org>
parents: 21721
diff changeset
201 min = octave::numeric_limits<R>::NaN ();
21275
85d8280c64f4 Return NaN for norm (..., +/-Inf) if input contains NaN (bug #32855).
Rik <rik@octave.org>
parents: 21273
diff changeset
202 else
85d8280c64f4 Return NaN for norm (..., +/-Inf) if input contains NaN (bug #32855).
Rik <rik@octave.org>
parents: 21273
diff changeset
203 min = std::min (min, std::abs (val));
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
204 }
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
205 operator R () { return min; }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
206 };
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
207
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
208 // norm accumulator for the 0-pseudonorm (hamming distance)
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
209 template <typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
210 class norm_accumulator_0
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
211 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
212 unsigned int num;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
213 public:
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
214 norm_accumulator_0 () : num (0) {}
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
215 template <typename U>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
216 void accum (U val)
19861
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
217 {
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
218 if (val != static_cast<U> (0)) ++num;
19755f4fc851 maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
219 }
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
220 operator R () { return num; }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
221 };
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 // OK, we're armed :) Now let's go for the fun
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
224
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
225 template <typename T, typename R, typename ACC>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
226 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
227 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
228 for (octave_idx_type i = 0; i < v.numel (); i++)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
229 acc.accum (v(i));
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
230
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
231 res = acc;
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
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
234 // dense versions
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
235 template <typename T, typename R, typename ACC>
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10145
diff changeset
236 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
237 {
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
238 res = MArray<R> (dim_vector (1, m.columns ()));
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
239 for (octave_idx_type j = 0; j < m.columns (); j++)
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 ACC accj = acc;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
242 for (octave_idx_type i = 0; i < m.rows (); i++)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
243 accj.accum (m(i, j));
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
244
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
245 res.xelem (j) = accj;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
246 }
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
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
249 template <typename T, typename R, typename ACC>
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10145
diff changeset
250 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
251 {
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
252 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
253 std::vector<ACC> acci (m.rows (), acc);
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
254 for (octave_idx_type j = 0; j < m.columns (); j++)
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 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
257 acci[i].accum (m(i, j));
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
258 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
259
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
260 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
261 res.xelem (i) = acci[i];
8303
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
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
264 // sparse versions
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
265 template <typename T, typename R, typename ACC>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
266 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
267 {
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
268 res = MArray<R> (dim_vector (1, m.columns ()));
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
269 for (octave_idx_type j = 0; j < m.columns (); j++)
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 ACC accj = acc;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
272 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
273 accj.accum (m.data (k));
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
274
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
275 res.xelem (j) = accj;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
276 }
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
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
279 template <typename T, typename R, typename ACC>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
280 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
281 {
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
282 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
283 std::vector<ACC> acci (m.rows (), acc);
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
284 for (octave_idx_type j = 0; j < m.columns (); j++)
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 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
287 acci[m.ridx (k)].accum (m.data (k));
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
288 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
289
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
290 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
291 res.xelem (i) = acci[i];
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
292 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
293
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
294 // now the dispatchers
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
295 #define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
296 template <typename T, typename R> \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
297 RES_TYPE FUNC_NAME (const ARG_TYPE& v, R p) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
298 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
299 RES_TYPE res; \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
300 if (p == 2) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
301 FUNC_NAME (v, res, norm_accumulator_2<R> ()); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
302 else if (p == 1) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
303 FUNC_NAME (v, res, norm_accumulator_1<R> ()); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
304 else if (lo_ieee_isinf (p)) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
305 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
306 if (p > 0) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
307 FUNC_NAME (v, res, norm_accumulator_inf<R> ()); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
308 else \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
309 FUNC_NAME (v, res, norm_accumulator_minf<R> ()); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
310 } \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
311 else if (p == 0) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
312 FUNC_NAME (v, res, norm_accumulator_0<R> ()); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
313 else if (p > 0) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
314 FUNC_NAME (v, res, norm_accumulator_p<R> (p)); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
315 else \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
316 FUNC_NAME (v, res, norm_accumulator_mp<R> (p)); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
317 return res; \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
318 }
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
319
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
320 DEFINE_DISPATCHER (vector_norm, MArray<T>, R)
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10145
diff changeset
321 DEFINE_DISPATCHER (column_norms, MArray<T>, MArray<R>)
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10145
diff changeset
322 DEFINE_DISPATCHER (row_norms, MArray<T>, MArray<R>)
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
323 DEFINE_DISPATCHER (column_norms, MSparse<T>, MArray<R>)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
324 DEFINE_DISPATCHER (row_norms, MSparse<T>, MArray<R>)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
325
21751
b571fc85953f maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents: 21724
diff changeset
326 // The approximate subproblem in Higham's method. Find lambda and mu such that
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
327 // norm ([lambda, mu], p) == 1 and norm (y*lambda + col*mu, p) is maximized.
21751
b571fc85953f maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents: 21724
diff changeset
328 // Real version. As in Higham's paper.
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
329 template <typename ColVectorT, typename R>
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
330 static void
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
331 higham_subp (const ColVectorT& y, const ColVectorT& col,
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
332 octave_idx_type nsamp, R p, R& lambda, R& mu)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
333 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
334 R nrm = 0;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
335 for (octave_idx_type i = 0; i < nsamp; i++)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
336 {
10145
fa01c1670b3e make p-norms breakable
Jaroslav Hajek <highegg@gmail.com>
parents: 10142
diff changeset
337 octave_quit ();
19924
277b12eed117 oct-norm.cc: Use static_cast<R> (M_PI) for cases where R is not double.
Rik <rik@octave.org>
parents: 19861
diff changeset
338 R fi = i * static_cast<R> (M_PI) / nsamp;
18084
8e056300994b Follow coding convention of defining and initializing only 1 variable per line in liboctave.
Rik <rik@octave.org>
parents: 17769
diff changeset
339 R lambda1 = cos (fi);
8e056300994b Follow coding convention of defining and initializing only 1 variable per line in liboctave.
Rik <rik@octave.org>
parents: 17769
diff changeset
340 R mu1 = sin (fi);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
341 R lmnr = std::pow (std::pow (std::abs (lambda1), p) +
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
342 std::pow (std::abs (mu1), p), 1/p);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
343 lambda1 /= lmnr; mu1 /= lmnr;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
344 R nrm1 = vector_norm (lambda1 * y + mu1 * col, p);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
345 if (nrm1 > nrm)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
346 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
347 lambda = lambda1;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
348 mu = mu1;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
349 nrm = nrm1;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
350 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
351 }
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
21751
b571fc85953f maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents: 21724
diff changeset
354 // Complex version. Higham's paper does not deal with complex case, so we use
b571fc85953f maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents: 21724
diff changeset
355 // a simple extension. First, guess the magnitudes as in real version, then
b571fc85953f maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents: 21724
diff changeset
356 // try to rotate lambda to improve further.
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
357 template <typename ColVectorT, typename R>
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
358 static void
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
359 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
360 octave_idx_type nsamp, R p,
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
361 std::complex<R>& lambda, std::complex<R>& mu)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
362 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
363 typedef std::complex<R> CR;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
364 R nrm = 0;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
365 lambda = 1.0;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
366 CR lamcu = lambda / std::abs (lambda);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
367 // Probe magnitudes
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
368 for (octave_idx_type i = 0; i < nsamp; i++)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
369 {
10145
fa01c1670b3e make p-norms breakable
Jaroslav Hajek <highegg@gmail.com>
parents: 10142
diff changeset
370 octave_quit ();
19924
277b12eed117 oct-norm.cc: Use static_cast<R> (M_PI) for cases where R is not double.
Rik <rik@octave.org>
parents: 19861
diff changeset
371 R fi = i * static_cast<R> (M_PI) / nsamp;
18084
8e056300994b Follow coding convention of defining and initializing only 1 variable per line in liboctave.
Rik <rik@octave.org>
parents: 17769
diff changeset
372 R lambda1 = cos (fi);
8e056300994b Follow coding convention of defining and initializing only 1 variable per line in liboctave.
Rik <rik@octave.org>
parents: 17769
diff changeset
373 R mu1 = sin (fi);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
374 R lmnr = std::pow (std::pow (std::abs (lambda1), p) +
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
375 std::pow (std::abs (mu1), p), 1/p);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
376 lambda1 /= lmnr; mu1 /= lmnr;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
377 R nrm1 = vector_norm (lambda1 * lamcu * y + mu1 * col, p);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
378 if (nrm1 > nrm)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
379 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
380 lambda = lambda1 * lamcu;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
381 mu = mu1;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
382 nrm = nrm1;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
383 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
384 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
385 R lama = std::abs (lambda);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
386 // Probe orientation
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
387 for (octave_idx_type i = 0; i < nsamp; i++)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
388 {
10145
fa01c1670b3e make p-norms breakable
Jaroslav Hajek <highegg@gmail.com>
parents: 10142
diff changeset
389 octave_quit ();
19924
277b12eed117 oct-norm.cc: Use static_cast<R> (M_PI) for cases where R is not double.
Rik <rik@octave.org>
parents: 19861
diff changeset
390 R fi = i * static_cast<R> (M_PI) / nsamp;
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
391 lamcu = CR (cos (fi), sin (fi));
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
392 R nrm1 = vector_norm (lama * lamcu * y + mu * col, p);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
393 if (nrm1 > nrm)
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 lambda = lama * lamcu;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
396 nrm = nrm1;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
397 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
398 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
399 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
400
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
401 // the p-dual element (should work for both real and complex)
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
402 template <typename T, typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
403 inline T elem_dual_p (T x, R p)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
404 {
21782
2aef506f3fec use namespace for lo-mappers.h functions
John W. Eaton <jwe@octave.org>
parents: 21751
diff changeset
405 return octave::math::signum (x) * std::pow (std::abs (x), p-1);
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
406 }
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 // the VectorT is used for vectors, but actually it has to be
21751
b571fc85953f maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents: 21724
diff changeset
409 // a Matrix type to allow all the operations. For instance SparseMatrix
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
410 // does not support multiplication with column/row vectors.
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
411 // the dual vector
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
412 template <typename VectorT, typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
413 VectorT dual_p (const VectorT& x, R p, R q)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
414 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
415 VectorT res (x.dims ());
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
416 for (octave_idx_type i = 0; i < x.numel (); i++)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
417 res.xelem (i) = elem_dual_p (x(i), p);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
418 return res / vector_norm (res, q);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
419 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
420
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
421 // Higham's hybrid method
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
422 template <typename MatrixT, typename VectorT, typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
423 R higham (const MatrixT& m, R p, R tol, int maxiter,
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
424 VectorT& x)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
425 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
426 x.resize (m.columns (), 1);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
427 // the OSE part
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
428 VectorT y(m.rows (), 1, 0), z(m.rows (), 1);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
429 typedef typename VectorT::element_type RR;
18084
8e056300994b Follow coding convention of defining and initializing only 1 variable per line in liboctave.
Rik <rik@octave.org>
parents: 17769
diff changeset
430 RR lambda = 0;
8e056300994b Follow coding convention of defining and initializing only 1 variable per line in liboctave.
Rik <rik@octave.org>
parents: 17769
diff changeset
431 RR mu = 1;
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
432 for (octave_idx_type k = 0; k < m.columns (); k++)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
433 {
10142
829e69ec3110 make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents: 9245
diff changeset
434 octave_quit ();
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
435 VectorT col (m.column (k));
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
436 if (k > 0)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
437 higham_subp (y, col, 4*k, p, lambda, mu);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
438 for (octave_idx_type i = 0; i < k; i++)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
439 x(i) *= lambda;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
440 x(k) = mu;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
441 y = lambda * y + mu * col;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
442 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
443
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
444 // the PM part
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
445 x = x / vector_norm (x, p);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
446 R q = p/(p-1);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
447
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
448 R gamma = 0, gamma1;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
449 int iter = 0;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
450 while (iter < maxiter)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
451 {
10142
829e69ec3110 make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents: 9245
diff changeset
452 octave_quit ();
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
453 y = m*x;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
454 gamma1 = gamma;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
455 gamma = vector_norm (y, p);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
456 z = dual_p (y, p, q);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
457 z = z.hermitian ();
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
458 z = z * m;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
459
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
460 if (iter > 0 && (vector_norm (z, q) <= gamma
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
461 || (gamma - gamma1) <= tol*gamma))
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
462 break;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
463
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
464 z = z.hermitian ();
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
465 x = dual_p (z, q, p);
20714
7b6d8c19dab0 Cuddle increment (++) and decrement (--) operators with variable in question.
Rik <rik@octave.org>
parents: 19924
diff changeset
466 iter++;
8303
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
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
469 return gamma;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
470 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
471
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
472 // derive column vector and SVD types
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
473
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
474 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
475
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
476 // 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
477 // 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
478 // constant from Octave.
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
479 static int max_norm_iter = 100;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
480
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
481 // version with SVD for dense matrices
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
482 template <typename MatrixT, typename VectorT, typename R>
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
483 R svd_matrix_norm (const MatrixT& m, R p, VectorT)
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
484 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
485 R res = 0;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
486 if (p == 2)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
487 {
22204
469c817eb256 svd: reduce code duplication with more use of template and macro.
Carnë Draug <carandraug@octave.org>
parents: 22197
diff changeset
488 svd<MatrixT> fact (m, svd<MatrixT>::Type::sigma_only);
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
489 res = fact.singular_values () (0,0);
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 else if (p == 1)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
492 res = xcolnorms (m, 1).max ();
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
493 else if (lo_ieee_isinf (p))
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
494 res = xrownorms (m, 1).max ();
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
495 else if (p > 1)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
496 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
497 VectorT x;
8995
1b097d86a61a remove a TODO in oct-norm.cc
Jaroslav Hajek <highegg@gmail.com>
parents: 8319
diff changeset
498 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
499 res = higham (m, p, sqrteps, max_norm_iter, x);
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
500 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
501 else
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
502 (*current_liboctave_error_handler) (p_less1_gripe);
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
503
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
504 return res;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
505 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
506
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
507 // SVD-free version for sparse matrices
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
508 template <typename MatrixT, typename VectorT, typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
509 R matrix_norm (const MatrixT& m, R p, VectorT)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
510 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
511 R res = 0;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
512 if (p == 1)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
513 res = xcolnorms (m, 1).max ();
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
514 else if (lo_ieee_isinf (p))
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
515 res = xrownorms (m, 1).max ();
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
516 else if (p > 1)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
517 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
518 VectorT x;
8995
1b097d86a61a remove a TODO in oct-norm.cc
Jaroslav Hajek <highegg@gmail.com>
parents: 8319
diff changeset
519 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
520 res = higham (m, p, sqrteps, max_norm_iter, x);
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
521 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
522 else
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
523 (*current_liboctave_error_handler) (p_less1_gripe);
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
524
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
525 return res;
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
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
528 // 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
529
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
530 #define DEFINE_XNORM_FUNCS(PREFIX, RTYPE) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
531 OCTAVE_API RTYPE xnorm (const PREFIX##ColumnVector& x, RTYPE p) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
532 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
533 return vector_norm (x, p); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
534 } \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
535 OCTAVE_API RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
536 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
537 return vector_norm (x, p); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
538 } \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
539 OCTAVE_API RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
540 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
541 return svd_matrix_norm (x, p, PREFIX##Matrix ()); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
542 } \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
543 OCTAVE_API RTYPE xfrobnorm (const PREFIX##Matrix& x) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
544 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
545 return vector_norm (x, static_cast<RTYPE> (2)); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
546 }
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 DEFINE_XNORM_FUNCS(, double)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
549 DEFINE_XNORM_FUNCS(Complex, double)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
550 DEFINE_XNORM_FUNCS(Float, float)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
551 DEFINE_XNORM_FUNCS(FloatComplex, float)
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 // this is needed to avoid copying the sparse matrix for xfrobnorm
21139
538b57866b90 consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents: 20714
diff changeset
554 template <typename T, typename R>
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
555 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
556 {
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
557 norm_accumulator_2<R> acc;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
558 for (octave_idx_type i = 0; i < n; i++)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
559 acc.accum (v[i]);
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
560
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
561 res = acc;
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
562 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
563
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
564 #define DEFINE_XNORM_SPARSE_FUNCS(PREFIX, RTYPE) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
565 OCTAVE_API RTYPE xnorm (const Sparse##PREFIX##Matrix& x, RTYPE p) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
566 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
567 return matrix_norm (x, p, PREFIX##Matrix ()); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
568 } \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
569 OCTAVE_API RTYPE xfrobnorm (const Sparse##PREFIX##Matrix& x) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
570 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
571 RTYPE res; \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
572 array_norm_2 (x.data (), x.nnz (), res); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
573 return res; \
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
574 }
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
575
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
576 DEFINE_XNORM_SPARSE_FUNCS(, double)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
577 DEFINE_XNORM_SPARSE_FUNCS(Complex, double)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
578
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
579 #define DEFINE_COLROW_NORM_FUNCS(PREFIX, RPREFIX, RTYPE) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
580 extern OCTAVE_API RPREFIX##RowVector \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
581 xcolnorms (const PREFIX##Matrix& m, RTYPE p) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
582 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
583 return column_norms (m, p); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
584 } \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
585 extern OCTAVE_API RPREFIX##ColumnVector \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
586 xrownorms (const PREFIX##Matrix& m, RTYPE p) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
587 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
588 return row_norms (m, p); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 22022
diff changeset
589 } \
8303
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
590
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
591 DEFINE_COLROW_NORM_FUNCS(, , double)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
592 DEFINE_COLROW_NORM_FUNCS(Complex, , double)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
593 DEFINE_COLROW_NORM_FUNCS(Float, Float, float)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
594 DEFINE_COLROW_NORM_FUNCS(FloatComplex, Float, float)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
595
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
596 DEFINE_COLROW_NORM_FUNCS(Sparse, , double)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
597 DEFINE_COLROW_NORM_FUNCS(SparseComplex, , double)
b11c31849b44 improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
598