|
3211
|
1 # Copyright (C) 1993, 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 [Uret,Ucols] = krylovb(A,V,k,eps1); |
|
|
20 # function [U,Ucols] = krylovb(A,V,k[,eps1]); |
|
|
21 # construct orthogonal basis U of block Krylov subspace; |
|
|
22 # [V AV A^2*V ... A^(k+1)*V]; |
|
|
23 # method used: householder reflections to guard against loss of |
|
|
24 # orthogonality |
|
|
25 # eps1: threshhold for 0 (default: 1e-12) |
|
|
26 # |
|
|
27 # outputs: |
|
|
28 # returned basis U is orthogonal matrix; due to "zeroed" |
|
|
29 # columns of product, may not satisfy A U = U H identity |
|
|
30 # Ucols: dimension of span of krylov subspace (based on eps1) |
|
|
31 # if k > m-1, krylov returns the Hessenberg decompostion of A. |
|
|
32 # |
|
|
33 # A. S. Hodel Feb 1993 |
|
|
34 # $Revision$ |
|
|
35 # $Log$ |
|
|
36 |
|
|
37 if(nargin == 3) |
|
|
38 eps1 = 1e-12; |
|
|
39 endif |
|
|
40 na = is_square(A); |
|
|
41 if( !na ) |
|
|
42 error("krylov: first argument must be a square matrix") |
|
|
43 endif |
|
|
44 |
|
|
45 [m,kb] = size(V); |
|
|
46 if(m != na); |
|
|
47 error("krylov: argument dimensions do not match") |
|
|
48 endif |
|
|
49 |
|
|
50 if( !is_scalar(k) ) |
|
|
51 error("krylov: third argument must be a scalar integer") |
|
|
52 endif |
|
|
53 |
|
|
54 Vnrm = norm(V,Inf); |
|
|
55 |
|
|
56 # compute factored QR |
|
|
57 [U,alpha,kb] = qrhouse(V,eps1*Vnrm); |
|
|
58 Q = krygetq(U,alpha,kb); |
|
|
59 Uret = Q; |
|
|
60 Ucols = kb; |
|
|
61 |
|
|
62 iter = 0; |
|
|
63 while (Ucols < na) & (kb > 0) & (iter < k) |
|
|
64 iter++; |
|
|
65 # multiply to get new vector; |
|
|
66 V = A*Q; |
|
|
67 |
|
|
68 # project off of previous vectors |
|
|
69 nzv = []; # set of reflection indices used so far |
|
|
70 nj = length(alpha); |
|
|
71 for i=1:nj |
|
|
72 hv = U(:,i); |
|
|
73 nzidx = find(hv != 0); # extract nonzero entries |
|
|
74 nzv = union(nzv,nzidx(1)); # save for call to qrhouse |
|
|
75 hv = hv(nzidx); |
|
|
76 vr = V(nzidx,:); |
|
|
77 av = alpha(i); |
|
|
78 V(nzidx,:) = vr - (av*hv)*(hv'*vr); |
|
|
79 end |
|
|
80 V(nzv,:) = 0; # clear entries covered by earlier reflections |
|
|
81 |
|
|
82 [hvk,alphk,kb] = qrhouse(V,eps1*Vnrm); # now get new QR factorization |
|
|
83 if(kb > 0) |
|
|
84 U = [U,hvk]; # append new data to Householder data structure |
|
|
85 alpha = [alpha;alphk]; |
|
|
86 Q2 = krygetq(U,alpha,kb); |
|
|
87 Uret = [Uret,Q2]; |
|
|
88 |
|
|
89 Ucols = Ucols + kb; |
|
|
90 Q = Q2; |
|
|
91 end |
|
|
92 end |
|
|
93 endfunction |
