Mercurial > octave
annotate liboctave/numeric/oct-norm.cc @ 21275:85d8280c64f4
Return NaN for norm (..., +/-Inf) if input contains NaN (bug #32855).
* oct-norm.cc (norm_accumulator_inf::accum, norm_accumulator_minf::accum):
Use xisnan to check for NaN value before continuing to accumulate.
author | Rik <rik@octave.org> |
---|---|
date | Tue, 16 Feb 2016 16:01:56 -0800 |
parents | cbced1c09916 |
children | 40de9f8f23a6 |
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
25 #ifdef HAVE_CONFIG_H |
21202
f7121e111991
maint: indent #ifdef blocks in liboctave and src directories.
Rik <rik@octave.org>
parents:
21139
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 |
7a756a7e145b
Add a citation to Higham's matrix norm paper
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11586
diff
changeset
|
65 // 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
|
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 { |
21275
85d8280c64f4
Return NaN for norm (..., +/-Inf) if input contains NaN (bug #32855).
Rik <rik@octave.org>
parents:
21273
diff
changeset
|
182 if (xisnan (val)) |
85d8280c64f4
Return NaN for norm (..., +/-Inf) if input contains NaN (bug #32855).
Rik <rik@octave.org>
parents:
21273
diff
changeset
|
183 max = octave_NaN; |
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: |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
196 norm_accumulator_minf () : min (octave_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 { |
21275
85d8280c64f4
Return NaN for norm (..., +/-Inf) if input contains NaN (bug #32855).
Rik <rik@octave.org>
parents:
21273
diff
changeset
|
200 if (xisnan (val)) |
85d8280c64f4
Return NaN for norm (..., +/-Inf) if input contains NaN (bug #32855).
Rik <rik@octave.org>
parents:
21273
diff
changeset
|
201 min = octave_NaN; |
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
224 // OK, we're armed :) Now let's go for the fun |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
225 |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
20714
diff
changeset
|
226 template <typename T, typename R, typename ACC> |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
227 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
|
228 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
229 for (octave_idx_type i = 0; i < v.numel (); i++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
230 acc.accum (v(i)); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
231 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
232 res = acc; |
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
235 // dense versions |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
20714
diff
changeset
|
236 template <typename T, typename R, typename ACC> |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10145
diff
changeset
|
237 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
|
238 { |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
239 res = MArray<R> (dim_vector (1, m.columns ())); |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
240 for (octave_idx_type j = 0; j < m.columns (); j++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
241 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
242 ACC accj = acc; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
243 for (octave_idx_type i = 0; i < m.rows (); i++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
244 accj.accum (m(i, j)); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
245 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
246 res.xelem (j) = accj; |
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 } |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
249 |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
20714
diff
changeset
|
250 template <typename T, typename R, typename ACC> |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10145
diff
changeset
|
251 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
|
252 { |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
253 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
|
254 std::vector<ACC> acci (m.rows (), acc); |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
255 for (octave_idx_type j = 0; j < m.columns (); j++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
256 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
257 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
|
258 acci[i].accum (m(i, j)); |
8303
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
261 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
|
262 res.xelem (i) = acci[i]; |
8303
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
265 // sparse versions |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
20714
diff
changeset
|
266 template <typename T, typename R, typename ACC> |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
267 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
|
268 { |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
269 res = MArray<R> (dim_vector (1, m.columns ())); |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
270 for (octave_idx_type j = 0; j < m.columns (); j++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
271 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
272 ACC accj = acc; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
273 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
|
274 accj.accum (m.data (k)); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
275 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
276 res.xelem (j) = accj; |
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 } |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
279 |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
20714
diff
changeset
|
280 template <typename T, typename R, typename ACC> |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
281 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
|
282 { |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
283 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
|
284 std::vector<ACC> acci (m.rows (), acc); |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
285 for (octave_idx_type j = 0; j < m.columns (); j++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
286 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
287 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
|
288 acci[m.ridx (k)].accum (m.data (k)); |
8303
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
291 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
|
292 res.xelem (i) = acci[i]; |
8303
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
295 // now the dispatchers |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
296 #define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE) \ |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
20714
diff
changeset
|
297 template <typename T, typename R> \ |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
298 RES_TYPE FUNC_NAME (const ARG_TYPE& v, R p) \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
299 { \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
300 RES_TYPE res; \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
301 if (p == 2) \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
302 FUNC_NAME (v, res, norm_accumulator_2<R> ()); \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
303 else if (p == 1) \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
304 FUNC_NAME (v, res, norm_accumulator_1<R> ()); \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
305 else if (lo_ieee_isinf (p)) \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
306 { \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
307 if (p > 0) \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
308 FUNC_NAME (v, res, norm_accumulator_inf<R> ()); \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
309 else \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
310 FUNC_NAME (v, res, norm_accumulator_minf<R> ()); \ |
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 else if (p == 0) \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
313 FUNC_NAME (v, res, norm_accumulator_0<R> ()); \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
314 else if (p > 0) \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
315 FUNC_NAME (v, res, norm_accumulator_p<R> (p)); \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
316 else \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
317 FUNC_NAME (v, res, norm_accumulator_mp<R> (p)); \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
318 return res; \ |
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
321 DEFINE_DISPATCHER (vector_norm, MArray<T>, R) |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10145
diff
changeset
|
322 DEFINE_DISPATCHER (column_norms, MArray<T>, MArray<R>) |
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10145
diff
changeset
|
323 DEFINE_DISPATCHER (row_norms, MArray<T>, MArray<R>) |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
324 DEFINE_DISPATCHER (column_norms, MSparse<T>, MArray<R>) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
325 DEFINE_DISPATCHER (row_norms, MSparse<T>, MArray<R>) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
326 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
327 // 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
|
328 // norm ([lambda, mu], p) == 1 and norm (y*lambda + col*mu, p) is maximized. |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
329 // 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
|
330 template <typename ColVectorT, typename R> |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
331 static void |
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
332 higham_subp (const ColVectorT& y, const ColVectorT& col, |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
333 octave_idx_type nsamp, R p, R& lambda, R& mu) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
334 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
335 R nrm = 0; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
336 for (octave_idx_type i = 0; i < nsamp; i++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
337 { |
10145
fa01c1670b3e
make p-norms breakable
Jaroslav Hajek <highegg@gmail.com>
parents:
10142
diff
changeset
|
338 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
|
339 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
|
340 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
|
341 R mu1 = sin (fi); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
342 R lmnr = std::pow (std::pow (std::abs (lambda1), p) + |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
343 std::pow (std::abs (mu1), p), 1/p); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
344 lambda1 /= lmnr; mu1 /= lmnr; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
345 R nrm1 = vector_norm (lambda1 * y + mu1 * col, p); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
346 if (nrm1 > nrm) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
347 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
348 lambda = lambda1; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
349 mu = mu1; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
350 nrm = nrm1; |
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 } |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
354 |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
355 // 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
|
356 // 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
|
357 // 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
|
358 template <typename ColVectorT, typename R> |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
359 static void |
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
360 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
|
361 octave_idx_type nsamp, R p, |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
362 std::complex<R>& lambda, std::complex<R>& mu) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
363 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
364 typedef std::complex<R> CR; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
365 R nrm = 0; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
366 lambda = 1.0; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
367 CR lamcu = lambda / std::abs (lambda); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
368 // Probe magnitudes |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
369 for (octave_idx_type i = 0; i < nsamp; i++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
370 { |
10145
fa01c1670b3e
make p-norms breakable
Jaroslav Hajek <highegg@gmail.com>
parents:
10142
diff
changeset
|
371 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
|
372 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
|
373 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
|
374 R mu1 = sin (fi); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
375 R lmnr = std::pow (std::pow (std::abs (lambda1), p) + |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
376 std::pow (std::abs (mu1), p), 1/p); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
377 lambda1 /= lmnr; mu1 /= lmnr; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
378 R nrm1 = vector_norm (lambda1 * lamcu * y + mu1 * col, p); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
379 if (nrm1 > nrm) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
380 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
381 lambda = lambda1 * lamcu; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
382 mu = mu1; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
383 nrm = nrm1; |
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 } |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
386 R lama = std::abs (lambda); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
387 // Probe orientation |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
388 for (octave_idx_type i = 0; i < nsamp; i++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
389 { |
10145
fa01c1670b3e
make p-norms breakable
Jaroslav Hajek <highegg@gmail.com>
parents:
10142
diff
changeset
|
390 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
|
391 R fi = i * static_cast<R> (M_PI) / nsamp; |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
392 lamcu = CR (cos (fi), sin (fi)); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
393 R nrm1 = vector_norm (lama * lamcu * y + mu * col, p); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
394 if (nrm1 > nrm) |
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 lambda = lama * lamcu; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
397 nrm = nrm1; |
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
402 // 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
|
403 template <typename T, typename R> |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
404 inline T elem_dual_p (T x, R p) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
405 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
406 return signum (x) * std::pow (std::abs (x), p-1); |
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 // 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
|
410 // a Matrix type to allow all the operations. For instance SparseMatrix |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
411 // does not support multiplication with column/row vectors. |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
412 // the dual vector |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
20714
diff
changeset
|
413 template <typename VectorT, typename R> |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
414 VectorT dual_p (const VectorT& x, R p, R q) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
415 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
416 VectorT res (x.dims ()); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
417 for (octave_idx_type i = 0; i < x.numel (); i++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
418 res.xelem (i) = elem_dual_p (x(i), p); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
419 return res / vector_norm (res, q); |
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
422 // 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
|
423 template <typename MatrixT, typename VectorT, typename R> |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
424 R higham (const MatrixT& m, R p, R tol, int maxiter, |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
425 VectorT& x) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
426 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
427 x.resize (m.columns (), 1); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
428 // the OSE part |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
429 VectorT y(m.rows (), 1, 0), z(m.rows (), 1); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
430 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
|
431 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
|
432 RR mu = 1; |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
433 for (octave_idx_type k = 0; k < m.columns (); k++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
434 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
9245
diff
changeset
|
435 octave_quit (); |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
436 VectorT col (m.column (k)); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
437 if (k > 0) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
438 higham_subp (y, col, 4*k, p, lambda, mu); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
439 for (octave_idx_type i = 0; i < k; i++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
440 x(i) *= lambda; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
441 x(k) = mu; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
442 y = lambda * y + mu * col; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
443 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
444 |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
445 // the PM part |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
446 x = x / vector_norm (x, p); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
447 R q = p/(p-1); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
448 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
449 R gamma = 0, gamma1; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
450 int iter = 0; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
451 while (iter < maxiter) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
452 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
9245
diff
changeset
|
453 octave_quit (); |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
454 y = m*x; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
455 gamma1 = gamma; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
456 gamma = vector_norm (y, p); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
457 z = dual_p (y, p, q); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
458 z = z.hermitian (); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
459 z = z * m; |
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 if (iter > 0 && (vector_norm (z, q) <= gamma |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
462 || (gamma - gamma1) <= tol*gamma)) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
463 break; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
464 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
465 z = z.hermitian (); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
466 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
|
467 iter++; |
8303
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
470 return gamma; |
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 |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
473 // derive column vector and SVD types |
8303
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 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
|
476 |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
477 // 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
|
478 // 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
|
479 // constant from Octave. |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
480 static int max_norm_iter = 100; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
481 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
482 // 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
|
483 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
|
484 R svd_matrix_norm (const MatrixT& m, R p, VectorT) |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
485 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
486 R res = 0; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
487 if (p == 2) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
488 { |
21273
cbced1c09916
better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents:
21202
diff
changeset
|
489 svd<MatrixT> fact (m, svd<MatrixT>::sigma_only); |
cbced1c09916
better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents:
21202
diff
changeset
|
490 res = fact.singular_values () (0,0); |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
491 } |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
492 else if (p == 1) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
493 res = xcolnorms (m, 1).max (); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
494 else if (lo_ieee_isinf (p)) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
495 res = xrownorms (m, 1).max (); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
496 else if (p > 1) |
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 VectorT x; |
8995
1b097d86a61a
remove a TODO in oct-norm.cc
Jaroslav Hajek <highegg@gmail.com>
parents:
8319
diff
changeset
|
499 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
|
500 res = higham (m, p, sqrteps, max_norm_iter, x); |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
501 } |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
502 else |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
503 (*current_liboctave_error_handler) (p_less1_gripe); |
8303
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 return res; |
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
508 // 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
|
509 template <typename MatrixT, typename VectorT, typename R> |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
510 R matrix_norm (const MatrixT& m, R p, VectorT) |
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 R res = 0; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
513 if (p == 1) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
514 res = xcolnorms (m, 1).max (); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
515 else if (lo_ieee_isinf (p)) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
516 res = xrownorms (m, 1).max (); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
517 else if (p > 1) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
518 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
519 VectorT x; |
8995
1b097d86a61a
remove a TODO in oct-norm.cc
Jaroslav Hajek <highegg@gmail.com>
parents:
8319
diff
changeset
|
520 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
|
521 res = higham (m, p, sqrteps, max_norm_iter, x); |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
522 } |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
523 else |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
524 (*current_liboctave_error_handler) (p_less1_gripe); |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
525 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
526 return res; |
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
529 // 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
|
530 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
531 #define DEFINE_XNORM_FUNCS(PREFIX, RTYPE) \ |
8319
c374691576f6
Fix for MSVC compilation
Michael Goffioul <michael.goffioul@gmail.com>
parents:
8309
diff
changeset
|
532 OCTAVE_API RTYPE xnorm (const PREFIX##ColumnVector& x, RTYPE p) \ |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
533 { return vector_norm (x, p); } \ |
8319
c374691576f6
Fix for MSVC compilation
Michael Goffioul <michael.goffioul@gmail.com>
parents:
8309
diff
changeset
|
534 OCTAVE_API RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \ |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
535 { return vector_norm (x, p); } \ |
8319
c374691576f6
Fix for MSVC compilation
Michael Goffioul <michael.goffioul@gmail.com>
parents:
8309
diff
changeset
|
536 OCTAVE_API RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \ |
21273
cbced1c09916
better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents:
21202
diff
changeset
|
537 { return svd_matrix_norm (x, p, PREFIX##Matrix ()); } \ |
8319
c374691576f6
Fix for MSVC compilation
Michael Goffioul <michael.goffioul@gmail.com>
parents:
8309
diff
changeset
|
538 OCTAVE_API RTYPE xfrobnorm (const PREFIX##Matrix& x) \ |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
539 { return vector_norm (x, static_cast<RTYPE> (2)); } |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
540 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
541 DEFINE_XNORM_FUNCS(, double) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
542 DEFINE_XNORM_FUNCS(Complex, double) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
543 DEFINE_XNORM_FUNCS(Float, float) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
544 DEFINE_XNORM_FUNCS(FloatComplex, float) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
545 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
546 // 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
|
547 template <typename T, typename R> |
8303
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
548 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
|
549 { |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
550 norm_accumulator_2<R> acc; |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
551 for (octave_idx_type i = 0; i < n; i++) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
552 acc.accum (v[i]); |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
553 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
554 res = acc; |
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 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
557 #define DEFINE_XNORM_SPARSE_FUNCS(PREFIX, RTYPE) \ |
8319
c374691576f6
Fix for MSVC compilation
Michael Goffioul <michael.goffioul@gmail.com>
parents:
8309
diff
changeset
|
558 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
|
559 { return matrix_norm (x, p, PREFIX##Matrix ()); } \ |
8319
c374691576f6
Fix for MSVC compilation
Michael Goffioul <michael.goffioul@gmail.com>
parents:
8309
diff
changeset
|
560 OCTAVE_API RTYPE xfrobnorm (const Sparse##PREFIX##Matrix& x) \ |
8303
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 RTYPE res; \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
563 array_norm_2 (x.data (), x.nnz (), res); \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
564 return res; \ |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
565 } |
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_XNORM_SPARSE_FUNCS(, double) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
568 DEFINE_XNORM_SPARSE_FUNCS(Complex, double) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
569 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
570 #define DEFINE_COLROW_NORM_FUNCS(PREFIX, RPREFIX, RTYPE) \ |
8319
c374691576f6
Fix for MSVC compilation
Michael Goffioul <michael.goffioul@gmail.com>
parents:
8309
diff
changeset
|
571 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
|
572 { return column_norms (m, p); } \ |
8319
c374691576f6
Fix for MSVC compilation
Michael Goffioul <michael.goffioul@gmail.com>
parents:
8309
diff
changeset
|
573 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
|
574 { return row_norms (m, p); } \ |
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_COLROW_NORM_FUNCS(, , double) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
577 DEFINE_COLROW_NORM_FUNCS(Complex, , double) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
578 DEFINE_COLROW_NORM_FUNCS(Float, Float, float) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
579 DEFINE_COLROW_NORM_FUNCS(FloatComplex, Float, float) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
580 |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
581 DEFINE_COLROW_NORM_FUNCS(Sparse, , double) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
582 DEFINE_COLROW_NORM_FUNCS(SparseComplex, , double) |
b11c31849b44
improve norm computation capabilities
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
583 |