GNU Octave
3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
Main Page
Namespaces
Classes
Files
File List
File Members
All
Classes
Namespaces
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Properties
Friends
Macros
Pages
liboctave
numeric
dbleQRP.cc
Go to the documentation of this file.
1
/*
2
3
Copyright (C) 1994-2013 John W. Eaton
4
Copyright (C) 2009 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 <cassert>
29
30
#include "
dbleQRP.h
"
31
#include "
f77-fcn.h
"
32
#include "
lo-error.h
"
33
#include "
oct-locbuf.h
"
34
35
extern
"C"
36
{
37
F77_RET_T
38
F77_FUNC
(dgeqp3, DGEQP3) (
const
octave_idx_type
&,
const
octave_idx_type&,
39
double
*,
const
octave_idx_type&,
40
octave_idx_type*,
double
*,
double
*,
41
const
octave_idx_type&, octave_idx_type&);
42
}
43
44
// It would be best to share some of this code with QR class...
45
46
QRP::QRP
(
const
Matrix
& a,
qr_type_t
qr_type)
47
:
QR
(), p ()
48
{
49
init
(a, qr_type);
50
}
51
52
void
53
QRP::init
(
const
Matrix
& a,
qr_type_t
qr_type)
54
{
55
assert (qr_type !=
qr_type_raw
);
56
57
octave_idx_type
m = a.
rows
();
58
octave_idx_type
n = a.
cols
();
59
60
octave_idx_type
min_mn = m < n ? m : n;
61
OCTAVE_LOCAL_BUFFER
(
double
, tau, min_mn);
62
63
octave_idx_type
info = 0;
64
65
Matrix
afact = a;
66
if
(m > n && qr_type ==
qr_type_std
)
67
afact.
resize
(m, m);
68
69
MArray<octave_idx_type>
jpvt (
dim_vector
(n, 1), 0);
70
71
if
(m > 0)
72
{
73
// workspace query.
74
double
rlwork;
75
F77_XFCN
(dgeqp3, DGEQP3, (m, n, afact.
fortran_vec
(),
76
m, jpvt.
fortran_vec
(), tau,
77
&rlwork, -1, info));
78
79
// allocate buffer and do the job.
80
octave_idx_type
lwork = rlwork;
81
lwork =
std::max
(lwork, static_cast<octave_idx_type> (1));
82
OCTAVE_LOCAL_BUFFER
(
double
, work, lwork);
83
F77_XFCN
(dgeqp3, DGEQP3, (m, n, afact.
fortran_vec
(),
84
m, jpvt.
fortran_vec
(), tau,
85
work, lwork, info));
86
}
87
else
88
for
(
octave_idx_type
i = 0; i < n; i++) jpvt(i) = i+1;
89
90
// Form Permutation matrix (if economy is requested, return the
91
// indices only!)
92
93
jpvt -=
static_cast<
octave_idx_type
>
(1);
94
p
=
PermMatrix
(jpvt,
true
);
95
96
97
form
(n, afact, tau, qr_type);
98
}
99
100
RowVector
101
QRP::Pvec
(
void
)
const
102
{
103
Array<double>
pa (
p
.
pvec
());
104
RowVector
pv (
MArray<double>
(pa) + 1.0);
105
return
pv;
106
}
Generated on Mon Dec 30 2013 03:04:48 for GNU Octave by
1.8.1.2