# HG changeset patch # User jwe # Date 910239394 0 # Node ID 440b2b28e74afc9c34172df769090d79e97d5ec4 # Parent 258d7b26c719c039fa842f1f6c1b616c73195b91 [project @ 1998-11-05 04:16:22 by jwe] diff --git a/scripts/ChangeLog b/scripts/ChangeLog --- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,15 @@ +Wed Nov 4 21:51:13 1998 John W. Eaton + + * linear-algebra/housh.m: New file from the OCST. + * linear-algebra/krygetq.m: Ditto. + * linear-algebra/krylov.m: Ditto. + * linear-algebra/krylovb.m: Ditto. + * linear-algebra/qrhouse.m: Ditto. + * general/is_duplicate_entry.m: Ditto. + + * general/is_symmetric.m: Call is_square instead of doing that + check in line. + Wed Oct 28 11:51:14 1998 John W. Eaton * general/is_square.m: diff --git a/scripts/linear-algebra/housh.m b/scripts/linear-algebra/housh.m new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/housh.m @@ -0,0 +1,64 @@ +# Copyright (C) 1995, 1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 2, or (at your option) any +# later version. +# +# Octave is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License +# along with Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +function [housv,beta,zer] = housh(x,j,z) + # function [housv,beta,zer] = housh(x,j,z) + # computes householder reflection vector housv to reflect x to be + # jth column of identity, i.e., (I - beta*housv*housv')x =e(j) + # inputs + # x: vector + # j: index into vector + # z: threshold for zero (usually should be the number 0) + # outputs: (see Golub and Van Loan) + # beta: If beta = 0, then no reflection need be applied (zer set to 0) + # housv: householder vector + # mar 6,1987 : rev dec 17,1988 + # rev sep 19,1991 (blas) + # translated from FORTRAN Aug 1995 + # A. S. Hodel + + # $Revision: 1.1.1.1 $ + # $Log$ + + # check for legal inputs + if( !is_vector(x) && !is_scalar(x)) + error("housh: first input must be a vector") + elseif( !is_scalar(j) ) + error("housh: second argment must be an integer scalar") + else + housv = x; + m = max(abs(housv)); + if (m ~= 0.0) + housv = housv/m; + alpha = norm(housv); + if (alpha > z) + beta = 1.0/(alpha*(alpha+abs(housv(j)))); + sg = sign(housv(j)); + if( sg == 0) + sg = 1; + endif + housv(j) = housv(j) + alpha*sg; + else + beta = 0.0; + endif + else + beta = 0.0; + endif + zer = (beta == 0); + endif +endfunction diff --git a/scripts/linear-algebra/krygetq.m b/scripts/linear-algebra/krygetq.m new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/krygetq.m @@ -0,0 +1,60 @@ +# Copyright (C) 1993, 1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 2, or (at your option) any +# later version. +# +# Octave is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License +# along with Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +function QQ = krygetq(HH,alpha,kb) +# function QQ = krygetq(HH,alpha,kb) +# used internally by krylovb +# construct last kb columns of orthogonal transform returned by qrhouse +# Inputs: +# HH,alpha: set of nh Householder reflections set returned by qrhouse +# kb: desired number of columns +# Outputs: +# QQ: n x kb orthogonal matrix; basis of orthogonal subspace provided +# by final kb reflections in HH. +# Note: due to details in qrhouse and krylovb, QQ is *not* necessarily +# the last kb columns defined by HH. span(QQ) is orthogonal to the +# subspace spanned by the first (nh-kb) reflections in HH. + +# Written by A. S. Hodel 1993--1995 +# $Revision$ +# $Log$ + +[n,m] = size(HH); +kb1 = m+1-kb; + +# construct permuted identity on the right rows: +# since qrhouse removes zero rows, we've got to check the rows to which +# the current householder reflections actually reflected their columns +Hs = HH(:,kb1:m); +if(is_vector(Hs)) + nzidx = find(Hs' != 0); +else + nzidx = find(max(abs(Hs')) != 0); +endif +nzidx = nzidx(1:kb); +QQ = zeros(n,kb); +QQ(nzidx,1:kb) = eye(kb); + +# back out QQ matrix +for i=m:-1:1 + hv = HH(:,i); + av = alpha(i); + QQ = QQ-av*hv*(hv'*QQ); +end; + +endfunction diff --git a/scripts/linear-algebra/krylov.m b/scripts/linear-algebra/krylov.m new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/krylov.m @@ -0,0 +1,130 @@ +# Copyright (C) 1996,1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 2, or (at your option) any +# later version. +# +# Octave is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License +# along with Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +function [U,H,k1] = krylov(A,v,k,eps1); + # function [U,H,k1] = krylov(A,v,k{,eps1}); + # construct orthogonal basis U of Krylov subspace; + # span([v Av A^2*v ... A^(k+1)*v]); + # via householder reflections; reflections are multiplied to obtain U + # eps1: threshhold for 0 relative to norm of current column (default: 1e-12) + # method used: householder reflections to guard against loss of + # orthogonality + # + # outputs: + # returned basis U is n x k+1; A*U(:,1:k) = U*H + # k1: dimension of span of krylov subspace (based on eps1) + # if k > m-1, krylov returns the Hessenberg decompostion of A. + + # Written by A. S. Hodel 1992 + # $Revision: 1.2 $ + # $Log$ + + if(nargin == 3) + eps1 = 1e-12; + endif + + if( !is_square(A) ) + error("first argument must be a square matrix") + else + [m,n] = size(v); + if(m != is_square(A)) + error("krylov: argument dimensions do not match") + elseif( !is_sample(k) ) + error("krylov: third argument must be a scalar integer") + elseif( k > m) + warning(["krylov: k is too large; reducing to ",num2str(m)]); + k = m-1; + endif + endif + + if(norm(v) == 0) + U = []; + H = []; + k1 = 0; + return + endif + + k1 = k+1; # Assume subspace basis has full rank + [m,n] = size(A); + [hv,alpha(1),z] = housh(v,1,0); + + # initial orthogonal vector + q = zeros(n,1); + q(1) = 1; + q = q - alpha*hv*hv'*q; + normq = norm(q); + normres = normq; + + U(:,1) = hv; + j = 1; + while(j <= k & normres > eps1*normq) + # multiply to get new vector; + q = A*q; + normq = norm(q); + + # multiply by householder reflections to obtain projected vector and the + # next column of H + for i=1:j + hv = U(i:n,i); + av = alpha(i); + q(i:n,1) = q(i:n,1)-av*hv*(hv'*q(i:n,1)); + endfor + + i = j+1; + # compute and apply next householder vector; + if(i <= n) + [hv,av,z] = housh(q(i:n,1),1,0); + alpha(i) = av; + q(i:n,1) = q(i:n,1)-av*hv*(hv'*q(i:n,1)); + U(i:n,i) = hv; + H(1:i,j) = q(1:i); + else + av = 0; + H(:,j) = q; # complete Hessenberg decomposition + endif + + # see if done yet + normres = norm(q(i:n)); + if(normres > eps1*normq) + j = j+1; + else + k1 = min(k1,j); # time to quit; norm of residual is small + endif + + # back out new q vector for next pass; + j1 = columns(U); + q = zeros(n,1); + q(j1) = 1; + for i=j1:-1:1; + hv = U(i:n,i); + av = alpha(i); + q(i:n,1) = q(i:n,1)-av*hv*(hv'*q(i:n,1)); + endfor + endwhile + + # back out U matrix ; + j1 = columns(U); + for i=j1:-1:1; + hv = U(i:n,i); + av = alpha(i); + U(:,i) = zeros(n,1); + U(i,i) = 1; + U(i:n,i:j1) = U(i:n,i:j1)-av*hv*(hv'*U(i:n,i:j1)); + endfor + +endfunction diff --git a/scripts/linear-algebra/krylovb.m b/scripts/linear-algebra/krylovb.m new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/krylovb.m @@ -0,0 +1,93 @@ +# Copyright (C) 1993, 1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 2, or (at your option) any +# later version. +# +# Octave is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License +# along with Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +function [Uret,Ucols] = krylovb(A,V,k,eps1); + # function [U,Ucols] = krylovb(A,V,k[,eps1]); + # construct orthogonal basis U of block Krylov subspace; + # [V AV A^2*V ... A^(k+1)*V]; + # method used: householder reflections to guard against loss of + # orthogonality + # eps1: threshhold for 0 (default: 1e-12) + # + # outputs: + # returned basis U is orthogonal matrix; due to "zeroed" + # columns of product, may not satisfy A U = U H identity + # Ucols: dimension of span of krylov subspace (based on eps1) + # if k > m-1, krylov returns the Hessenberg decompostion of A. +# +# A. S. Hodel Feb 1993 +# $Revision$ +# $Log$ + + if(nargin == 3) + eps1 = 1e-12; + endif + na = is_square(A); + if( !na ) + error("krylov: first argument must be a square matrix") + endif + + [m,kb] = size(V); + if(m != na); + error("krylov: argument dimensions do not match") + endif + + if( !is_scalar(k) ) + error("krylov: third argument must be a scalar integer") + endif + + Vnrm = norm(V,Inf); + + # compute factored QR + [U,alpha,kb] = qrhouse(V,eps1*Vnrm); + Q = krygetq(U,alpha,kb); + Uret = Q; + Ucols = kb; + + iter = 0; + while (Ucols < na) & (kb > 0) & (iter < k) + iter++; + # multiply to get new vector; + V = A*Q; + + # project off of previous vectors + nzv = []; # set of reflection indices used so far + nj = length(alpha); + for i=1:nj + hv = U(:,i); + nzidx = find(hv != 0); # extract nonzero entries + nzv = union(nzv,nzidx(1)); # save for call to qrhouse + hv = hv(nzidx); + vr = V(nzidx,:); + av = alpha(i); + V(nzidx,:) = vr - (av*hv)*(hv'*vr); + end + V(nzv,:) = 0; # clear entries covered by earlier reflections + + [hvk,alphk,kb] = qrhouse(V,eps1*Vnrm); # now get new QR factorization + if(kb > 0) + U = [U,hvk]; # append new data to Householder data structure + alpha = [alpha;alphk]; + Q2 = krygetq(U,alpha,kb); + Uret = [Uret,Q2]; + + Ucols = Ucols + kb; + Q = Q2; + end + end +endfunction diff --git a/scripts/linear-algebra/qrhouse.m b/scripts/linear-algebra/qrhouse.m new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/qrhouse.m @@ -0,0 +1,87 @@ +# Copyright (C) 1992, 1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 2, or (at your option) any +# later version. +# +# Octave is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License +# along with Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +function [hv,alph,kb] = qrhouse(VV,eps1) +# function [hv,alph,kb] = qrhouse(VV,eps1) +# construct orthogonal basis of span(VV) with Householder vectors +# Q R = VV; Q may be obtained via routine krygetq; R is upper triangular +# if all rows of VV are nonzero; otherwise it's a permuted uppert +# triangular matrix with zero rows matching those of VV. +# Inputs: +# VV: matrix +# eps1: zero threshhold; default: 0. Used to check if a column +# of reduced R is already upper triangular; entries with magnitude +# <= eps1 are considered to be 0. +# Outputs: +# hv: matrix of householder reflection vectors as returned by housh +# alpha: vector of householder reflection values as returned by housh +# kb: computed rank of matrix +# qrhouse is used in krylovb for block Arnoldi iteration +# +# Reference: Golub and Van Loan, MATRIX COMPUTATIONS, 2nd ed. + +# Written by A. S. Hodel, 1992 +# $Revision: 1.2 $ +# $Log: qrhouse.m,v $ +# Revision 1.2 1998/08/26 21:08:29 hodelas +# Fixed bug in controllability analysis; code isolates zero rows of +# input matrices; correctly checks zero threshhold +# + +if(nargin < 2) + usage("[hv,alph,kb] = qrhouse(VV,eps1)"); +endif + +# extract only those rows of VV that are nonzero +if(is_vector(VV)) nzidx = find(abs(VV') > 0); +else nzidx = find(max(abs(VV')) > 0); +endif + +VVlen = rows(VV); # remember original size + +if(is_vector(VV)) VV = VV(nzidx); +else VV = VV(nzidx,:); +endif + +[Vr,Vc] = size(VV); nits = min(Vr,Vc); +kb = 0; kbnext = kb+1; +for ii = 1:nits; + kb = kbnext; # select column number of hv to build + hh = VV(:,ii); # extract next column of VV; ignore items 1:(ii-1). + idx = kb-1; # for zeroing + if(max(abs(hh(kb:Vr))) > eps1) + [hv(kb:Vr,kb),alph(kb),z] = housh(hh(kb:Vr),1,0); + if(kb>1) + hv(1:idx,kb) = 0; # zero top of hv for safety + endif + VV = VV - alph(kb)*hv(:,kb)*(hv(:,kb)'*VV); + kbnext = kb+1; + else + kb = kb-1; + end +endfor +if(kb <=0) + hv = []; + alph = []; +else + hvs = hv(:,1:kb); # remove extraneous vectors, expand to original size + hv = zeros(VVlen,kb); + hv(nzidx,:) = hvs; + alph = alph(1:kb); +end +endfunction