#!/bin/sh # This is a shell archive (produced by GNU sharutils 4.2). # To extract the files from this archive, save it to some FILE, remove # everything before the `!/bin/sh' line above, then type `sh FILE'. # # Made on 2004-09-23 16:15 CDT by . # Source directory was `/www/caam551/public_html/MatlabCode'. # # Existing files will *not* be overwritten unless `-c' is specified. # # This shar contains: # length mode name # ------ ---------- ------------------------------------------ # 1349 -rw-rw-r-- Arnold.m # 785 -rw-rw-r-- Arnoldi.m # 1000 -rw-rw-r-- ArnoldiC.m # 1450 -rw-r--r-- CGStest.m # 1882 -rw-r--r-- CGStestN.m # 950 -rw-r--r-- Chol.m # 2040 -rw-r--r-- DemoAngle.m # 1107 -rw-r--r-- EQNcondemo.m # 430 -rw-r--r-- Givens.m # 1984 -rw-r--r-- Gmres.m # 944 -rw-r--r-- Gmres0.m # 1997 -rw-rw-r-- GmresIR.m # 2175 -rw-rw-r-- GmresIR_plt.m # 2176 -rw-r--r-- GmresP.m # 727 -rw-r--r-- HHbad.m # 733 -rw-r--r-- HHgood.m # 6773 -rw-r--r-- Hump.m # 1899 -rw-rw-r-- Iram.m # 1833 -rw-rw-r-- Iram1.m # 3062 -rw-rw-r-- Iram_Cnv.m # 4352 -rw-rw-r-- Iram_Plt.m # 4345 -rw-rw-r-- Iram_Plt1.m # 2013 -rw-r--r-- LScondemo.m # 1094 -rw-r--r-- LUfac.m # 1098 -rw-r--r-- LUfacC.m # 1819 -rw-r--r-- LUfacP.m # 1068 -rw-r--r-- LUfacR.m # 2481 -rw-rw-r-- Lanconv_gui.m # 4049 -rw-rw-r-- Lanconv_update.m # 942 -rw-rw-r-- Lanczos.m # 782 -rw-r--r-- QRcgs.m # 1655 -rw-r--r-- QRcgsF.m # 1220 -rw-r--r-- QRfac.m # 1277 -rw-r--r-- QRfac2.m # 1583 -rw-r--r-- QRgivens.m # 858 -rw-r--r-- QRmgs.m # 1755 -rw-rw-r-- QRpiv.m # 2083 -rw-rw-r-- QRpivD.m # 2084 -rw-r--r-- SparseQRex.m # 2775 -rw-rw-r-- Sym_ritz.m # 3130 -rw-r--r-- Threshtest.m # 4319 -rw-r--r-- UNsparsedem.m # 417 -rw-rw-r-- bluemap.m # 1970 -rw-rw-r-- compLinSlvrs.m # 371 -rw-rw-r-- driveA.m # 289 -rw-rw-r-- driveA2.m # 290 -rw-rw-r-- driveAC.m # 466 -rw-rw-r-- driveGM.m # 5108 -rw-r--r-- driveGM0.m # 1312 -rw-r--r-- driveGMLU.m # 1216 -rw-r--r-- driveHHB.m # 1370 -rw-r--r-- driveHHG.m # 1221 -rw-rw-r-- driveIRA.m # 1156 -rw-rw-r-- driveIRA_Cnv.m # 1153 -rw-rw-r-- driveIRA_Plt.m # 2293 -rw-rw-r-- driveIRA_SI.m # 411 -rw-rw-r-- driveL.m # 315 -rw-rw-r-- driveL2.m # 189 -rw-rw-r-- driveP.m # 175 -rw-rw-r-- driveQR.m # 120 -rw-rw-r-- driveQRR.m # 227 -rw-rw-r-- driveQRRS.m # 175 -rw-rw-r-- driveQRS.m # 600 -rw-rw-r-- eigconv.m # 2954 -rw-rw-r-- eigen_update.m # 5004 -rw-rw-r-- eigenui.m # 619 -rw-r--r-- evalB.m # 768 -rw-rw-r-- krylov.m # 7224 -rw-rw-rw- mmread.m # 170 -rw-rw-r-- orthdemo.m # 2135 -rw-rw-r-- polyR.m # 331 -rw-rw-r-- pow.m # 338 -rw-rw-r-- pow_inf.m # 295 -rw-rw-r-- pow_rq.m # 2056 -rw-r--r-- powersR.m # 2704 -rw-rw-r-- qriter.m # 4331 -rw-rw-r-- qriterR.m # 4350 -rw-rw-r-- qriterRS.m # 2898 -rw-rw-r-- qriterS.m # 1148 -rw-rw-r-- qrstep.m # 64 -rw-rw-r-- readsher.m # 704 -rw-rw-r-- ritzest.m # 1169 -rw-rw-r-- select_shifts.m # 3623 -rw-r--r-- sparsdem.m # save_IFS="${IFS}" IFS="${IFS}:" gettext_dir=FAILED locale_dir=FAILED first_param="$1" for dir in $PATH do if test "$gettext_dir" = FAILED && test -f $dir/gettext \ && ($dir/gettext --version >/dev/null 2>&1) then set `$dir/gettext --version 2>&1` if test "$3" = GNU then gettext_dir=$dir fi fi if test "$locale_dir" = FAILED && test -f $dir/shar \ && ($dir/shar --print-text-domain-dir >/dev/null 2>&1) then locale_dir=`$dir/shar --print-text-domain-dir` fi done IFS="$save_IFS" if test "$locale_dir" = FAILED || test "$gettext_dir" = FAILED then echo=echo else TEXTDOMAINDIR=$locale_dir export TEXTDOMAINDIR TEXTDOMAIN=sharutils export TEXTDOMAIN echo="$gettext_dir/gettext -s" fi touch -am 1231235999 $$.touch >/dev/null 2>&1 if test ! -f 1231235999 && test -f $$.touch; then shar_touch=touch else shar_touch=: echo $echo 'WARNING: not restoring timestamps. Consider getting and' $echo "installing GNU \`touch', distributed in GNU File Utilities..." echo fi rm -f 1231235999 $$.touch # if mkdir _sh24386; then $echo 'x -' 'creating lock directory' else $echo 'failed to create lock directory' exit 1 fi # ============= Arnold.m ============== if test -f 'Arnold.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Arnold.m' '(file already exists)' else $echo 'x -' extracting 'Arnold.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Arnold.m' && function [V,H,f] = Arnold(A,V,H,f,k,m); % % Input: A -- an n by n matrix % V -- an n by k orthogonal matrix % H -- a k by k upper Hessenberg matrix % f -- an nonzero n vector % % with AV = VH + fe_k' (if k > 1) % % k -- a positive integer (k << n assumed) % % m -- a positive integer (k < m << n assumed) % % % Output: V -- an n by m orthogonal matrix % H -- a m by m upper Hessenberg matrix % f -- an n vector % % % with AV = VH + fe_k' % % Leading k columns of V agree with input V % % assures V'V = I_m with DGKS correction % % % D.C. Sorensen % 2 March 2000 % X n = length(f); X if (k == 1), X X beta = norm(f); X if ( beta == 0) , f_is_zero_error = beta, return; end X X v = f/beta; X w = A*v; X alpha = v'*w; X X X f = w - v*alpha; X c = v'*f; X f = f - v*c; X alpha = alpha + c; X X V(:,1) = v; H(1,1) = alpha; X end X X for j = k+1:m, X X beta = norm(f); X v = f/beta; X X H(j,j-1) = beta; X V(:,j) = v; X X w = A*v; X h = V(:,1:j)'*w; X f = w - V(:,1:j)*h; X c = V(:,1:j)'*f; X f = f - V(:,1:j)*c; X h = h + c; X X H(1:j,j) = h; X X end X SHAR_EOF $shar_touch -am 08270935104 'Arnold.m' && chmod 0664 'Arnold.m' || $echo 'restore of' 'Arnold.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Arnold.m:' 'MD5 check failed' 3f68071c0c4a5d833d4d9a986eba1b66 Arnold.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Arnold.m'`" test 1349 -eq "$shar_count" || $echo 'Arnold.m:' 'original size' '1349,' 'current size' "$shar_count!" fi fi # ============= Arnoldi.m ============== if test -f 'Arnoldi.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Arnoldi.m' '(file already exists)' else $echo 'x -' extracting 'Arnoldi.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Arnoldi.m' && function [V,H,f] = Arnoldi(A,k,v); % % Input: A -- an n by n matrix % k -- a positive integer (k << n assumed) % v -- an n vector (v .ne. 0 assumed) % % Output: V -- an n by k orthogonal matrix % H -- a k by k upper Hessenberg matrix % f -- an n vector % % % with AV = VH + fe_k' % % % D.C. Sorensen % 10 Feb 00 % X n = length(v); X H = zeros(k); X V = zeros(n,k); X X v = v/norm(v); X w = A*v; X alpha = v'*w; X X X V(:,1) = v; H(1,1) = alpha; X f = w - v*alpha; X X for j = 2:k, X X beta = norm(f); X v = f/beta; X X H(j,j-1) = beta; X V(:,j) = v; X X w = A*v; X h = V(:,1:j)'*w; X f = w - V(:,1:j)*h; X X H(1:j,j) = h; X X end X X SHAR_EOF $shar_touch -am 02241523100 'Arnoldi.m' && chmod 0664 'Arnoldi.m' || $echo 'restore of' 'Arnoldi.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Arnoldi.m:' 'MD5 check failed' f2c2a57a1ecf8c1bf78214dcd9f25fc1 Arnoldi.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Arnoldi.m'`" test 785 -eq "$shar_count" || $echo 'Arnoldi.m:' 'original size' '785,' 'current size' "$shar_count!" fi fi # ============= ArnoldiC.m ============== if test -f 'ArnoldiC.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'ArnoldiC.m' '(file already exists)' else $echo 'x -' extracting 'ArnoldiC.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'ArnoldiC.m' && function [V,H,f] = ArnoldiC(A,k,v); % % Input: A -- an n by n matrix % k -- a positive integer (k << n assumed) % v -- an n vector (v .ne. 0 assumed) % % Output: V -- an n by k orthogonal matrix % H -- a k by k upper Hessenberg matrix % f -- an n vector % % % with AV = VH + fe_k' % % assures V'V = I_k with DGKS correction % % % D.C. Sorensen % 10 Feb 00 % X n = length(v); X H = zeros(k); X V = zeros(n,k); X X v = v/norm(v); X w = A*v; X alpha = v'*w; X X X f = w - v*alpha; X c = v'*f; X f = f - v*c; X alpha = alpha + c; X X V(:,1) = v; H(1,1) = alpha; X X X for j = 2:k, X X beta = norm(f); X v = f/beta; X X H(j,j-1) = beta; X V(:,j) = v; X X w = A*v; X h = V(:,1:j)'*w; X f = w - V(:,1:j)*h; X c = V(:,1:j)'*f; X f = f - V(:,1:j)*c; X h = h + c; X X H(1:j,j) = h; X X end X X SHAR_EOF $shar_touch -am 02241526100 'ArnoldiC.m' && chmod 0664 'ArnoldiC.m' || $echo 'restore of' 'ArnoldiC.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'ArnoldiC.m:' 'MD5 check failed' f37ff20d0182c697a34f81a2984876d7 ArnoldiC.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'ArnoldiC.m'`" test 1000 -eq "$shar_count" || $echo 'ArnoldiC.m:' 'original size' '1000,' 'current size' "$shar_count!" fi fi # ============= CGStest.m ============== if test -f 'CGStest.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'CGStest.m' '(file already exists)' else $echo 'x -' extracting 'CGStest.m' '(binary)' sed 's/^X//' << 'SHAR_EOF' | uudecode && begin 600 CGStest.m M#0HE#0HE("!4:&ES('-CF%T:6]N+@T*)0T*)2`@ M("!#1U,@+2!#;&%S`T* M#0H@("!;52Q272`]('%R*')A;F1N*#$U,"PX,"DL,"D[#0H@("!;5BQ272`] M('%R*')A;F1N*#@P*2PP*3L-"B`@(`T*("`@"!!#0HE("!W:71H('-I;F=U;&%R('9A;'5EF4H02D[#0HE#0HE M("!$:7-P;&%Y(&-O;F1I=&EO;B!N=6UB97(@;V8@00T*)0T*("`@2V]N9"`] M(',H,2DOC(@/2!E<',J;VYE2AR,2PG2AR M,RPG9RHG*3L-"B`@('-E;6EL;V=Y*'HQ+"C(L)RTG*3L-"B`@('-E;6EL;V=Y*',L)VLG*0T*("`@;&5G96YD*"=#1U,G M+"=-1U,G+"=#1U-F)RPG4W%R="A%<',I)RP@)T5P&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'CGStest.m:' 'MD5 check failed' a97b5d49bf78aae7b8004428da28990d CGStest.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'CGStest.m'`" test 1450 -eq "$shar_count" || $echo 'CGStest.m:' 'original size' '1450,' 'current size' "$shar_count!" fi fi # ============= CGStestN.m ============== if test -f 'CGStestN.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'CGStestN.m' '(file already exists)' else $echo 'x -' extracting 'CGStestN.m' '(binary)' sed 's/^X//' << 'SHAR_EOF' | uudecode && begin 600 CGStestN.m M)0T*)2`@5&AI2!&:7@-"@T* M)0T*)2`@($-O;7!U=&4@=&AE('1HC(@ M/2!E<',J;VYE2AR,2PG2AZ,2PG+2XG*3L-"B`@('-E;6EL M;V=Y*'HR+"2AS+"=K)RD-"B`@(&QE9V5N9"@G M0T=3)RPG34=3)RPG0T=31B&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'CGStestN.m:' 'MD5 check failed' 9088bb059fc33af213f8fcf1d78435af CGStestN.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'CGStestN.m'`" test 1882 -eq "$shar_count" || $echo 'CGStestN.m:' 'original size' '1882,' 'current size' "$shar_count!" fi fi # ============= Chol.m ============== if test -f 'Chol.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Chol.m' '(file already exists)' else $echo 'x -' extracting 'Chol.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Chol.m' && X function [L,ierr] = Chol(A); X % X % X % Input: A an symmetric n by n matrix. X % X % Output: L a lower triangular matrix with A = LL' X % X % ierr an error flag. ierr = k >0 means failure at step k X % = 0 means normal exit X % X % X % X % D.C. Sorensen X % 25 Aug 03 X %----------------------------------------------------------------- X % X [n,n] = size(A); X ierr = 0; X % X for k = 1:n, X % X % exit if A is not positive definite X % X if (A(k,k) <= 0), ierr = k; return; end X % X % Compute main diagonal elt. and then scale the k-th column X % X X A(k,k) = sqrt(A(k,k)); X A(k+1:n,k) = A(k+1:n,k)/A(k,k); X % X % Update lower triangle of the trailing (n-k) by (n-k) block X % X X for j = k+1:n, X X A(j:n,j) = A(j:n,j) - A(j:n,k)*A(j,k); X end X end X L = tril(A); SHAR_EOF $shar_touch -am 08270936104 'Chol.m' && chmod 0644 'Chol.m' || $echo 'restore of' 'Chol.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Chol.m:' 'MD5 check failed' 3a25885dcf6b71666c0d743114b4380e Chol.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Chol.m'`" test 950 -eq "$shar_count" || $echo 'Chol.m:' 'original size' '950,' 'current size' "$shar_count!" fi fi # ============= DemoAngle.m ============== if test -f 'DemoAngle.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'DemoAngle.m' '(file already exists)' else $echo 'x -' extracting 'DemoAngle.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'DemoAngle.m' && % % This script demonstrates the increasing relative % projections in the Least Sqares problem for a % sequence of RHS b_j of increasing norm that % all have the same projections onto Range(A) % % To get a good aspect ratio on the plot % stretch the window by hand along the z - axis % after the first plot appears % % D.C. Sorensen % 4 Oct 01 %---------------------------------------------------- % X X for tau1 = [10 50 100 150], X X m = 200; X t = linspace(0,(6.5*pi/4),m); X X x = cos(t); X y = sin(t); X z = t/10; X X x1 = [0 x(m)]; X x2 = [0 y(m)]; X x3 = [0 z(m)]; X X z1 = [x(m) x(m)]; X z2 = [y(m) y(m)]; X z3 = 0*x3; X X plot3(x1,x2,z3,'b') X hold % % Plot four vectors of increasing length % that all have the same projection % onto the Range(A) % X X for tau = [10 50 100 150], X plot3(x1,x2,tau*x3,'k') X plot3(x1(2),x2(2),tau*z(m),'k+') X end X X plot3(z1,z2,tau*x3,'m') X X t = linspace(0,2*pi,m); X x0 = cos(t); X y0 = sin(t); X w = 0*t; X mu = 2; X x = x0*mu; X y = y0*mu; X plot3(x,y,w,'c'); X x = x0*mu; X X tau = tau1 % % Construct a perturbation to the RHS b % % of size mu*norm(b) % % where mu = .005 % X X b = [x1(2) x2(2) tau*z(m)]; X bnorm = norm(b); X mu = .005*bnorm; % % Plot the projection of the circle of such perturbations % in the Range(A) - in this case the x,y plane % X x = x1(2) + x0*mu; X y = x2(2) + y0*mu; X plot3(x,y,w,'r'); % % Plot the circle of such perturbations about b % parallel to the x,y plane % X w = tau*z(m)*ones(m,1); X plot3(x,y,w,'k'); X X grid X X hold X if (tau1 < 20), X disp(' To get a good aspect ratio on the plot ') X disp(' stretch the window by mouse along the z - axis ') X disp(' after the first plot appears ') X end X disp(' ') X disp('strike any key to continue') X pause X end X SHAR_EOF $shar_touch -am 09240814103 'DemoAngle.m' && chmod 0644 'DemoAngle.m' || $echo 'restore of' 'DemoAngle.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'DemoAngle.m:' 'MD5 check failed' 934d5d9b1144353f00cdc8746e6a1e48 DemoAngle.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'DemoAngle.m'`" test 2040 -eq "$shar_count" || $echo 'DemoAngle.m:' 'original size' '2040,' 'current size' "$shar_count!" fi fi # ============= EQNcondemo.m ============== if test -f 'EQNcondemo.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'EQNcondemo.m' '(file already exists)' else $echo 'x -' extracting 'EQNcondemo.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'EQNcondemo.m' && % % This matlab script % is designed to demonstrate the validity of the bound. % % norm(x1 -x)/norm(x) .le. cond(A)*norm(b1 - b)/norm(b) % % where Ax = b , Ax1 = b1. % % Also demonstrates small residual does not imply accurate answer. % You MUST know condition number. % % D.C. Sorensen 9/27/99 % % modified 10 Oct 01 (DCS) X X format short e X n = 100; X m = 90; X k = 13; X Xout = zeros(k,5); X A = randn(n); X [U,S,V] = svd(A); % % prepare rhs. b % X b = U(:,1:m)*randn(m,1); % % X z = randn(n,1); X tau = 100*eps; X x = A\b; X kk = cond(A); X s = diag(S); X b1 = b + tau*z; X for j = 1:k, % % increase condition by factor of 10 % X s(m+1:100) = s(m+1:100)/10; X S = diag(s); % % form new matrix and solve for x1 and x % X A1 = U*S*V'; X x1 = A1\b1; X x = A1\b; X resid = norm(b - A1*x1); X condA1 = cond(A1); X xerr = norm(x1 - x)/norm(x); X berr = norm(b1 - b)/norm(b); X Xout(j,:) = [ xerr condA1*berr round(-log10(xerr)) condA1 resid]; X end X X disp(' ') X disp(' x-err estimate digits cond(A1) resid ') X disp(Xout) X X X X SHAR_EOF $shar_touch -am 09240813103 'EQNcondemo.m' && chmod 0644 'EQNcondemo.m' || $echo 'restore of' 'EQNcondemo.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'EQNcondemo.m:' 'MD5 check failed' d66cd4c689bba8e25e732bd5f7026242 EQNcondemo.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'EQNcondemo.m'`" test 1107 -eq "$shar_count" || $echo 'EQNcondemo.m:' 'original size' '1107,' 'current size' "$shar_count!" fi fi # ============= Givens.m ============== if test -f 'Givens.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Givens.m' '(file already exists)' else $echo 'x -' extracting 'Givens.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Givens.m' && function [c,s,r] = Givens(x,y); % % computes c,s,r such that % % [c s] [x] = [r] % [-s c] [y] [0] % % with c*c + s*s = 1; % X if (y == 0), X c = 1; s = 0; r = x; X else X if (abs(x) >= abs(y)) X t = y/x; r = sqrt(1 + t*t); X c = 1/r; X s = t*c; X r = x*r; X else X t = x/y; r = sqrt(1 + t*t); X s = 1/r; X c = t*s; X r = y*r; X end X end X SHAR_EOF $shar_touch -am 09161400103 'Givens.m' && chmod 0644 'Givens.m' || $echo 'restore of' 'Givens.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Givens.m:' 'MD5 check failed' 0a6779e6474e5841998e1b3bc7242ff7 Givens.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Givens.m'`" test 430 -eq "$shar_count" || $echo 'Givens.m:' 'original size' '430,' 'current size' "$shar_count!" fi fi # ============= Gmres.m ============== if test -f 'Gmres.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Gmres.m' '(file already exists)' else $echo 'x -' extracting 'Gmres.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Gmres.m' && function [x,res] = Gmres(A,b,kmax,tol); % % Usage [x,rho] = Gmres(A,b,tol); % % Solves Ax = b % % Input: A -- an n by n matrix % % % kmax -- max no. GMRES steps % % tol -- stopping tolerance % % % Output: x -- an n vector... approx solution Ax = b % % res -- Residual history res(j) = (norm(b -Ax_j) % % % D.C. Sorensen % 13 March 2001 % X k = kmax; X q = eye(k+1,1); % % The rhs b is scaled to unit length % size will be restored after convergence % X theta = norm(b); X v = b/theta; X X % X n = length(v); X R = zeros(k); X V = zeros(n,k); X X v = v/norm(v); X w = A*v; X alpha = v'*w; X X f = w - v*alpha; X c = v'*f; X f = f - v*c; X alpha = alpha + c; X X R(1,1) = alpha; X V(:,1) = v; X rho = 1; X X Q = eye(k+1); X res = []; X X X X j = 1; X while (rho > tol & j < k), X X j = j+1; X beta = norm(f); X v = f/beta; % % Compute the Givens transformation G to zero % last row of updated H (H is not stored) % to get H = QR % X [G,t] = qr([R(j-1,j-1) ; beta]); X Q(:,j-1:j) = Q(:,j-1:j)*G; X % % Update rhs q = Q'*e_1 % % Q should not be stored, it should be % represented by the Givens transformations % X q(j-1:j) = G'*q(j-1:j); X R(j-1,j-1) = t(1); X rho = abs(q(j)); X X res = [res; rho]; X X V(:,j) = v; X X if (rho > tol), X w = A*v; X h = V(:,1:j)'*w; X f = w - V(:,1:j)*h; X c = V(:,1:j)'*f; X f = f - V(:,1:j)*c; X h = h + c; X X R(1:j,j) = Q(1:j,1:j)'*h; X end X X end % % After convergence (or max iter exceeded) % compute the solution x. Must multiply by % theta to restore rhs b to original length % X % res = res*theta; X res = res; % relative residual X X y = R(1:j-1,1:j-1)\q(1:j-1); X x = V(:,1:j-1)*y*theta; X X SHAR_EOF $shar_touch -am 09271136102 'Gmres.m' && chmod 0644 'Gmres.m' || $echo 'restore of' 'Gmres.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Gmres.m:' 'MD5 check failed' e022328a15ef587b182935aa113bb8a6 Gmres.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Gmres.m'`" test 1984 -eq "$shar_count" || $echo 'Gmres.m:' 'original size' '1984,' 'current size' "$shar_count!" fi fi # ============= Gmres0.m ============== if test -f 'Gmres0.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Gmres0.m' '(file already exists)' else $echo 'x -' extracting 'Gmres0.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Gmres0.m' && X function [x,rho,ritz,hritz] = Gmres0(A,b,k); % % Usage [x,rho] = GmresIR(A,b,k); % % Solves Ax = b % % Input: A -- an n by n matrix % % % k -- a positive integer % % % Output: x -- an n vector... approx solution Ax = b % % rho -- norm(b - Ax) % % ritz -- Ritz values = roots of FOM polynomial % % hritz -- Harmonic Ritz values = roots of GMres polynomial % % % D.C. Sorensen % 13 March 2001 % X ek = zeros(k,1); ek(k) = 1; X theta = norm(b); X v = b/theta; X n = length(v); X X [V,H,f] = ArnoldiC(A,k,v); X X Rnorm = norm(A*V - V*H - f*ek') X X beta = norm(f); X if (k < n), X [Q,R] = qr([H ; beta*ek']); X rho = abs(Q(1,k+1)*theta); X else X [Q,R] = qr(H); X rho = 0; X end X q = Q(1,1:k)'; X y = R(1:k,1:k)\q; X X x = V*y*theta; X X g = ((H')\ek)*(beta^2); X hritz = eig( H + g*ek'); X ritz = eig(H); X SHAR_EOF $shar_touch -am 09271127102 'Gmres0.m' && chmod 0644 'Gmres0.m' || $echo 'restore of' 'Gmres0.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Gmres0.m:' 'MD5 check failed' 1bf99580d95a62969cb44021eb698ab7 Gmres0.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Gmres0.m'`" test 944 -eq "$shar_count" || $echo 'Gmres0.m:' 'original size' '944,' 'current size' "$shar_count!" fi fi # ============= GmresIR.m ============== if test -f 'GmresIR.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'GmresIR.m' '(file already exists)' else $echo 'x -' extracting 'GmresIR.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'GmresIR.m' && function [x,rho,Resid,nitr] = GmresIR(A,b,k,m); % % Usage [x,rho,Resid,nitr] = GmresIR(A,b,k,m); % % Solves Ax = b % % Input: A -- an n by n matrix % % % k -- a positive integer % % m -- a positive integer k < m << n % % Recommended: m > 2k % % b -- % % Output: x -- an n vector... approx solution Ax = b % % rho -- norm(b - Ax) % % % D.C. Sorensen % 2 March 2000 % X [n,n] = size(A); X Maxitr = 200; X Resid = zeros(Maxitr,1); X X V = zeros(n,m); H = zeros(m,m); X which = 'SM'; X tol = sqrt(eps); ritz = 1; nitr = 0; X rho = norm(b); X e1 = eye(m,1); X em = zeros(m,1); em(m) = 1; X X x = zeros(n,1); X ko = k; X nitr = 0; X X X f = b; X [V,H,f] = Arnold(A,V,H,f,1,m); % check this X e1 = eye(m+1,1); X beta = norm(f); X y = [H ; beta*em']\e1*rho; X X x = V*y; % r = b - A*x; X r = V(:,1)*rho - V*(H*y) - f*y(m); X rho = norm(r); X X while (rho > tol & nitr < Maxitr), X X nitr = nitr + 1; X Resid(nitr) = rho; % g = (H')\em*beta*beta; % [mu, ritz]= select_shifts((H + g*em'),which); X X G = ([H ; beta*em']'*[H ; beta*em'])\H'; X [mu, ritz]= select_shifts(G,which); X X Q = eye(m); X j = m; X while(j > k) % % Apply unwanted eigenvalues as shifts with QR steps % X [Q,H] = qrstep(Q,H,mu(j),1,m); % % Decrease j by 2 if m_j imaginary and by 1 if real % X if (abs(imag(mu(j))) > 0), X j = j-2; X else X j = j-1; X end X end X ko = j; X X f = V*Q(:,ko+1)*H(ko+1,ko) + f*Q(m,ko); X V(:,1:ko) = V*Q(:,1:ko); X X [V,H,f] = Arnold(A,V,H,f,ko,m); X X t = V'*r; X X beta = norm(f); X t = [t ; f'*r/beta]; X y = [H ; beta*em']\t; X X x = x + V*y; % r = r - A*V*y; X r = r - V*(H*y) - f*y(m); X rho = norm(r); X end X X SHAR_EOF $shar_touch -am 04051446101 'GmresIR.m' && chmod 0664 'GmresIR.m' || $echo 'restore of' 'GmresIR.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'GmresIR.m:' 'MD5 check failed' edb659c0f4298ddd7a67dff868e077a0 GmresIR.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'GmresIR.m'`" test 1997 -eq "$shar_count" || $echo 'GmresIR.m:' 'original size' '1997,' 'current size' "$shar_count!" fi fi # ============= GmresIR_plt.m ============== if test -f 'GmresIR_plt.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'GmresIR_plt.m' '(file already exists)' else $echo 'x -' extracting 'GmresIR_plt.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'GmresIR_plt.m' && function [x,rho,Resid,nitr] = GmresIR_plt(A,b,k,m); % % Solves Ax = b % % Input: A -- an n by n matrix % % % k -- a positive integer % % m -- a positive integer k < m << n % % Recommended: m > 2k % % b -- % % Output: x -- an n vector... approx solution Ax = b % % rho -- norm(b - Ax) % % % D.C. Sorensen % 2 March 2000 % X [n,n] = size(A); X Maxitr = 200; X Resid = zeros(Maxitr,1); X X avals = eig(A); X X V = zeros(n,m); H = zeros(m,m); X which = 'SM'; X tol = sqrt(eps); ritz = 1; nitr = 0; X rho = norm(b); X e1 = eye(m,1); X em = zeros(m,1); em(m) = 1; X X x = zeros(n,1); X ko = k; X nitr = 0; X X X f = b; X [V,H,f] = Arnold(A,V,H,f,1,m); % check this X e1 = eye(m+1,1); X beta = norm(f); X y = [H ; beta*em']\e1*rho; X X x = V*y; % r = b - A*x; X r = V(:,1)*rho - V*(H*y) - f*y(m); X rho = norm(r); X X while (rho > tol & nitr < Maxitr), X X nitr = nitr + 1; X Resid(nitr) = rho; % g = (H')\em*beta*beta; % [mu, ritz]= select_shifts((H + g*em'),which); X X G = ([H ; beta*em']'*[H ; beta*em'])\H'; X G = inv(G); X [mu, ritz]= select_shifts(G,which); X X hvals = eig(H); X X plot(real(avals), imag(avals), 'k+') X hold X plot(real(mu), imag(mu), 'ro') X plot(real(hvals), imag(hvals), 'b*') X hold X pause X X Q = eye(m); X j = m; X while(j > k) % % Apply unwanted eigenvalues as shifts with QR steps % X [Q,H] = qrstep(Q,H,mu(j),1,m); % % Decrease j by 2 if m_j imaginary and by 1 if real % X if (abs(imag(mu(j))) > 0), X j = j-2; X else X j = j-1; X end X end X ko = j; X X f = V*Q(:,ko+1)*H(ko+1,ko) + f*Q(m,ko); X V(:,1:ko) = V*Q(:,1:ko); X X [V,H,f] = Arnold(A,V,H,f,ko,m); X X t = V'*r; X X beta = norm(f); X t = [t ; f'*r/beta]; X y = [H ; beta*em']\t; X X x = x + V*y; % r = r - A*V*y; X r = r - V*(H*y) - f*y(m); X rho = norm(r); X end X X SHAR_EOF $shar_touch -am 04051446101 'GmresIR_plt.m' && chmod 0664 'GmresIR_plt.m' || $echo 'restore of' 'GmresIR_plt.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'GmresIR_plt.m:' 'MD5 check failed' 8cab33ad055b2afd5d59fb075edd9b26 GmresIR_plt.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'GmresIR_plt.m'`" test 2175 -eq "$shar_count" || $echo 'GmresIR_plt.m:' 'original size' '2175,' 'current size' "$shar_count!" fi fi # ============= GmresP.m ============== if test -f 'GmresP.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'GmresP.m' '(file already exists)' else $echo 'x -' extracting 'GmresP.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'GmresP.m' && function [x,res] = GmresP(A,L,U,b,kmax,tol); % % Usage [x,rho] = GmresP(A,L,U,b,tol); % % Solves Ax = b % % Input: A -- an n by n matrix % % L,U -- The LU decompostion of a preconditioner % M = LU. Solves M\A x = M\b % % % kmax -- max no. GMRES steps % % tol -- stopping tolerance % % % Output: x -- an n vector... approx solution Ax = b % % res -- Residual history res(j) = (norm(b -Ax_j) % % % D.C. Sorensen % 13 March 2001 % X k = kmax; X q = eye(k+1,1); % % The rhs b is scaled to unit length % size will be restored after convergence X X v = U\(L\b); X theta = norm(v); X v = v/theta; X X % X n = length(v); X R = zeros(k); X V = zeros(n,k); X X v = v/norm(v); X w = U\(L\(A*v)); X X alpha = v'*w; X X f = w - v*alpha; X c = v'*f; X f = f - v*c; X alpha = alpha + c; X X R(1,1) = alpha; X V(:,1) = v; X rho = 1; X X Q = eye(k+1); X res = [norm(b)]; X X X X j = 1; X while (rho > tol & j < k), X X j = j+1; X beta = norm(f); X v = f/beta; % % Compute the Givens transformation G to zero % last row of updated H (H is not stored) % to get H = QR % X [G,t] = qr([R(j-1,j-1) ; beta]); X Q(:,j-1:j) = Q(:,j-1:j)*G; X % % Update rhs q = Q'*e_1 % % Q should not be stored, it should be % represented by the Givens transformations % X q(j-1:j) = G'*q(j-1:j); X R(j-1,j-1) = t(1); X rho = abs(q(j)); X % res = [res; rho]; X X V(:,j) = v; X X if (rho > tol), X w = U\(L\(A*v)); X h = V(:,1:j)'*w; X f = w - V(:,1:j)*h; X c = V(:,1:j)'*f; X f = f - V(:,1:j)*c; X h = h + c; X X R(1:j,j) = Q(1:j,1:j)'*h; X end X % % After convergence (or max iter exceeded) % compute the solution x. Must multiply by % theta to restore rhs b to original length % X % res = res*theta; % res = res; X y = R(1:j-1,1:j-1)\q(1:j-1); X x = V(:,1:j-1)*y*theta; X rho = norm(b - A*x); X res = [res; rho]; X end X X SHAR_EOF $shar_touch -am 10291507103 'GmresP.m' && chmod 0644 'GmresP.m' || $echo 'restore of' 'GmresP.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'GmresP.m:' 'MD5 check failed' 51c98a7bc0eaa47808c3301dda889e00 GmresP.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'GmresP.m'`" test 2176 -eq "$shar_count" || $echo 'GmresP.m:' 'original size' '2176,' 'current size' "$shar_count!" fi fi # ============= HHbad.m ============== if test -f 'HHbad.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'HHbad.m' '(file already exists)' else $echo 'x -' extracting 'HHbad.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'HHbad.m' && function [v,tau] = HHbad(a); % % Input: a - an m vector % % Output: v,tau - % The Householder Transformatin is I - tau*v*v' % % Demonstrates results of wrong sign choice in HH transformation % % D.C. Sorensen % 15 Sep 04 %----------------------------------------------------------------- X X m = length(a); X e1 = eye(m,1); X % % Compute the Householder vector v % % I - tau(k)*v*v' is the transformation % X v = a; % % Wrong Choice of Sign !!! % X rho = sign(v(1))*norm(v); X X if (abs(rho) > 0), X v(1) = v(1) - rho; X tau = -v(1)/rho; X v = v/v(1); X end X X w = a - v*(tau*v'*a); X X disp('Diff = norm(e1*rho - w)') X Diff = norm(e1*rho - w) SHAR_EOF $shar_touch -am 09160908104 'HHbad.m' && chmod 0644 'HHbad.m' || $echo 'restore of' 'HHbad.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'HHbad.m:' 'MD5 check failed' a8bbc08df4097baaf7a32bc4cc66ec89 HHbad.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'HHbad.m'`" test 727 -eq "$shar_count" || $echo 'HHbad.m:' 'original size' '727,' 'current size' "$shar_count!" fi fi # ============= HHgood.m ============== if test -f 'HHgood.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'HHgood.m' '(file already exists)' else $echo 'x -' extracting 'HHgood.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'HHgood.m' && function [v,tau] = HHgood(a); % % Input: a - an m vector % % Output: v,tau - % The Householder Transformatin is I - tau*v*v' % % Demonstrates results of correct sign choice in HH transformation % % D.C. Sorensen % 15 Sep 04 %----------------------------------------------------------------- X X m = length(a); X e1 = eye(m,1); X % % Compute the Householder vector v % % I - tau(k)*v*v' is the transformation % X v = a; % % Correct Choice of Sign !!! % X rho = -sign(v(1))*norm(v); X X if (abs(rho) > 0), X v(1) = v(1) - rho; X tau = -v(1)/rho; X v = v/v(1); X end X X w = a - v*(tau*v'*a); X X disp('Diff = norm(e1*rho - w)') X Diff = norm(e1*rho - w) SHAR_EOF $shar_touch -am 09160908104 'HHgood.m' && chmod 0644 'HHgood.m' || $echo 'restore of' 'HHgood.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'HHgood.m:' 'MD5 check failed' 978618020fecdf957d6a69d553b74e5d HHgood.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'HHgood.m'`" test 733 -eq "$shar_count" || $echo 'HHgood.m:' 'original size' '733,' 'current size' "$shar_count!" fi fi # ============= Hump.m ============== if test -f 'Hump.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Hump.m' '(file already exists)' else $echo 'x -' extracting 'Hump.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Hump.m' && % % This script contains a code that illustrates convergence % behavior for Jacobi, Gauss-Seidel, and SOR iterations. % % Example 0 illustrates behavior of all three on a well % behaved problem - asympotic convergence rate % is predictive % % Example 1 illustrates behavior of SOR on a highly % non-normal problem - asympotic convergence rate % is not predictive at all % % Example 2 illustrates behavior of Gauss Seidel on a flow % problem - with vs. against the direction of the % flow. Asympotic convergence rate is predictive % with the flow but not against it % %---------------------------------------------------------- % % D.C. Sorensen % 11 Sep 02 % % Slight modification of code from M. Embree (10 Sep 02) % X %%%%%%%%%%%%%%%%%%%%%% Example 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % PDE type matrix for 1-d diffusion reaction equation % % X alpha = 4; X X n = 100; X A = alpha*eye(n) - diag(ones(n-1,1),-1) - diag(ones(n-1,1),1); X X D = diag(diag(A)); E = -tril(A,-1); F = -triu(A,1); X X X % % Take b = 0, so that Ax=b has solution x=0. % Start with x_0 = 1e-2*ones(N,1); X X X % % Apply Jacobi iteration % X G = D\(E + F); X x0 = 1e-2*ones(n,1); X xk = x0; X maxit = 400; X resvec = [norm(A*x0)]; X for j=1:maxit X xk = G*xk; X resvec = [resvec;norm(A*xk)]; X if resvec(end) < 1e-15, break, end X end X figure(1), clf X semilogy([0:length(resvec)-1],resvec, 'g-', 'linewidth',2) X rho = max(abs(eig(G))); hold X SpecRadJac = max(abs(eig(G))) % % Apply Gauss-Seidel iteration % X G = (D - E)\F; X x0 = 1e-2*ones(n,1); X xk = x0; X maxit = 400; X resvec = [norm(A*x0)]; X for j=1:maxit X xk = G*xk; X resvec = [resvec;norm(A*xk)]; X if resvec(end) < 1e-15, break, end X end X semilogy([0:length(resvec)-1],resvec, 'k-', 'linewidth',2) X X SpecRadGS = max(abs(eig(G))) X % % Apply SOR iteration with w = 1.5 % X w = 1.5; X G = (D - w*E)\((1-w)*D + w*F); X x0 = 1e-2*ones(n,1); X xk = x0; X maxit = 400; X resvec = [norm(A*x0)]; X for j=1:maxit X xk = G*xk; X resvec = [resvec;norm(A*xk)]; X if resvec(end) < 1e-15, break, end X end X semilogy([0:length(resvec)-1],resvec, 'b-', 'linewidth',2) X SpecRadSOR = max(abs(eig(G))) X % % Apply SOR iteration with w = optimal % X w = 2/(1 + sqrt(1 - rho*rho)); X G = (D - w*E)\((1-w)*D + w*F); X x0 = 1e-2*ones(n,1); X xk = x0; X maxit = 400; X resvec = [norm(A*x0)]; X for j=1:maxit X xk = G*xk; X resvec = [resvec;norm(A*xk)]; X if resvec(end) < 1e-15, break, end X end X semilogy([0:length(resvec)-1],resvec, 'r-', 'linewidth',2) X X X xlabel('iteration, k', 'fontsize', 18) X ylabel('residual norm, ||r_k||', 'fontsize', 18) X legend('Jacobi','G-S','SOR','SOR-opt') hold X SpecRadOpt = max(abs(eig(G))) X disp('any key to continue') pause X %%%%%%%%%%%%%%%%%%%%%% Example 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % (from Quarteroni, Sacco, Saleri; matrix credited to an % NPL tech report by Hammarling & Wilkinson, 1976) % N.B. spectral radius of iteration matrix is 0.5. % % Apply SOR iteration with w = 1.5 % X n = 100; X A = 1.5*eye(n) + diag(ones(n-1,1),-1); X w = 1.5; X D = diag(diag(A)); E = -tril(A,-1); F = -triu(A,1); X G = (D - w*E)\((1-w)*D + w*F); X r = max(abs(eig(G))); t = 1; X % % Take b = 0, so that Ax=b has solution x=0. % Start with x_0 = 1e-2*ones(N,1); X X x0 = 1e-2*ones(n,1); X xk = x0; X maxit = 400; X resvec = [norm(A*x0)]; X atote = [t]; X for j=1:maxit X xk = G*xk; X resvec = [resvec;norm(A*xk)]; X t = t*r; X atote = [atote; t]; X if resvec(end) < 1e-15, break, end X end X figure(2), clf X semilogy([0:length(resvec)-1],resvec, 'b-', 'linewidth',2) X xlabel('iteration, k', 'fontsize', 18) X ylabel('residual norm, ||r_k||', 'fontsize', 18) X X adjust = resvec(end)/atote(end); X atote = atote*adjust; X X hold; X semilogy([310:length(resvec)-1],atote(311:end), 'g:', 'linewidth',2) X legend('actual','predicted') X hold; X X specradG = r X Con_ratio = resvec(382)/resvec(381) X X %break disp('any key to continue') pause %%%%%%%%%%%%%%%%%%%%%% Example 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % (from Elman and Chernesky, SINUM 30 (1993), 1268-1290) % 1-d convection diffusion problem % -nu * u''(x) + gamma*u'(x) = f(x) % for x \in [0,1] with u(0) = alpha, u(1) = beta, nu>0, gamma>0. % % We will solve this system using Gauss-Seidel. We will choose gamma >0 % and in this case the "flow" of G-S is in the direction of the % fluid flow (with the wind). % % We will then choose gamma < 0 and agian solve the system using Gauss-Seidel. % In this case the "flow" of G-S is opposite to the direction of the % fluid flow (against the wind). % % Another way to look at this is to keep gamma > 0 but order the nodes % (grid points) in standard order (with the wind) and then in reverse % order (against the wind). % The spectral radius is the same if we label the nodes "with the wind" % (i.e., from left to right) or "against the wind", but there is a % transient plateau if we order against the wind. (For the usual left-to-right % ordering of the unkowns, "with the wind" corresponds to normal Gauss-Seidel, % while "against the wind" is Gauss-Seidel with the roles of upper and lower % triangular parts swapped. X X n = 64; X nu = 1; X gamma = 3*n/2; X w = 1; X % % With Wind: Take gamma positive ... flow is to right % X a = -nu-gamma/(2*n); b = 2*nu; c = -nu+gamma/(2*n); X A = a*diag(ones(n-1,1),-1) + ... X b*eye(n) + ... X c*diag(ones(n-1,1),1); X Gdown = tril(A)\(-triu(A,1)); % downwind (with wind) G-S iteration matrix X % % Against Wind: Take gamma negative ... flow is to left % X X a = -nu+gamma/(2*n); b = 2*nu; c = -nu-gamma/(2*n); X A = a*diag(ones(n-1,1),-1) + ... X b*eye(n) + ... X c*diag(ones(n-1,1),1); X Gup = tril(A)\(-triu(A,1)); % upwind (against wind) G-S iteration matrix X % Solve Ax=0 X X b = zeros(n,1); X x0 = ones(n,1); X maxit = 150; X resdown = [norm(A*x0)]; X r = max(abs(eig(Gdown))); t = 1; X atote = [t]; X X xk = x0; X for j=1:maxit X xk = Gdown*xk; X resdown = [resdown;norm(A*xk)]; X t = t*r; X atote = [atote; t]; X if resdown(end) < 1e-15, break, end X end X adjust = resdown(end)/atote(end); X atote = atote*adjust; X X resup = [norm(A*x0)]; X xk = x0; X for j=1:maxit X xk = Gup*xk; X hold; X resup = [resup;norm(A*xk)]; X if resup(end) < 1e-15, break, end X end X hold; X X figure(3), clf X semilogy([0:length(resdown)-1],resdown, 'b-', 'linewidth',2); hold on X semilogy([40:length(resdown)-1],atote(41:end), 'g:', 'linewidth',2) X semilogy([0:length(resup)-1],resup, 'r-', 'linewidth',2 ); X hold; X axis([0 140 1e-15 10]) X legend('downwind','predicted','upwind') X xlabel('iteration, k', 'fontsize', 18) X hold; X ylabel('residual norm, ||r_k||', 'fontsize', 18) X X specrad_Gup = max(abs(eig(Gup))) X specrad_Gdown = max(abs(eig(Gdown))) X SHAR_EOF $shar_touch -am 10101449103 'Hump.m' && chmod 0644 'Hump.m' || $echo 'restore of' 'Hump.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Hump.m:' 'MD5 check failed' 868251b94f2650b528f3127542d1d373 Hump.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Hump.m'`" test 6773 -eq "$shar_count" || $echo 'Hump.m:' 'original size' '6773,' 'current size' "$shar_count!" fi fi # ============= Iram.m ============== if test -f 'Iram.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Iram.m' '(file already exists)' else $echo 'x -' extracting 'Iram.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Iram.m' && function [V,R,ritz] = Iram(A,which,k,m,f); % % % Input: A -- an n by n matrix % % % which -- a string indicating which eigenvalues are % wanted. Options are: % % 'LR', 'SR', 'LM','LA','SA' % % k -- a positive integer % % m -- a positive integer k < m << n % % Recommended: m > 2k % % f -- an nonzero n vector % % If f does not appear, f <- randn(n,1) will occur % % % Output: V -- an n by k orthogonal matrix % % R -- a k by k quasi upper triangular matrix % % ritz -- a vector containing Ritz estimates for the % k Ritz values (eigenvalues of R); % % % with norm(AV - VR) ~ norm(ritz) % % Hence, a partial real Schur form of order k % % % D.C. Sorensen % 2 March 2000 % X [n,n] = size(A); X if (nargin < 5), f = randn(n,1); end X X V = zeros(n,m); H = zeros(m,m); X X tol = sqrt(eps); ritz = 1; nitr = 0; X X ko = 1; X X msave = m; ksave = k; X while (norm(ritz) > tol & nitr < 200), X m = msave + 4; k = ksave + 4; X nitr = nitr+1; X X [V,H,f] = Arnold(A,V,H,f,ko,m); X X [mu, ritz]= select_shifts(H,which); X X Q = eye(m); X j = m; X while(j > k) % % Apply unwanted eigenvalues as shifts with QR steps % X [Q,H] = qrstep(Q,H,mu(j),1,m); % % Decrease j by 2 if m_j imaginary and by 1 if real % X if (abs(imag(mu(j))) > 0), X j = j-2; X else X j = j-1; X end X end X ko = j; X X ritz = norm(f)*ritz(1:ksave); X f = V*Q(:,ko+1)*H(ko+1,ko) + f*Q(m,ko); X V(:,1:ko) = V*Q(:,1:ko); X end % % Change basis to partial Schur form % X X [Q,R] = schur(H(1:ko,1:ko)); X V = V(:,1:ko)*Q; SHAR_EOF $shar_touch -am 03141529100 'Iram.m' && chmod 0664 'Iram.m' || $echo 'restore of' 'Iram.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Iram.m:' 'MD5 check failed' 36b57f438143b40e3447a18b6507a73d Iram.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Iram.m'`" test 1899 -eq "$shar_count" || $echo 'Iram.m:' 'original size' '1899,' 'current size' "$shar_count!" fi fi # ============= Iram1.m ============== if test -f 'Iram1.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Iram1.m' '(file already exists)' else $echo 'x -' extracting 'Iram1.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Iram1.m' && function [V,R,ritz] = Iram(A,which,k,m,f); % % % Input: A -- an n by n matrix % % % which -- a string indicating which eigenvalues are % wanted. Options are: % % 'LR', 'SR', 'LM','LA','SA' % % k -- a positive integer % % m -- a positive integer k < m << n % % Recommended: m > 2k % % f -- an nonzero n vector % % If f does not appear, f <- randn(n,1) will occur % % % Output: V -- an n by k orthogonal matrix % % R -- a k by k quasi upper triangular matrix % % ritz -- a vector containing Ritz estimates for the % k Ritz values (eigenvalues of R); % % % with norm(AV - VR) ~ norm(ritz) % % Hence, a partial real Schur form of order k % % % D.C. Sorensen % 2 March 2000 % X [n,n] = size(A); X if (nargin < 5), f = randn(n,1); end X X V = zeros(n,m); H = zeros(m,m); X X tol = sqrt(eps); ritz = 1; nitr = 0; X X ko = 1; X X while (norm(ritz) > tol & nitr < 200), X nitr = nitr+1; X X [V,H,f] = Arnold(A,V,H,f,ko,m); X X [mu, ritz]= select_shifts(H,which); X X Q = eye(m); X j = m; X while(j > k) % % Apply unwanted eigenvalues as shifts with QR steps % X [Q,H] = qrstep(Q,H,mu(j),1,m); % % Decrease j by 2 if m_j imaginary and by 1 if real % X if (abs(imag(mu(j))) > 0), X j = j-2; X else X j = j-1; X end X end X ko = j; X X ritz = norm(f)*ritz(1:ko); X f = V*Q(:,ko+1)*H(ko+1,ko) + f*Q(m,ko); X V(:,1:ko) = V*Q(:,1:ko); X end % % Change basis to partial Schur form % X X [Q,R] = schur(H(1:ko,1:ko)); X V = V(:,1:ko)*Q; SHAR_EOF $shar_touch -am 03141529100 'Iram1.m' && chmod 0664 'Iram1.m' || $echo 'restore of' 'Iram1.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Iram1.m:' 'MD5 check failed' 9d3f088e6cfae00cc3e04187998220cc Iram1.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Iram1.m'`" test 1833 -eq "$shar_count" || $echo 'Iram1.m:' 'original size' '1833,' 'current size' "$shar_count!" fi fi # ============= Iram_Cnv.m ============== if test -f 'Iram_Cnv.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Iram_Cnv.m' '(file already exists)' else $echo 'x -' extracting 'Iram_Cnv.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Iram_Cnv.m' && function [V,R,ritz] = Iram_Cnv(A,eigA,which,k,m,f); % % % Input: A -- an n by n matrix % % eigA -- eigenvalues of A % % which -- a string indicating which eigenvalues are % wanted. Options are: % % 'LR', 'SR', 'LM','LA','SA' % % k -- a positive integer % % m -- a positive integer k < m << n % % Recommended: m > 2k % % f -- an nonzero n vector % % If f does not appear, f <- randn(n,1) will occur % % % Output: V -- an n by k orthogonal matrix % % R -- a k by k quasi upper triangular matrix % % ritz -- a vector containing Ritz estimates for the % k Ritz values (eigenvalues of R); % % % with norm(AV - VR) ~ norm(ritz) % % Hence, a partial real Schur form of order k % % % D.C. Sorensen % 2 March 2000 % X nmax = 30; X ph1 = zeros(nmax,1); X ph2 = zeros(nmax,1); X X plot(real(eigA),imag(eigA),'k+'); X hold on; X X [n,n] = size(A); X if (nargin < 5), f = randn(n,1); end X X V = zeros(n,m); H = zeros(m,m); X X tol = sqrt(eps); ritz = 1; nitr = 0; X X ko = 1; X X while (norm(ritz) > tol & nitr < nmax), X nitr = nitr+1; X X [V,H,f] = Arnold(A,V,H,f,ko,m); X X [mu, ritz]= select_shifts(H,which); X X Q = eye(m); X j = m; X while(j > k) % % Apply unwanted eigenvalues as shifts with QR steps % X [Q,H] = qrstep(Q,H,mu(j),1,m); % % Decrease j by 2 if m_j imaginary and by 1 if real % X if (abs(imag(mu(j))) > 0), X j = j-2; X else X j = j-1; X end X end X ko = j; X X yvec = norm(f)*ritz; X ritz = norm(f)*ritz(1:ko); X f = V*Q(:,ko+1)*H(ko+1,ko) + f*Q(m,ko); X V(:,1:ko) = V*Q(:,1:ko); X X Ritz_est = ritz' X Nitrs = nitr X X X % % ---------------------Graphics to Display Convergence----------------- % X X X ph1(nitr) = plot(real(mu(1:ko)), imag(mu(1:ko)),'o'); X set(ph1(nitr),'MarkeredgeColor',[1 1 1]) X X ph2(nitr) = plot(real(mu(ko+1:m)), imag(mu(ko+1:m)),'o'); X set(ph2(nitr),'MarkeredgeColor',[1 1 1]) X end X X mlt = 4; X colors = bluemap(mlt*nitr); X X for j = 1:nitr X set(ph1(j),'Color',colors((mlt-1)*nitr+j,:)) % set(ph1(j),'Color',colors(mlt*nitr-2,:)) X set(ph1(j),'MarkerSize',7) X set(ph1(j),'MarkeredgeColor',[1 1 1]) X set(ph1(j),'MarkerFaceColor', colors((mlt-1)*nitr+j,:)); % set(ph1(j),'MarkerFaceColor', colors(mlt*nitr-2,:)); X X set(ph2(j),'Color',colors(j,:)) X set(ph2(j),'MarkerSize',7) X set(ph2(j),'MarkeredgeColor',[1 1 1]) X set(ph2(j),'MarkerFaceColor', colors(j,:)); X pause(.5) X end % % Change basis to partial Schur form % X X [Q,R] = schur(H(1:ko,1:ko)); X V = V(:,1:ko)*Q; X mu = eig(R); X ph = plot(real(mu), imag(mu),'ro'); X set(ph,'MarkerSize',7) X X hold off SHAR_EOF $shar_touch -am 03141529100 'Iram_Cnv.m' && chmod 0664 'Iram_Cnv.m' || $echo 'restore of' 'Iram_Cnv.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Iram_Cnv.m:' 'MD5 check failed' 6583c59d052ad1ceae05014d781353e7 Iram_Cnv.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Iram_Cnv.m'`" test 3062 -eq "$shar_count" || $echo 'Iram_Cnv.m:' 'original size' '3062,' 'current size' "$shar_count!" fi fi # ============= Iram_Plt.m ============== if test -f 'Iram_Plt.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Iram_Plt.m' '(file already exists)' else $echo 'x -' extracting 'Iram_Plt.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Iram_Plt.m' && function [V,R,ritz] = Iram_Plt(A,which,k,m,f); % % % Input: A -- an n by n matrix % % % which -- a string indicating which eigenvalues are % wanted. Options are: % % 'LR', 'SR', 'LM','LA','SA' % % k -- a positive integer % % m -- a positive integer k < m << n % % Recommended: m > 2k % % f -- an nonzero n vector % % If f does not appear, f <- randn(n,1) will occur % % % Output: V -- an n by k orthogonal matrix % % R -- a k by k quasi upper triangular matrix % % ritz -- a vector containing Ritz estimates for the % k Ritz values (eigenvalues of R); % % % with norm(AV - VR) ~ norm(ritz) % % Hence, a partial real Schur form of order k % % % D.C. Sorensen % 2 March 2000 % X [n,n] = size(A); X if (nargin < 5), f = randn(n,1); end X X V = zeros(n,m); H = zeros(m,m); X X tol = sqrt(eps); ritz = 1; nitr = 0; X X ko = 1; X msave = m; ksave = k; X while (norm(ritz) > tol & nitr < 200), X nitr = nitr+1; X m = msave+4; k = ksave+4; X X [V,H,f] = Arnold(A,V,H,f,ko,m); X X [mu, ritz]= select_shifts(H,which); X X Q = eye(m); X j = m; X while(j > k) % % Apply unwanted eigenvalues as shifts with QR steps % X [Q,H] = qrstep(Q,H,mu(j),1,m); % % Decrease j by 2 if m_j imaginary and by 1 if real % X if (abs(imag(mu(j))) > 0), X j = j-2; X else X j = j-1; X end X end X ko = j; X X yvec = norm(f)*ritz; X ritz = norm(f)*ritz(1:ksave); X f = V*Q(:,ko+1)*H(ko+1,ko) + f*Q(m,ko); X V(:,1:ko) = V*Q(:,1:ko); X X Ritz_est = ritz' X Nitrs = nitr X X X % % ---------------------Graphics to Display Convergence----------------- % X n1 = m; X x = ones(n1*(n1+1)/2 , 1); y = x; X phd = zeros(n1-1,1); X phy = zeros(n1-1,1); X X m = 16; X colors = bluemap(m); X x = [-1 -1 n1+2 n1+2]; X y = [-1 n1+2 -1 n1+2]; X X ph1 = plot(x,y,'o'); X axis tight; X set(ph1,'MarkeredgeColor',[1 1 1]) X hold on; X indx = 1; X Hmax = max(abs(diag(H,-1))); X for j = 1:n1, X for i = 1:j, X ii = n1-i+1; jj = n1-j+1; X x(indx) = ii; y(indx) = j; X indx = indx + 1; X end X jj = min(n1-1,j); X i = min(n1,jj+1); X gval = max(eps, abs(H(i,jj)/Hmax)); X gval = -log(gval); X gval = max(1,ceil(gval)); X gval = min(16,gval); X gval = 16 - gval + 1; X X yval = max(eps, abs(yvec(jj))); X yval = -log(yval); X yval = max(1,ceil(yval)); X yval = min(16,yval); X yval = 16 - yval + 1; X X ii = n1-i+1; jj = j; X phy(j) = plot(jj,0,'o'); X phd(j) = plot(jj,ii,'o'); X axis tight; X set(phd(j),'Color',colors(m,:)) X set(phd(j),'MarkerSize',9) X set(phd(j),'MarkeredgeColor',[1 1 1]) X X set(phy(j),'Color',colors(m,:)) X set(phy(j),'MarkerSize',9) X set(phy(j),'MarkeredgeColor',[1 1 1]) X X if(abs(H(i,j)) > 0), X set(phd(j),'MarkerFaceColor', colors(gval,:)); X end X if(abs(yvec(j)) > 100*eps), X set(phy(j),'MarkerFaceColor', colors(yval,:)); X end X end X ph = plot(x,y,'o'); X axis tight; X set(ph,'Color',colors(m,:)) X set(ph,'MarkerSize',9) X set(ph,'MarkeredgeColor',[1 1 1]) X set(ph,'MarkerFaceColor',colors(m,:)) X x1 = [ .7 k+.3]; X y1 = [n1-k-.3 n1+.3]; X x2 = [ .7 .7]; X x3 = [ k+.3 k+.3]; X y2 = [n1-k-.3 n1-k-.3]; X y3 = [n1+.3 n1+.3]; X plot(x1,y2,'m--',x1,y3,'m--',x2,y1,'m--',x3,y1,'m--'); X axis tight; X pause(.2) X hold off X X end % % Change basis to partial Schur form % X X [Q,R] = schur(H(1:ko,1:ko)); X V = V(:,1:ko)*Q; X SHAR_EOF $shar_touch -am 03141529100 'Iram_Plt.m' && chmod 0664 'Iram_Plt.m' || $echo 'restore of' 'Iram_Plt.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Iram_Plt.m:' 'MD5 check failed' d1a6298924bb14ebf23e3e7287c65e88 Iram_Plt.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Iram_Plt.m'`" test 4352 -eq "$shar_count" || $echo 'Iram_Plt.m:' 'original size' '4352,' 'current size' "$shar_count!" fi fi # ============= Iram_Plt1.m ============== if test -f 'Iram_Plt1.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Iram_Plt1.m' '(file already exists)' else $echo 'x -' extracting 'Iram_Plt1.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Iram_Plt1.m' && function [V,R,ritz] = Iram_Plt(A,which,k,m,f); % % % Input: A -- an n by n matrix % % % which -- a string indicating which eigenvalues are % wanted. Options are: % % 'LR', 'SR', 'LM','LA','SA' % % k -- a positive integer % % m -- a positive integer k < m << n % % Recommended: m > 2k % % f -- an nonzero n vector % % If f does not appear, f <- randn(n,1) will occur % % % Output: V -- an n by k orthogonal matrix % % R -- a k by k quasi upper triangular matrix % % ritz -- a vector containing Ritz estimates for the % k Ritz values (eigenvalues of R); % % % with norm(AV - VR) ~ norm(ritz) % % Hence, a partial real Schur form of order k % % % D.C. Sorensen % 2 March 2000 % X [n,n] = size(A); X if (nargin < 5), f = randn(n,1); end X X V = zeros(n,m); H = zeros(m,m); X X tol = sqrt(eps); ritz = 1; nitr = 0; X X ko = 1; X msave = m; ksave = k; X while (norm(ritz) > tol & nitr < 200), X nitr = nitr+1; X m = msave; k = ksave; X X [V,H,f] = Arnold(A,V,H,f,ko,m); X X [mu, ritz]= select_shifts(H,which); X X Q = eye(m); X j = m; X while(j > k) % % Apply unwanted eigenvalues as shifts with QR steps % X [Q,H] = qrstep(Q,H,mu(j),1,m); % % Decrease j by 2 if m_j imaginary and by 1 if real % X if (abs(imag(mu(j))) > 0), X j = j-2; X else X j = j-1; X end X end X ko = j; X X yvec = norm(f)*ritz; X ritz = norm(f)*ritz(1:ko); X f = V*Q(:,ko+1)*H(ko+1,ko) + f*Q(m,ko); X V(:,1:ko) = V*Q(:,1:ko); X X Ritz_est = ritz' X Nitrs = nitr X X X % % ---------------------Graphics to Display Convergence----------------- % X n1 = m; X x = ones(n1*(n1+1)/2 , 1); y = x; X phd = zeros(n1-1,1); X phy = zeros(n1-1,1); X X m = 16; X colors = bluemap(m); X x = [-1 -1 n1+2 n1+2]; X y = [-1 n1+2 -1 n1+2]; X X ph1 = plot(x,y,'o'); X axis tight; X set(ph1,'MarkeredgeColor',[1 1 1]) X hold on; X indx = 1; X Hmax = max(abs(diag(H,-1))); X for j = 1:n1, X for i = 1:j, X ii = n1-i+1; jj = n1-j+1; X x(indx) = ii; y(indx) = j; X indx = indx + 1; X end X jj = min(n1-1,j); X i = min(n1,jj+1); X gval = max(eps, abs(H(i,jj)/Hmax)); X gval = -log(gval); X gval = max(1,ceil(gval)); X gval = min(16,gval); X gval = 16 - gval + 1; X X yval = max(eps, abs(yvec(jj))); X yval = -log(yval); X yval = max(1,ceil(yval)); X yval = min(16,yval); X yval = 16 - yval + 1; X X ii = n1-i+1; jj = j; X phy(j) = plot(jj,0,'o'); X phd(j) = plot(jj,ii,'o'); X axis tight; X set(phd(j),'Color',colors(m,:)) X set(phd(j),'MarkerSize',9) X set(phd(j),'MarkeredgeColor',[1 1 1]) X X set(phy(j),'Color',colors(m,:)) X set(phy(j),'MarkerSize',9) X set(phy(j),'MarkeredgeColor',[1 1 1]) X X if(abs(H(i,j)) > 0), X set(phd(j),'MarkerFaceColor', colors(gval,:)); X end X if(abs(yvec(j)) > 100*eps), X set(phy(j),'MarkerFaceColor', colors(yval,:)); X end X end X ph = plot(x,y,'o'); X axis tight; X set(ph,'Color',colors(m,:)) X set(ph,'MarkerSize',9) X set(ph,'MarkeredgeColor',[1 1 1]) X set(ph,'MarkerFaceColor',colors(m,:)) X x1 = [ .7 k+.3]; X y1 = [n1-k-.3 n1+.3]; X x2 = [ .7 .7]; X x3 = [ k+.3 k+.3]; X y2 = [n1-k-.3 n1-k-.3]; X y3 = [n1+.3 n1+.3]; X plot(x1,y2,'m--',x1,y3,'m--',x2,y1,'m--',x3,y1,'m--'); X axis tight; X pause(.2) X hold off X X end % % Change basis to partial Schur form % X X [Q,R] = schur(H(1:ko,1:ko)); X V = V(:,1:ko)*Q; X SHAR_EOF $shar_touch -am 03141529100 'Iram_Plt1.m' && chmod 0664 'Iram_Plt1.m' || $echo 'restore of' 'Iram_Plt1.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Iram_Plt1.m:' 'MD5 check failed' 55c555bdb0012a528ede6a64d929c8c5 Iram_Plt1.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Iram_Plt1.m'`" test 4345 -eq "$shar_count" || $echo 'Iram_Plt1.m:' 'original size' '4345,' 'current size' "$shar_count!" fi fi # ============= LScondemo.m ============== if test -f 'LScondemo.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'LScondemo.m' '(file already exists)' else $echo 'x -' extracting 'LScondemo.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'LScondemo.m' && X % % This matlab script % is designed to demonstrate the validity of the bound. % % norm(x1 -x)/norm(x) .le. (cond(A)/cos(theta))*norm(b1 - b)/norm(b) % % where Ax = b , Ax1 = b1. % 2 2 % % i.e. x and x1 solve least squares problem for rhs b and rhs b1 % respectively % % Here the angle theta is increased while Pb (Orth. projection of b) % remains constant. % % b1 = b + tau*w , where tau = eps*norm(b) % % and norm(w) = 1 % % D.C. Sorensen 10 Oct 01 % X format short e X m = 100; X n = 30; X d = 20; X k = 10; X k1 = 3; X X Xout = zeros(k1,5); X A = randn(m,n); X [U,S,V] = svd(A,0); % % prepare rhs. b % X b = U(:,1:d)*randn(d,1); % % X z = randn(m,1); X z = z - U*(U'*z); X z = z - U*(U'*z); X z = z/norm(z); X X b0 = b/norm(b); % % Note bo is Pb for all of the constructed right hand sides % % while z is a unit vector orthogonal to Range(A) % X S0 = S; X fac = 10; X for itr = 1:k, X X b = b0 + fac*z; % % construct RHS b % % Note Pz = 0 so Pb = b0 % % and norm(b) = sqrt(1 + fac^2) % % since b0'*z = 0 and both are unit vectors. % X fac = fac*10; X X tau = norm(b)*eps; X w = randn(m,1); w = w/norm(w); X X s = diag(S0); X b1 = b + tau*w; X for j = 1:k1, % % increase condition by factor of 10 % X s(d+1:n) = s(d+1:n)/10; X S = diag(s); % % form new matrix and solve for x1 and x % X A1 = U*S*V'; X x1 = A1\b1; X x = A1\b; X condA1 = cond(A1); X xerr = norm(x1 - x)/norm(x); X berr1 = norm(b1 - b)/norm(b); X berr2 = norm(b1 - b)/norm(b0); X Xout(j,:) = [ xerr condA1*berr2 condA1*norm(b)/norm(b0) condA1*berr1 condA1 ]; X end X X disp(' ') X disp(' x-err estimate2 cond(A1)/ estimate1 cond(A1) ') X disp(' cos(theta) ') X disp(Xout) X disp(' ') X disp(' norm(b) 1/cos(theta) ') X disp([ norm(b) norm(b)/norm(b0)]) X pause X end X X X SHAR_EOF $shar_touch -am 09240813103 'LScondemo.m' && chmod 0644 'LScondemo.m' || $echo 'restore of' 'LScondemo.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'LScondemo.m:' 'MD5 check failed' 41ca089c4fbe50825b0cd0493d3d9b94 LScondemo.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'LScondemo.m'`" test 2013 -eq "$shar_count" || $echo 'LScondemo.m:' 'original size' '2013,' 'current size' "$shar_count!" fi fi # ============= LUfac.m ============== if test -f 'LUfac.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'LUfac.m' '(file already exists)' else $echo 'x -' extracting 'LUfac.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'LUfac.m' && X function [A,p] = LUfac(A); X % X % Input: A an n by n matrix. X % X % Output: A overwritten with the LU decompositon X % L in strict lower triangle, U in upper triangle X % X % p an integer array of length n containing pivot X % information. X % X % X % D.C. Sorensen X % 3 Oct 00 X %----------------------------------------------------------------- X % X [n,n] = size(A); X p = zeros(n,1); X for k = 1:n-1, X % X % Find index of max elt. below diagonal, adjust for offset X % and record in p(k) X % X [tk,ik] = max(abs(A(k:n,k))); p(k) = ik+k-1; X % X % Interchange row k with row p(k) X % X t = A(k,k:n); A(k,k:n) = A(p(k),k:n); A(p(k),k:n) = t; X % X % Scale by 1/A(k,k) to determine the kth column of L X % X A(k+1:n,k) = A(k+1:n,k)/A(k,k); X % X % Update lower (n-k) by (n-k) block X % X X for j = k+1:n, X for i = k+1:n, X X A(i,j) = A(i,j) - A(i,k)*A(k,j); X end X end X end X SHAR_EOF $shar_touch -am 08280936103 'LUfac.m' && chmod 0644 'LUfac.m' || $echo 'restore of' 'LUfac.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'LUfac.m:' 'MD5 check failed' 0b1e6ab49dfbfef2a2da3bf0311309af LUfac.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'LUfac.m'`" test 1094 -eq "$shar_count" || $echo 'LUfac.m:' 'original size' '1094,' 'current size' "$shar_count!" fi fi # ============= LUfacC.m ============== if test -f 'LUfacC.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'LUfacC.m' '(file already exists)' else $echo 'x -' extracting 'LUfacC.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'LUfacC.m' && X function [A,p] = LUfacC(A); X % X % LUfac organized by columns X % X % Input: A an n by n matrix. X % X % Output: A overwritten with the LU decompositon X % L in strict lower triangle, U in upper triangle X % X % p an integer array of length n containing pivot X % information. X % X % X % D.C. Sorensen X % 25 Aug 03 X %----------------------------------------------------------------- X % X [n,n] = size(A); X p = zeros(n,1); X for k = 1:n-1, X % X % Find index of max elt. below diagonal, adjust for offset X % and record in p(k) X % X [tk,ik] = max(abs(A(k:n,k))); p(k) = ik+k-1; X % X % Interchange row k with row p(k) X % X t = A(k,k:n); A(k,k:n) = A(p(k),k:n); A(p(k),k:n) = t; X % X % Scale by 1/A(k,k) to determine the kth column of L X % X A(k+1:n,k) = A(k+1:n,k)/A(k,k); X % X % Update lower (n-k) by (n-k) block X % X X for j = k+1:n, X A(k+1:n,j) = A(k+1:n,j) - A(k+1:n,k)*A(k,j); X end X end X SHAR_EOF $shar_touch -am 08280936103 'LUfacC.m' && chmod 0644 'LUfacC.m' || $echo 'restore of' 'LUfacC.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'LUfacC.m:' 'MD5 check failed' 1aac721ca4b259bb0ee7f47e934bd845 LUfacC.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'LUfacC.m'`" test 1098 -eq "$shar_count" || $echo 'LUfacC.m:' 'original size' '1098,' 'current size' "$shar_count!" fi fi # ============= LUfacP.m ============== if test -f 'LUfacP.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'LUfacP.m' '(file already exists)' else $echo 'x -' extracting 'LUfacP.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'LUfacP.m' && X function [L,U,P] = LUfacP(A); X % X % Input: A an n by n matrix. (A is unaltered on return) X % X % A is first overwritten with the LU decompositon X % L in strict lower triangle, U in upper triangle X % and then L,U are explicitly constructed such that X % X % PA = LU X % X % Output: L a unit lower triangular n by n matrix X % X % U an n by n upper triangular matrix X % X % P an n by n permutation matrix X % X % X % X % X % D.C. Sorensen X % 3 Oct 00 X %----------------------------------------------------------------- X % X [n,n] = size(A); X p = zeros(n,1); X for k = 1:n-1, X % X % Find index of max elt. below diagonal, adjust for offset X % and record in p(k) X % X [tk,ik] = max(abs(A(k:n,k))); p(k) = ik+k-1; X % X % Interchange row k with row p(k) X % X t = A(k,k:n); A(k,k:n) = A(p(k),k:n); A(p(k),k:n) = t; X % X % Scale by 1/A(k,k) to determine the kth column of L X % X A(k+1:n,k) = A(k+1:n,k)/A(k,k); X % X % Update lower (n-k) by (n-k) block X % X X for j = k+1:n, X for i = k+1:n, X A(i,j) = A(i,j) - A(i,k)*A(k,j); X end X end X end X % X % Apply permutations successively to the identity to get P X % and to the multipliers (lower triangle of overwritten A) X % Note: the k-th permutation is applied to columns 1 : k-1. X % X P = eye(n); X t = P(1,:); P(1,:) = P(p(1),:); P(p(1),:) = t; X X for k = 2:n-1; X t = A(k,1:k-1); A(k,1:k-1) = A(p(k),1:k-1); A(p(k),1:k-1) = t; X t = P(k,:); P(k,:) = P(p(k),:); P(p(k),:) = t; X end X L = eye(n) + tril(A,-1); X U = triu(A); SHAR_EOF $shar_touch -am 08280936103 'LUfacP.m' && chmod 0644 'LUfacP.m' || $echo 'restore of' 'LUfacP.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'LUfacP.m:' 'MD5 check failed' 2839f1207d31c88417640792f25f79a1 LUfacP.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'LUfacP.m'`" test 1819 -eq "$shar_count" || $echo 'LUfacP.m:' 'original size' '1819,' 'current size' "$shar_count!" fi fi # ============= LUfacR.m ============== if test -f 'LUfacR.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'LUfacR.m' '(file already exists)' else $echo 'x -' extracting 'LUfacR.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'LUfacR.m' && function [A,p] = LUfacR(A); X % X % LUfac organized by rank one update X % X % Input: A an n by n matrix. X % X % Output: A overwritten with the LU decompositon X % L in strict lower triangle, U in upper triangle X % X % p an integer array of length n containing pivot X % information. X % X % X % D.C. Sorensen X % 25 Aug 03 X %----------------------------------------------------------------- X % X [n,n] = size(A); X p = zeros(n,1); X for k = 1:n-1, X % X % Find index of max elt. below diagonal, adjust for offset X % and record in p(k) X % X [tk,ik] = max(abs(A(k:n,k))); p(k) = ik+k-1; X % X % Interchange row k with row p(k) X % X t = A(k,k:n); A(k,k:n) = A(p(k),k:n); A(p(k),k:n) = t; X % X % Scale by 1/A(k,k) to determine the kth column of L X % X A(k+1:n,k) = A(k+1:n,k)/A(k,k); X % X % Update lower (n-k) by (n-k) block X % X X A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n); X end X SHAR_EOF $shar_touch -am 08280937103 'LUfacR.m' && chmod 0644 'LUfacR.m' || $echo 'restore of' 'LUfacR.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'LUfacR.m:' 'MD5 check failed' afb64bec6338c5d9f2bbda146e58dbff LUfacR.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'LUfacR.m'`" test 1068 -eq "$shar_count" || $echo 'LUfacR.m:' 'original size' '1068,' 'current size' "$shar_count!" fi fi # ============= Lanconv_gui.m ============== if test -f 'Lanconv_gui.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Lanconv_gui.m' '(file already exists)' else $echo 'x -' extracting 'Lanconv_gui.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Lanconv_gui.m' && function Lanconv_update() X global EIG font_size = 10; X EIG.fig = figure('Color',[0.8 0.8 0.8], ... X 'Name','Eigenvalue Calculation', ... X 'Units','normalized', ... X 'Position',[.01 .05 .98 .90], ... X 'Tag','Fig1'); X EIG.axes(1) = axes('Parent',EIG.fig, ... X 'CameraUpVector',[0 1 0], ... X 'Color',[1 1 1], ... X 'FontSize',font_size, ... X 'FontWeight','bold', ... X 'LineWidth',1, ... X 'NextPlot','add', ... X 'Position',[0.01 0.09 0.98 0.40], ... X 'Tag','Axes1', ... X 'XColor',[0 0 0], ... X 'YColor',[0 0 0], ... X 'ZColor',[0 0 0]); X EIG.axes(2) = axes('Parent',EIG.fig, ... X 'CameraUpVector',[0 1 0], ... X 'Color',[1 1 1], ... X 'FontSize',font_size, ... X 'FontWeight','bold', ... X 'LineWidth',1, ... X 'NextPlot','add', ... X 'Position',[0.01 0.55 0.98 0.40], ... X 'Tag','Axes2', ... X 'XColor',[0 0 0], ... X 'YColor',[0 0 0], ... X 'ZColor',[0 0 0]); X EIG.str_ortho = uicontrol('Parent',EIG.fig, ... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',font_size, ... X 'FontWeight','bold', ... X 'Position',[0.01 0.00 0.53 0.05], ... X 'String','Ortho', ... X 'Style','text', ... X 'Tag','StaticText3'); X EIG.push_calc = uicontrol('Parent',EIG.fig, ... X 'Callback','Lanconv_update(''Calc'');',... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'ForegroundColor',[1 0 0], ... X 'FontSize',font_size, ... X 'FontWeight','bold', ... X 'Position',[0.55 0.00 0.09 0.05], ... X 'String','Calculate', ... X 'Tag','PushButton0' ... X ); X EIG.push_step = uicontrol('Parent',EIG.fig, ... X 'Callback','Lanconv_update(''Step'');',... X 'Interruptible','off',... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',font_size, ... X 'FontWeight','bold', ... X 'Position',[0.65 0.00 0.09 0.05], ... X 'String','Step', ... X 'Tag','PushButton2' ... X ); X EIG.push_continue = uicontrol('Parent',EIG.fig, ... X 'Callback','Lanconv_update(''Continue'');',... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',font_size, ... X 'FontWeight','bold', ... X 'Position',[0.75 0.00 0.09 0.05], ... X 'String','Continue', ... X 'Tag','PushButton3' ... X ); X EIG.push_exit = uicontrol('Parent',EIG.fig, ... X 'Callback','Lanconv_update(''Exit'');',... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',font_size, ... X 'FontWeight','bold', ... X 'Position',[0.90 0.00 0.09 0.05], ... X 'String','Exit', ... X 'Tag','PushButton4' ... X ); X %------------------------------------------- Lanconv_update('Init'); SHAR_EOF $shar_touch -am 02261255100 'Lanconv_gui.m' && chmod 0664 'Lanconv_gui.m' || $echo 'restore of' 'Lanconv_gui.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Lanconv_gui.m:' 'MD5 check failed' 2ac768e3fd69abda566391e1d4b6b958 Lanconv_gui.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Lanconv_gui.m'`" test 2481 -eq "$shar_count" || $echo 'Lanconv_gui.m:' 'original size' '2481,' 'current size' "$shar_count!" fi fi # ============= Lanconv_update.m ============== if test -f 'Lanconv_update.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Lanconv_update.m' '(file already exists)' else $echo 'x -' extracting 'Lanconv_update.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Lanconv_update.m' && function Lanconv_update(who) format short e X global EIG global n k T t m m0 fmax X true = 1; false= 0; X switch who X % --- initialize info --- X case {'Init'} X fprintf('%s\n',who); X X EIG.current_step = 0; X set(EIG.str_ortho,'String','Orthogonality Results...') X figure(EIG.fig) X hold off X EIG.calc = false; X set(EIG.push_calc,'ForegroundColor',[1 0 0]); X % --- calculate info --- X case {'Calc'} X fprintf('%s\n',who); X EIG.current_step = 0; X figure(EIG.fig) X % X % set up matrix X % X inc = 5; X n = 200; k = 100; X % Graph Resolution X X m = 5000; m0 = 500; X X A = randn(n); A = A + A'; X % X % X v = randn(n,1); %[V,T,f] = Lanczos(A,k,v); X [V,T,f] = ArnoldiC(A,k,v); X X ek = zeros(k,1); ek(k)= 1; X Resid = norm(A*V - V*T - f*ek') X orthtestV = norm(eye(k) - V'*V) X orthtestf = norm(V'*f) X set(EIG.str_ortho,'String',['res:' num2str(Resid) ', V:' num2str(orthtestV) ', f:' num2str(orthtestf)]); X X t = eig(A); X t = sort(t); X X X k0 = k; X fmax = 50; X X Nsteps = 10; X EIG.steps = 10:inc:min(k0-1, 10+Nsteps*inc); X EIG.calc = true; X X % --- step, next iteration --- X case {'Step'} X fprintf('%s\n',who); X if (EIG.calc == true & length(EIG.steps)>0) X EIG.iters = EIG.steps(1); X EIG.steps = EIG.steps(2:end); X Lanconv_update('Iter'); X end X % --- continue through all iterations --- X case {'Continue'} X fprintf('%s\n',who); X if (EIG.calc == true) X EIG.iters = EIG.steps; X EIG.steps = []; X Lanconv_update('Iter'); X end X % --- iteration --- X case {'Iter'} X fprintf('%s\n',who); X figure(EIG.fig); iterations=EIG.iters X for k = EIG.iters X X [Y,D] = eig(T(1:k,1:k)); X y = T(k+1,k)*Y(k,:)'; X [d,ii] = sort(diag(D)); X y = y(ii); X y = y.*y; X X Ynorm = norm(eye(k) - Y'*Y) X X k1 = 6; X t1 = [t(1:k1);t(n-k1+1:n)]; X d1 = [d(1:k1);d(k-k1+1:k)]; X t__d__diff = [t1 d1 abs(t1-d1)] X X X X h = (d(k) - d(k-k1+1))/(m-2*m0); X e = ones(m,1); X x = zeros(m,1); X x(1) = d(k-k1+1) - m0*h; X for j = 2:m, X x(j) = x(j-1) + h; X end X X f = zeros(m,1); X for j = 1:k, X f = f + y(j)*e./(x - e*d(j)); X end X f = -f; X g = e*T(k+1,k+1) - x; X X for j = 1:m, X if (abs(f(j)) > fmax), X f(j) = sign(f(j))*fmax; X end X end X vert = [-fmax:1:fmax]'; X dvert = ones(length(vert),1); X horiz = zeros(m,1); X X axes(EIG.axes(1)),plot(x,f,'k-') X hold on X axes(EIG.axes(1)),plot(x,g,'b-') X for j = k-k1+1:k X axes(EIG.axes(1)),plot(d(j)*dvert,vert,'c--') X end X axes(EIG.axes(1)),plot(x,horiz,'k-') X dd = sort(eig(T(1:k+1,1:k+1))); X axes(EIG.axes(1)),plot(dd(k-k1+2:k+1),zeros(k1,1),'r.') X X for j = k-k1+2:k+1 X [xx,jj] = min(abs(x - dd(j))); X fj = [0 f(jj)]; X xj = [dd(j) dd(j)]; X axes(EIG.axes(1)),plot(xj,fj,'r--') X end X X axes(EIG.axes(1)),title('Lanczos Convergence Right End') X hold off X X h = (d(k1) - d(1))/(m-2*m0); X e = ones(m,1); X x = zeros(m,1); X x(1) = d(1) - m0*h; X for j = 2:m, X x(j) = x(j-1) + h; X end X X f = zeros(m,1); X for j = 1:k, X f = f + y(j)*e./(x - e*d(j)); X end X f = -f; X g = e*T(k+1,k+1) - x; X X for j = 1:m, X if (abs(f(j)) > fmax), X f(j) = sign(f(j))*fmax; X end X end X vert = [-fmax:1:fmax]'; X dvert = ones(length(vert),1); X horiz = zeros(m,1); X X axes(EIG.axes(2)),plot(x,f,'k-') X hold on X plot(x,g,'b-') X for j = 1:k1 X axes(EIG.axes(2)),plot(d(j)*dvert,vert,'c--') X end X axes(EIG.axes(2)),plot(x,horiz,'k-') X dd = sort(eig(T(1:k+1,1:k+1))); X axes(EIG.axes(2)),plot(dd(1:k1),zeros(k1,1),'r.') X X for j = 1:k1 X [xx,jj] = min(abs(x - dd(j))); X fj = [0 f(jj)]; X xj = [dd(j) dd(j)]; X axes(EIG.axes(2)),plot(xj,fj,'r--') X end X X axes(EIG.axes(2)),title('Lanczos Convergence Left End') X hold off X pause(0.5) X end X % --- exit the program --- X case {'Exit'} X fprintf('%s\n',who); X close(EIG.fig); X % --- unknown request --- X otherwise X fprintf('\n\nUNKNOWN option for Lanconv_update -> %s\n',who); X EIG X end X SHAR_EOF $shar_touch -am 02261256100 'Lanconv_update.m' && chmod 0664 'Lanconv_update.m' || $echo 'restore of' 'Lanconv_update.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Lanconv_update.m:' 'MD5 check failed' 41ccb707383ff0f0047d7ce9551aa08c Lanconv_update.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Lanconv_update.m'`" test 4049 -eq "$shar_count" || $echo 'Lanconv_update.m:' 'original size' '4049,' 'current size' "$shar_count!" fi fi # ============= Lanczos.m ============== if test -f 'Lanczos.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Lanczos.m' '(file already exists)' else $echo 'x -' extracting 'Lanczos.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Lanczos.m' && X function [V,T,f] = Lanczos(A,k,v); % % Input: A -- an n by n matrix (A = A' assumed) % k -- a positive integer (k << n assumed) % v -- an n vector (v .ne. 0 assumed) % % Output: V -- an n by k orthogonal matrix % T -- a k by k symmetric tridiagonal matrix % f -- an n vector % % % with AV = VT + fe_k' % % In real life you would not store V, and you would store T as two % vectors a = diag(T) , b = diag(T,-1) % % D.C. Sorensen % 21 Feb 00 % X n = length(v); X T = zeros(k); X V = zeros(n,k); X X v1 = v/norm(v); X X f = A*v1; X alpha = v1'*f; X f = f - v1*alpha; X X V(:,1) = v1; T(1,1) = alpha; X X for j = 2:k, X X beta = norm(f); X v0 = v1; v1 = f/beta; X X f = A*v1 - v0*beta; X alpha = v1'*f; X f = f - v1*alpha; X X T(j,j-1) = beta; T(j-1,j) = beta; T(j,j) = alpha; X V(:,j) = v1; X X end X X SHAR_EOF $shar_touch -am 02241526100 'Lanczos.m' && chmod 0664 'Lanczos.m' || $echo 'restore of' 'Lanczos.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Lanczos.m:' 'MD5 check failed' ec0159a04988bf039bfc87af0eb0e689 Lanczos.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Lanczos.m'`" test 942 -eq "$shar_count" || $echo 'Lanczos.m:' 'original size' '942,' 'current size' "$shar_count!" fi fi # ============= QRcgs.m ============== if test -f 'QRcgs.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'QRcgs.m' '(file already exists)' else $echo 'x -' extracting 'QRcgs.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'QRcgs.m' && X function [Q,R,ierr] = QRcgs(A); % % Input: A - an m by n matrix with m .ge. n % % Output: Q,R - The Classical Gram Schmidt QR factorization % % ierr = 0 successful % j A not full rank detected at step j % % Factors A = QR , where A is m by n % Q'Q = I_n, R is n by n upper triangular % % Q is m by n % % D.C. Sorensen % 26 Oct 00 %----------------------------------------------------------------- X X [m,n] = size(A); X ierr = 0; X X rho = norm(A(:,1)); if ( rho == 0), ierr = 1; return; end X Q = [A(:,1)/rho]; R = [rho]; X for j = 2:n, X r = Q'*A(:,j); X q = A(:,j) - Q*r; X rho = norm(q); if ( rho == 0), ierr = j; return; end X q = q/rho; X Q = [Q,q]; R = [R r ; zeros(1,j-1) rho]; X end X X X X SHAR_EOF $shar_touch -am 08270937104 'QRcgs.m' && chmod 0644 'QRcgs.m' || $echo 'restore of' 'QRcgs.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'QRcgs.m:' 'MD5 check failed' 91030a53be76e5823bb020ed2f174350 QRcgs.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'QRcgs.m'`" test 782 -eq "$shar_count" || $echo 'QRcgs.m:' 'original size' '782,' 'current size' "$shar_count!" fi fi # ============= QRcgsF.m ============== if test -f 'QRcgsF.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'QRcgsF.m' '(file already exists)' else $echo 'x -' extracting 'QRcgsF.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'QRcgsF.m' && X function [Q,R,ierr] = QRcgsF(A); % % Input: A - an m by n matrix with m .ge. n % % Output: Q,R - The Classical Gram Schmidt QR factorization % % ierr = 0 successful % j A not full rank detected at step j % % Factors A = QR , where A is m by n % Q'Q = I_n, R is n by n upper triangular % % Maintains orthogonality using DGKS fix % % Q is m by n % % D.C. Sorensen % 26 Oct 00 % Modified 28 Sep 01 by DCS % included additional test % for very pathological cases % of rank deficient matrices. % %----------------------------------------------------------------- X X [m,n] = size(A); X ierr = 0; X X rho = norm(A(:,1)); if ( rho == 0), ierr = 1; break; end X Q = [A(:,1)/rho]; R = [rho]; X for j = 2:n, X r = Q'*A(:,j); X q = A(:,j) - Q*r; % % Apply DGKS correction % X c = Q'*q; X q = q - Q*c; X r = r + c; X X rho = norm(q); if ( rho == 0), ierr = j; break; end X q = q/rho; % % The following test asks if A(:,j) is numerically in Range(Q). % If the test is passed then rho is indistinguishable from 0 % when compared to norm(r) and the component q is meaningless % X if ( rho < 10*eps*norm(r)), % % Construct a random q and orthogonalize against Q % Since rho is indistinguishable from 0, q*rho will be of % negligible size for any unit vector q % X q = randn(m,1); X c = Q'*q; X q = q - Q*c; X c = Q'*q; X q = q - Q*c; X q = q/norm(q); X end X X Q = [Q,q]; R = [R r ; zeros(1,j-1) rho]; X X end X X X X X SHAR_EOF $shar_touch -am 09161400103 'QRcgsF.m' && chmod 0644 'QRcgsF.m' || $echo 'restore of' 'QRcgsF.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'QRcgsF.m:' 'MD5 check failed' d42bcb523ec15240df035741fbede703 QRcgsF.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'QRcgsF.m'`" test 1655 -eq "$shar_count" || $echo 'QRcgsF.m:' 'original size' '1655,' 'current size' "$shar_count!" fi fi # ============= QRfac.m ============== if test -f 'QRfac.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'QRfac.m' '(file already exists)' else $echo 'x -' extracting 'QRfac.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'QRfac.m' && function [QR,tau] = QRfac(A); % % Input: A - an m by n matrix with m .ge. n % % Output: QR - The Householder QR factorization % in product form. QR can overwrite A. % % tau - A vector of length n containing % "pivot" information % % Factors A = QR , where A is m by n % Q'Q = I_m, R is m by n upper triangular % % The leading n by n block of R is an % n by n upper triangular matrix R_n % % In practic A is overwritten by the % Q in product form and R_n in the % upper triangle % % D.C. Sorensen % 3 Oct 00 %----------------------------------------------------------------- X X [m,n] = size(A); X tau = zeros(n,1); X X for k = 1:n, % % Compute the Householder vector v % % I - tau(k)*v*v' is the transformation % X v = A(k:m,k); X rho = sign(v(1))*norm(v); X if (abs(rho) > 0), X v(1) = v(1) + rho; X tau(k) = v(1)/rho; X v = v/v(1); X X for j = k+1:n % % Apply the Householder transformation to the j-th column of A % X alpha = tau(k)*(v'*A(k:m,j)); X A(k:m,j) = A(k:m,j) - v*alpha; X X end X A(k:m,k) = v; X A(k,k) = -rho; X X end X end X QR = A; X X X X X SHAR_EOF $shar_touch -am 09240814103 'QRfac.m' && chmod 0644 'QRfac.m' || $echo 'restore of' 'QRfac.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'QRfac.m:' 'MD5 check failed' 52218eb6f61cd2d9a7532d81306f7679 QRfac.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'QRfac.m'`" test 1220 -eq "$shar_count" || $echo 'QRfac.m:' 'original size' '1220,' 'current size' "$shar_count!" fi fi # ============= QRfac2.m ============== if test -f 'QRfac2.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'QRfac2.m' '(file already exists)' else $echo 'x -' extracting 'QRfac2.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'QRfac2.m' && X function [QR,tau] = QRfac2(A); % % Input: A - an m by n matrix with m .ge. n % % Output: QR - The Householder QR factorization % in product form. QR can overwrite A. % % tau - A vector of length n containing % "pivot" information % % Factors A = QR , where A is m by n % Q'Q = I_m, R is m by n upper triangular % % The leading n by n block of R is an % n by n upper triangular matrix R_n % % In practic A is overwritten by the % Q in product form and R_n in the % upper triangle % % D.C. Sorensen % 3 Oct 00 %----------------------------------------------------------------- % X [m,n] = size(A); X tau = zeros(n,1); X X for k = 1:n, % % Compute the Householder vector v % % I - tau(k)*v*v' is the transformation % v = [1 ; A(k+1:m,k)]; % X rho = sign(A(k,k))*norm(A(k:m,k)); X if (abs(rho) > 0), X A(k,k) = A(k,k) + rho; X tau(k) = A(k,k)/rho; X A(k:m,k) = A(k:m,k)/A(k,k); X X for j = k+1:n % % Apply the Householder transformation to the j-th column of A % X alpha = tau(k)*(A(k:m,k)'*A(k:m,j)); X A(k:m,j) = A(k:m,j) - A(k:m,k)*alpha; X X end X A(k,k) = -rho; X X end X end X QR = A; X X X X X SHAR_EOF $shar_touch -am 09240814103 'QRfac2.m' && chmod 0644 'QRfac2.m' || $echo 'restore of' 'QRfac2.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'QRfac2.m:' 'MD5 check failed' b1939c724e44b5cde3618a235ea6ac66 QRfac2.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'QRfac2.m'`" test 1277 -eq "$shar_count" || $echo 'QRfac2.m:' 'original size' '1277,' 'current size' "$shar_count!" fi fi # ============= QRgivens.m ============== if test -f 'QRgivens.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'QRgivens.m' '(file already exists)' else $echo 'x -' extracting 'QRgivens.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'QRgivens.m' && X function [x,rhoLS] = QRgivens(A,b); % % Input: A - an m by n matrix with m .ge. n % % b - an m vector % % Output: x - The solution to min norm(b - Ax) % % rhoLS = min norm(b - Ax) % % Factors A = QR , where A is m by n % Q'Q = I_m, R is m by n upper triangular % by Givens method and reduces b at the % same time. % % Then the solution to Rx = b1 is obtained % % % D.C. Sorensen % 30 Oct 00 %----------------------------------------------------------------- X X [m,n] = size(A); X tau = zeros(n,1); X if (m == n), x = A\b; rhoLS = 0; break; end % % Reduce the initial n+1 by n+1 submatrix of [A b] % to upper triangular form % X R = [A(1:n+1,:) b(1:n+1)] ; X for j = 2:n+1, X for i = 1:j-1, X [c,s,r] = Givens(R(i,i),R(j,i)); X R(i,i) = r; R(j,i) = 0; X t = c*R(i, i+1:n+1) + s*R(j,i+1:n+1); X R(j,i+1:n+1) = -s*R(i, i+1:n+1) + c*R(j,i+1:n+1); X R(i, i+1:n+1) = t; X end X end % % Process the remaining m - (n+1) rows of [A b] % X for k = n+2:m, X a = [A(k,:) b(k)]; X for i = 1:n+1, X [c,s,r] = Givens(R(i,i),a(i)); X R(i,i) = r; a(i) = 0; X t = c*R(i, i+1:n+1) + s*a(i+1:n+1); X a(i+1:n+1) = -s*R(i, i+1:n+1) + c*a(i+1:n+1); X R(i, i+1:n+1) = t; X end X end % % Solve the upper triangular system Rx = b1 in place % b1 resides in R(1:n,n+1); % X x = R(1:n,n+1); X for j = n:-1:2 X x(j) = x(j)/R(j,j); X x(1:j-1) = x(1:j-1) - R(1:j-1,j)*x(j); X end X x(1) = x(1)/R(1,1); X rhoLS = abs(R(n+1,n+1)); SHAR_EOF $shar_touch -am 09171120103 'QRgivens.m' && chmod 0644 'QRgivens.m' || $echo 'restore of' 'QRgivens.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'QRgivens.m:' 'MD5 check failed' 95b3196e913073ee63d1426848762d4c QRgivens.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'QRgivens.m'`" test 1583 -eq "$shar_count" || $echo 'QRgivens.m:' 'original size' '1583,' 'current size' "$shar_count!" fi fi # ============= QRmgs.m ============== if test -f 'QRmgs.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'QRmgs.m' '(file already exists)' else $echo 'x -' extracting 'QRmgs.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'QRmgs.m' && X function [Q,R,ierr] = QRcgs(A); % % Input: A - an m by n matrix with m .ge. n % % Output: Q,R - The Modified Gram Schmidt QR factorization % % ierr = 0 successful % j A not full rank detected at step j % % Factors A = QR , where A is m by n % Q'Q = I_n, R is n by n upper triangular % % Q is m by n % % D.C. Sorensen % 26 Oct 00 %----------------------------------------------------------------- X X [m,n] = size(A); X ierr = 0; X R = zeros(n); Q = zeros(m,n); X X rho = norm(A(:,1)); if ( rho == 0), ierr = 1; return; end X Q(:,1) = [A(:,1)/rho]; R(1,1) = rho; X for j = 2:n, X q = A(:,j); X for i = 1:j-1 X R(i,j) = Q(:,i)'*q; X q = q - Q(:,i)*R(i,j); X end X rho = norm(q); if ( rho == 0), ierr = j; return; end X Q(:,j) = q/rho; R(j,j) = rho; X end X X X X SHAR_EOF $shar_touch -am 08270937104 'QRmgs.m' && chmod 0644 'QRmgs.m' || $echo 'restore of' 'QRmgs.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'QRmgs.m:' 'MD5 check failed' d5d8fcb633897582e529a15ce6ca34d7 QRmgs.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'QRmgs.m'`" test 858 -eq "$shar_count" || $echo 'QRmgs.m:' 'original size' '858,' 'current size' "$shar_count!" fi fi # ============= QRpiv.m ============== if test -f 'QRpiv.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'QRpiv.m' '(file already exists)' else $echo 'x -' extracting 'QRpiv.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'QRpiv.m' && X function [QR,tau,p] = QRpiv(A); % % Input: A - an m by n matrix with m .ge. n % % Output: QR - The Householder QR factorization % in product form. QR can overwrite A. % % tau - A vector of length n containing % "pivot" information % % p - an integer vector with column pivoting % % Factors A = QRP , where A is m by n % Q'Q = I_m, R is m by n upper triangular % % P is a permutation represented by p % % The leading n by n block of R is an % n by n upper triangular matrix R_n % % In practic A is overwritten by the % Q in product form and R_n in the % upper triangle % % D.C. Sorensen % 3 Oct 00 %----------------------------------------------------------------- X X [m,n] = size(A); X tau = zeros(n,1); X cnorm = zeros(n,1); X p = zeros(n,1); X for j = 1:n, X cnorm(j) = norm(A(:,j)); X end X for k = 1:n, % % compute pivot index and interchange columns % X [temp,jk] = max(cnorm(k:n)); p(k) = jk + k - 1; X t = A(:,k); A(:,k) = A(:,p(k)); A(:,p(k)) = t; X temp = cnorm(k); cnorm(k) = cnorm(p(k)); cnorm(p(k)) = temp; X % % Compute the Householder vector v % % I - tau(k)*v*v' is the transformation % X v = A(k:m,k); X rho = sign(v(1))*norm(v); X if (abs(rho) > 0), X v(1) = v(1) + rho; X tau(k) = v(1)/rho; X v = v/v(1); X X for j = k+1:n % % Apply the Householder transformation to the j-th column of A % X alpha = tau(k)*(v'*A(k:m,j)); X A(k:m,j) = A(k:m,j) - v*alpha; X X end X A(k:m,k) = v; X A(k,k) = -rho; X end % % update cnorm % X for j = k+1:n, X cnorm(j) = cnorm(j)*sqrt(1.0 - (A(k,j)/cnorm(j))^2); X end X end X QR = A; X X X X X SHAR_EOF $shar_touch -am 11071618100 'QRpiv.m' && chmod 0664 'QRpiv.m' || $echo 'restore of' 'QRpiv.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'QRpiv.m:' 'MD5 check failed' d32df6e0f74dbbced60b514d887f8f13 QRpiv.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'QRpiv.m'`" test 1755 -eq "$shar_count" || $echo 'QRpiv.m:' 'original size' '1755,' 'current size' "$shar_count!" fi fi # ============= QRpivD.m ============== if test -f 'QRpivD.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'QRpivD.m' '(file already exists)' else $echo 'x -' extracting 'QRpivD.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'QRpivD.m' && function [QR, Q, R, tau,p] = QRpivD(A); % % Input: A - an m by n matrix with m .ge. n % % Output: QR - The Householder QR factorization % in product form. QR can overwrite A. % % Q - Full m by m orthogonal Q % % R - m by n upper triangular R % % AP = QR (P represented by p -- see below) % % tau - A vector of length n containing % "pivot" information % % p - an integer vector with column pivoting % Note: AP does not equal A(:,p) % % Factors AP = QR , where A is m by n % Q'Q = I_m, R is m by n upper triangular % % P is a permutation represented by p % % The leading n by n block of R is an % n by n upper triangular matrix R_n % % In practic A is overwritten by the % Q in product form and R_n in the % upper triangle % % D.C. Sorensen % 3 Oct 00 %----------------------------------------------------------------- X X [m,n] = size(A); X tau = zeros(n,1); X R = zeros(m,n); X Q = eye(m); X X cnorm = zeros(n,1); X p = zeros(n,1); X for j = 1:n, X cnorm(j) = norm(A(:,j)); X end X for k = 1:n, % % compute pivot index and interchange columns % X [temp,jk] = max(cnorm(k:n)); p(k) = jk + k - 1; X t = A(:,k); A(:,k) = A(:,p(k)); A(:,p(k)) = t; X temp = cnorm(k); cnorm(k) = cnorm(p(k)); cnorm(p(k)) = temp; X % % Compute the Householder vector v % % I - tau(k)*v*v' is the transformation % X v = A(k:m,k); X rho = sign(v(1))*norm(v); X if (abs(rho) > 0), X v(1) = v(1) + rho; X tau(k) = v(1)/rho; X v = v/v(1); X X for j = k+1:n % % Apply the Householder transformation to the j-th column of A % X alpha = tau(k)*(v'*A(k:m,j)); X A(k:m,j) = A(k:m,j) - v*alpha; X X end X A(k:m,k) = v; X A(k,k) = -rho; X end % % update cnorm % X for j = k+1:n, X cnorm(j) = cnorm(j)*sqrt(1.0 - (A(k,j)/cnorm(j))^2); X end X Q(:,k:m) = Q(:,k:m) - tau(k)*(Q(:,k:m)*v)*v'; X R(1:k,k) = A(1:k,k); X end X QR = A; X X X X X SHAR_EOF $shar_touch -am 11071618100 'QRpivD.m' && chmod 0664 'QRpivD.m' || $echo 'restore of' 'QRpivD.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'QRpivD.m:' 'MD5 check failed' dc23e7537c62fbeba070e90299bede28 QRpivD.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'QRpivD.m'`" test 2083 -eq "$shar_count" || $echo 'QRpivD.m:' 'original size' '2083,' 'current size' "$shar_count!" fi fi # ============= SparseQRex.m ============== if test -f 'SparseQRex.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'SparseQRex.m' '(file already exists)' else $echo 'x -' extracting 'SparseQRex.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'SparseQRex.m' && X % This demo illustrates QR with column pivoting for % sparsity. A sparse matrix is generated with % columns permuted randomly to give an irregular % sparsity pattern to A'A. The symrcm ordering % is applied to this sparsity pattern to get % a permutation P such that AP has a good sparsity % pattern for obtaining a cholesky factor R that % is sparse (in this case banded). The rows % of AP can be processed one at a time via the % Row-wise Givens factorization to produce the % QR factorization of AP with the resulting sparsity % pattern. % % % See George and Heath, "Solution of Sparse Linear % Least Squares Problems Using Givens Rotations," % LAA 77, 165-187, (1980). % % D.C. Sorensen % 22 Sep 03 % % Set up Matrix % X n = 20; X a = 0; b = 2*pi; X h = (b-a)/(n+1); X X x = [a:h:b]'; X X m = 100; X h1 = (b-a)/(m+1); X t = [a:h1:b]'; X m = length(t); X % % set up a sparse A using % a B-spline basis for the interval [a,b] % partitioned by x % X m = length(t); n = length(x); X A = zeros(m,n); X for i = 1:m, X [j, a,b,c,d, ierr] = evalB(t(i),x); X if(j > 1), A(i,j-1) = a; end X A(i,j) = b; X if(j < n), A(i,j+1) = c; end X if (j < n-1), A(i,j+2) = d; end X end X % % Now give a random permutation to the columns % of A to give an irregular sparsity pattern % X p = randperm(n); X A = A(:,p); % % Demonstrate that the sparsity of A'A and the % resulting factorization A = QR are undesirable % X X X spy(A) X title('Sparsity of A','fontsize',20) X pause X X echo on X M = A'*A; X spy(M) X title('Sparsity of A''*A','fontsize',20) X pause X X [Q,R] = qr(A,0); X spy(R) X title('Sparsity of R','fontsize',20) X pause X % % Apply SYMRCM to re-order the columns of A % for a more favorable sparsity pattern % X X p1 = symrcm(M); X X A1 = A(:,p1); X spy(A1); X pause X X % % Show the modified sparsity pattern % X X M1 = A1'*A1; X spy(M1) X title('Sparsity of M1 = (AP)''*AP','fontsize',20) X pause X X [Q1,R1] = qr(A1,0); X spy(R1) X title('Sparsity of R1','fontsize',20) X X echo off SHAR_EOF $shar_touch -am 09240949103 'SparseQRex.m' && chmod 0644 'SparseQRex.m' || $echo 'restore of' 'SparseQRex.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'SparseQRex.m:' 'MD5 check failed' 818df0b4bf977a812c95ac60a6561827 SparseQRex.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'SparseQRex.m'`" test 2084 -eq "$shar_count" || $echo 'SparseQRex.m:' 'original size' '2084,' 'current size' "$shar_count!" fi fi # ============= Sym_ritz.m ============== if test -f 'Sym_ritz.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Sym_ritz.m' '(file already exists)' else $echo 'x -' extracting 'Sym_ritz.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'Sym_ritz.m' && X X % D.C. Sorensen X % 26 Feb, 2000 X % X X X format short e X X % Graph Resolution h = 1/m X X m = 5000; m0 = 500; X X % Increment between successive k values X X inc = 5; X X % X % set up matrix X % X n = 200; k = 100; X A = randn(n); A = A + A'; X % X % X v = randn(n,1); %[V,T,f] = Lanczos(A,k,v); X [V,T,f] = ArnoldiC(A,k,v); X X ek = zeros(k,1); ek(k)= 1; X Resid = norm(A*V - V*T - f*ek') X X orthtestV = norm(eye(k) - V'*V) X ortestf = norm(V'*f) X X t = eig(A); X t = sort(t); X X X k0 = k; X fmax = 50; X for k = 10:inc:min(k0-1, 10+5*inc), X X [Y,D] = eig(T(1:k,1:k)); X y = T(k+1,k)*Y(k,:)'; X [d,ii] = sort(diag(D)); X y = y(ii); X y = y.*y; X X Ynorm = norm(eye(k) - Y'*Y) X X k1 = 6; X t1 = [t(1:k1);t(n-k1+1:n)]; X d1 = [d(1:k1);d(k-k1+1:k)]; X t__d__diff = [t1 d1 abs(t1-d1)] X X X X h = (d(k) - d(k-k1+1))/(m-2*m0); X e = ones(m,1); X x = zeros(m,1); X x(1) = d(k-k1+1) - m0*h; X for j = 2:m, X x(j) = x(j-1) + h; X end X X f = zeros(m,1); X for j = 1:k, X f = f + y(j)*e./(x - e*d(j)); X end X f = -f; X g = e*T(k+1,k+1) - x; X X for j = 1:m, X if (abs(f(j)) > fmax), X f(j) = sign(f(j))*fmax; X end X end X vert = [-fmax:1:fmax]'; X dvert = ones(length(vert),1); X horiz = zeros(m,1); X X subplot(2,1,2),plot(x,f,'k-') X hold on X subplot(2,1,2),plot(x,g,'b-') X for j = k-k1+1:k X subplot(2,1,2),plot(d(j)*dvert,vert,'c--') X end X subplot(2,1,2),plot(x,horiz,'k-') X dd = sort(eig(T(1:k+1,1:k+1))); X subplot(2,1,2),plot(dd(k-k1+2:k+1),zeros(k1,1),'r.') X X for j = k-k1+2:k+1 X [xx,jj] = min(abs(x - dd(j))); X fj = [0 f(jj)]; X xj = [dd(j) dd(j)]; X subplot(2,1,2),plot(xj,fj,'r--') X end X X subplot(2,1,2),title('Lanczos Convergence Right End') X hold off X X h = (d(k1) - d(1))/(m-2*m0); X e = ones(m,1); X x = zeros(m,1); X x(1) = d(1) - m0*h; X for j = 2:m, X x(j) = x(j-1) + h; X end X X f = zeros(m,1); X for j = 1:k, X f = f + y(j)*e./(x - e*d(j)); X end X f = -f; X g = e*T(k+1,k+1) - x; X X for j = 1:m, X if (abs(f(j)) > fmax), X f(j) = sign(f(j))*fmax; X end X end X vert = [-fmax:1:fmax]'; X dvert = ones(length(vert),1); X horiz = zeros(m,1); X X subplot(2,1,1),plot(x,f,'k-') X hold on X plot(x,g,'b-') X for j = 1:k1 X subplot(2,1,1),plot(d(j)*dvert,vert,'c--') X end X subplot(2,1,1),plot(x,horiz,'k-') X dd = sort(eig(T(1:k+1,1:k+1))); X subplot(2,1,1),plot(dd(1:k1),zeros(k1,1),'r.') X X for j = 1:k1 X [xx,jj] = min(abs(x - dd(j))); X fj = [0 f(jj)]; X xj = [dd(j) dd(j)]; X subplot(2,1,1),plot(xj,fj,'r--') X end X X subplot(2,1,1),title('Lanczos Convergence Left End') X hold off X pause X end SHAR_EOF $shar_touch -am 02260847100 'Sym_ritz.m' && chmod 0664 'Sym_ritz.m' || $echo 'restore of' 'Sym_ritz.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Sym_ritz.m:' 'MD5 check failed' 421867a5da4d51a62698b8987632147e Sym_ritz.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Sym_ritz.m'`" test 2775 -eq "$shar_count" || $echo 'Sym_ritz.m:' 'original size' '2775,' 'current size' "$shar_count!" fi fi # ============= Threshtest.m ============== if test -f 'Threshtest.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'Threshtest.m' '(file already exists)' else $echo 'x -' extracting 'Threshtest.m' '(binary)' sed 's/^X//' << 'SHAR_EOF' | uudecode && begin 600 Threshtest.m M)0T*)2`@5&AI2`](%M=.PT*)0T*)2`@"!F"!M87)K M970-"B4-"@T*("`@96-H;R!O;CL-"B`@($$@/2!M;7)E860H)W=E""2A!*0T*("`@=&ET M;&4H)TYA='5R86P@3W)D97)I;FF4G+#$U*0T*("`@:&]L9`T*("`@9&ES<"AS<')I;G1F*"<@54U&('=I M=&@@5&AR97-H;VQD('!I=F]T:6YG+"`@=&AR97-H(#T@)60G+'1H2`](%M.6FAI2A,+"=R)RD-"B`@ M('1I=&QE*"=,($9A8W1OF5R;W-,52`](&YN>BA, M*2`K(&YN>BA5*0T*("`@3EIH:7-T;W)Y(#T@6TY::&ES=&]R>2`L($YO;GIE M2`G*0T*("`@9&ES<"@G("`G*0T*("`@<&%U2A,+"=R)RD-"B`@('1I=&QE*"=,($9A8W1OF5R;W-,52`](&YN>BA,*2`K(&YN>BA5*0T*("`@3EIH:7-T;W)Y M(#T@6TY::&ES=&]R>2`L($YO;GIE2`G*0T*("`@9&ES<"@G("`G*0T*("`@ M<&%UF4G+#$U*0T*("`@:&]L9`T*("`@9&ES M<"AS<')I;G1F*"<@54U&('=I=&@@5&AR97-H;VQD('!I=F]T:6YG+"`@=&AR M97-H(#T@)60G+'1H2`](%M.6FAI2!A(&-O;7!A2PN,RD-"B`@('AL86)E;"@G("!4:')E M&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'Threshtest.m:' 'MD5 check failed' dc990de8da4c48b69511cac3c1c761d0 Threshtest.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'Threshtest.m'`" test 3130 -eq "$shar_count" || $echo 'Threshtest.m:' 'original size' '3130,' 'current size' "$shar_count!" fi fi # ============= UNsparsedem.m ============== if test -f 'UNsparsedem.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'UNsparsedem.m' '(file already exists)' else $echo 'x -' extracting 'UNsparsedem.m' '(binary)' sed 's/^X//' << 'SHAR_EOF' | uudecode && begin 600 UNsparsedem.m M)0T*)2`@5&AI""2`G*0T*("`@9&ES<"@G("`G*0T*("`@<&%U2A!*2`G*0T*("`@9&ES<"@G("`G M*0T*("`@9&ES<"@G(%-T2`G*0T*("`@9&ES<"@G("`G*0T* M("`@<&%UF4G+#$U*0T*("`@:&]L M9`T*#0H@("!.;VYZ97)OF5R;W-,55T[#0H@("!D:7-P*"<@ M("2A!*#HL<"DI#0H@("!T:71L92@G M8V]L;6UD($]R9&5R:6YG)RPG1F]N='-I>F4G+#$U*0T*("`@9&ES<"@G('-P M>2A!*#HL<"DI("2A,+"=R)RD-"B`@('1I=&QE*"=, M($9A8W1O2`](%M.6FAI2`G M*0T*("`@9&ES<"@G("`G*0T*("`@<&%U2A5+"=G)RD-"B`@('1I=&QE*"=5($9A8W1O2`](%M.6FAI2A5+"=G)RD-"B`@('1I=&QE*"=5($9A8W1O MF5R;W-,55T[#0H@("!D:7-P*"<@("2A5+"=G)RD-"B`@('1I=&QE*"=5($9A8W1OF5R;W-,55T[#0H@("!D M:7-P*"<@("F4G+#$U M*0T*("`@:&]L9`T*("`@9&ES<"AS<')I;G1F*"<@54U&('=I=&@@5&AR97-H M;VQD('!I=F]T:6YG+"`@=&AR97-H(#T@)60G+'1H2`](%M. M6FAI2!A(&-O;7!A2PN,BD- M"B`@('AL86)E;"@G("!N871U&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'UNsparsedem.m:' 'MD5 check failed' a90b830d90106c549ac7d837f25e8d1f UNsparsedem.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'UNsparsedem.m'`" test 4319 -eq "$shar_count" || $echo 'UNsparsedem.m:' 'original size' '4319,' 'current size' "$shar_count!" fi fi # ============= bluemap.m ============== if test -f 'bluemap.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'bluemap.m' '(file already exists)' else $echo 'x -' extracting 'bluemap.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'bluemap.m' && function [colors] = bluemap(m) if (length(m)>1) X m=length(m); end %tail = [0.95; 0.8; 0.65; 0.0]; %tail = [0.95; 0.8; 0.65; .59]; tail = [0.95; 0.8; 0.7]; n=m-length(tail); X % red = .1+[[n:-1:1] [-0.1:1:n-0.1]]'/(4*n); % green = [n:-.5:0]'/n; % blue = ones(2*n+1,1); X red = [log([n:-2:1 2:2:n-2 1])]'/(4*log(n)); green = [log(n:-1:1)/log(n)]'; blue = ones(n,1); colors = [ red green blue ; tail*[0 0 1]]; SHAR_EOF $shar_touch -am 03141614100 'bluemap.m' && chmod 0664 'bluemap.m' || $echo 'restore of' 'bluemap.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'bluemap.m:' 'MD5 check failed' 63e6f6978d8738c71c498dd08e3ac903 bluemap.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'bluemap.m'`" test 417 -eq "$shar_count" || $echo 'bluemap.m:' 'original size' '417,' 'current size' "$shar_count!" fi fi # ============= compLinSlvrs.m ============== if test -f 'compLinSlvrs.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'compLinSlvrs.m' '(file already exists)' else $echo 'x -' extracting 'compLinSlvrs.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'compLinSlvrs.m' && % Code: compLinSlvrs % % This code compares the nonsymmetric iterative solvers from % Matlab: gmres, gmres5, bicgstab, bicg, qmr, cgs % % It is set up to solve the Sherman 2 problem from Harwell Bowing % collection (available from Matrix Market). The banded % matrix consiting of bands [-10:10] is used as a pre-conditioner. % % Assumes files sherman2.mtx sherman2_rhs1.mtx % % from Matrix Market exist % %------------------- % % D.C. Sorensen % 30 Mar 01 % X X A = mmread('sherman2.mtx'); X b = mmread('sherman2_rhs1.mtx'); X % % Construct a preconditioner from the diagonals [-10:10] % X d = [-10:10]; X B = spdiags(A,d); X [n,n] = size(A); X M = spdiags(B,d,n,n); X spy(A) X hold X spy(M,'r') X hold X X pause X X [L,U] = lu(M); X X spy(U) X hold X spy(L,'r') X hold X X pause X X spy(L,'r') X hold X spy(U) X hold X X pause X flpcnt = []; X X flops(0); X [x,flag,relres,iter,resvec] = gmres(A,b,n,1e-8,1,L,U); X flpcnt = [flpcnt, flops]; X X resvec = resvec/norm(b); X semilogy(resvec,'k') X hold X pause X X flops(0); X [x,flag,relres,iter,resvec] = gmres(A,b,5,1e-8,100,L,U); X flpcnt = [flpcnt, flops]; X X resvec = resvec/norm(b); X semilogy(resvec,'m') X pause X X X flops(0); X [x,flag,relres,iter,resvec] = bicgstab(A,b,1e-8,100,L,U); X flpcnt = [flpcnt, flops]; X resvec = resvec/norm(b); X semilogy(resvec) X pause X X flops(0); X [x,flag,relres,iter,resvec] = qmr(A,b,1e-8,100,L,U); X flpcnt = [flpcnt, flops]; X resvec = resvec/norm(b); X semilogy(resvec, 'c') X pause X X flops(0); X [x,flag,relres,iter,resvec] = bicg(A,b,1e-8,100,L,U); X flpcnt = [flpcnt, flops]; X resvec = resvec/norm(b); X semilogy(resvec,'g') X pause X X flops(0); X [x,flag,relres,iter,resvec] = cgs(A,b,1e-8,100,L,U); X flpcnt = [flpcnt, flops]; X resvec = resvec/norm(b); X semilogy(resvec,'r') X pause X hold X legend('gmres','gmres5','bicgstab','qmr','bicg','cgs') X X method = [' gmres ','gmres5 ', 'bicgstab ',' qmr ',' bicg ',' cgs '] X flpcnt = flpcnt X X sort(flpcnt') X SHAR_EOF $shar_touch -am 04051458101 'compLinSlvrs.m' && chmod 0664 'compLinSlvrs.m' || $echo 'restore of' 'compLinSlvrs.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'compLinSlvrs.m:' 'MD5 check failed' 552b4df5da25c31a5f30a495aaad1c59 compLinSlvrs.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'compLinSlvrs.m'`" test 1970 -eq "$shar_count" || $echo 'compLinSlvrs.m:' 'original size' '1970,' 'current size' "$shar_count!" fi fi # ============= driveA.m ============== if test -f 'driveA.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveA.m' '(file already exists)' else $echo 'x -' extracting 'driveA.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveA.m' && X % X % set up matrix X % X n = 100; k = 10; X X [Q,R] = qr(randn(n)); X %t = max(abs(diag(R))); X %R(n,n) = t*n; X A = Q*R*Q'; X % X % X v = randn(n,1); X [V,H,f] = Arnoldi(A,k,v); X X ek = zeros(k,1); ek(k)= 1; X Resid = norm(A*V - V*H - f*ek') X X X orthtestV = norm(eye(k) - V'*V) X ortestf = norm(V'*f) X X t = eig(A); X t = sort(abs(t)); X X s = sort(abs(eig(H))); X [t(n-k+1:n) s] SHAR_EOF $shar_touch -am 02101456100 'driveA.m' && chmod 0664 'driveA.m' || $echo 'restore of' 'driveA.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveA.m:' 'MD5 check failed' ff2835493852539621e0a4baa6eb7689 driveA.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveA.m'`" test 371 -eq "$shar_count" || $echo 'driveA.m:' 'original size' '371,' 'current size' "$shar_count!" fi fi # ============= driveA2.m ============== if test -f 'driveA2.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveA2.m' '(file already exists)' else $echo 'x -' extracting 'driveA2.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveA2.m' && X X [V,H,f] = Arnoldi(A,k,v); X X ek = zeros(k,1); ek(k)= 1; X Resid = norm(A*V - V*H - f*ek') X X X orthtestV = norm(eye(k) - V'*V) X ortestf = norm(V'*f) X X t = eig(A); X t = sort(abs(t)); X X s = sort(abs(eig(H))); X k1 = 10; X t__s__diff = [t(n-k1+1:n) s(k-k1+1:k) abs(t(n-k1+1:n) - s(k-k1+1:k)) ] SHAR_EOF $shar_touch -am 02101456100 'driveA2.m' && chmod 0664 'driveA2.m' || $echo 'restore of' 'driveA2.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveA2.m:' 'MD5 check failed' 3e9a9bd9f45392a81cf420246540d905 driveA2.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveA2.m'`" test 289 -eq "$shar_count" || $echo 'driveA2.m:' 'original size' '289,' 'current size' "$shar_count!" fi fi # ============= driveAC.m ============== if test -f 'driveAC.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveAC.m' '(file already exists)' else $echo 'x -' extracting 'driveAC.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveAC.m' && X X [V,H,f] = ArnoldiC(A,k,v); X X ek = zeros(k,1); ek(k)= 1; X Resid = norm(A*V - V*H - f*ek') X X X orthtestV = norm(eye(k) - V'*V) X ortestf = norm(V'*f) X X t = eig(A); X t = sort(abs(t)); X X s = sort(abs(eig(H))); X k1 = 10; X t__s__diff = [t(n-k1+1:n) s(k-k1+1:k) abs(t(n-k1+1:n) - s(k-k1+1:k)) ] SHAR_EOF $shar_touch -am 02101456100 'driveAC.m' && chmod 0664 'driveAC.m' || $echo 'restore of' 'driveAC.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveAC.m:' 'MD5 check failed' e52a0f1bd50dc083efe09a8b7de8bbad driveAC.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveAC.m'`" test 290 -eq "$shar_count" || $echo 'driveAC.m:' 'original size' '290,' 'current size' "$shar_count!" fi fi # ============= driveGM.m ============== if test -f 'driveGM.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveGM.m' '(file already exists)' else $echo 'x -' extracting 'driveGM.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveGM.m' && X % n = 100; % A = randn(n); b = randn(n,1); X d = [ones(40,1) ; 2*ones(40,1) ; 3*ones(40,1)]; X n = length(d); % d = d + .0001*randn(n,1); X b = randn(n,1); X [Q,R] = qr(randn(n)); X X tau = 1; X R = R*tau; X for j = 1:n, R(j,j) = d(j); end % R = diag(d); % R = diag(d) + .001*randn(n); X A = Q*R*Q'; X X X X k = 70; X tol = 100*sqrt(eps); X [x,r] = Gmres(A,b,k,tol); X X X X X semilogy(r) X hold X semilogy(r,'*') X hold X X X SHAR_EOF $shar_touch -am 04051446101 'driveGM.m' && chmod 0664 'driveGM.m' || $echo 'restore of' 'driveGM.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveGM.m:' 'MD5 check failed' 6b7de077d47918c62d3fca0101a1f1d7 driveGM.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveGM.m'`" test 466 -eq "$shar_count" || $echo 'driveGM.m:' 'original size' '466,' 'current size' "$shar_count!" fi fi # ============= driveGM0.m ============== if test -f 'driveGM0.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveGM0.m' '(file already exists)' else $echo 'x -' extracting 'driveGM0.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveGM0.m' && X function [outstring] = driveGM0(instring); % % Usage disp(driveGM0(instring)); % % Demonstrates GMRES convergence on various types of matrices % % for the linear system Ax = b % % Graphical displays are given of % % norm(b - Ax) convergence history % % placement of roots of GMRES and FOM polynomials % % After each 10 iterations a surface and contour % plot of log(abs(p(z))) where p(z) is the GMRES % polynomial. % % Cases are specified with % % Input: A character string % % instring = 'random' % % specifies order 50 random system % shows slow convergence % % instring = 'non-normal' % % specifies an order 120 non-normal system % with eigenvalues clustered at lambda = 1,2,3 % % instring = 'normal' % % specifies an order 120 normal system % with eigenvalues clustered at lambda = 1,2,3 % having the same distribution as the non-normal % system % % instring = 'symmetric_3' % % specifies an order 120 symmetric system % with eigenvalues at lambda = 1,2,3 % % instring = 'symmetric_cluster' % % specifies an order 120 symmetric system % with clustered eigenvalues at lambda = 1,2,3 % % % % D.C. Sorensen % % Fall 2002 % User interface added Oct. 2003 % X X MaxIter = 5; Itr = 1; X X switch lower(instring) X case 'random' X n = 50; X A = randn(n); b = randn(n,1); X s = eig(A); X outstring = 'Random order 50 system'; X X case {'normal', 'non-normal'} X d = [ones(20,1) ; 2*ones(60,1) ; 3*ones(40,1)]; X n = length(d); X b = randn(n,1); X [Q,R] = qr(randn(n)); % % change value of tau here to vary non-normality % larger tau gives higher non-normality % % tau = 0; % tau = .05; X tau = .5; X X R = R*tau; X for j = 1:n, R(j,j) = d(j); end X X A = Q*R*Q'; X s = eig(A); X X switch lower(instring) X case 'normal' X A = Q*diag(s)*Q'; X outstring = 'Normal system with 3 clusters of eigenvalues'; X otherwise X outstring = 'Non-normal system with 3 clusters of eigenvalues'; X end X X case 'symmetric_3' X X d = [ones(20,1) ; 2*ones(60,1) ; 3*ones(40,1)]; X n = length(d); X b = randn(n,1); X [Q,R] = qr(randn(n)); X A = Q*diag(d)*Q'; s = d; X outstring = 'Symmetric system with 3 eigenvalues'; X X case 'symmetric_cluster' X X d = [ones(20,1) ; 2*ones(60,1) ; 3*ones(40,1)]; X n = length(d); X b = randn(n,1); X [Q,R] = qr(randn(n)); X d = d + .05*randn(n,1); X A = Q*diag(d)*Q'; s = d; X outstring = 'Symmetric system with 3 clusters of eigenvalues'; X X otherwise X outstring = 'Unknown case'; return; X end % % Begin computing and displaying results % X r = []; X h1 = 1; X h2 = 2; X h3 = 3; X t1 = [1:12]; X t2 = (1.0e-8)*ones(12,1); X t2(1) = 10; t2(12) = 10; X for k = 1:1:n, X X [x,rho,ritz,hritz] = Gmres0(A,b,k); X X r = [r, norm(b - A*x)]; X X figure(h1) X clf X semilogy(t1,t2,'k.') X hold X semilogy(r) X semilogy(r,'*') X hold X X figure(h2) X clf X plot(real(s),imag(s),'k+') X hold X t = hritz; X plot(real(t),imag(t),'o') X t = ritz; X plot(real(t),imag(t),'r*') X legend('eigs','hritz','ritz') X hold X pause(.7) X if k == 1, X disp(' Arrange Graphs and Strike a Key') X pause X end X X if (mod(k,10) == 0 ), X disp('Strike a Key') X pause X logP = 0; X t = hritz; % % to see FOM poly put t = ritz % % t = ritz; % X figure(h2) X x = axis; X hx = (x(2)-x(1))/101; X hy = (x(4)-x(3))/101; X xx = x(1)+hx:hx:x(2)-hx; X yy = x(3)+hx:hy:x(4)-hy; X yy = x(3)+hy:hy:x(4)-hy; X m = length(yy); X n = length(xx); X e = ones(n,1); X X = e*xx + i*yy'*e'; X P = zeros(m,n); X p0 = sum(log(abs(t))); X for jj = 1:length(t), X P = P + log(abs(X - t(jj))); X end X P = P - p0*ones(m,n); % % P should have the values of log(abs(p(z)/p(0))) % for z in region of axis where p is GMRES poly % X hold X [cc,hh] = contour(xx,yy,P); % clabel(cc,hh); X hold X figure(h3) X clf X surfc(P); X colorbar X disp('Strike a Key') X Itr = Itr + 1; X pause X if (Itr > MaxIter), return, end; X end X end X X X X X X X SHAR_EOF $shar_touch -am 08270937104 'driveGM0.m' && chmod 0644 'driveGM0.m' || $echo 'restore of' 'driveGM0.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveGM0.m:' 'MD5 check failed' 65c03ee3063146da2cd3a10ba536fb7c driveGM0.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveGM0.m'`" test 5108 -eq "$shar_count" || $echo 'driveGM0.m:' 'original size' '5108,' 'current size' "$shar_count!" fi fi # ============= driveGMLU.m ============== if test -f 'driveGMLU.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveGMLU.m' '(file already exists)' else $echo 'x -' extracting 'driveGMLU.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveGMLU.m' && X % Code: compLinSlvrs % % This code compares GMRES with Precond GMRES % % It is set up to solve the Sherman 2 problem from Harwell Bowing % collection (available from Matrix Market). The banded % matrix consiting of bands [-10:10] is used as a pre-conditioner. % % Assumes files sherman2.mtx sherman2_rhs1.mtx % % from Matrix Market exist % %------------------- % % D.C. Sorensen % 30 Mar 01 % X X A = mmread('sherman2.mtx'); X b = mmread('sherman2_rhs1.mtx'); X X h1 = 1; X h2 = 2; X h3 = 3; X % % Construct a preconditioner from the diagonals [-10:10] % X d = [-10:10]; X B = spdiags(A,d); X [n,n] = size(A); X M = spdiags(B,d,n,n); X X figure(h1) X clf X spy(A) X hold X X disp('Strike Key') X pause X k = 150; X tol = 100*sqrt(eps); X [x,r] = Gmres(A,b,k,tol); X X X figure(h2) X clf X semilogy(r) X hold X semilogy(r,'*') X hold X X X X X X X disp('Strike Key') X pause X X figure(h1) X clf X spy(A) X hold X spy(M,'r') X hold X X disp('Strike Key') X pause X X dtol = 1e-6; X [L,U] = lu(M); X %opts = struct('droptol',dtol,'milu',1,'udiag',1); %[L,U] = luinc(A,opts); X X X spy(L,'r') X hold X spy(U) X hold X X disp('Strike Key') X pause X X X k = 150; X tol = 100*sqrt(eps); X [x,r] = GmresP(A,L,U,b,k,tol); X X X X X semilogy(r) X hold X semilogy(r,'*') X hold X X X SHAR_EOF $shar_touch -am 09271138102 'driveGMLU.m' && chmod 0644 'driveGMLU.m' || $echo 'restore of' 'driveGMLU.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveGMLU.m:' 'MD5 check failed' f8c11fd1337507a54ebf5d58c95ef1fc driveGMLU.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveGMLU.m'`" test 1312 -eq "$shar_count" || $echo 'driveGMLU.m:' 'original size' '1312,' 'current size' "$shar_count!" fi fi # ============= driveHHB.m ============== if test -f 'driveHHB.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveHHB.m' '(file already exists)' else $echo 'x -' extracting 'driveHHB.m' '(binary)' sed 's/^X//' << 'SHAR_EOF' | uudecode && begin 600 driveHHB.m M)0T*)2`@(%1H:7,@64H-BD@+2!1)RI1 M*0T*#0H-"G!A=7-E#0H-"@T*)2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM M+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM M+0T*#0IT:&5T82`]('-Q64H-BPQ*2`K('1H M971A*G)A;F1N*#8L,2D-"EMV+'1A=5T@/2!(2&)A9"AA*3L-"E$@/2!E>64H M-BD@+2!T874J=BIV)SL-"FYO64H-BPQ*2`K('1H971A*G)A;F1N M*#8L,2D-"EMV+'1A=5T@/2!(2&)A9"AA*3L-"E$@/2!E>64H-BD@+2!T874J M=BIV)SL-"FYO&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveHHB.m:' 'MD5 check failed' b78d869b99394fbf8d2d35a016a4a36c driveHHB.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveHHB.m'`" test 1216 -eq "$shar_count" || $echo 'driveHHB.m:' 'original size' '1216,' 'current size' "$shar_count!" fi fi # ============= driveHHG.m ============== if test -f 'driveHHG.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveHHG.m' '(file already exists)' else $echo 'x -' extracting 'driveHHG.m' '(binary)' sed 's/^X//' << 'SHAR_EOF' | uudecode && begin 600 driveHHG.m M)0T*)2`@(%1H:7,@2!I64H-BD@+2!1)RI1*0T*#0H-"G!A=7-E M#0H-"B4M+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM M+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2TM+2T-"@T*=&AE=&$@/2!S M<7)T*&5P64H-BPQ*2`K('1H971A*G)A;F1N*#8L,2D-"EMV M+'1A=5T@/2!(2&=O;V0H82D[#0I1(#T@97EE*#8I("T@=&%U*G8J=B<[#0IN M;W)M*&5Y92@V*2`M(%$G*E$I#0H-"@T*<&%U64H-BD@+2!T874J=BIV)SL-"FYO64H M-BPQ*2`K('1H971A*G)A;F1N*#8L,2D-"EMV+'1A=5T@/2!(2&=O;V0H82D[ M#0I1(#T@97EE*#8I("T@=&%U*G8J=B<[#0IN;W)M*&5Y92@V*2`M(%$G*E$I 4#0H-"@T*96-H;R!O9F8-"@T*#0HJ ` end SHAR_EOF $shar_touch -am 09160908104 'driveHHG.m' && chmod 0644 'driveHHG.m' || $echo 'restore of' 'driveHHG.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveHHG.m:' 'MD5 check failed' 56633996e5ac4d696137d8ab80943805 driveHHG.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveHHG.m'`" test 1370 -eq "$shar_count" || $echo 'driveHHG.m:' 'original size' '1370,' 'current size' "$shar_count!" fi fi # ============= driveIRA.m ============== if test -f 'driveIRA.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveIRA.m' '(file already exists)' else $echo 'x -' extracting 'driveIRA.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveIRA.m' && X format short e X n = 500; X X A = randn(n); f = randn(n,1); X X m = 20; k = 6; X X % % uncomment one and only one of the following selections % X X which = 'LR'; % select largest real part X % which = 'SR'; % select smallest real part X % which = 'LM'; % select largest magnitude X % A = A+A'; which = 'LA'; % select largest algebraic X % A = A+A'; which = 'SA'; % select smallest algebraic X X tic, X t = eig(A); X timeA = toc X X tic %[V,R,ritz] = Iram1(A,which,k,m,f); X [V,R,ritz] = Iram(A,which,k,m,f); X timeR = toc X X Resid__Ritz = [norm(A*V - V*R) norm(ritz)] X tr = eig(R); k1 = length(tr); X % X switch which X X case {'LR', 'LA'} % sort for largest real part % X [s,ii] = sort(-real(t)); t = t(ii); X [s,ii] = sort(-real(tr)); tr = tr(ii); X X X case {'SR', 'SA'} % sort for smallest real part % X [s,ii] = sort(real(t)); t = t(ii); X [s,ii] = sort(real(tr)); tr = tr(ii); X X X case {'LM'} % sort for largest magnitude % X [s,ii] = sort(-abs(t)); t = t(ii); X [s,ii] = sort(-abs(tr)); tr = tr(ii); X X end X X eigA__eigR__diff = [ t(1:k1) tr abs(abs(t(1:k1)) - abs(tr))] X X SHAR_EOF $shar_touch -am 03141530100 'driveIRA.m' && chmod 0664 'driveIRA.m' || $echo 'restore of' 'driveIRA.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveIRA.m:' 'MD5 check failed' 0a49bb91443fd7bee284e9237a9635fb driveIRA.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveIRA.m'`" test 1221 -eq "$shar_count" || $echo 'driveIRA.m:' 'original size' '1221,' 'current size' "$shar_count!" fi fi # ============= driveIRA_Cnv.m ============== if test -f 'driveIRA_Cnv.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveIRA_Cnv.m' '(file already exists)' else $echo 'x -' extracting 'driveIRA_Cnv.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveIRA_Cnv.m' && X X X format short e X n = 200; X X A = randn(n); f = randn(n,1); X X m = 20; k = 6; X X % % uncomment one and only one of the following selections % X X which = 'LR'; % select largest real part X % which = 'SR'; % select smallest real part X % which = 'LM'; % select largest magnitude X % A = A+A'; which = 'LA'; % select largest algebraic X % A = A+A'; which = 'SA'; % select smallest algebraic X X X t = eig(A); X [V,R,ritz] = Iram_Cnv(A,t,which,k,m,f); X X Resid__Ritz = [norm(A*V - V*R) norm(ritz)] X tr = eig(R); k1 = length(tr); X % X switch which X X case {'LR', 'LA'} % sort for largest real part % X [s,ii] = sort(-real(t)); t = t(ii); X [s,ii] = sort(-real(tr)); tr = tr(ii); X X X case {'SR', 'SA'} % sort for smallest real part % X [s,ii] = sort(real(t)); t = t(ii); X [s,ii] = sort(real(tr)); tr = tr(ii); X X X case {'LM'} % sort for largest magnitude % X [s,ii] = sort(-abs(t)); t = t(ii); X [s,ii] = sort(-abs(tr)); tr = tr(ii); X X end X X eigA__eigR__diff = [ t(1:k1) tr abs(abs(t(1:k1)) - abs(tr))] X X SHAR_EOF $shar_touch -am 03141530100 'driveIRA_Cnv.m' && chmod 0664 'driveIRA_Cnv.m' || $echo 'restore of' 'driveIRA_Cnv.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveIRA_Cnv.m:' 'MD5 check failed' 4346b054068b3fd587f298e371c8c270 driveIRA_Cnv.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveIRA_Cnv.m'`" test 1156 -eq "$shar_count" || $echo 'driveIRA_Cnv.m:' 'original size' '1156,' 'current size' "$shar_count!" fi fi # ============= driveIRA_Plt.m ============== if test -f 'driveIRA_Plt.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveIRA_Plt.m' '(file already exists)' else $echo 'x -' extracting 'driveIRA_Plt.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveIRA_Plt.m' && X X format short e X n = 500; X X A = randn(n); f = randn(n,1); X X m = 20; k = 5; X X % % uncomment one and only one of the following selections % X X which = 'LR'; % select largest real part X % which = 'SR'; % select smallest real part X % which = 'LM'; % select largest magnitude X % A = A+A'; which = 'LA'; % select largest algebraic X % A = A+A'; which = 'SA'; % select smallest algebraic X X X [V,R,ritz] = Iram_Plt(A,which,k,m,f); X X Resid__Ritz = [norm(A*V - V*R) norm(ritz)] X tr = eig(R); k1 = length(tr); X t = eig(A); X % X switch which X X case {'LR', 'LA'} % sort for largest real part % X [s,ii] = sort(-real(t)); t = t(ii); X [s,ii] = sort(-real(tr)); tr = tr(ii); X X X case {'SR', 'SA'} % sort for smallest real part % X [s,ii] = sort(real(t)); t = t(ii); X [s,ii] = sort(real(tr)); tr = tr(ii); X X X case {'LM'} % sort for largest magnitude % X [s,ii] = sort(-abs(t)); t = t(ii); X [s,ii] = sort(-abs(tr)); tr = tr(ii); X X end X X eigA__eigR__diff = [ t(1:k1) tr abs(abs(t(1:k1)) - abs(tr))] X X SHAR_EOF $shar_touch -am 03141530100 'driveIRA_Plt.m' && chmod 0664 'driveIRA_Plt.m' || $echo 'restore of' 'driveIRA_Plt.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveIRA_Plt.m:' 'MD5 check failed' 174752d0880749910c15b26ca60031a6 driveIRA_Plt.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveIRA_Plt.m'`" test 1153 -eq "$shar_count" || $echo 'driveIRA_Plt.m:' 'original size' '1153,' 'current size' "$shar_count!" fi fi # ============= driveIRA_SI.m ============== if test -f 'driveIRA_SI.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveIRA_SI.m' '(file already exists)' else $echo 'x -' extracting 'driveIRA_SI.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveIRA_SI.m' && % % This script illustrates convergence of IRAM with Shift-Invert % % % Recommended: Invoke % % figure(1); figure(2); % % and resize to taste prior to executing this script % % D.C. Sorensen % 14 Mar 2000 % X format short e X n = 500; X X A = randn(n); f = randn(n,1); X X m = 20; k = 5; X % % Do m steps of Arnoldi to get an idea of the location of the spectrum of A % X X H = zeros(m); V = zeros(n,m); [V,R,f] = Arnold(A,V,H,f,1,m); X % % Compute the eigenvalues of the projected matrix H % X X tr = eig(R); k1 = length(tr); X [s,ii] = sort(-real(tr)); tr = tr(ii); X % % Uncomment this statment if you want to see eigenvalues near the % % RIGHT END of the spectrum % X X sigma = real(tr(1)) + .1*abs(real(tr(1)) - real(tr(k1))); %--------------------------------------------------------- X % % Uncomment this statment if you want to see eigenvalues near the % % CENTER of the spectrum % X %sigma = sum(real(tr))/k1; %-------------------------- X X SIGMA = sigma X X f = randn(n,1); X X % % Compute A1 = (A - sigma I)^{-1} % X A1 = inv(A - sigma*eye(n)); X % % This should NOT BE DONE IN PRACTICE !!!! % % Instead do [L,U] = lu(A - sigma I) ONCE % % and solve w = U\(L\v) when a matrix-vector product % % is called for in Arnoldi % X X X which1 = 'LM'; X figure(1) X [V,R,ritz] = Iram_Plt(A1,which1,k,m,f); X X R = V'*(A*V); X X Resid__Ritz = [norm(A*V - V*R) norm(ritz)] X tr = eig(R); k1 = length(tr); X t = eig(A); X % % sort eigenvalues of A and R by order of distance from the % shift sigma % X [s,ii] = sort(abs(t-sigma)); t = t(ii); X X [s,ii] = sort(abs(tr-sigma)); tr = tr(ii); X % % Compute points on a circle centered at sigma % and just enclosing the converged Ritz values % X rho = abs(sigma - tr(k1)) X kpts = 50; X tau = 2*pi*sqrt(-1)/kpts; X s = zeros(kpts+1,1); X for j = 0:kpts, X s(j+1) = sigma + rho*exp(j*tau); X end X X X eigA__eigR__diff = [ t(1:k1) tr abs(abs(t(1:k1)) - abs(tr))] X % % Display eigenvalues of A with black +'s % % and Ritz values of R with red o's % % Draw a cyan dashed circle centered at sigma and enclosing % the computed Ritz values. % X X figure(2) X plot(real(t),imag(t),'+',real(tr),imag(tr),'ro',real(sigma),imag(sigma),'g*',real(s),imag(s),'c--') X SHAR_EOF $shar_touch -am 03141530100 'driveIRA_SI.m' && chmod 0664 'driveIRA_SI.m' || $echo 'restore of' 'driveIRA_SI.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveIRA_SI.m:' 'MD5 check failed' 111288893a91a23a28bc2fd00397b7ea driveIRA_SI.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveIRA_SI.m'`" test 2293 -eq "$shar_count" || $echo 'driveIRA_SI.m:' 'original size' '2293,' 'current size' "$shar_count!" fi fi # ============= driveL.m ============== if test -f 'driveL.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveL.m' '(file already exists)' else $echo 'x -' extracting 'driveL.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveL.m' && X X format short e X % X % set up matrix X % X n = 200; k = 10; X X A = randn(n); A = A + A'; X % X % X v = randn(n,1); X [V,T,f] = Lanczos(A,k,v); X X ek = zeros(k,1); ek(k)= 1; X Resid = norm(A*V - V*T - f*ek') X X X orthtestV = norm(eye(k) - V'*V) X ortestf = norm(V'*f) X X t = eig(A); X t = sort(t); X X s = sort(eig(T)); X k1 = 8; X t1 = [t(1:k1);t(n-k1+1:n)]; X s1 = [s(1:k1);s(k-k1+1:k)]; X t__s__diff = [t1 s1 abs(t1-s1)] SHAR_EOF $shar_touch -am 02241527100 'driveL.m' && chmod 0664 'driveL.m' || $echo 'restore of' 'driveL.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveL.m:' 'MD5 check failed' da98206e5f723f210aadc4615f79c1b1 driveL.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveL.m'`" test 411 -eq "$shar_count" || $echo 'driveL.m:' 'original size' '411,' 'current size' "$shar_count!" fi fi # ============= driveL2.m ============== if test -f 'driveL2.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveL2.m' '(file already exists)' else $echo 'x -' extracting 'driveL2.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveL2.m' && X X X format short e X X [V,T,f] = Lanczos(A,k,v); X X ek = zeros(k,1); ek(k)= 1; X Resid = norm(A*V - V*T - f*ek') X X X orthtestV = norm(eye(k) - V'*V) X ortestf = norm(V'*f) X X X t = eig(A); X t = sort(t); X X s = sort(eig(T)); X k1 = 8; X t1 = [t(1:k1);t(n-k1+1:n)]; X s1 = [s(1:k1);s(k-k1+1:k)]; X t__s__diff = [t1 s1 abs(t1-s1)] SHAR_EOF $shar_touch -am 02241527100 'driveL2.m' && chmod 0664 'driveL2.m' || $echo 'restore of' 'driveL2.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveL2.m:' 'MD5 check failed' e5131f4b3cd418412c9d008d25a63eec driveL2.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveL2.m'`" test 315 -eq "$shar_count" || $echo 'driveL2.m:' 'original size' '315,' 'current size' "$shar_count!" fi fi # ============= driveP.m ============== if test -f 'driveP.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveP.m' '(file already exists)' else $echo 'x -' extracting 'driveP.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveP.m' && X % X % set up matrix X % X n = 10; X X [Q,R] = qr(randn(n)); X d = randn(n,1); X [t,ii] = sort(abs(d)); X d(ii(n)) = 10*d(ii(n-1)); X A = Q*diag(d)*Q'; X % X % X pow_rq(A) X pause X X pow_inf(A) SHAR_EOF $shar_touch -am 02101456100 'driveP.m' && chmod 0664 'driveP.m' || $echo 'restore of' 'driveP.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveP.m:' 'MD5 check failed' 50d07bd07d9190eb6514b60e418bb4f2 driveP.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveP.m'`" test 189 -eq "$shar_count" || $echo 'driveP.m:' 'original size' '189,' 'current size' "$shar_count!" fi fi # ============= driveQR.m ============== if test -f 'driveQR.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveQR.m' '(file already exists)' else $echo 'x -' extracting 'driveQR.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveQR.m' && % % Drives the shifted QR iteration qriter % format short e n = 25; A = randn(n); X [V,H] = qriter(A); X s = eig(H); t = eig(A); s = sort(abs(s)); t = sort(abs(t)); X [t s t-s] SHAR_EOF $shar_touch -am 03021402100 'driveQR.m' && chmod 0664 'driveQR.m' || $echo 'restore of' 'driveQR.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveQR.m:' 'MD5 check failed' d868dba180a89686ad7642f4f6644b98 driveQR.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveQR.m'`" test 175 -eq "$shar_count" || $echo 'driveQR.m:' 'original size' '175,' 'current size' "$shar_count!" fi fi # ============= driveQRR.m ============== if test -f 'driveQRR.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveQRR.m' '(file already exists)' else $echo 'x -' extracting 'driveQRR.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveQRR.m' && % % Drives the implicit polynomial restart QR method qriterR % n = 100; A = randn(n); %A = A+A'; X [V,H] = qriterR(A); X SHAR_EOF $shar_touch -am 03021402100 'driveQRR.m' && chmod 0664 'driveQRR.m' || $echo 'restore of' 'driveQRR.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveQRR.m:' 'MD5 check failed' 715caa8b02935e86c2efea6490f6912f driveQRR.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveQRR.m'`" test 120 -eq "$shar_count" || $echo 'driveQRR.m:' 'original size' '120,' 'current size' "$shar_count!" fi fi # ============= driveQRRS.m ============== if test -f 'driveQRRS.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveQRRS.m' '(file already exists)' else $echo 'x -' extracting 'driveQRRS.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveQRRS.m' && % % Drives the implicit polynomial restart QR iteration qriterRS % for symmetric matrix A % n = 100; A = randn(n); A = A+A'; X [V,H] = qriterRS(A); X %s = eig(H); t = eig(A); %s = sort(abs(s)); t = sort(abs(t)); % %[t s t-s] SHAR_EOF $shar_touch -am 03021402100 'driveQRRS.m' && chmod 0664 'driveQRRS.m' || $echo 'restore of' 'driveQRRS.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveQRRS.m:' 'MD5 check failed' 98d517f6083cfa97d34f5ba66c8cd2f3 driveQRRS.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveQRRS.m'`" test 227 -eq "$shar_count" || $echo 'driveQRRS.m:' 'original size' '227,' 'current size' "$shar_count!" fi fi # ============= driveQRS.m ============== if test -f 'driveQRS.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'driveQRS.m' '(file already exists)' else $echo 'x -' extracting 'driveQRS.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'driveQRS.m' && % % drives the symmetric QR iteration qriterS % n = 25; A = randn(n); A = A+A'; X [V,H] = qriterS(A); X s = eig(H); t = eig(A); s = sort(abs(s)); t = sort(abs(t)); X [t s t-s] SHAR_EOF $shar_touch -am 03021402100 'driveQRS.m' && chmod 0664 'driveQRS.m' || $echo 'restore of' 'driveQRS.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'driveQRS.m:' 'MD5 check failed' 9955d5dff45e7a12d71770dc82378472 driveQRS.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'driveQRS.m'`" test 175 -eq "$shar_count" || $echo 'driveQRS.m:' 'original size' '175,' 'current size' "$shar_count!" fi fi # ============= eigconv.m ============== if test -f 'eigconv.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'eigconv.m' '(file already exists)' else $echo 'x -' extracting 'eigconv.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'eigconv.m' && X X % X % set up matrix X % X n = 40; k = 40; X X A = randn(n); X %A = A + A'; X % X % X v = randn(n,1); X [V,H,f] = ArnoldiC(A,k,v); X X ek = zeros(k,1); ek(k)= 1; X Resid = norm(A*V - V*H - f*ek') X X orthtestV = norm(eye(k) - V'*V) X X t = eig(A); X plot(real(t),imag(t),'k+'); X hold X %for k = 4:4:n, X h = n/2; X for k = 1:h, X s = eig(H(1:k,1:k)); X plot(real(s),imag(s),'c.'); X nvals = k X pause X end X for k = h+1:2*h, X s = eig(H(1:k,1:k)); X plot(real(s),imag(s),'m.'); X nvals = k X pause X end X plot(real(t),imag(t),'bo'); X plot(real(t),imag(t),'bd'); X plot(real(t),imag(t),'b+'); X hold X SHAR_EOF $shar_touch -am 02171349100 'eigconv.m' && chmod 0664 'eigconv.m' || $echo 'restore of' 'eigconv.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'eigconv.m:' 'MD5 check failed' 0edbd815695c27c4a3dd7d93d4ccef31 eigconv.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'eigconv.m'`" test 600 -eq "$shar_count" || $echo 'eigconv.m:' 'original size' '600,' 'current size' "$shar_count!" fi fi # ============= eigen_update.m ============== if test -f 'eigen_update.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'eigen_update.m' '(file already exists)' else $echo 'x -' extracting 'eigen_update.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'eigen_update.m' && function eigen_update(who) X global EIG global n k H t symmetric X true = 1; false= 0; X switch who X % --- initialize info --- X case {'Init' 'Function' 'Dim' 'Sym' 'Steps'} X fprintf('%s\n',who); X EIG.nvals = 0; X set(EIG.edit_res,'String','') X set(EIG.edit_ortho,'String','') X figure(EIG.fig) X hold off X EIG.calc = false; X set(EIG.push_calc,'ForegroundColor',[1 0 0]); X symmetric = false; X switch get(EIG.pop_sym,'Value') X case 1 X symmetric = false; X case 2 X symmetric = true; X case 3 X symmetric = true; X otherwise X symmetric = false; X end X switch get(EIG.pop_fcn,'Value') X case 1 X fprintf('ArnoldiC\n'); X otherwise X fprintf('Unknown\n'); X end X % --- calculate info --- X case {'Calc'} X fprintf('%s\n',who); X EIG.nvals = 0; X figure(EIG.fig) X hold off X % X % set up matrix X % X n = str2num(get(EIG.edit_dim,'String')); X k = str2num(get(EIG.edit_steps,'String'));; X X A = randn(n); X if (symmetric == true) X A = A + A'; X end X X v = randn(n,1); X [V,H,f] = ArnoldiC(A,k,v); X X ek = zeros(k,1); ek(k)= 1; X Resid = norm(A*V - V*H - f*ek'); X set(EIG.edit_res,'String',num2str(Resid)) X X orthtestV = norm(eye(k) - V'*V); X set(EIG.edit_ortho,'String',num2str(orthtestV)) X X t = eig(A); X X EIG.calc = true; X set(EIG.push_calc,'ForegroundColor',[0 0 0]); X set(EIG.push_step,'String','Step') X % --- plot eigenvalues --- X case {'PlotEigenValues'} X fprintf('%s\n',who); X X if (EIG.calc == true) X figure(EIG.fig); X X ph = plot(real(t),imag(t),'k+'); X hold on X if (get(EIG.pop_sym,'Value')==3) X ax=axis; ax = [ax(1) ax(2) -1 n]; axis(ax); X ph=plot([t t NaN*t]', [zeros(size(t)) n*ones(size(t)) NaN*ones(size(t))]','r-'); X set(ph,'Color',[.9 .9 .9]); X ph = plot(real(t),imag(t),'k+'); X end X end X % --- step, next iteration --- X case {'Step'} X fprintf('%s\n',who); X if (EIG.calc == true) X EIG.end_iter = EIG.nvals+1; X eigen_update('Iter'); X end X % --- continue through all iterations --- X case {'Continue'} X fprintf('%s\n',who); X if (EIG.calc == true) X EIG.end_iter = k; X eigen_update('Iter'); X end X % --- iteration --- X case {'Iter'} X fprintf('%s\n',who); X figure(EIG.fig); X if (EIG.calc == true) X m = ceil(n/2); X colors = [ .1+[[m:-1:1] [-0.1:1:m-0.1]]'/4 [m:-.5:0]' m*ones(2*m+1,1) ]/m; X for j = EIG.nvals+1:EIG.end_iter X X s = eig(H(1:j,1:j)); X X if (get(EIG.pop_sym,'Value')==3) X ph = plot(real(s),n-j+1,'o'); X else X ph = plot(real(s),imag(s),'o'); X end X set(ph,'Color',colors(j,:)) X set(ph,'MarkerSize',7) X set(ph,'MarkeredgeColor',[1 1 1]) X set(ph,'MarkerFaceColor',colors(j,:)) X EIG.nvals = j; X pause(.5) X end X fprintf('\tstep: %d\n',j); X set(EIG.push_step,'String',['Step (' num2str(j) ')']) X end % --- exit the program --- X case {'Exit'} X fprintf('%s\n',who); X close(EIG.fig); X % --- unknown request --- X otherwise X fprintf('\n\nUNKNOWN option for eigen_update -> %s\n',who); X EIG X end X SHAR_EOF $shar_touch -am 02231621100 'eigen_update.m' && chmod 0664 'eigen_update.m' || $echo 'restore of' 'eigen_update.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'eigen_update.m:' 'MD5 check failed' e0d9a512483746d81e0bf5b58c2c971c eigen_update.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'eigen_update.m'`" test 2954 -eq "$shar_count" || $echo 'eigen_update.m:' 'original size' '2954,' 'current size' "$shar_count!" fi fi # ============= eigenui.m ============== if test -f 'eigenui.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'eigenui.m' '(file already exists)' else $echo 'x -' extracting 'eigenui.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'eigenui.m' && function eigenui() %load eigenui global EIG X EIG.fig = figure('Color',[0.8 0.8 0.8], ... X 'Name','Eigenvalues', ... X 'Units','normalized', ... X 'Position',[.01 .05 .98 .9], ... X 'Tag','Fig1'); X EIG.axes = axes('Parent',EIG.fig, ... X 'CameraUpVector',[0 1 0], ... X 'Color',[1 1 1], ... X 'FontSize',10, ... X 'FontWeight','bold', ... X 'LineWidth',1, ... X 'NextPlot','add', ... X 'Position',[0.05 0.05 0.64 0.93], ... X 'Tag','Axes1', ... X 'XColor',[0 0 0], ... X 'YColor',[0 0 0], ... X 'ZColor',[0 0 0]); X EIG.pop_fcn = uicontrol('Parent',EIG.fig, ... X 'Callback','eigen_update(''Function'');',... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.70 0.91 0.30 0.08], ... X 'String',{'ArnoldiC'}, ... X 'Style','popupmenu', ... X 'Tag','PopupMenu1', ... X 'Value',1); X EIG.edit_dim = uicontrol('Parent',EIG.fig, ... X 'Callback','eigen_update(''Dim'');',... X 'Units','normalized', ... X 'BackgroundColor',[1 1 1], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.90 0.84 0.10 0.08], ... X 'String','40', ... X 'Style','edit', ... X 'Tag','EditText1'); EIG.str_dim = uicontrol('Parent',EIG.fig, ... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.70 0.82 0.20 0.08], ... X 'String','Dimension', ... X 'Style','text', ... X 'Tag','StaticText1'); X EIG.pop_sym = uicontrol('Parent',EIG.fig, ... X 'Callback','eigen_update(''Sym'');',... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.70 0.73 0.30 0.08], ... X 'String',{'nonsymmetric' 'symmetric' 'symmetric w/f'}, ... X 'Style','popupmenu', ... X 'Tag','PopupMenu2', ... X 'Value',1); X EIG.edit_steps = uicontrol('Parent',EIG.fig, ... X 'Callback','eigen_update(''Steps'');',... X 'Units','normalized', ... X 'BackgroundColor',[1 1 1], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.90 0.66 0.10 0.08], ... X 'String','40', ... X 'Style','edit', ... X 'Tag','EditText2'); EIG.str_steps = uicontrol('Parent',EIG.fig, ... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.70 0.64 0.20 0.08], ... X 'String','Steps', ... X 'Style','text', ... X 'Tag','StaticText2'); X EIG.edit_ortho = uicontrol('Parent',EIG.fig, ... X 'Units','normalized', ... X 'BackgroundColor',[1 1 1], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.80 0.57 0.20 0.08], ... X 'String','', ... X 'Style','edit', ... X 'Tag','EditText3'); EIG.str_ortho = uicontrol('Parent',EIG.fig, ... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.70 0.55 0.10 0.08], ... X 'String','Ortho', ... X 'Style','text', ... X 'Tag','StaticText3'); X EIG.edit_res = uicontrol('Parent',EIG.fig, ... X 'Units','normalized', ... X 'BackgroundColor',[1 1 1], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.80 0.48 0.20 0.08], ... X 'String','', ... X 'Style','edit', ... X 'Tag','EditText4'); EIG.str_res = uicontrol('Parent',EIG.fig, ... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.70 0.46 0.10 0.08], ... X 'String','Res', ... X 'Style','text', ... X 'Tag','StaticText4'); X EIG.push_calc = uicontrol('Parent',EIG.fig, ... X 'Callback','eigen_update(''Calc'');',... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'ForegroundColor',[1 0 0], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.70 0.37 0.30 0.08], ... X 'String','Calculate', ... X 'Tag','PushButton0' ... X ); X EIG.push_ploteig = uicontrol('Parent',EIG.fig, ... X 'Callback','eigen_update(''PlotEigenValues'');',... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.70 0.28 0.30 0.08], ... X 'String','Plot Eigenvalues', ... X 'Tag','PushButton1' ... X ); X EIG.push_step = uicontrol('Parent',EIG.fig, ... X 'Callback','eigen_update(''Step'');',... X 'Interruptible','off',... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.70 0.19 0.30 0.08], ... X 'String','Step', ... X 'Tag','PushButton2' ... X ); X EIG.push_continue = uicontrol('Parent',EIG.fig, ... X 'Callback','eigen_update(''Continue'');',... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.70 0.10 0.30 0.08], ... X 'String','Continue', ... X 'Tag','PushButton3' ... X ); X EIG.push_exit = uicontrol('Parent',EIG.fig, ... X 'Callback','eigen_update(''Exit'');',... X 'Units','normalized', ... X 'BackgroundColor',[0.8 0.8 0.8], ... X 'FontSize',14, ... X 'FontWeight','bold', ... X 'Position',[0.70 0.01 0.30 0.08], ... X 'String','Exit', ... X 'Tag','PushButton4' ... X ); X %------------------------------------------- eigen_update('Init'); X SHAR_EOF $shar_touch -am 02231621100 'eigenui.m' && chmod 0664 'eigenui.m' || $echo 'restore of' 'eigenui.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'eigenui.m:' 'MD5 check failed' 14588e665eca7c67ac911ff42d3b1592 eigenui.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'eigenui.m'`" test 5004 -eq "$shar_count" || $echo 'eigenui.m:' 'original size' '5004,' 'current size' "$shar_count!" fi fi # ============= evalB.m ============== if test -f 'evalB.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'evalB.m' '(file already exists)' else $echo 'x -' extracting 'evalB.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'evalB.m' && X X function [j, a,b,c,d, ierr] = evalB(t,x); % % Evaluates Bsplines at value t in the interval % [x(1) : x(n)] with equally spaced points x(j) % X n = length(x); X ierr = 0; X X if (t < x(1) | t > x(n)), ierr = 1; break; end X X h = x(2) - x(1); X X j = 1; X while (t > x(j)), j = j+1; end X if(j > 1), j = j - 1; end % % t is in [x(j),x(j+1)) % % evaluate the four B_k(t) values that % are nonzero on this interval % X X z = (x(j+1) - t)/h; X a = z^3; X b = 1 + 3*z + 3*z^2 - 3*z^3; X X z = (t - x(j))/h; X c = 1 + 3*z + 3*z^2 - 3*z^3; X d = z^3; X SHAR_EOF $shar_touch -am 09240953103 'evalB.m' && chmod 0644 'evalB.m' || $echo 'restore of' 'evalB.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'evalB.m:' 'MD5 check failed' 4a3a68435bbcdd2be02c53a8db2d8095 evalB.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'evalB.m'`" test 619 -eq "$shar_count" || $echo 'evalB.m:' 'original size' '619,' 'current size' "$shar_count!" fi fi # ============= krylov.m ============== if test -f 'krylov.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'krylov.m' '(file already exists)' else $echo 'x -' extracting 'krylov.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'krylov.m' && % % set up matrix % X n = 100; k = 10; X X [Q,R] = qr(randn(n)); X t = max(abs(diag(R))); X R(n,n) = t*n; X A = Q*R*Q'; % % X v = randn(n,1); v = v/norm(v); X K = zeros(n,k); X K(:,1) = v; X resid = zeros(k,1); % % Form the Krylov matrix from the power sequence % X for j = 1:k, X w = A*v; X v = w/norm(w); X K(:,j) = v; X end X format long e X [V,R] = qr(K,0); % % compute the projected matrix H and th residual F % X H = V'*(A*V); X F = A*V - V*H; X format short e X X Hval = H X X for j = 1:k, X resid(j) = norm(F(:,j)); X end % % plot norms of residual columns % X semilogy(resid,'*') X hold X semilogy((resid)) X hold X % % contour and mesh plots of log abs H % X rev = k:-1:1; X figure X contour(log(abs(H(rev,:)))) X colorbar X mesh(log(abs(H(rev,:)))) X colorbar SHAR_EOF $shar_touch -am 02101456100 'krylov.m' && chmod 0664 'krylov.m' || $echo 'restore of' 'krylov.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'krylov.m:' 'MD5 check failed' c462fce4ce9c11957797c3416ae85940 krylov.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'krylov.m'`" test 768 -eq "$shar_count" || $echo 'krylov.m:' 'original size' '768,' 'current size' "$shar_count!" fi fi # ============= mmread.m ============== if test -f 'mmread.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'mmread.m' '(file already exists)' else $echo 'x -' extracting 'mmread.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'mmread.m' && function [A,rows,cols,entries,rep,field,symm] = mmread(filename) % % function [A] = mmread(filename) % % function [A,rows,cols,entries,rep,field,symm] = mmread(filename) % % Reads the contents of the Matrix Market file 'filename' % into the matrix 'A'. 'A' will be either sparse or full, % depending on the Matrix Market format indicated by % 'coordinate' (coordinate sparse storage), or % 'array' (dense array storage). The data will be duplicated % as appropriate if symmetry is indicated in the header. % % Optionally, size information about the matrix can be % obtained by using the return values rows, cols, and % entries, where entries is the number of nonzero entries % in the final matrix. Type information can also be retrieved % using the optional return values rep (representation), field, % and symm (symmetry). % X mmfile = fopen(filename,'r'); if ( mmfile == -1 ) X disp(filename); X error('File not found'); end; X header = fgets(mmfile); if (header == -1 ) X error('Empty file.') end X % NOTE: If using a version of Matlab for which strtok is not % defined, substitute 'gettok' for 'strtok' in the % following lines, and download gettok.m from the % Matrix Market site. [head0,header] = strtok(header); % see note above [head1,header] = strtok(header); [rep,header] = strtok(header); [field,header] = strtok(header); [symm,header] = strtok(header); head1 = lower(head1); rep = lower(rep); field = lower(field); symm = lower(symm); if ( length(symm) == 0 ) X disp(['Not enough words in header line of file ',filename]) X disp('Recognized format: ') X disp('%%MatrixMarket matrix representation field symmetry') X error('Check header line.') end if ( ~ strcmp(head0,'%%MatrixMarket') ) X error('Not a valid MatrixMarket header.') end if ( ~ strcmp(head1,'matrix') ) X disp(['This seems to be a MatrixMarket ',head1,' file.']); X disp('This function only knows how to read MatrixMarket matrix files.'); X disp(' '); X error(' '); end X % Read through comments, ignoring them X commentline = fgets(mmfile); while length(commentline) > 0 & commentline(1) == '%', X commentline = fgets(mmfile); end X % Read size information, then branch according to % sparse or dense format X if ( strcmp(rep,'coordinate')) % read matrix given in sparse X % coordinate matrix format X X [sizeinfo,count] = sscanf(commentline,'%d%d%d'); X while ( count == 0 ) X commentline = fgets(mmfile); X if (commentline == -1 ) X error('End-of-file reached before size information was found.') X end X [sizeinfo,count] = sscanf(commentline,'%d%d%d'); X if ( count > 0 & count ~= 3 ) X error('Invalid size specification line.') X end X end X rows = sizeinfo(1); X cols = sizeinfo(2); X entries = sizeinfo(3); X X if ( strcmp(field,'real') ) % real valued entries: X X [T,count] = fscanf(mmfile,'%f',3); X T = [T; fscanf(mmfile,'%f')]; X if ( size(T) ~= 3*entries ) X message = ... X str2mat('Data file does not contain expected amount of data.',... X 'Check that number of data lines matches nonzero count.'); X disp(message); X error('Invalid data.'); X end X T = reshape(T,3,entries)'; X A = sparse(T(:,1), T(:,2), T(:,3), rows , cols); X X elseif ( strcmp(field,'complex')) % complex valued entries: X X T = fscanf(mmfile,'%f',4); X T = [T; fscanf(mmfile,'%f')]; X if ( size(T) ~= 4*entries ) X message = ... X str2mat('Data file does not contain expected amount of data.',... X 'Check that number of data lines matches nonzero count.'); X disp(message); X error('Invalid data.'); X end X T = reshape(T,4,entries)'; X A = sparse(T(:,1), T(:,2), T(:,3) + T(:,4)*sqrt(-1), rows , cols); X X elseif ( strcmp(field,'pattern')) % pattern matrix (no values given): X X T = fscanf(mmfile,'%f',2); X T = [T; fscanf(mmfile,'%f')]; X if ( size(T) ~= 2*entries ) X message = ... X str2mat('Data file does not contain expected amount of data.',... X 'Check that number of data lines matches nonzero count.'); X disp(message); X error('Invalid data.'); X end X T = reshape(T,2,entries)'; X A = sparse(T(:,1), T(:,2), ones(entries,1) , rows , cols); X X end X elseif ( strcmp(rep,'array') ) % read matrix given in dense X % array (column major) format X X [sizeinfo,count] = sscanf(commentline,'%d%d'); X while ( count == 0 ) X commentline = fgets(mmfile); X if (commentline == -1 ) X error('End-of-file reached before size information was found.') X end X [sizeinfo,count] = sscanf(commentline,'%d%d'); X if ( count > 0 & count ~= 2 ) X error('Invalid size specification line.') X end X end X rows = sizeinfo(1); X cols = sizeinfo(2); X entries = rows*cols; X if ( strcmp(field,'real') ) % real valued entries: X A = fscanf(mmfile,'%f',1); X A = [A; fscanf(mmfile,'%f')]; X if ( strcmp(symm,'symmetric') | strcmp(symm,'hermitian') | strcmp(symm,'skew-symmetric') ) X for j=1:cols-1, X currenti = j*rows; X A = [A(1:currenti); zeros(j,1);A(currenti+1:length(A))]; X end X elseif ( ~ strcmp(symm,'general') ) X disp('Unrecognized symmetry') X disp(symm) X disp('Recognized choices:') X disp(' symmetric') X disp(' hermitian') X disp(' skew-symmetric') X disp(' general') X error('Check symmetry specification in header.'); X end X A = reshape(A,rows,cols); X elseif ( strcmp(field,'complex')) % complx valued entries: X tmpr = fscanf(mmfile,'%f',1); X tmpi = fscanf(mmfile,'%f',1); X A = tmpr+tmpi*i; X for j=1:entries-1 X tmpr = fscanf(mmfile,'%f',1); X tmpi = fscanf(mmfile,'%f',1); X A = [A; tmpr + tmpi*i]; X end X if ( strcmp(symm,'symmetric') | strcmp(symm,'hermitian') | strcmp(symm,'skew-symmetric') ) X for j=1:cols-1, X currenti = j*rows; X A = [A(1:currenti); zeros(j,1);A(currenti+1:length(A))]; X end X elseif ( ~ strcmp(symm,'general') ) X disp('Unrecognized symmetry') X disp(symm) X disp('Recognized choices:') X disp(' symmetric') X disp(' hermitian') X disp(' skew-symmetric') X disp(' general') X error('Check symmetry specification in header.'); X end X A = reshape(A,rows,cols); X elseif ( strcmp(field,'pattern')) % pattern (makes no sense for dense) X disp('Matrix type:',field) X error('Pattern matrix type invalid for array storage format.'); X else % Unknown matrix type X disp('Matrix type:',field) X error('Invalid matrix type specification. Check header against MM documentation.'); X end end X % % If symmetric, skew-symmetric or Hermitian, duplicate lower % triangular part and modify entries as appropriate: % X if ( strcmp(symm,'symmetric') ) X A = A + A.' - diag(diag(A)); X entries = nnz(A); elseif ( strcmp(symm,'hermitian') ) X A = A + A' - diag(diag(A)); X entries = nnz(A); elseif ( strcmp(symm,'skew-symmetric') ) X A = A - A'; X entries = nnz(A); end X fclose(mmfile); % Done. X SHAR_EOF $shar_touch -am 12050854101 'mmread.m' && chmod 0666 'mmread.m' || $echo 'restore of' 'mmread.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'mmread.m:' 'MD5 check failed' ed540e073957b876c5f04dc17d2f034c mmread.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'mmread.m'`" test 7224 -eq "$shar_count" || $echo 'mmread.m:' 'original size' '7224,' 'current size' "$shar_count!" fi fi # ============= orthdemo.m ============== if test -f 'orthdemo.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'orthdemo.m' '(file already exists)' else $echo 'x -' extracting 'orthdemo.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'orthdemo.m' && X [X,D] = eig(A); %norm(imag(X)) [d,ii] = sort(-abs(diag(D))); X v = X(:,ii(1:4))*ones(4,1); X driveA2 pause X v = v + sqrt(eps)*randn(n,1)*norm(v); X driveA2 pause X driveAC X SHAR_EOF $shar_touch -am 02101456100 'orthdemo.m' && chmod 0664 'orthdemo.m' || $echo 'restore of' 'orthdemo.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'orthdemo.m:' 'MD5 check failed' 13c7ce0ed5ab2462d2606f4fd62bf41e orthdemo.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'orthdemo.m'`" test 170 -eq "$shar_count" || $echo 'orthdemo.m:' 'original size' '170,' 'current size' "$shar_count!" fi fi # ============= polyR.m ============== if test -f 'polyR.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'polyR.m' '(file already exists)' else $echo 'x -' extracting 'polyR.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'polyR.m' && % % % Demonstrates polynomial restart action on a vector with respect % to a symmetric positive definite matrix A % % Displays Chebyshev filter polynomial T constructed to damp unwanted % coefficients and then displays how expansion coefficients decay % % as v <- T(A)v is done over and over % % D.C. Sorensen % 8 Feb 01 % X X n = 20; X A = randn(n); X A = A'*A; X X [Q,D] = eig(A); X [d,ii] = sort(diag(D)); X X k = 5; % % Uncomment to make a gap between first k eigenvalues and the rest % X d(k+1:n) = d(k+1:n) + 10; X X Q = Q(:,ii); D = diag(d); X Diff = norm(A - Q*D*Q') X X % % Construct an interval [a,b] containing unwanted eigenvalues d(k+1:n) % X X a = d(k) + .99*(d(k+1) - d(k)); b = d(n)*(1.1); X X c = (a + b)/2; X r = (b - a)/2; X % % Construct the roots of Chebyshev polynomial T on [a,b] % X p = 10; X t = zeros(p,1); X tau = pi/(p+1); X for j = 1:p, X t(j) = c + r*cos(j*tau); X end % % get coefficients and normalize so T(0) = 1 % X g = poly(t); X g = g/g(p+1); % % Plot the polynomial and the eigenvalues of A % X X h = b/(201); X x = zeros(200); X for j = 1:200, x(j) = j*h; end X pval = polyval(g,x); X plot(x,pval); X hold X plot(d(1:k), zeros(k,1), 'g*') X plot(d(k+1:n), zeros(n-k,1), 'r*') X hold X pause X X figure % % Repeatedly apply v <- T(A)v and normalize % Display decaying coefficients of the unwanted eigenvectors % in the expansion of v in the eigenvector basis Q % % Note: You would not construct and apply the Chebychev polynomials in % this way in a production setting. You should use the 3-term % recursion to apply T(A) to a vector, where the recursion has % been adjusted to the specified interval [a,b]. % X X y = randn(n,1); X v = Q*y; X v = v/norm(v); X semilogy(abs(y),'o') X hold X while(norm(y(k+1:n)) > sqrt(eps)), X v = polyvalm(g,A)*v; X v = v/norm(v); X y = Q'*v; X semilogy(abs(y),'o') X end X semilogy(abs(y),'ro') X hold X X X X X SHAR_EOF $shar_touch -am 02081300101 'polyR.m' && chmod 0664 'polyR.m' || $echo 'restore of' 'polyR.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'polyR.m:' 'MD5 check failed' 8ed45fa0d338caeff04f5b72c7c235e4 polyR.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'polyR.m'`" test 2135 -eq "$shar_count" || $echo 'polyR.m:' 'original size' '2135,' 'current size' "$shar_count!" fi fi # ============= pow.m ============== if test -f 'pow.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'pow.m' '(file already exists)' else $echo 'x -' extracting 'pow.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'pow.m' && X v = randn(n,1); v = v/norm(v); X K = zeros(n,k); X K(:,1) = v; X resid = zeros(k,1); X for j = 1:k, X w = A*v; X v = w/norm(w); X K(:,j) = v; X end X format long e X [V,R] = qr(K); X H = V'*(A*V); X F = A*V - V*H; X format short e X Hval = H X X for j = 1:k, X resid(j) = norm(F(:,j)); X end X plot(resid,'*') X hold X plot(resid) X hold SHAR_EOF $shar_touch -am 02101456100 'pow.m' && chmod 0664 'pow.m' || $echo 'restore of' 'pow.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'pow.m:' 'MD5 check failed' 1a21551502e0014463d4ff26092f18c6 pow.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'pow.m'`" test 331 -eq "$shar_count" || $echo 'pow.m:' 'original size' '331,' 'current size' "$shar_count!" fi fi # ============= pow_inf.m ============== if test -f 'pow_inf.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'pow_inf.m' '(file already exists)' else $echo 'x -' extracting 'pow_inf.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'pow_inf.m' && X function pow_inf(A); X X [n,n] = size(A); X nval = n X s = eig(A); X [d,ii] = max(abs(s)); X lam0 = s(ii); X t = zeros(20,1); X s = lam0*ones(20,1); X X v = randn(n,1); v = v/norm(v,'inf'); X for j = 1:20, X w = A*v; X [lam, i ] = max(abs(w)); X lam = w(i); X v = w/lam; X t(j) = lam; X end X format long e X lambda_vals_INF = [t abs(t-s)] SHAR_EOF $shar_touch -am 02101456100 'pow_inf.m' && chmod 0664 'pow_inf.m' || $echo 'restore of' 'pow_inf.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'pow_inf.m:' 'MD5 check failed' 0e4e0626ce609bd32647c764174adfd7 pow_inf.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'pow_inf.m'`" test 338 -eq "$shar_count" || $echo 'pow_inf.m:' 'original size' '338,' 'current size' "$shar_count!" fi fi # ============= pow_rq.m ============== if test -f 'pow_rq.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'pow_rq.m' '(file already exists)' else $echo 'x -' extracting 'pow_rq.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'pow_rq.m' && X function pow_rq(A); X X [n,n] = size(A); X s = eig(A); X [d,ii] = max(abs(s)); X lam0 = s(ii); X t = zeros(20,1); X s = lam0*ones(20,1); X X v = randn(n,1); v = v/norm(v); X for j = 1:20, X w = A*v; X lam = v'*w; X v = w/norm(w); X t(j) = lam; X end X format long e X lambda_vals_RQ = [t abs(t-s)] SHAR_EOF $shar_touch -am 02101456100 'pow_rq.m' && chmod 0664 'pow_rq.m' || $echo 'restore of' 'pow_rq.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'pow_rq.m:' 'MD5 check failed' a011c555d993f82d64939f3e1b638b04 pow_rq.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'pow_rq.m'`" test 295 -eq "$shar_count" || $echo 'pow_rq.m:' 'original size' '295,' 'current size' "$shar_count!" fi fi # ============= powersR.m ============== if test -f 'powersR.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'powersR.m' '(file already exists)' else $echo 'x -' extracting 'powersR.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'powersR.m' && function powersR(rho,tau); % % This function illustrates the behaviour of norm(R^k) % for k = 1:300. % % Inputs: % rho - A scalar specifying the intial spectral radius. % It is assumed that 0 < rho < 1. % % tau - A non-negative number used to control the degree of % non-normality. % % An upper triangular matrix R = rho*D + tau*U is generated % and then norm(R^k,1) is recorded for each power k and the % sequence is plotted. The diagonal matrix D is generated % to have spectral radius 1, and U is generated to be % strictly upper triangular. So the spectral radius of R for % the three plots will be rho, rho/10, rho/100. % % % Recommended Values rho = .1, tau = .001, .1, 1, 10 % % %---------------------------------------------------------- % % D.C. Sorensen % 10 Oct 03 X X nplots = 3; X n = 100; m = 300; X X D = diag(randn(n,1)); X [Q,R] = qr(randn(n)); X X for j = 1:n, R(j,j) = 0; end X U = R; % % U is strict upper triangular. % X X dmax = max(abs(diag(D))); X D = D/dmax; % % D is diagonal with spec. radius 1 % X X X t = zeros(m,nplots); X for j = 1:nplots, X R = rho*D + tau*U; % % R is upper triangular with spec. radius rho % X Rk = eye(n); X for k = 1:m, X Rk = R*Rk; X t(k,j) = norm(Rk,1); X end X SpecRad = rho X rho = rho/10; X end % % Plot the three convergence histories % to compare % X figure(1), clf X semilogy([1:m],t(:,1), 'r-', 'linewidth',2) X ylabel('power norm, ||R^k||', 'fontsize', 18) X xlabel('iteration k', 'fontsize', 18) X hold X disp('Strike a Key') X pause X semilogy([1:m],t(:,2), 'k-', 'linewidth',2) X ylabel('power norm, ||R^k||', 'fontsize', 18) X xlabel('iteration k', 'fontsize', 18) X disp('Strike a Key') X pause X semilogy([1:m],t(:,3), 'b-', 'linewidth',2) X xlabel('iteration k', 'fontsize', 18) X ylabel('power norm, ||R^k||', 'fontsize', 18) X legend('\rho_1 = rho ','\rho_2 = rho/10','\rho_3 = rho/100') X hold X X X X X X X SHAR_EOF $shar_touch -am 10101449103 'powersR.m' && chmod 0644 'powersR.m' || $echo 'restore of' 'powersR.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'powersR.m:' 'MD5 check failed' 246c7159a634e112bf9918a74fa44ff1 powersR.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'powersR.m'`" test 2056 -eq "$shar_count" || $echo 'powersR.m:' 'original size' '2056,' 'current size' "$shar_count!" fi fi # ============= qriter.m ============== if test -f 'qriter.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'qriter.m' '(file already exists)' else $echo 'x -' extracting 'qriter.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'qriter.m' && function [V,H] = qriter(A); % % Performs QR iteration in real arithmetic % % Input: A -- a real n by n matrix % % Output: V -- a real n by n orthogonal matrix % % H -- a real n by n upper triangular matrix % % AV = VH % % D.C. Sorensen % 2 Mar 2000 % X [n,n] = size(A); % % Reduce A to Hessenberg form AV = VH % X v = randn(n,1); X [V,H] = ArnoldiC(A,n,v); X X k = n; X while(k > 2) % % Select the Wilkinson shift % X mu = eig(H(k-1:k,k-1:k)); X if (abs(mu(1) - H(k,k)) < abs(mu(2) - H(k,k)) ), X X [V,H] = qrstep(V,H,mu(1),1,k); X X else X X [V,H] = qrstep(V,H,mu(2),1,k); X X end % % Deflate: Set subdiagonal to zero if it is small compared to % the two adjacent diagonal elements X X if (abs(H(k-1,k-2)) < 10*eps*max(abs(H(k-2,k-2)),abs(H(k-1,k-1)))), X X H(k-1,k-2) = 0; X k = k - 2; X X elseif (abs(H(k,k-1)) < 10*eps*max(abs(H(k-1,k-1)),abs(H(k,k)))), X X H(k,k-1) = 0; X k = k - 1; X X end X % % ---------------Graphics Display of Convergence-------------------- % X X x = ones(n*(n+1)/2 , 1); y = x; X phd = zeros(n-1,1); X X m = 16; X colors = bluemap(m); X axis([0 n+1 0 n+1]); X x = [-1 -1 n+2 n+2]; X y = [-1 n+2 -1 n+2]; X ph1 = plot(x,y,'o'); X axis tight; X set(ph1,'MarkeredgeColor',[1 1 1]) X hold on; X indx = 1; X Hmax = max(abs(diag(H,-1))); X for j = 1:n, X for i = 1:j, X ii = n-i+1; jj = n-j+1; X x(indx) = ii; y(indx) = j; X indx = indx + 1; X end X jj = min(n-1,j); X i = min(n,jj+1); X gval = max(eps, abs(H(i,jj)/Hmax)); X gval = -log(gval); X gval = max(1,ceil(gval)); X gval = min(16,gval); X gval = 16 - gval + 1; X ii = n-i+1; jj = j; % ii = i; jj = n - jj -1 ; % phd(j) = plot(ii,jj,'o'); X phd(j) = plot(jj,ii,'o'); X axis tight; X set(phd(j),'Color',colors(m,:)) X set(phd(j),'MarkerSize',9) X set(phd(j),'MarkeredgeColor',[1 1 1]) % Hdiag = abs(H(i,jj)/Hmax) % Hdiag = abs(diag(H,-1)/Hmax) % Hdiag2= Hdiag*colors(1,:) X if(abs(H(i,j)) > 0), X set(phd(j),'MarkerFaceColor', colors(gval,:)); X end X end X ph = plot(x,y,'o'); X axis tight; X set(ph,'Color',colors(m,:)) X set(ph,'MarkerSize',9) X set(ph,'MarkeredgeColor',[1 1 1]) X set(ph,'MarkerFaceColor',colors(m,:)) X pause(.1) % pause X hold off X axis tight; X end X X X X SHAR_EOF $shar_touch -am 03021402100 'qriter.m' && chmod 0664 'qriter.m' || $echo 'restore of' 'qriter.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'qriter.m:' 'MD5 check failed' f04137eeaf50d8f35be16ea16390b356 qriter.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'qriter.m'`" test 2704 -eq "$shar_count" || $echo 'qriter.m:' 'original size' '2704,' 'current size' "$shar_count!" fi fi # ============= qriterR.m ============== if test -f 'qriterR.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'qriterR.m' '(file already exists)' else $echo 'x -' extracting 'qriterR.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'qriterR.m' && function [V,H] = qriterR(A); % % Performs QR iteration in real arithmetic % % Input: A -- a real n by n matrix % % Output: V -- a real n by n orthogonal matrix % % H -- a real n by n upper triangular matrix % % AV = VH % % leading k by k block of H has eigenvalues of % largest real part % % % D.C. Sorensen % 2 Mar 2000 % X format short e X [n1,n1] = size(A); % % Reduce A to Hessenberg form AV = VH % X v = randn(n1,1); X [V,H] = ArnoldiC(A,n1,v); X t1 = eig(A); X [t2,ii] = sort(-real(t1)); X t1 = t1(ii); X X m = 20; k = 6; X t1 = t1(1:m); % % Find eigenvalues of leading m by m block % and sort acording to largest real part % X mu = eig(H(1:m,1:m)); X [t,ii] = sort(-real(mu)); X mu = mu(ii); X X iter = 1; ynorm = 1; X while( iter < 25 & ynorm > 1000*eps) X X iter = iter + 1; X m = 20; k = 6; X j = m; X while(j > k) % % Apply unwanted eigenvalues as shifts with QR steps % X [V,H] = qrstep(V,H,mu(j),1,n1); % % Decrease j by 2 if m_j imaginary and by 1 if real % X if (abs(imag(mu(j))) > 0), X j = j-2; X else X j = j-1; X end X end X % % Display some convergence information % X Hvals = H(1:j,1:j); X [Y,D] = eig(Hvals); X y = abs(Y(j,:)*H(j+1,j)); X y = y' X ynorm = norm(y(1:j-1)); X Niter = iter X % % Compute eigenvalues of leading m by m block X X Hvals = H(1:m,1:m); X [Y,D] = eig(Hvals); X mu = diag(D); % % Sort the eigenvalues according to largest real part % X [t,ii] = sort(-real(mu)); X mu = mu(ii); X % % Display some convergence information % X X yvec = abs(Y(m,:)*H(m+1,m)); X yvec = yvec(ii)'; X y_diff_mu_lam = [yvec abs(abs(mu) - abs(t1)) abs(mu) abs(t1)] X X % % ---------------------Graphics to Display Convergence----------------- % X n = m; X x = ones(n*(n+1)/2 , 1); y = x; X phd = zeros(n-1,1); X phy = zeros(n-1,1); X X m = 16; X colors = bluemap(m); X x = [-1 -1 n+2 n+2]; X y = [-1 n+2 -1 n+2]; X X ph1 = plot(x,y,'o'); X axis tight; X set(ph1,'MarkeredgeColor',[1 1 1]) X hold on; X indx = 1; X Hmax = max(abs(diag(Hvals,-1))); X for j = 1:n, X for i = 1:j, X ii = n-i+1; jj = n-j+1; X x(indx) = ii; y(indx) = j; X indx = indx + 1; X end X jj = min(n-1,j); X i = min(n,jj+1); X gval = max(eps, abs(Hvals(i,jj)/Hmax)); X gval = -log(gval); X gval = max(1,ceil(gval)); X gval = min(16,gval); X gval = 16 - gval + 1; X X yval = max(eps, abs(yvec(jj))); X yval = -log(yval); X yval = max(1,ceil(yval)); X yval = min(16,yval); X yval = 16 - yval + 1; X X ii = n-i+1; jj = j; X phy(j) = plot(jj,0,'o'); X phd(j) = plot(jj,ii,'o'); X axis tight; X set(phd(j),'Color',colors(m,:)) X set(phd(j),'MarkerSize',9) X set(phd(j),'MarkeredgeColor',[1 1 1]) X X set(phy(j),'Color',colors(m,:)) X set(phy(j),'MarkerSize',9) X set(phy(j),'MarkeredgeColor',[1 1 1]) X X if(abs(Hvals(i,j)) > 0), X set(phd(j),'MarkerFaceColor', colors(gval,:)); X end X if(abs(yvec(j)) > 100*eps), X set(phy(j),'MarkerFaceColor', colors(yval,:)); X end X end X ph = plot(x,y,'o'); X axis tight; X set(ph,'Color',colors(m,:)) X set(ph,'MarkerSize',9) X set(ph,'MarkeredgeColor',[1 1 1]) X set(ph,'MarkerFaceColor',colors(m,:)) X x1 = [ .7 k+.3]; X y1 = [n-k-.3 n+.3]; X x2 = [ .7 .7]; X x3 = [ k+.3 k+.3]; X y2 = [n-k-.3 n-k-.3]; X y3 = [n+.3 n+.3]; X plot(x1,y2,'m--',x1,y3,'m--',x2,y1,'m--',x3,y1,'m--'); X axis tight; X pause(.2) X hold off X end X X X X X X SHAR_EOF $shar_touch -am 03021402100 'qriterR.m' && chmod 0664 'qriterR.m' || $echo 'restore of' 'qriterR.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'qriterR.m:' 'MD5 check failed' 2f86f67a36380319d5066da8b3f84017 qriterR.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'qriterR.m'`" test 4331 -eq "$shar_count" || $echo 'qriterR.m:' 'original size' '4331,' 'current size' "$shar_count!" fi fi # ============= qriterRS.m ============== if test -f 'qriterRS.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'qriterRS.m' '(file already exists)' else $echo 'x -' extracting 'qriterRS.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'qriterRS.m' && function [V,H] = qriterRS(A); % % % % Performs QR iteration in real arithmetic % % Shifts are selected from "unwanted Ritz values % from the leading m by m block. % % Input: A -- a real n by n symmetric matrix % % Output: V -- a real n by n orthogonal matrix % % H -- a real n by n tridiagonal matrix % % AV = VH % % The leading k by k block contains the algebraically % largest eigenvalues of A. % % Here m = 20, k = 6 % % Note: In real life, you would take advantage of symmetry. % This code does not do that. % % D.C. Sorensen % 2 Mar 2000 % % X format short e X [n1,n1] = size(A); X v = randn(n1,1); % % Reduce A to tridiagonal form with full Arnoldi process % X X [V,H] = ArnoldiC(A,n1,v); X t1 = eig(A); X [t2,ii] = sort(-real(t1)); X t1 = t1(ii); X X m = 20; k = 6; X t1 = t1(1:m); % % Find eigenvalues of leading m by m block % and sort acording to largest real part % X mu = eig(H(1:m,1:m)); X [t,ii] = sort(-real(mu)); X mu = mu(ii); X iter = 1; ynorm = 1; X while( iter < 25 & ynorm > 1000*eps) X iter = iter + 1; X m = 20; k = 6; X j = m; X while(j > k), % % Apply unwanted eigenvalues as shifts with QR steps % X X [V,H] = qrstep(V,H,mu(j),1,n1); X j = j-1; X end X % % Display some convergence information % X Hvals = H(1:j,1:j); X [Y,D] = eig(Hvals); X y = abs(Y(j,:)*H(j+1,j)); X y = y' X ynorm = norm(y(1:j-1)); X Niter = iter % % Compute eigenvalues of leading m by m block % X Hvals = H(1:m,1:m); X [Y,D] = eig(Hvals); % % Sort the eigenvalues according to largest real part % X X mu = diag(D); X [t,ii] = sort(-real(mu)); X mu = mu(ii); % % Display some convergence information % X X yvec = abs(Y(m,:)*H(m+1,m)); X yvec = yvec(ii)'; X y_diff_mu_lam = [yvec abs(abs(mu) - abs(t1)) abs(mu) abs(t1)] % % ---------------------Graphics to Display Convergence----------------- % X n = m; X x = ones(n*(n+1)/2 , 1); y = x; X phd = zeros(n-1,1); X phy = zeros(n-1,1); X X m = 16; X colors = bluemap(m); X x = [-1 -1 n+2 n+2]; X y = [-1 n+2 -1 n+2]; X X ph1 = plot(x,y,'o'); X axis tight; X set(ph1,'MarkeredgeColor',[1 1 1]) X hold on; X indx = 1; X Hmax = max(abs(diag(Hvals,-1))); X X x = [1:n]; X y = [n:-1:1]; X X for j = 1:n-1, X jj = min(n-1,j); X i = min(n,jj+1); X gval = max(eps, abs(Hvals(i,jj)/Hmax)); X gval = -log(gval); X gval = max(1,ceil(gval)); X gval = min(16,gval); X gval = 16 - gval + 1; X X yval = max(eps, abs(yvec(jj))); X yval = -log(yval); X yval = max(1,ceil(yval)); X yval = min(16,yval); X yval = 16 - yval + 1; X X ii = n-i+1; jj = j; X phy(j) = plot(jj,0,'o'); X phd(j) = plot([jj,jj+1],[ii, ii+1],'o'); X axis tight; X set(phd(j),'Color',colors(m,:)) X set(phd(j),'MarkerSize',9) X set(phd(j),'MarkeredgeColor',[1 1 1]) X X set(phy(j),'Color',colors(m,:)) X set(phy(j),'MarkerSize',9) X set(phy(j),'MarkeredgeColor',[1 1 1]) X X if(abs(Hvals(i,j)) > 0), X set(phd(j),'MarkerFaceColor', colors(gval,:)); X end X if(abs(yvec(j)) > 100*eps), X set(phy(j),'MarkerFaceColor', colors(yval,:)); X end X end X x = [x n]; y = [y 0]; X ph = plot(x,y,'o'); X axis tight; X set(ph,'Color',colors(m,:)) X set(ph,'MarkerSize',9) X set(ph,'MarkeredgeColor',[1 1 1]) X set(ph,'MarkerFaceColor',colors(m,:)) X x1 = [ .7 k+.3]; X y1 = [n-k-.3 n+.3]; X x2 = [ .7 .7]; X x3 = [ k+.3 k+.3]; X y2 = [n-k-.3 n-k-.3]; X y3 = [n+.3 n+.3]; X plot(x1,y2,'m--',x1,y3,'m--',x2,y1,'m--',x3,y1,'m--'); X axis tight; X pause(.2) X hold off X end X X X X X X SHAR_EOF $shar_touch -am 03021402100 'qriterRS.m' && chmod 0664 'qriterRS.m' || $echo 'restore of' 'qriterRS.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'qriterRS.m:' 'MD5 check failed' 7ac4e096997a9cbed8f350891dbce463 qriterRS.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'qriterRS.m'`" test 4350 -eq "$shar_count" || $echo 'qriterRS.m:' 'original size' '4350,' 'current size' "$shar_count!" fi fi # ============= qriterS.m ============== if test -f 'qriterS.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'qriterS.m' '(file already exists)' else $echo 'x -' extracting 'qriterS.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'qriterS.m' && X function [V,H] = qriterS(A); % % % Performs QR iteration in real arithmetic % % Input: A -- a real n by n symmetric matrix % % Output: V -- a real n by n orthogonal matrix % % H -- a real n by n diagonal matrix % % AV = VH % % Note: In real life, you would take advantage of symmetry. % This code does not do that. % % D.C. Sorensen % 2 Mar 2000 % X % X [n,n] = size(A); % % Reduce A to tridiagonal form with full Arnoldi process % X v = randn(n,1); X [V,H] = ArnoldiC(A,n,v); X X k1 = 1; k2 = n; X while(k2 > 1) X if (k2 == 2), % % Diagonalize the leading 2 by 2 block after all % the rest of the eigenvalues have converged. % X [Q,D] = eig(H(1:2,1:2)); X X H(1:2,:) = Q'*H(1:2,:); X H(:,1:2) = H(:,1:2)*Q; X V(:,1:2) = V(:,1:2)*Q; X H(2,1) = 0; X k2 = k2 - 1; X else % % Select the Wilkinson shift % X mu = eig(H(k2-1:k2,k2-1:k2)); X if (abs(mu(1) - H(k2,k2)) < abs(mu(2) - H(k2,k2)) ), X [V,H] = qrstep(V,H,mu(1),k1,k2); X else X [V,H] = qrstep(V,H,mu(2),k1,k2); X end % % Deflate: Set off-diagonal to zero if it is small compared to % the two adjacent diagonal elements X X if (abs(H(k2,k2-1)) < 10*eps*max(abs(H(k2-1,k2-1)),abs(H(k2,k2)))), X H(k2,k2-1) = 0; X H(k2-1,k2) = 0; X k2 = k2 - 1; X end X end X % % ---------------Graphics Display of Convergence-------------------- % X X X x = ones(n*(n+1)/2 , 1); y = x; X phd = zeros(n-1,1); X X m = 16; X colors = bluemap(m); X axis([0 n+1 0 n+1]); X x = [-1 -1 n+2 n+2]; X y = [-1 n+2 -1 n+2]; X ph1 = plot(x,y,'o'); X axis tight; X set(ph1,'MarkeredgeColor',[1 1 1]) X hold on; X indx = 1; X Hmax = max(abs(diag(H,-1))); X Hmax = max(Hmax,1); X x = [1:n]; X y = [n:-1:1]; X for j = 1:n-1, X jj = min(n-1,j); X i = min(n,jj+1); X gval = max(eps, abs(H(i,jj)/Hmax)); X gval = -log(gval); X gval = max(1,ceil(gval)); X gval = min(16,gval); X gval = 16 - gval + 1; X ii = n-i+1; jj = j; X X phd(j) = plot([jj,jj+1],[ii, ii+1],'o'); X axis tight; X set(phd(j),'Color',colors(m,:)) X set(phd(j),'MarkerSize',9) X set(phd(j),'MarkeredgeColor',[1 1 1]) X X if(abs(H(i,j)) > 0), X set(phd(j),'MarkerFaceColor', colors(gval,:)); X end X end X ph = plot(x,y,'o'); X axis tight; X set(ph,'Color',colors(m,:)) X set(ph,'MarkerSize',9) X set(ph,'MarkeredgeColor',[1 1 1]) X set(ph,'MarkerFaceColor',colors(m,:)) X pause(.1) % pause X hold off X axis tight; X end X X X X SHAR_EOF $shar_touch -am 03021402100 'qriterS.m' && chmod 0664 'qriterS.m' || $echo 'restore of' 'qriterS.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'qriterS.m:' 'MD5 check failed' e89dc3cddd90e5b656dd6ef083d63433 qriterS.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'qriterS.m'`" test 2898 -eq "$shar_count" || $echo 'qriterS.m:' 'original size' '2898,' 'current size' "$shar_count!" fi fi # ============= qrstep.m ============== if test -f 'qrstep.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'qrstep.m' '(file already exists)' else $echo 'x -' extracting 'qrstep.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'qrstep.m' && function [V,H] = qrstep(V,H,mu,k1,k2); % % Input: V -- a real square orthogonal matrix % % H -- a real square upper Hessenberg matrix % % mu -- a real or complex shift % % k1,k2 -- pointers to submatrix in k1:k2 block % % Output: V -- a real orthogonal matrix % % H -- a real square upper Hessenberg matrix % % V <- VQ; H <- Q'HQ; % % Q corresponds to a single real shift or % a double complex shift depending on mu % % D.C. Sorensen % 2 Mar 2000 % X k = k2-k1+1; X kr = k1:k2; X n = length(V(:,1)); X X eta = imag(mu); X X if (abs(eta) > 0), % mu is imaginary -- apply double shift % X xi = real(mu); X [Q,R] = qr((H(kr,kr) - xi*eye(k))^2 + eta^2*eye(k)); X X else % mu is real -- apply single shift % X [Q,R] = qr(H(kr,kr) - mu*eye(k)); X X end X X H(kr,:) = Q'*H(kr,:); X H(:,kr) = H(:,kr)*Q; X V(:,kr) = V(:,kr)*Q; X % % clean up rounding error noise below first subdiagonal % X for j = k1:k2, X H(j+2:n,j) = H(j+2:n,j)*0; X end SHAR_EOF $shar_touch -am 03021402100 'qrstep.m' && chmod 0664 'qrstep.m' || $echo 'restore of' 'qrstep.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'qrstep.m:' 'MD5 check failed' 2c0797a3c0a479b779a32ec97082ded2 qrstep.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'qrstep.m'`" test 1148 -eq "$shar_count" || $echo 'qrstep.m:' 'original size' '1148,' 'current size' "$shar_count!" fi fi # ============= readsher.m ============== if test -f 'readsher.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'readsher.m' '(file already exists)' else $echo 'x -' extracting 'readsher.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'readsher.m' && A = mmread('sherman2.mtx'); b = mmread('sherman2_rhs1.mtx'); X X SHAR_EOF $shar_touch -am 04051446101 'readsher.m' && chmod 0664 'readsher.m' || $echo 'restore of' 'readsher.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'readsher.m:' 'MD5 check failed' c608aa57c1a7f6b9d20938ff55365f6d readsher.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'readsher.m'`" test 64 -eq "$shar_count" || $echo 'readsher.m:' 'original size' '64,' 'current size' "$shar_count!" fi fi # ============= ritzest.m ============== if test -f 'ritzest.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'ritzest.m' '(file already exists)' else $echo 'x -' extracting 'ritzest.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'ritzest.m' && X % X % set up matrix X % X format short e X n = 50; k = 10; X % [Q,R] = qr(randn(n)); % A = Q*R*Q'; X A = randn(n); A = A + A'; X % X % X v = randn(n,1); X X [V,H,f] = ArnoldiC(A,n,v); X X ek = zeros(n,1); ek(n)= 1; X Resid = norm(A*V - V*H - f*ek') X X orthtestV = norm(eye(n) - V'*V) X ortestf = norm(V'*f) X X t = eig(A); X t = sort(abs(t)); X X k1 = 10; X for j = k1:n-1, X [Y,D] = eig(H(1:j,1:j)); X d = diag(D); X [s,ii] = sort(abs(d)); X X ritz = abs(H(j+1,j))*abs(Y(j,ii(j-k1+1:j))); ritz = ritz'; X X t__s__diff__ritz = [ritz abs(t(n-k1+1:n) - s(j-k1+1:j)) t(n-k1+1:n) s(j-k1+1:j)] X X X y = Y(:,ii(j)); x = V(:,1:j)*y; theta = d(ii(j)); X X normRTZ = norm(A*x - x*theta) X pause X end SHAR_EOF $shar_touch -am 02171349100 'ritzest.m' && chmod 0664 'ritzest.m' || $echo 'restore of' 'ritzest.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'ritzest.m:' 'MD5 check failed' 63cb48d83944f196d9dbcd3d1d0d78c4 ritzest.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'ritzest.m'`" test 704 -eq "$shar_count" || $echo 'ritzest.m:' 'original size' '704,' 'current size' "$shar_count!" fi fi # ============= select_shifts.m ============== if test -f 'select_shifts.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'select_shifts.m' '(file already exists)' else $echo 'x -' extracting 'select_shifts.m' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'select_shifts.m' && function [w, q]= select_shifts(H,which); % % Compute the eigenvalues of H and sort according to % % which = 'LR', 'SR', 'LM','LA','SA' % % Input: H -- a real square upper Hessenberg matrix % % which -- a string indicating which eigenvalues are % wanted % % Output: w -- eigenvalues of H sorted according to which % most wanted to least wanted order % % q -- absolute value of last component of eigenvectors % of H in same order as the eigenvalues in w % % X [q,w] = eig(H); X w = diag(w); X X m = length(w); X q = q(m,:); % % select filter mechanism by activating appropriate choice below % X switch which X X case {'LR', 'LA'} % sort for largest real part % X [s,ir] = sort(-real(w)); X % shifts are smallest real part X X case {'SR', 'SA'} % sort for smallest real part % X [s,ir] = sort(real(w)); X % shifts are largest real part X X case {'LM'} % sort for largest magnitude % X [s,ir] = sort(-abs(w)); X X end X w = w(ir); X q = abs(q(ir)); X X SHAR_EOF $shar_touch -am 03031313100 'select_shifts.m' && chmod 0664 'select_shifts.m' || $echo 'restore of' 'select_shifts.m' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'select_shifts.m:' 'MD5 check failed' d8abd46e3a0e76cae2800224ce96ea1a select_shifts.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'select_shifts.m'`" test 1169 -eq "$shar_count" || $echo 'select_shifts.m:' 'original size' '1169,' 'current size' "$shar_count!" fi fi # ============= sparsdem.m ============== if test -f 'sparsdem.m' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'sparsdem.m' '(file already exists)' else $echo 'x -' extracting 'sparsdem.m' '(binary)' sed 's/^X//' << 'SHAR_EOF' | uudecode && begin 600 sparsdem.m M)0T*)2`@5&AI2`](%M=.PT*)0T*)2`@"!M87)K970-"B4-"@T*("`@96-H;R!O;CL-"B`@ M($$@/2!M;7)E860H)V)L8VMH;VQE+FUT>"2`G*0T*("`@9&ES<"@G("`G*0T*("`@ M<&%U64H;BD[#0H@("!D:7-P*"<@02`]($$@*R`Q M,#`P*G-P97EE*&XI.R2A!*0T*("`@=&ET;&4H)TYA='5R86P@ M3W)D97)I;FBA,*0T*("`@3EIH:7-T;W)Y(#T@6TY: M:&ES=&]R>2`L($YO;GIE6UR;6,@;W)D97)I;F6UR8VTH02D[("2`G*0T*("`@9&ES<"@G("`G*0T*("`@<&%U M2`](%M.6FAI2`G*0T*("`@9&ES<"@G("`G*0T*("`@<&%U6UM;60H02D[#0H@("!D:7-P*"<@<"`]('-Y;6UM9"A!*3L@)RD- M"@T*("`@6UM971R:6,@36EN($1E M9W)E92`G+"=&;VYT2A,+"=R)RD-"B`@(&1IF5R;W-,(#T@;FYZ*$PI#0H@("!. M6FAIF5R;W-,73L-"B`@(&1I6UA;60- M"B4@(&UO6UA;60H02D[#0H@("!D:7-P*"<@ M<"`]('-Y;6%M9"A!*3L@)RD-"@T*("`@2`G*0T*("`@9&ES<"@G("`G*0T*("`@<&%U2`](%M.6FAI2`G*0T*("`@9&ES<"@G("`G*0T*("`@<&%U6UR8VT@("`@("!S>6UM;60@("`@("!S>6UA M;60@("`G+"=&;VYT&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'sparsdem.m:' 'MD5 check failed' 622a9d2b51bc017cddc198ef10abcf16 sparsdem.m SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'sparsdem.m'`" test 3623 -eq "$shar_count" || $echo 'sparsdem.m:' 'original size' '3623,' 'current size' "$shar_count!" fi fi rm -fr _sh24386 exit 0