Mercurial > hg > octave-avbm

changeset 3211:440b2b28e74a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1998-11-05 04:16:22 by jwe]
author jwe
date Thu, 05 Nov 1998 04:16:34 +0000
parents 258d7b26c719
children bf61c443a366
files scripts/ChangeLog scripts/linear-algebra/housh.m scripts/linear-algebra/krygetq.m scripts/linear-algebra/krylov.m scripts/linear-algebra/krylovb.m scripts/linear-algebra/qrhouse.m
diffstat 6 files changed, 446 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,15 @@
+Wed Nov  4 21:51:13 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* 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  <jwe@bevo.che.wisc.edu>
 
 	* general/is_square.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
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
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
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
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