annotate src/svd.cc @ 620:8e4e7e5f307e

[project @ 1994-08-16 04:36:32 by jwe]
author jwe
date Tue, 16 Aug 1994 04:36:32 +0000
parents 20fbad23ae51
children 6636e6198f01
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7
d68036bcad4c [project @ 1993-08-08 03:53:43 by jwe]
jwe
parents: 1
diff changeset
1 // f-svd.cc -*- C++ -*-
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
2 /*
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
3
453
393e95f46b51 [project @ 1994-06-06 00:05:20 by jwe]
jwe
parents: 240
diff changeset
4 Copyright (C) 1993, 1994 John W. Eaton
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
5
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
6 This file is part of Octave.
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
7
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
8 Octave is free software; you can redistribute it and/or modify it
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
9 under the terms of the GNU General Public License as published by the
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
10 Free Software Foundation; either version 2, or (at your option) any
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
11 later version.
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
12
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
13 Octave is distributed in the hope that it will be useful, but WITHOUT
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
16 for more details.
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
17
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
18 You should have received a copy of the GNU General Public License
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
19 along with Octave; see the file COPYING. If not, write to the Free
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
21
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
22 */
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
23
240
a99f28f5e351 [project @ 1993-11-30 20:24:36 by jwe]
jwe
parents: 162
diff changeset
24 #ifdef HAVE_CONFIG_H
a99f28f5e351 [project @ 1993-11-30 20:24:36 by jwe]
jwe
parents: 162
diff changeset
25 #include "config.h"
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
26 #endif
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
27
453
393e95f46b51 [project @ 1994-06-06 00:05:20 by jwe]
jwe
parents: 240
diff changeset
28 #include "dbleSVD.h"
393e95f46b51 [project @ 1994-06-06 00:05:20 by jwe]
jwe
parents: 240
diff changeset
29 #include "CmplxSVD.h"
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
30
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
31 #include "tree-const.h"
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
32 #include "user-prefs.h"
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
33 #include "gripes.h"
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
34 #include "error.h"
544
20fbad23ae51 [project @ 1994-07-22 05:04:44 by jwe]
jwe
parents: 539
diff changeset
35 #include "help.h"
519
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
36 #include "defun-dld.h"
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
37
519
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
38 DEFUN_DLD ("svd", Fsvd, Ssvd, 2, 3,
539
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
39 "S = svd (X) or [U, S, V] = svd (X [, 0])\n\
519
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
40 \n\
539
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
41 Compute the singular value decomposition of X. Given a second input\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
42 argument, an `economy' sized factorization is computed that omits\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
43 unnecessary rows and columns of U and V")
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
44 {
497
88614b380d6e [project @ 1994-07-08 02:00:57 by jwe]
jwe
parents: 453
diff changeset
45 Octave_object retval;
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
46
519
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
47 int nargin = args.length ();
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
48
539
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
49 if (nargin < 2 || nargin > 3 || nargout == 2 || nargout > 3)
519
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
50 {
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
51 print_usage ("svd");
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
52 return retval;
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
53 }
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
54
497
88614b380d6e [project @ 1994-07-08 02:00:57 by jwe]
jwe
parents: 453
diff changeset
55 tree_constant arg = args(1).make_numeric ();
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
56
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
57 if (error_state)
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
58 return retval;
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
59
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
60 if (arg.rows () == 0 || arg.columns () == 0)
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
61 {
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
62 int flag = user_pref.propagate_empty_matrices;
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
63 if (flag != 0)
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
64 {
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
65 if (flag < 0)
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
66 gripe_empty_arg ("svd", 0);
516
309fc59f66ee [project @ 1994-07-13 02:31:31 by jwe]
jwe
parents: 506
diff changeset
67
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
68 Matrix m;
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
69 retval.resize (3, m);
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
70 }
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
71 else
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
72 gripe_empty_arg ("svd", 1);
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
73
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
74 return retval;
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
75 }
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
76
539
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
77 SVD::type type = (nargin == 3) ? SVD::economy : SVD::std;
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
78
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
79 if (arg.is_real_type ())
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
80 {
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
81 Matrix tmp = arg.matrix_value ();
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
82
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
83 SVD result (tmp, type);
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
84
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
85 DiagMatrix sigma = result.singular_values ();
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
86
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
87 if (nargout == 0 || nargout == 1)
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
88 {
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
89 retval(0) = tree_constant (sigma.diag (), 1);
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
90 }
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
91 else
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
92 {
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
93 retval(2) = result.right_singular_matrix ();
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
94 retval(1) = sigma;
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
95 retval(0) = result.left_singular_matrix ();
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
96 }
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
97 }
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
98 else if (arg.is_complex_type ())
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
99 {
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
100 ComplexMatrix ctmp = arg.complex_matrix_value ();
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
101
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
102 ComplexSVD result (ctmp, type);
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
103
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
104 DiagMatrix sigma = result.singular_values ();
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
105
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
106 if (nargout == 0 || nargout == 1)
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
107 {
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
108 retval(0) = tree_constant (sigma.diag (), 1);
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
109 }
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
110 else
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
111 {
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
112 retval(2) = result.right_singular_matrix ();
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
113 retval(1) = sigma;
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
114 retval(0) = result.left_singular_matrix ();
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
115 }
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
116 }
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
117 else
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
118 {
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
119 gripe_wrong_type_arg ("svd", arg);
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
120 return retval;
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
121 }
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
122
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
123 return retval;
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
124 }
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
125
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
126 /*
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
127 ;;; Local Variables: ***
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
128 ;;; mode: C++ ***
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
129 ;;; page-delimiter: "^/\\*" ***
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
130 ;;; End: ***
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
131 */