|
3211
|
1 # Copyright (C) 1992, 1998 A. Scottedward Hodel |
|
|
2 # |
|
|
3 # This file is part of Octave. |
|
|
4 # |
|
|
5 # Octave is free software; you can redistribute it and/or modify it |
|
|
6 # under the terms of the GNU General Public License as published by the |
|
|
7 # Free Software Foundation; either version 2, or (at your option) any |
|
|
8 # later version. |
|
|
9 # |
|
|
10 # Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
11 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
12 # FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
13 # for more details. |
|
|
14 # |
|
|
15 # You should have received a copy of the GNU General Public License |
|
|
16 # along with Octave; see the file COPYING. If not, write to the Free |
|
|
17 # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. |
|
|
18 |
|
|
19 function [hv,alph,kb] = qrhouse(VV,eps1) |
|
3238
|
20 # function [hv,alph,kb] = qrhouse(VV{,eps1}) |
|
3211
|
21 # construct orthogonal basis of span(VV) with Householder vectors |
|
|
22 # Q R = VV; Q may be obtained via routine krygetq; R is upper triangular |
|
|
23 # if all rows of VV are nonzero; otherwise it's a permuted uppert |
|
|
24 # triangular matrix with zero rows matching those of VV. |
|
|
25 # Inputs: |
|
|
26 # VV: matrix |
|
|
27 # eps1: zero threshhold; default: 0. Used to check if a column |
|
|
28 # of reduced R is already upper triangular; entries with magnitude |
|
|
29 # <= eps1 are considered to be 0. |
|
|
30 # Outputs: |
|
|
31 # hv: matrix of householder reflection vectors as returned by housh |
|
|
32 # alpha: vector of householder reflection values as returned by housh |
|
|
33 # kb: computed rank of matrix |
|
|
34 # qrhouse is used in krylovb for block Arnoldi iteration |
|
|
35 # |
|
|
36 # Reference: Golub and Van Loan, MATRIX COMPUTATIONS, 2nd ed. |
|
|
37 |
|
|
38 # Written by A. S. Hodel, 1992 |
|
|
39 |
|
3238
|
40 if(nargin < 1 | nargin > 2) |
|
|
41 usage("[hv,alph,kb] = qrhouse(VV{,eps1})"); |
|
|
42 elseif(nargin == 1) # default value for eps set to 0 |
|
|
43 eps1 = 0; |
|
3211
|
44 endif |
|
|
45 |
|
3237
|
46 |
|
3211
|
47 # extract only those rows of VV that are nonzero |
|
|
48 if(is_vector(VV)) nzidx = find(abs(VV') > 0); |
|
3237
|
49 else nzidx = find(max(abs(VV')) > 0); endif |
|
3211
|
50 VVlen = rows(VV); # remember original size |
|
|
51 |
|
|
52 if(is_vector(VV)) VV = VV(nzidx); |
|
3237
|
53 else VV = VV(nzidx,:); endif |
|
3211
|
54 |
|
|
55 [Vr,Vc] = size(VV); nits = min(Vr,Vc); |
|
|
56 for ii = 1:nits; |
|
3237
|
57 # permute maximum row entry to (ii,ii) position |
|
|
58 Vrowi = VV(ii,1:Vc); # pivot maximum entry in this row to lead position |
|
|
59 Vrm = max(abs(Vrowi)); |
|
|
60 Vmidx = min(find(abs(Vrowi) == Vrm)); |
|
|
61 if(Vmidx > eps1) |
|
|
62 kb = ii; # update computed rank |
|
|
63 idx = kb-1; |
|
|
64 if(Vmidx != ii) |
|
|
65 [VV(:,kb),VV(:,Vmidx)] = swap(VV(:,kb),VV(:,Vmidx)); |
|
|
66 endif |
|
|
67 hh = VV(:,ii); # extract next column of VV; ignore items 1:(ii-1). |
|
3211
|
68 [hv(kb:Vr,kb),alph(kb),z] = housh(hh(kb:Vr),1,0); |
|
|
69 if(kb>1) |
|
3237
|
70 hv(1:idx,kb) = 0; # zero top of hv for safety |
|
3211
|
71 endif |
|
3237
|
72 # project off of current Householder vector |
|
3211
|
73 VV = VV - alph(kb)*hv(:,kb)*(hv(:,kb)'*VV); |
|
|
74 else |
|
3237
|
75 break; |
|
|
76 endif |
|
3211
|
77 endfor |
|
|
78 if(kb <=0) |
|
|
79 hv = []; |
|
|
80 alph = []; |
|
|
81 else |
|
|
82 hvs = hv(:,1:kb); # remove extraneous vectors, expand to original size |
|
|
83 hv = zeros(VVlen,kb); |
|
|
84 hv(nzidx,:) = hvs; |
|
|
85 alph = alph(1:kb); |
|
|
86 end |
|
|
87 endfunction |
