annotate src/qr.cc @ 1315:611d403c7f3d

[project @ 1995-06-25 19:56:32 by jwe]
author jwe
date Sun, 25 Jun 1995 19:56:32 +0000
parents b6360f2d4fa6
children 19c10b8657d5
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-qr.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
1009
dfe01093f657 [project @ 1995-01-04 04:05:12 by jwe]
jwe
parents: 718
diff changeset
4 Copyright (C) 1993, 1994, 1995 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
1315
611d403c7f3d [project @ 1995-06-25 19:56:32 by jwe]
jwe
parents: 1192
diff changeset
20 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
1
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
1192
b6360f2d4fa6 [project @ 1995-03-30 21:38:35 by jwe]
jwe
parents: 1009
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 "dbleQR.h"
393e95f46b51 [project @ 1994-06-06 00:05:20 by jwe]
jwe
parents: 240
diff changeset
29 #include "CmplxQR.h"
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
30
539
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
31 #include "dbleQRP.h"
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
32 #include "CmplxQRP.h"
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
33
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
34 #include "tree-const.h"
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
35 #include "user-prefs.h"
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
36 #include "gripes.h"
636
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
37 #include "utils.h"
544
20fbad23ae51 [project @ 1994-07-22 05:04:44 by jwe]
jwe
parents: 539
diff changeset
38 #include "help.h"
519
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
39 #include "defun-dld.h"
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
40
701
0a81458ef677 [project @ 1994-09-15 02:23:24 by jwe]
jwe
parents: 636
diff changeset
41 DEFUN_DLD_BUILTIN ("qr", Fqr, Sqr, 2, 2,
539
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
42 "[Q, R] = qr (X): form Q unitary and R upper triangular such\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
43 that Q * R = X\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
44 \n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
45 [Q, R] = qr (X, 0): form the economy decomposition such that if X is\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
46 if X is m by n then only the first n columns of Q\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
47 are computed.\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
48 \n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
49 [Q, R, P] = qr (X): form QRP factorization of X where\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
50 P is a permutation matrix such that\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
51 A * P = Q * R\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
52 \n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
53 [Q, R, P] = qr (X, 0): form the economy decomposition with \n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
54 permutation vector P such that Q * R = X (:, P)\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
55 \n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
56 qr (X) alone returns the output of the LAPACK routine dgeqrf, such\n\
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
57 that R = triu (qr (X))")
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
58 {
519
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
59 Octave_object retval;
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
60
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
61 int nargin = args.length ();
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
62
712
36ba0576bd1b [project @ 1994-09-19 14:18:15 by jwe]
jwe
parents: 701
diff changeset
63 if (nargin != 1 && nargin != 2 || nargout > 3)
519
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
64 {
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
65 print_usage ("qr");
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
66 return retval;
b9284136189a [project @ 1994-07-19 14:40:20 by jwe]
jwe
parents: 516
diff changeset
67 }
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
68
712
36ba0576bd1b [project @ 1994-09-19 14:18:15 by jwe]
jwe
parents: 701
diff changeset
69 tree_constant arg = args(0);
718
e81d3a66725e [project @ 1994-09-21 14:58:18 by jwe]
jwe
parents: 712
diff changeset
70
e81d3a66725e [project @ 1994-09-21 14:58:18 by jwe]
jwe
parents: 712
diff changeset
71 int arg_is_empty = empty_arg ("qr", arg.rows (), arg.columns ());
e81d3a66725e [project @ 1994-09-21 14:58:18 by jwe]
jwe
parents: 712
diff changeset
72
e81d3a66725e [project @ 1994-09-21 14:58:18 by jwe]
jwe
parents: 712
diff changeset
73 if (arg_is_empty < 0)
636
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
74 return retval;
718
e81d3a66725e [project @ 1994-09-21 14:58:18 by jwe]
jwe
parents: 712
diff changeset
75 else if (arg_is_empty > 0)
e81d3a66725e [project @ 1994-09-21 14:58:18 by jwe]
jwe
parents: 712
diff changeset
76 return Octave_object (3, Matrix ());
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
77
539
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
78 QR::type type = nargout == 1 ? QR::raw
712
36ba0576bd1b [project @ 1994-09-19 14:18:15 by jwe]
jwe
parents: 701
diff changeset
79 : (nargin == 2 ? QR::economy : QR::std);
539
5ec10a984241 [project @ 1994-07-21 22:40:04 by jwe]
jwe
parents: 519
diff changeset
80
636
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
81 if (arg.is_real_type ())
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
82 {
636
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
83 Matrix m = arg.matrix_value ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
84
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
85 if (! error_state)
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
86 {
636
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
87 if (nargout < 3)
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
88 {
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
89 QR fact (m, type);
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
90 retval(1) = fact.R ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
91 retval(0) = fact.Q ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
92 }
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
93 else
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
94 {
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
95 QRP fact (m, type);
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
96 retval(2) = fact.P ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
97 retval(1) = fact.R ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
98 retval(0) = fact.Q ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
99 }
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
100 }
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
101 }
636
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
102 else if (arg.is_complex_type ())
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
103 {
636
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
104 ComplexMatrix m = arg.complex_matrix_value ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
105
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
106 if (! error_state)
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
107 {
636
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
108 if (nargout < 3)
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
109 {
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
110 ComplexQR fact (m, type);
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
111 retval(1) = fact.R ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
112 retval(0) = fact.Q ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
113 }
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
114 else
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
115 {
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
116 ComplexQRP fact (m, type);
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
117 retval(2) = fact.P ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
118 retval(1) = fact.R ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
119 retval(0) = fact.Q ();
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
120 }
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
121 }
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
122 }
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
123 else
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
124 {
636
fae2bd91c027 [project @ 1994-08-23 18:39:50 by jwe]
jwe
parents: 620
diff changeset
125 gripe_wrong_type_arg ("qr", arg);
620
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
126 }
8e4e7e5f307e [project @ 1994-08-16 04:36:32 by jwe]
jwe
parents: 544
diff changeset
127
1
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
128 return retval;
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
129 }
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
130
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
131 /*
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
132 ;;; Local Variables: ***
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
133 ;;; mode: C++ ***
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
134 ;;; page-delimiter: "^/\\*" ***
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
135 ;;; End: ***
78fd87e624cb [project @ 1993-08-08 01:13:40 by jwe]
jwe
parents:
diff changeset
136 */