GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
mgorth.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2009-2015 Carlo de Falco
4 Copyright (C) 2010 VZLU Prague
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include "oct-norm.h"
29 #include "defun.h"
30 #include "error.h"
31 #include "gripes.h"
32 
33 template <class ColumnVector, class Matrix, class RowVector>
34 static void
36 {
37  octave_idx_type Vc = V.columns ();
38  h = RowVector (Vc + 1);
39  for (octave_idx_type j = 0; j < Vc; j++)
40  {
41  ColumnVector Vcj = V.column (j);
42  h(j) = RowVector (Vcj.hermitian ()) * x;
43  x -= h(j) * Vcj;
44  }
45 
46  h(Vc) = xnorm (x);
47  if (real (h(Vc)) > 0)
48  x = x / h(Vc);
49 }
50 
51 DEFUN (mgorth, args, nargout,
52  "-*- texinfo -*-\n\
53 @deftypefn {Built-in Function} {[@var{y}, @var{h}] =} mgorth (@var{x}, @var{v})\n\
54 Orthogonalize a given column vector @var{x} with respect to a set of\n\
55 orthonormal vectors comprising the columns of @var{v} using the modified\n\
56 Gram-Schmidt method.\n\
57 \n\
58 On exit, @var{y} is a unit vector such that:\n\
59 \n\
60 @example\n\
61 @group\n\
62  norm (@var{y}) = 1\n\
63  @var{v}' * @var{y} = 0\n\
64  @var{x} = [@var{v}, @var{y}]*@var{h}'\n\
65 @end group\n\
66 @end example\n\
67 \n\
68 @end deftypefn")
69 {
70  octave_value_list retval;
71 
72  int nargin = args.length ();
73 
74  if (nargin != 2 || nargout > 2)
75  {
76  print_usage ();
77  return retval;
78  }
79 
80  octave_value arg_x = args(0);
81  octave_value arg_v = args(1);
82 
83  if (arg_v.ndims () != 2 || arg_x.ndims () != 2 || arg_x.columns () != 1
84  || arg_v.rows () != arg_x.rows ())
85  {
86  error ("mgorth: V should be a matrix, and X a column vector with"
87  " the same number of rows as V.");
88  return retval;
89  }
90 
91  if (! arg_x.is_numeric_type () && ! arg_v.is_numeric_type ())
92  {
93  error ("mgorth: X and V must be numeric");
94  }
95 
96  bool iscomplex = (arg_x.is_complex_type () || arg_v.is_complex_type ());
97  if (arg_x.is_single_type () || arg_v.is_single_type ())
98  {
99  if (iscomplex)
100  {
105  do_mgorth (x, V, h);
106  retval(1) = h;
107  retval(0) = x;
108  }
109  else
110  {
112  FloatMatrix V = arg_v.float_matrix_value ();
113  FloatRowVector h;
114  do_mgorth (x, V, h);
115  retval(1) = h;
116  retval(0) = x;
117  }
118  }
119  else
120  {
121  if (iscomplex)
122  {
126  do_mgorth (x, V, h);
127  retval(1) = h;
128  retval(0) = x;
129  }
130  else
131  {
133  Matrix V = arg_v.matrix_value ();
134  RowVector h;
135  do_mgorth (x, V, h);
136  retval(1) = h;
137  retval(0) = x;
138  }
139  }
140 
141  return retval;
142 }
143 
144 /*
145 %!test
146 %! for ii=1:100
147 %! assert (abs (mgorth (randn (5, 1), eye (5, 4))), [0 0 0 0 1]', eps);
148 %! endfor
149 
150 %!test
151 %! a = hilb (5);
152 %! a(:, 1) /= norm (a(:, 1));
153 %! for ii = 1:5
154 %! a(:, ii) = mgorth (a(:, ii), a(:, 1:ii-1));
155 %! endfor
156 %! assert (a' * a, eye (5), 1e10);
157 */
ColumnVector column_vector_value(bool frc_str_conv=false, bool frc_vec_conv=false) const
Definition: ov.cc:1658
int ndims(void) const
Definition: ov.h:479
octave_idx_type rows(void) const
Definition: ov.h:473
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:795
FloatColumnVector float_column_vector_value(bool frc_str_conv=false, bool frc_vec_conv=false) const
Definition: ov.cc:1869
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type length(void) const
Definition: oct-obj.h:89
bool is_numeric_type(void) const
Definition: ov.h:663
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
static void do_mgorth(ColumnVector &x, const Matrix &V, RowVector &h)
Definition: mgorth.cc:35
octave_idx_type columns(void) const
Definition: ov.h:475
bool is_complex_type(void) const
Definition: ov.h:654
OCTAVE_API double xnorm(const ColumnVector &x, double p)
Definition: oct-norm.cc:536
Definition: dMatrix.h:35
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:773
MArray< T > hermitian(T(*fcn)(const T &)=0) const
Definition: MArray.h:86
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type const double const double octave_idx_type double * V
Definition: CmplxGEPBAL.cc:49
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:791
FloatComplexColumnVector float_complex_column_vector_value(bool frc_str_conv=false, bool frc_vec_conv=false) const
Definition: ov.cc:1877
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:776
ColumnVector column(octave_idx_type i) const
Definition: dMatrix.cc:645
ComplexColumnVector complex_column_vector_value(bool frc_str_conv=false, bool frc_vec_conv=false) const
Definition: ov.cc:1666
bool is_single_type(void) const
Definition: ov.h:611
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:156
octave_idx_type columns(void) const
Definition: Array.h:322
F77_RET_T const double * x