#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'chap3/backsub.m' <<'END_OF_FILE' Xfunction [y] = backsub(T,b); X%BACKSUB Back substitution X%y = backsub(T,b) computes the solution y of a nonsingular upper X%triangular system Ty = b using back substitution process. X%This program implements Algorithm 3.1.3 of the book. X%Input : Matrix T and vector b X%output : vector y X X X [m,n] = size(T); X if m~=n X error('matrix T is not square') X end; X y = zeros(n,1); X for i = n:-1:1 X sum = 0; X if (i ~= n) X sum = T(i,i+1:n)*y(i+1:n); X end; X if (T(i,i) ==0) X error('matrix T is singular') X end; X y(i) = (b(i) - sum ) / T(i,i); X end; END_OF_FILE if test 645 -ne `wc -c <'chap3/backsub.m'`; then echo shar: \"'chap3/backsub.m'\" unpacked with wrong size! fi # end of 'chap3/backsub.m' fi if test -f 'chap3/testing.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap3/testing.m'\" else echo shar: Extracting \"'chap3/testing.m'\" \(1558 characters\) sed "s/^X//" >'chap3/testing.m' <<'END_OF_FILE' X!rm ho Xdiary ho Xchdir ../chap3 Xn=7;a=rand(n,n);t=triu(a);b=rand(n,1); Xy=backsub(t,b);y1=backsubcopy(t,b);y-y1 Xchdir ../chap4 Xa=rand(n,n);b=rand(n,1); Xx=gauss(a,b);x1=gauss(a,b);x-x1 Xa1=a; a=housmulp(a,b);a1=housmulpcopy(a1,b);a-a1 Xa=rand(n,n);t=triu(a); Xt1=t; t=invuptr(t);t1=invuptrcopy(t1);t-t1 Xa=rand(n,n);a1=a; Xa=phousmul(a,b);a1=phousmulcopy(a1,b);a-a1 Xchdir ../chap5 Xa=rand(n,n);a1=a; X[a,u,m,q]=compiv(a);[a1,u1,m1,q1]=compivcopy(a1);a-a1,u-u1,m-m1,q-q1 Xa=rand(n,n);a1=a; X[h,p,a,v]=houshess(a);[h1,p1,a1,v1]=houshesscopy(a1);h-h1,p-p1,a-a1,v-v1 Xa=rand(n,n);a1=a; X[q,r,a,v]=housqr(a);[q1,r1,a1,v1]=housqrcopy(a1);q-q1,r-r1,a-a1,v-v1 Xa=rand(n,n);a1=a; X[q,r,a,v]=housqrn(a);[q1,r1,a1,v1]=housqrncopy(a1);q-q1,r-r1,a-a1,v-v1 Xa=rand(n,n);a1=a; X[a,l,u]=lugsel(a);[a1,l1,u1]=lugselcopy(a1); Xa=rand(n,n);a1=a; X[a,u,m]=parpiv(a);[a1,u1,m1]=parpivcopy(a1); Xchdir ../chap6 Xa=rand(n,n);a=a*a'+eye(n,n);a1=a; X[a,h]=choles(a);[a1,h1]=cholescopy(a1);a-a1,h-h1 Xa=rand(n,n); l=tril(a);b=rand(n,1); Xy=forelm(l,b);y1=forelmcopy(l,b);y-y1 Xa=rand(n,n);x0=rand(n,1);b=rand(n,1);ep=0.0001;num=20; X[x,it]=gaused(a,x0,b,ep,num);[x1,it1]=gausedcopy(a,x0,b,ep,num); Xx-x1,it-it1 Xa=rand(n,n);b=rand(b,1); Xx=gausswf(a,b);x1=gausswfcopy(a,b);x-x1 Xa=rand(n,n);x0=rand(n,1);b=rand(n,1);ep=0.0001;num=20; X[x,it]=sucov(a,x0,b,w,ep,num);[x1,it1]=sucovcopy(a,x0,b,w,ep,num); Xx-x1,it-it1 Xchdir ../chap7 Xa=rand(n,n); X[q,r]=clgrsch(a);[q1,r1]=clgrschcopy(a);q-q1,r-r1 Xa=rand(n,n); X[q,r]=mdgrsch(a);[q1,r1]=mdgrschcopy(a);q-q1,r-r1 Xa=rand(n,n); Xx=varcovar(a);x1=varcovarcopy(a);x-x1 Xdiary off END_OF_FILE if test 1558 -ne `wc -c <'chap3/testing.m'`; then echo shar: \"'chap3/testing.m'\" unpacked with wrong size! fi # end of 'chap3/testing.m' fi if test ! -d 'chap4' ; then echo shar: Creating directory \"'chap4'\" mkdir 'chap4' fi if test -f 'chap4/gauss.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap4/gauss.m'\" else echo shar: Extracting \"'chap4/gauss.m'\" \(1028 characters\) sed "s/^X//" >'chap4/gauss.m' <<'END_OF_FILE' Xfunction x = gauss(A,b) X%GAUSS Linear system solution using Gaussian elimination with economy in storage X%x = gauss(A,b) computes the solution x of the linear X%system Ax = b using Gaussian Elimination. The upper triangular X%part of A is overwritten by the triangular matrix produced at the end of the X%(n-1)th step and the multipliers are stored in the lower triangular part of A. X%On output b contains the final transformed vector. X%This program calls the MATCOM program BACKSUB. X%This program implements Algorithm 4.2.3 of the book. X%input : Matrix A and vector b X%output : vector x X X [m,n] = size(A); X for k = 1 : n-1 X if (A(k,k) == 0) X disp('the algorithm has encountered a zero pivot') X x=[]; X return; X end; X A(k+1:n,k) = -A(k+1:n,k)/ A(k,k); X b(k+1:n) = b(k+1:n) + A(k+1:n,k) * b(k); 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 U = triu(A); X x = backsub(U,b); END_OF_FILE if test 1028 -ne `wc -c <'chap4/gauss.m'`; then echo shar: \"'chap4/gauss.m'\" unpacked with wrong size! fi # end of 'chap4/gauss.m' fi if test -f 'chap4/housmulp.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap4/housmulp.m'\" else echo shar: Extracting \"'chap4/housmulp.m'\" \(580 characters\) sed "s/^X//" >'chap4/housmulp.m' <<'END_OF_FILE' Xfunction A = housmulp(A,u) X%HOUSMULP Post Multiplication By a Householder Matrix X%A = housmulp(A,u) computes the post-multiplication of X%a matrix A by the Householder matrix H generated by a X%vector u. The output matrix A contains the product AH. X%See Section 5.4 of the book. X%input : Matrix A and vector u X%output : Matrix A X X [m1,n] = size(A); X beta = 2/(u'*u); X for i = 1 : m1 X alpha = 0; X alpha = alpha + u(1:n) * A(i,1:n); X alpha = beta * alpha; X A(i,1:n) = A(i,1:n) - (alpha*u(1:n))'; X end; X end; X END_OF_FILE if test 580 -ne `wc -c <'chap4/housmulp.m'`; then echo shar: \"'chap4/housmulp.m'\" unpacked with wrong size! fi # end of 'chap4/housmulp.m' fi if test -f 'chap4/invuptr.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap4/invuptr.m'\" else echo shar: Extracting \"'chap4/invuptr.m'\" \(695 characters\) sed "s/^X//" >'chap4/invuptr.m' <<'END_OF_FILE' Xfunction T = invuptr(T); X%INVUPTR Inverse of an upper triangular matrix X%T = invuptr(T) computes the inverse of a nonsingular upper triangular X%matrix T. The output matrix T contains the inverse of T. X%This program implements Algorithm 4.2.2 of the book. X%Input : Matrix T X%output : Matrix T X X [m,n] = size(T); X if m~=n X disp('matrix T is not square') X return; X end; X s = eye(n,n); X for k = n:-1:1 X if ( T(k,k) == 0) X disp('matrix T is singular') X return; X end; X T(k,k) = 1/T(k,k); X for i = k-1 :-1 :1 X sum = 0; X sum = sum + T(i,i+1:k)*T(i+1:k,k); X T(i,k) = -sum/T(i,i); X end; X end; END_OF_FILE if test 695 -ne `wc -c <'chap4/invuptr.m'`; then echo shar: \"'chap4/invuptr.m'\" unpacked with wrong size! fi # end of 'chap4/invuptr.m' fi if test -f 'chap4/phousmul.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap4/phousmul.m'\" else echo shar: Extracting \"'chap4/phousmul.m'\" \(581 characters\) sed "s/^X//" >'chap4/phousmul.m' <<'END_OF_FILE' Xfunction A = phousmul(A,u) X%PHOUSMUL Pre Multiplication By Householder Matrix X%A = phousmul(A,u) computes the pre-multiplication X%of a matrix A by the Householder matrix generated X%by the vector u. The output matrix A contains the product HA. X%This program implements Algorithm 4.2.1 of the book. X%input : Matrix A and vector u X%output : Matrix A X X [m,n] = size(A); X beta = 2/(u'*u); X for j = 1 : n X alpha = 0; X alpha = alpha + u(1:m)'*A(1:m,j); X alpha = beta * alpha; X A(1:m,j) = A(1:m,j) - alpha * u(1:m); X end; END_OF_FILE if test 581 -ne `wc -c <'chap4/phousmul.m'`; then echo shar: \"'chap4/phousmul.m'\" unpacked with wrong size! fi # end of 'chap4/phousmul.m' fi if test ! -d 'chap5' ; then echo shar: Creating directory \"'chap5'\" mkdir 'chap5' fi if test -f 'chap5/compiv.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/compiv.m'\" else echo shar: Extracting \"'chap5/compiv.m'\" \(1701 characters\) sed "s/^X//" >'chap5/compiv.m' <<'END_OF_FILE' Xfunction [A,U,M,Q] = compiv(A); X%COMPIV Triangularization using Gaussian elimination with complete pivoting X%[A,U,M,Q] = compiv(A) produces an upper triangular matrix U, X%a permuted lower triangular matrix M, and a permutation matrix Q, X%using complete pivoting so that MAQ = U. The lower triangular part X%of the output matrix A contains the multipliers and the upper X%triangular part of A contains U. X%This program implements Algorithm 5.2.3 of the book. X%input : Matrix A X%output : Matrices A, U, M and Q. X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') X return; X end; X M = eye(n,n); X Q = eye(n,n); X U = zeros(n,n); X for k = 1:n-1 X M1 = eye(n,n); X d = k ; X e = k; X s = A(k,k); X for i = k:n X for j = k : n X if abs(A(i,j)) > abs(s) X s = A(i,j); X d = i; X e = j; X end; X end; X end; X if s == 0 X disp('the algorithm has encountered a zero pivot') X A=[]; X U=[]; X M=[]; X Q=[]; X return; X end X p(k) = d; X if (d ~= k | e ~=k) X [A(k,k:n),A(d,k:n)] = inter(A(k,k:n),A(d,k:n)); X [Q(:,k),Q(:,e)] = inter(Q(:,k),Q(:,e)); X [A(:,k),A(:,e)] = inter(A(:,k),A(:,e)); X [M(k,:),M(d,:)] = inter(M(k,:),M(d,:)); X end; X M1(k+1:n,k) = -A(k+1:n,k)/A(k,k); X A(k+1:n,k) = M1(k+1:n,k); X A(k+1:n,k+1:n) = A(k+1:n,k+1:n) +M1(k+1:n,k) * A(k,k+1:n); X%-updating M--- X M(k+1:n,1:n) = M(k+1:n,1:n) +M1(k+1:n,k) * M(k,1:n); X%---- X end; X U = triu(A); X end; X END_OF_FILE if test 1701 -ne `wc -c <'chap5/compiv.m'`; then echo shar: \"'chap5/compiv.m'\" unpacked with wrong size! fi # end of 'chap5/compiv.m' fi if test -f 'chap5/givhess.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/givhess.m'\" else echo shar: Extracting \"'chap5/givhess.m'\" \(1010 characters\) sed "s/^X//" >'chap5/givhess.m' <<'END_OF_FILE' Xfunction [H,P] = givhess(A) X%GIVHS Givens Hessenberg reduction X%[H,P] = givhess(A) produces an upper Hessenberg matrix H X%and an orthogonal matrix P using the Givens X%method so that H is orthogonally similar to A. PAP' = H. X%P is the product of (n-i+1) Givens rotations. X%This program implements Algorithm 5.5.4 of the book. X%input : Matrix A X%output : Matrices H and P X X [m,n] = size(A); X P = eye(m,m); X for p= 1: n-2 X for q = p+2:n X x=[A(p+1,p) A(q,p)]'; X [c,s] = givzero(x); X k = p+1; X l = q; X p1 = c * P(k,:) + s*P(l,:); X p2 = -s * P(k,:) + c*P(l,:); X P(k,:) = p1; X P(l,:) = p2; X a1 = c * A(k,:) + s*A(l,:); X a2 =-s * A(k,:) + c *A(l,:); X A(k,:) = a1; X A(l,:) = a2; X a3 = c * A(:,k) + s*A(:,l); X a4 = -s * A(:,k) + c *A(:,l); X A(:,k) = a3; X A(:,l) = a4; X end X end X H = A; X end; X X X END_OF_FILE if test 1010 -ne `wc -c <'chap5/givhess.m'`; then echo shar: \"'chap5/givhess.m'\" unpacked with wrong size! fi # end of 'chap5/givhess.m' fi if test -f 'chap5/givqr.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/givqr.m'\" else echo shar: Extracting \"'chap5/givqr.m'\" \(889 characters\) sed "s/^X//" >'chap5/givqr.m' <<'END_OF_FILE' Xfunction [Q,R] = givqr(A) X%GIVQR QR factorization by Givens rotation X%[Q,R] = givqr(A) produces an orthogonal matrix Q X%and a matrix R of the same size as A X%with zeros below the diagonal, such that A = QR, X%using Givens rotations. X%This program calls MATCOM program GIVZERO. X%This program implements Algorithm 5.5.4 of the book and X%also explicitly computes Q and R. X%input : Matrix A X%output : Matrices Q and R X X [m,n] = size(A); X Q = eye(m,m); X for k= 1:min(n,m-1) X for l = k+1:m X x=[A(k,k) A(l,k)]'; X [c,s] = givzero(x); X A1 = c * A(k,:) + s*A(l,:); X A2 =-s * A(k,:) + c *A(l,:); X A(k,:) = A1; X A(l,:) = A2; X A3 = c * Q(:,k) + s*Q(:,l); X A4 = -s * Q(:,k) + c *Q(:,l); X Q(:,k) = A3; X Q(:,l) = A4; X end X end X R = A; X end; X END_OF_FILE if test 889 -ne `wc -c <'chap5/givqr.m'`; then echo shar: \"'chap5/givqr.m'\" unpacked with wrong size! fi # end of 'chap5/givqr.m' fi if test -f 'chap5/givrot.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/givrot.m'\" else echo shar: Extracting \"'chap5/givrot.m'\" \(492 characters\) sed "s/^X//" >'chap5/givrot.m' <<'END_OF_FILE' Xfunction J = givrot(i,j,c,s,n) X%GIVROT Givens rotation matrix from Givens paramters c and s. X%J = givrot(i,j,c,s,n) forms the n x n Givens rotation matrix J X%from the Givens paramters c and s. A(i,i) = c, A(i,j) = s, X%A(j,i) = -s, A(j,j) = c. The other entries of J are the X%same as those of an n x n Identity matrix. X%input : Integers i, j, n and scalars c and s X%output : Matrix J X X X J = eye(n,n); X J(i,i) = c; X J(i,j) = s; X J(j,i) =-s ; X J(j,j) = c; END_OF_FILE if test 492 -ne `wc -c <'chap5/givrot.m'`; then echo shar: \"'chap5/givrot.m'\" unpacked with wrong size! fi # end of 'chap5/givrot.m' fi if test -f 'chap5/givszero.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/givszero.m'\" else echo shar: Extracting \"'chap5/givszero.m'\" \(587 characters\) sed "s/^X//" >'chap5/givszero.m' <<'END_OF_FILE' Xfunction A = givszero(i,j,A) X%GIVSZERO Givens specific zeroing in a vector. X%A = givszero(i,j,A) creates a zero in the (j,i)th position of a X%matrix A using Givens rotation matrix J. The output matrix A X%contains the product JA such that its (j,i)th entry is X%zero. This program calls the MATCOM program PGIVMUL and X%GIVZERO. X%This program implements Algorithm 5.5.3 of the book. X%input : Integers i, j and Matrix A X%output : Matrix A X X [m,n] = size(A); X x = zeros(2,1); X x(1) = A(i,i); X x(2) = A(j,i); X [c,s] = givzero(x); X A = pgivmul(A,i,j,c,s); END_OF_FILE if test 587 -ne `wc -c <'chap5/givszero.m'`; then echo shar: \"'chap5/givszero.m'\" unpacked with wrong size! fi # end of 'chap5/givszero.m' fi if test -f 'chap5/givzero.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/givzero.m'\" else echo shar: Extracting \"'chap5/givzero.m'\" \(592 characters\) sed "s/^X//" >'chap5/givzero.m' <<'END_OF_FILE' X function [c,s] = givzero(x) X%GIVZERO Givens zeroing in a vector x. X%[c,s] = givzero(x) produces the Givens parameters c and s X%for a 2-vector x such that J(1,2,theta)x has a zero X%in the second place. c = cos(theta), s = sin(theta). X%This program implements Algorithm 5.5.1 of the book. X%input : Vector x X%output : Scalars c and s X X X if abs(x(2)) > abs(x(1)) X t = x(1)/x(2); X s = 1/((1+t*t)^.5); X c = s*t ; X else X t = x(2)/x(1); X c = 1/((1+t*t)^.5); X s = c*t ; X end X X X end; X END_OF_FILE if test 592 -ne `wc -c <'chap5/givzero.m'`; then echo shar: \"'chap5/givzero.m'\" unpacked with wrong size! fi # end of 'chap5/givzero.m' fi if test -f 'chap5/houshess.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/houshess.m'\" else echo shar: Extracting \"'chap5/houshess.m'\" \(1518 characters\) sed "s/^X//" >'chap5/houshess.m' <<'END_OF_FILE' Xfunction [H,P,A,v] = houshess(A) X%HOUSHESS Householder Hessenberg reduction X%[H,P,A,v] = houshess(A) produces an upper Hessenberg matrix H from A X%using the Householder method such that H is orthogonally similar to A. X%PAP' = H. X%The upper Hessenberg part of the output matrix A X%contains the upper Hessenberg matrix H, and the lower Hessenberg X%part contains the components u_k+2,k, through u_n,k of the X%vectors u_n-k = (u_k+1,k, ..., u_n,k)', k = 1, ... n-2. X%The output vector v contains the first components X%u_k+1,k of the vector u_n-k. This program calls X%the MATCOM programs HOUSZERO and HOUSMULP. X%This program implements Algorithm 5.4.4 of the book. X%The matrix P is also explicitly formed. X%input : Matrix A X%output : Matrices H, P, A and vector v X X [m,n] = size(A); X P = eye(m,m); X for k = 1 : n-2 X [u,sigma] = houszero(A(k+1:n,k)); X A(k+1,k) = sigma; X s1 = size(u); X A(k+2:n,k) = u(2:s1); X v(k) = u(1); X beta = 2/(u'*u); X for j = k+1:n X s = 0; X s = s + u(1:n-k)' * A(k+1:n,j); X s = s * beta; X A(k+1:n,j) = A(k+1:n,j) - s * u(1:n-k); X end; X for i = 1:n X X s = 0; X s = s + u(1:n-k) * A(i,k+1:n); X s = s * beta; X A(i,k+1:n) = A(i,k+1:n) - (s * u(1:n-k))'; X end; X A; X P(1:n,k+1:n) = housmulp(P(1:n,k+1:n),u); X end; X H = triu(A,-1); X P = P'; X end; X END_OF_FILE if test 1518 -ne `wc -c <'chap5/houshess.m'`; then echo shar: \"'chap5/houshess.m'\" unpacked with wrong size! fi # end of 'chap5/houshess.m' fi if test -f 'chap5/housmulp.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/housmulp.m'\" else echo shar: Extracting \"'chap5/housmulp.m'\" \(580 characters\) sed "s/^X//" >'chap5/housmulp.m' <<'END_OF_FILE' Xfunction A = housmulp(A,u) X%HOUSMULP Post Multiplication By a Householder Matrix X%A = housmulp(A,u) computes the post-multiplication of X%a matrix A by the Householder matrix H generated by a X%vector u. The output matrix A contains the product AH. X%See Section 5.4 of the book. X%input : Matrix A and vector u X%output : Matrix A X X [m1,n] = size(A); X beta = 2/(u'*u); X for i = 1 : m1 X alpha = 0; X alpha = alpha + u(1:n) * A(i,1:n); X alpha = beta * alpha; X A(i,1:n) = A(i,1:n) - (alpha*u(1:n))'; X end; X end; X END_OF_FILE if test 580 -ne `wc -c <'chap5/housmulp.m'`; then echo shar: \"'chap5/housmulp.m'\" unpacked with wrong size! fi # end of 'chap5/housmulp.m' fi if test -f 'chap5/housqr.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/housqr.m'\" else echo shar: Extracting \"'chap5/housqr.m'\" \(1177 characters\) sed "s/^X//" >'chap5/housqr.m' <<'END_OF_FILE' Xfunction [Q,R,A,v] = housqr(A) X%HOUSQR Householder QR factorization X%[Q,R,A,v] = housqr(A) produces an orthogonal matrix Q X%and an upper triangular matrix R so that A = QR, using X%Householder matrices. The upper triangular part of X%the output matrix A contains R, and the lower triangular part X%contains the components u_k+1,k through u_nk of the X%vectors u_n-k=1 = (u_kk, .... u_nk)'. The output vector v X%contains first components u_kk of the vector u_n-k+1. X%This program calls MATCOM programs HOUSZERO and HOUSMULP. X%This program implements Algorithm 5.4.3 of the book. X%input : A square matrix A X%output : Matrices Q, R, A and vector v X X [m,n] = size(A); X Q = eye(m,m); X for k = 1 : n-1 X [x,sigma] = houszero(A(k:n,k)); X Q(1:n,k:n) = housmulp(Q(1:n,k:n),x); X A(k,k) = sigma; X s1 = size(x); X A(k+1:n,k) = x(2:s1) ; X v(k) = x(1); X beta = 2/(x' * x); X for j = k + 1 : n X s = 0; X s = s + x(1:n-k+1)' * A(k:n,j); X s = beta * s; X A(k:n,j) = A(k:n,j) - s * x(1:n-k+1); X end; X end; X R = triu(A); X end; END_OF_FILE if test 1177 -ne `wc -c <'chap5/housqr.m'`; then echo shar: \"'chap5/housqr.m'\" unpacked with wrong size! fi # end of 'chap5/housqr.m' fi if test -f 'chap5/housqrn.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/housqrn.m'\" else echo shar: Extracting \"'chap5/housqrn.m'\" \(1249 characters\) sed "s/^X//" >'chap5/housqrn.m' <<'END_OF_FILE' Xfunction [Q,R,A,v] = housqrn(A) X%HOUSQRN Householder QR factorization of a nonsquare matrix X%[Q,R,A,v] = housqrn(A) produces an orthogonal matrix Q X%and an upper triangular matrix R of the same size as X%A with zeros below the diagonal, so that A = QR, using X%Householder matrices. The upper triangular part of X%the output matrix A contains R, and the lower triangular part X%contains the components u_k+1,k through u_nk of the X%vectors u_n-k=1 = (u_kk, .... u_nk)'. The output vector v X%contains first components u_kk of the vector u_n-k+1. X%This program calls MATCOM programs HOUSZERO and HOUSMULP. X%see Section 5.4.2 of the book. X%input : Matrix A X%output : Matrices Q, R, A and vector v X X [m,n] = size(A); X S= min(n,m-1); X Q = eye(m,m); X for k = 1 : S X [x,sigma] = houszero(A(k:m,k)); X Q(1:m,k:m) = housmulp(Q(1:m,k:m),x); X A(k,k) = sigma ; X s1 = size(x); X A(k+1:m,k) = x(2:s1); X v(k) = x(1); X beta = 2/(x'*x); X for j = k+1:n X s = 0; X s = s + x(1:m-k+1)' * A(k:m,j); X s = beta * s; X A(k:m,j) = A(k:m,j) - s * x(1:m-k+1); X end; X end; X R = triu(A); X end; END_OF_FILE if test 1249 -ne `wc -c <'chap5/housqrn.m'`; then echo shar: \"'chap5/housqrn.m'\" unpacked with wrong size! fi # end of 'chap5/housqrn.m' fi if test -f 'chap5/houszero.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/houszero.m'\" else echo shar: Extracting \"'chap5/houszero.m'\" \(578 characters\) sed "s/^X//" >'chap5/houszero.m' <<'END_OF_FILE' Xfunction [u,sigma] = houszero(x) X%HOUSZERO Zeros in a vector using a Householder matrix. X%[u,sigma] = houszero(x) produces a vector u X%defining a Householder matrix H, and a scalar sigma X%such that Hx = [sigma, 0, ...., 0]'. X%This program implements Algorithm 5.4.1 of the book. X%input : vector x X%output : vector u, and scalar sigma X X [m,n] = size(x); X mm = max(abs(x)); X x = x/mm; X s = sign(x(1)); X if s == 0 X s = 1; X end; X sigma = s * norm(x,2); X u = x + sigma * eye(m,1); X sigma = -mm * sigma; END_OF_FILE if test 578 -ne `wc -c <'chap5/houszero.m'`; then echo shar: \"'chap5/houszero.m'\" unpacked with wrong size! fi # end of 'chap5/houszero.m' fi if test -f 'chap5/lugsel.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/lugsel.m'\" else echo shar: Extracting \"'chap5/lugsel.m'\" \(1043 characters\) sed "s/^X//" >'chap5/lugsel.m' <<'END_OF_FILE' Xfunction [A,L,U] = lugsel(A) X%LUGSEL LU factorization using Gaussian elimination without pivoting X%[A,L,U] = lugsel(A) produces the LU factorization of a matrix A using X%Gaussian elimination without pivoting : A = LU. L is unit lower triangular X%and U is upper triangular. This program implements X%Algorithm 5.2.1 of the book. The matrices L and U are X% also explicitly formed. The upper triangular part of the X%matrix A contains the matrix U and the lower triangular X%part of the matrix A contains the multipliers. X%input : Matrix A X%output : Matrices A, L and U X X [m,n] = size(A); X for k = 1 : n-1 X if (A(k,k) == 0) X disp('the algorithm has encountered a zero pivot') X L=[]; X U=[]; X return; X end; X A(k+1:n,k) = A(k+1:n,k)/ A(k,k); 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 U = triu(A); X L = tril(A,-1); X for i = 1:n X L(i,i) = 1; X end; X end; X END_OF_FILE if test 1043 -ne `wc -c <'chap5/lugsel.m'`; then echo shar: \"'chap5/lugsel.m'\" unpacked with wrong size! fi # end of 'chap5/lugsel.m' fi if test -f 'chap5/parpiv.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/parpiv.m'\" else echo shar: Extracting \"'chap5/parpiv.m'\" \(1459 characters\) sed "s/^X//" >'chap5/parpiv.m' <<'END_OF_FILE' Xfunction [A,U,M] = parpiv(A); X%PARPIV Triangularization using Gaussian Elimination with partial pivoting X%[A,U,M] = parpiv(A) produces an upper triangular matrix U X%and a permuted lower triangular matrix M using X%partial pivoting, so that MA = U. The lower triangular part X%of the output matrix A contains the multiplers and the upper triangular part X%contains U. X%This program implements Algorithm 5.2.2 of the book. X%input : Matrix A X%output : Matrices A, U and M X X [m,n] = size(A); X if m~=n X disp('matrix T is not square') X return; X end; X M = eye(n,n); X U = zeros(n,n); X for k = 1:n-1 X M1 = eye(n,n); X d = k ; X s = A(k,k); X for i = k+1:n X if abs(A(i,k)) > abs(s) X s = A(i,k); X d = i; X end; X end; X if (s == 0) X disp('the algorithm has encountered a zero pivot') X A =[]; X U =[]; X M =[]; X return; X end; X p(k) = d; X if d ~= k X [A(k,k:n),A(d,k:n)] = inter(A(k,k:n),A(d,k:n)); X [M(k,:),M(d,:)] = inter(M(k,:),M(d,:)); X end; X M1(k+1:n,k) = -A(k+1:n,k)/A(k,k); X A(k+1:n,k) = M1(k+1:n,k); X A(k+1:n,k+1:n) = A(k+1:n,k+1:n) +M1(k+1:n,k) * A(k,k+1:n); X%------------updating M--------- X M(k+1:n,1:n) = M(k+1:n,1:n) +M1(k+1:n,k) * M(k,1:n); X%-end updating X end; X U = triu(A); X end; END_OF_FILE if test 1459 -ne `wc -c <'chap5/parpiv.m'`; then echo shar: \"'chap5/parpiv.m'\" unpacked with wrong size! fi # end of 'chap5/parpiv.m' fi if test -f 'chap5/pgivmul.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap5/pgivmul.m'\" else echo shar: Extracting \"'chap5/pgivmul.m'\" \(566 characters\) sed "s/^X//" >'chap5/pgivmul.m' <<'END_OF_FILE' Xfunction A = pgivmul(A,i,j,c,s) X%PGIVMUL Premultiplication by a Givens Matrix. X%A = pgivmul(A,i,j,c,s) computes the premultiplication X%of the matrix A by the Givens matrix, J(i,j,theta), where X%1 <= i <=j <= m, and c = cos(theta), s = sin(theta). X%The output matrix A contains the product JA. X%This program implements Algorithm 5.5.2 of the book. X%Input : Matrix A, Givens parameters c and s, indices i and j X%Output : Matrix A X X [m,n] = size(A); X a1 = A(i,:); X a2 = A(j,:); X A(i,:) = c * a1 + s * a2; X A(j,:) = -s * a1 + c * a2; END_OF_FILE if test 566 -ne `wc -c <'chap5/pgivmul.m'`; then echo shar: \"'chap5/pgivmul.m'\" unpacked with wrong size! fi # end of 'chap5/pgivmul.m' fi if test ! -d 'chap6' ; then echo shar: Creating directory \"'chap6'\" mkdir 'chap6' fi if test -f 'chap6/choles.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/choles.m'\" else echo shar: Extracting \"'chap6/choles.m'\" \(894 characters\) sed "s/^X//" >'chap6/choles.m' <<'END_OF_FILE' Xfunction [A,H] = choles(A); X%CHOLES Cholesky factorization X%[A,H] = choles(A) computes the Cholesky factorization of a X%symmetric positive matrix A : A = HH'. H is a lower triangular X%matrix with positive diagonal entries. On output, the lower X%triangular part ofthe output matrix A contains H. Note that MATLAB CHOL(A) X%produces an upper triangular matrix R such that R'*R = A. This X%program implements Algorithm 6.4.4 of the book. X%input : Matrix A X%output : Matrices A and H X X [m,n] = size(A); X if m~=n X disp('matrix A is not square'); X return; X end; X for k = 1:n X for i = 1:k-1 X sum = 0; X if (i ~= 1) X sum =A(i,1:i-1)*A(k,1:i-1)'; X end; X A(k,i) = (A(k,i) - sum)/A(i,i); X end X sum = 0; X if (k ~= 1) X sum =A(k,1:k-1)*A(k,1:k-1)'; X end; X A(k,k) = (A(k,k) - sum) ^ 0.5; X end X H = tril(A); X end; END_OF_FILE if test 894 -ne `wc -c <'chap6/choles.m'`; then echo shar: \"'chap6/choles.m'\" unpacked with wrong size! fi # end of 'chap6/choles.m' fi if test -f 'chap6/congrad.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/congrad.m'\" else echo shar: Extracting \"'chap6/congrad.m'\" \(1169 characters\) sed "s/^X//" >'chap6/congrad.m' <<'END_OF_FILE' Xfunction [x,iter] = congrad(A,x0,b,ep,numitr) X%CONGRAD Classical conjugate gradient method X%[x,iter] = congrad(A,x0,b,ep,numitr) computes the solution x X%of a symmetric positive definite linear system Ax = b using X%the Conjugate Gradient method. X%x0 is the initial approximation, ep is the tolerance, X%and numitr is the user supplied number of iterations. X%If the conjugate-gradient method converged, iter contains X%the iteration number needed to converge. If the conjugate-gradient method X%did not converge, iter contains numitr. X%This program implements Algorithm 6.10.4 of the book. X%input : Matrix A, vectors x0 and b, scalar ep, and integer numitr X%output : Vector x, and integer iter. X X [m,n]=size(A); X p = b - A * x0; X r = p; X for k = 1 : numitr X iter = k; X w = A * p; X conol = (norm(r))^2 ; X alpha = conol/ (p' * w); X x = x0 + alpha * p; X rnew = r - alpha * w; X conv = (norm(rnew))^2; X if conv < ep X return; X end ; X beta = conv / conol; X p = rnew + beta * p; X r = rnew; X end; X end; END_OF_FILE if test 1169 -ne `wc -c <'chap6/congrad.m'`; then echo shar: \"'chap6/congrad.m'\" unpacked with wrong size! fi # end of 'chap6/congrad.m' fi if test -f 'chap6/forelm.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/forelm.m'\" else echo shar: Extracting \"'chap6/forelm.m'\" \(665 characters\) sed "s/^X//" >'chap6/forelm.m' <<'END_OF_FILE' Xfunction y = forelm(L,b); X%FORELM Forward elimination X%y = forelm(L,b) computes the solution y of a nonsingular lower X%triangular system Ly = b using forward elimination process. X%This program implements Algorithm 6.4.1 of the book. X%Input : Matrix L and vector b X%output : vector y X X [m,n] = size(L); X if m~=n X disp('matrix L is not square') X return; X end; X y = zeros(n,1); X y(1) = b(1) / L(1,1); X for i = 2:n X sum =L(i,1:i-1)*y(1:i-1); X if ( L(i,i) == 0) X disp('matrix L is singular') X return; X end; X y(i) = (b(i) - sum ) / L(i,i); X end; X end; END_OF_FILE if test 665 -ne `wc -c <'chap6/forelm.m'`; then echo shar: \"'chap6/forelm.m'\" unpacked with wrong size! fi # end of 'chap6/forelm.m' fi if test -f 'chap6/gaused.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/gaused.m'\" else echo shar: Extracting \"'chap6/gaused.m'\" \(1228 characters\) sed "s/^X//" >'chap6/gaused.m' <<'END_OF_FILE' Xfunction [x,iter] = gaused(A,x0,b,ep,numitr); X%GAUSED Gauss-Seidel method X%[x,iter] = gaused(A,x0,b,ep,numitr) computes the solution x of X%the linear system Ax = b, using the Gauss-Seidel method. X%x0 is the initial approximation, ep is the tolerance, numitr is the X%user supplied number of iterations. If the Gauss-Seidel method X%converged, iter contains the number of iterations needed to converge. X%If the Gauss-Seidel method did not converge, iter contains numitr. X%This program implements Algorithm 6.10.2 of the book. X%input : Matrix A, vectors x0 and b, scalar ep, and integer numitr X%output : vector x and integer iter X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') X return; X end; X xnew = zeros(m,1); X x = x0; X for k = 1 : numitr X iter = k; X for i = 1:n X s1 = 0; X s2 = 0; X if (i ~= 1) X s1 = A(i,1:i-1) * xnew(1:i-1); X end; X if (i ~= n) X s2 = A(i,i+1:n) * x(i+1:n); X end X xnew(i) = (b(i) - s1 - s2) / A(i,i); X end; X if (norm(xnew - x) / norm(x) ) < ep X x = xnew ; X return; X end X x = xnew; X end; END_OF_FILE if test 1228 -ne `wc -c <'chap6/gaused.m'`; then echo shar: \"'chap6/gaused.m'\" unpacked with wrong size! fi # end of 'chap6/gaused.m' fi if test -f 'chap6/gausswf.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/gausswf.m'\" else echo shar: Extracting \"'chap6/gausswf.m'\" \(1339 characters\) sed "s/^X//" >'chap6/gausswf.m' <<'END_OF_FILE' Xfunction x = gausswf(A,b); X%GAUSSWF Linear system solution without explicit factorization X%x = gausswf(A,b) computes the solution x of the system Ax = b X%using Gaussian elimination with partial pivoting. The factorization X%MA = U is never explicitly formed. X%The program calls the MATCOM program BACKSUB and INTER. X%The program implements Algorithm 6.4.2 of the book. X%input : Matrix A and vector b X%output : vector x X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') X return; X end; X s1 = b ; X U = zeros(n,n); X for k = 1:n-1 X M1 = eye(n,n); X d = k ; X s = A(k,k); X for i = k+1:n X if abs(A(i,k)) > abs(s) X s = A(i,k); X d = i; X end; X end; X if (s == 0) X disp('the algorithm has encountered a zero pivot'); X x =[]; X return; X end; X p(k) = d; X if d ~= k X [A(k,k:n),A(d,k:n)] = inter(A(k,k:n),A(d,k:n)); X [b(k),b(d)] = inter(b(k),b(d)); X end; X M1(k+1:n,k) = -A(k+1:n,k)/A(k,k); X A(k+1:n,k) = M1(k+1:n,k); X A(k+1:n,k+1:n) = A(k+1:n,k+1:n) +M1(k+1:n,k) * A(k,k+1:n); X b(k+1:n) = b(k+1:n) + M1(k+1:n,k) * b(k); X end; X U = triu(A); X x = backsub(U,b); X end; END_OF_FILE if test 1339 -ne `wc -c <'chap6/gausswf.m'`; then echo shar: \"'chap6/gausswf.m'\" unpacked with wrong size! fi # end of 'chap6/gausswf.m' fi if test -f 'chap6/hagcond1.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/hagcond1.m'\" else echo shar: Extracting \"'chap6/hagcond1.m'\" \(1404 characters\) sed "s/^X//" >'chap6/hagcond1.m' <<'END_OF_FILE' Xfunction [cnd,iter] = hagcond1(A); X%HAGCOND1 Hager's norm-1 condition number estimator X%[CND,ITER] = hagcond1(A) produces CND, the norm-1 condition number estimate X%using Hager's Algorithm. Integer iter is the number of X%iterations needed for the algorithm to converge. X%This program calls the MATCOM programs BACKSUB and PARPIV. X%The program implements Algorithm 6.7.1 of the book. X%input : Matrix A X%output : Scalar cnd and integer iter X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X rho = 0; X b = zeros(n,1); X y = zeros(n,1); X for i = 1 : n X b(i) = 1/n; X end; X flag = 0; X iter = 0; X count = 0; X while flag == 0 X%solve Ax = b using parpiv X [storea,U,M] = parpiv(A); X bprime = M*b; X x = backsub(U,bprime) ; X if norm(x,1) <= rho X flag = 1; X return; X end; X iter = iter + 1; X rho = norm(x,1); X for i = 1:n X if x(i) >= 0 X y(i) = 1; X else X y(i) = -1; X end; X end; X%Solve A'z = y using parpiv X c1 = A'; X [storea,U,M] = parpiv(c1); X bprime = M * y; X z = backsub(U,bprime); X [j] = absmax(z); X if abs(z(j)) > z'*b X b = zeros(n,1); X b(j) = 1; X else X flag = 1; X end; X cnd = rho * norm(A,1); X end; X end; X END_OF_FILE if test 1404 -ne `wc -c <'chap6/hagcond1.m'`; then echo shar: \"'chap6/hagcond1.m'\" unpacked with wrong size! fi # end of 'chap6/hagcond1.m' fi if test -f 'chap6/icholes.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/icholes.m'\" else echo shar: Extracting \"'chap6/icholes.m'\" \(833 characters\) sed "s/^X//" >'chap6/icholes.m' <<'END_OF_FILE' Xfunction [A,L] = icholes(A); X%ICHOLES Incomplete Cholesky factorization X%[A,L] = icholes(A) produces the incomplete Cholesky factor A of X% a sparse symmetric positive definite matrix A. The program implements X% Algorithm 6.10.6 of the book. X%input : Matrix A X%output : Matrices A and L X X [m,n] = size(A); X if m~=n X disp('matrix A is not square'); X return; X end; X L = eye(n,n); X L(1,1) = A(1,1) ^ 0.5; X for i = 2:n X for j = 1:i-1 X if A(i,j) == 0; X L(i,j) = 0; X else X sum = 0; X for k = 1:j-1 X sum = sum + L(i,k)*L(j,k); X end; X L(i,j) = (A(i,j) - sum)/L(j,j); X end; X end; X sum = 0; X for k = 1:i-1 X sum = sum + L(i,k) * L(i,k); X end X L(i,i) = (A(i,i) - sum) ^ 0.5; X end; X end; END_OF_FILE if test 833 -ne `wc -c <'chap6/icholes.m'`; then echo shar: \"'chap6/icholes.m'\" unpacked with wrong size! fi # end of 'chap6/icholes.m' fi if test -f 'chap6/incholes.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/incholes.m'\" else echo shar: Extracting \"'chap6/incholes.m'\" \(544 characters\) sed "s/^X//" >'chap6/incholes.m' <<'END_OF_FILE' Xfunction T = incholes(A); X%INCHOLES Inverse of a symmetric positive definite matrix using Cholesky factorization X%T = incholes(A) computes the inverse of a symmetric positive X%definite matrix A using its Cholesky factor H. X%inv(A) = inv(H')inv(H). X%This program calls the MATCOM program CHOLES. X%See section 6.5.3 of the book. X%input : Matrix A X%output : Matrix T X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X [C,H] = choles(A); X Hinv = inv(H); X T = Hinv'*Hinv; X end; X END_OF_FILE if test 544 -ne `wc -c <'chap6/incholes.m'`; then echo shar: \"'chap6/incholes.m'\" unpacked with wrong size! fi # end of 'chap6/incholes.m' fi if test -f 'chap6/incompiv.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/incompiv.m'\" else echo shar: Extracting \"'chap6/incompiv.m'\" \(376 characters\) sed "s/^X//" >'chap6/incompiv.m' <<'END_OF_FILE' Xfunction T = incompiv(A); X%INCOMPIV Inverse by complete pivoting X%T = incompiv(A) computes the inverse of a matrix A using complete X%pivoting : inv(A) = Q*inv(U)*M. X%input : Matrix A X%output : Matrix T X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X [C,U,M,Q] = compiv(A); X T = Q * (inv(U)) * M; X end; END_OF_FILE if test 376 -ne `wc -c <'chap6/incompiv.m'`; then echo shar: \"'chap6/incompiv.m'\" unpacked with wrong size! fi # end of 'chap6/incompiv.m' fi if test -f 'chap6/inlu.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/inlu.m'\" else echo shar: Extracting \"'chap6/inlu.m'\" \(394 characters\) sed "s/^X//" >'chap6/inlu.m' <<'END_OF_FILE' Xfunction T = inlu(A); X%INLU Inverse using the LU factorization X%T = inlu(A) computes the inverse T of A using LU factorization X%obtained by Gaussian elimination without pivoting. X%input : Matrix A X%output : Matrix T X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X [L,U] = lu(A); X T = inv(U) * inv(L); X end; END_OF_FILE if test 394 -ne `wc -c <'chap6/inlu.m'`; then echo shar: \"'chap6/inlu.m'\" unpacked with wrong size! fi # end of 'chap6/inlu.m' fi if test -f 'chap6/inparpiv.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/inparpiv.m'\" else echo shar: Extracting \"'chap6/inparpiv.m'\" \(429 characters\) sed "s/^X//" >'chap6/inparpiv.m' <<'END_OF_FILE' Xfunction T = inparpiv(A); X%INPARPIV Inverse by partial pivoting X%T = inparpiv(A) computes the inverse T of a matrix A using Gaussian X%elimination with partial pivoting : inv(A) = inv(U) * M X%This program uses MATCOM program PARPIV. X%input : Matrix A X%output : Matrix T X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X [C,U,M] = parpiv(A); X T = (U) \ M; X end; END_OF_FILE if test 429 -ne `wc -c <'chap6/inparpiv.m'`; then echo shar: \"'chap6/inparpiv.m'\" unpacked with wrong size! fi # end of 'chap6/inparpiv.m' fi if test -f 'chap6/iterref.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/iterref.m'\" else echo shar: Extracting \"'chap6/iterref.m'\" \(961 characters\) sed "s/^X//" >'chap6/iterref.m' <<'END_OF_FILE' Xfunction [x,iter] = iterref(A,b,x0,ep,numitr); X%ITERREF Iterative refinement X%[x,iter] = iterref(A,b,x0,ep,numitr) iteratively computes succesive refinements X%of an initial solution x0 of the system Ax = b. ep is X%the tolerance. numitr is the user supplied number of iterations. X%If the process converged to the desired accuracy, iter contains X%the iteration number needed to converge. X%If the process did not converge, iter contains numitr. X%This program implements Algorithm 6.9.1 of the book. X%input : Matrix A, vectors b and x0, scalar ep and integer numitr X%output : vector x and integer iter X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X x = x0; X for k1 = 1:numitr X iter = k1; X d = A; X r = b - A * x; X y = d\r; X xnew = x + y ; X if (norm(y)) / norm(xnew) < ep X x = xnew; X return; X end; X x = xnew; X end; X end; END_OF_FILE if test 961 -ne `wc -c <'chap6/iterref.m'`; then echo shar: \"'chap6/iterref.m'\" unpacked with wrong size! fi # end of 'chap6/iterref.m' fi if test -f 'chap6/jacobi.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/jacobi.m'\" else echo shar: Extracting \"'chap6/jacobi.m'\" \(1208 characters\) sed "s/^X//" >'chap6/jacobi.m' <<'END_OF_FILE' Xfunction [x,iter] = jacobi(A,x0,b,ep,numitr); X%JACOBI Jacobi method X%[x,iter] = jacobi(A,x0,b,ep,numitr) computes the solution x of Ax = b using X%the jacobi iterative method. ep is the tolerance. numitr is the X%user supplied number of iteration. x0 is the initial X%approximate to the solution. X%If the Jacobi method converged, iter contains the iteration number X%needed to converge. X%If the Jacobi method did not converge, iter contains numitr. X%This program implements Algorithm 6.10.1 of the book. X%input : Matrix A and vectors x0 and b, scalar ep and integer numitr X%output : vector x and integer iter X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X xnew = zeros(m,1); X x = x0; X for k = 1:numitr X iter = k; X for i = 1:n X s = 0; X for j = 1:n X if (j ~= i) X s = A(i,j) * x(j) + s; X end; X end; X xnew(i) = (b(i) - s) / A(i,i); X end; X if (norm(xnew - x) / norm(x) ) < ep X iter = k; X x = xnew; X return; X end X x = xnew; X end; X end; END_OF_FILE if test 1208 -ne `wc -c <'chap6/jacobi.m'`; then echo shar: \"'chap6/jacobi.m'\" unpacked with wrong size! fi # end of 'chap6/jacobi.m' fi if test -f 'chap6/nichol.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/nichol.m'\" else echo shar: Extracting \"'chap6/nichol.m'\" \(941 characters\) sed "s/^X//" >'chap6/nichol.m' <<'END_OF_FILE' Xfunction [L,D] = nichol(A); X%NICHOL No-fill Incomplete LDL'. X%[L,D] = nichol(A) computes the no-fill incomplete Cholesky X%factorization of A, without any square roots : A = LDL'. X%This program implements Algorithm 6.10.7 of the book. X%input : Matrix A X%output : Matrices L and D X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X L = eye(n,n); X D = eye(n,n); X D(1,1) = A(1,1); X for i = 2:n X for j = 1:i-1 X if A(i,j) == 0 X L(i,j) = 0; X else X sum = 0 ; X for k = 1 : j-1 X sum = sum + L(i,k) * D(k,k) * L(j,k); X end; X L(i,j) = (A(i,j) - sum) / D(j,j); X end; X end; X sum = 0; X for k = 1 : i-1 X sum = sum + L(i,k) * L(i,k) * D(k,k); X end; X D(i,i) = A(i,i) - sum; X end; X end; END_OF_FILE if test 941 -ne `wc -c <'chap6/nichol.m'`; then echo shar: \"'chap6/nichol.m'\" unpacked with wrong size! fi # end of 'chap6/nichol.m' fi if test -f 'chap6/shermor.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/shermor.m'\" else echo shar: Extracting \"'chap6/shermor.m'\" \(531 characters\) sed "s/^X//" >'chap6/shermor.m' <<'END_OF_FILE' Xfunction [H] = shermor(A,u,v); X%SHERMOR Sherman Morrison formula X%H = shermor(A,u,v) computes the inverse of a matrix obtained X%by rank-one change in a matrix A, using the Sherman-Morrison X%formula : inv(A - uv'). u and v are column vectors. X%see section 6.5.2 of the book X%input : Matrix A and vectors u and v X%output : Matrix H X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') X return; X end; X ainv = inv(A); X alpha = 1/(1-v'*ainv*u); X H = ainv + alpha*ainv*u*v'*ainv; X end; END_OF_FILE if test 531 -ne `wc -c <'chap6/shermor.m'`; then echo shar: \"'chap6/shermor.m'\" unpacked with wrong size! fi # end of 'chap6/shermor.m' fi if test -f 'chap6/sucov.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap6/sucov.m'\" else echo shar: Extracting \"'chap6/sucov.m'\" \(1481 characters\) sed "s/^X//" >'chap6/sucov.m' <<'END_OF_FILE' Xfunction [x,iter] = sucov(A,x0,b,w,ep,numitr); X%SUCOV Successive overrelaxation X%[x,iter] = sucov(A,x0,b,w,ep,numitr) computes the solution x X%of the linear system Ax = b using the successive X%overrelaxation iterative method. x0 is the initial X%solution, numitr is the number of iterations to be performed, X%specified by the user and w is the relaxation X%parameter. (w > 1). If w = 1 then the successive X%overrelaxation iterative method reduces to the X%Gauss-Seidel iterative method. ep is the tolerance X%If the successive overrelaxation method converged, X%iter contains the iteration number needed to converge. If the successive X%overrelaxation method did not converge, iter contains X%numitr. X%This program implements Algorithm 6.10.3 of the book. X%input : Matrix A, vectors x0 and b, scalars w and ep and integer numitr X%output : vector x and integer iter X X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X xnew = zeros(m,1); X x = x0; X for k = 1 : numitr X iter = k; X for i = 1:n X s1 = 0; X if i ~= 1 X s1 = A(i,1:i-1) * xnew(1:i-1); X end; X s2 = 0; X if (i~=n) X s2 = A(i,i+1:n) * x(i+1:n) ; X end; X xnew(i) = (w/A(i,i))*(b(i) - s1 - s2) + (1 - w) *x(i); X end; X if (norm(xnew - x) / norm(x) ) < ep X x = xnew ; X return; X end X x = xnew; X end; END_OF_FILE if test 1481 -ne `wc -c <'chap6/sucov.m'`; then echo shar: \"'chap6/sucov.m'\" unpacked with wrong size! fi # end of 'chap6/sucov.m' fi if test ! -d 'chap7' ; then echo shar: Creating directory \"'chap7'\" mkdir 'chap7' fi if test -f 'chap7/clgrsch.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap7/clgrsch.m'\" else echo shar: Extracting \"'chap7/clgrsch.m'\" \(664 characters\) sed "s/^X//" >'chap7/clgrsch.m' <<'END_OF_FILE' Xfunction [Q,R] = clgrsch(A); X%CLGRSCH Classical Gram-Schmidt for QR factorization X%[Q,R] = clgrsch(A) computes the QR factorization X%of an m x n matrix A using the classical Gram-Schmidt method : X%A = QR, R is n x n upper triangular, Q is m x n and has orthonormal columns. X%This program implements Algorithm 7.8.3 of the book. X%input : Matrix A X%output : Matrices Q and R X X [m,n] = size(A); X for k = 1 :n X R(1:k-1,k) = Q(:,1:k-1)'* A(:,k); X sum = 0; X for i = 1:k-1 X sum = sum+ R(i,k) * Q(:,i); X end X Q(:,k) = A(:,k) - sum; X R(k,k) = norm(Q(:,k)); X Q(:,k) = Q(:,k) / R(k,k); X end; X end; END_OF_FILE if test 664 -ne `wc -c <'chap7/clgrsch.m'`; then echo shar: \"'chap7/clgrsch.m'\" unpacked with wrong size! fi # end of 'chap7/clgrsch.m' fi if test -f 'chap7/lsfrmgs.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap7/lsfrmgs.m'\" else echo shar: Extracting \"'chap7/lsfrmgs.m'\" \(569 characters\) sed "s/^X//" >'chap7/lsfrmgs.m' <<'END_OF_FILE' Xfunction x = lsfrmgs(A,b) X%LSFRMGS Least-squares solution by MGS. X%x = lsfrmgs(A,b) computes the least squares solution X%x of the full-rank overdetermined system Ax = b using modified X%Gram-Schmidt. X%This program calls the MATCOM programs BACKSUB and MDGRSCH. X%This program implements Algorithm 7.8.5 of the book. X%input : Matrix A and vector b X%output : vector x X X b1 = b; X [m,n] = size(A); X [Q,R] = mdgrsch(A); X for k = 1:n X del(k) = Q(:,k)'*b1; X b1 = b1 - del(k) * Q(:,k); X end; X x = backsub(R,del); X end; END_OF_FILE if test 569 -ne `wc -c <'chap7/lsfrmgs.m'`; then echo shar: \"'chap7/lsfrmgs.m'\" unpacked with wrong size! fi # end of 'chap7/lsfrmgs.m' fi if test -f 'chap7/lsfrnme.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap7/lsfrnme.m'\" else echo shar: Extracting \"'chap7/lsfrnme.m'\" \(568 characters\) sed "s/^X//" >'chap7/lsfrnme.m' <<'END_OF_FILE' Xfunction [x] = lsfrnme(A,b); X%LSFRNME Least-squares solution using normal equations X%x = lsfrnme(A,b) computes the least squares solution x to the X%full-rank overdetermined system Ax = b using the X%Normal equations. X%This program calls the MATCOM programs CHOLES, FORELM and X%BACKSUB. X%This program implements Algorithm 7.8.1 of the book. X%input : Matrix A and vector b X%output : vector x X X [m,n] = size(A); X y = zeros(n,1); X c = A'*b; X HHtr = A'*A; X [HHtr,H] = choles(HHtr); X y = forelm(H,c) ; X U = H'; X x = backsub(U,y); X end; END_OF_FILE if test 568 -ne `wc -c <'chap7/lsfrnme.m'`; then echo shar: \"'chap7/lsfrnme.m'\" unpacked with wrong size! fi # end of 'chap7/lsfrnme.m' fi if test -f 'chap7/lsfrqrh.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap7/lsfrqrh.m'\" else echo shar: Extracting \"'chap7/lsfrqrh.m'\" \(754 characters\) sed "s/^X//" >'chap7/lsfrqrh.m' <<'END_OF_FILE' Xfunction x = lsfrqrh(A,b); X%LSFRQRH Least Squares solutions using Householder QR X%x = lsfrqrh(A,b) computes the least squares solution x to the full X%rank overdetermined system Ax = b using the Householder-Golub method. X%This program calls the MATCOM programs BACKSUB, HOUSZERO and PHOUSMUL. X%This program implements Algorithm 7.8.2 of the book. X%input : Matrix A and vector b X%output : vector x X X [m,n] = size(A); X y = zeros(n,1); X s= min(n,m-1); X for k = 1 : s X [u,sigma] = houszero(A(k:m,k)); X A(k:m,k:n) = phousmul(A(k:m,k:n),u); X b(k:m) = phousmul(b(k:m),u); X end; X r1 = A; X c = b; X ran = rank(r1); X R = r1(1:ran,:); X x = backsub(R,c); X end; END_OF_FILE if test 754 -ne `wc -c <'chap7/lsfrqrh.m'`; then echo shar: \"'chap7/lsfrqrh.m'\" unpacked with wrong size! fi # end of 'chap7/lsfrqrh.m' fi if test -f 'chap7/lsitrn1.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap7/lsitrn1.m'\" else echo shar: Extracting \"'chap7/lsitrn1.m'\" \(774 characters\) sed "s/^X//" >'chap7/lsitrn1.m' <<'END_OF_FILE' Xfunction x = lsitrn1(A,b,numitr) X%LSITRN1 Linear system analog least-squares iterative refinement. X%x = lsitrn1(A,b,numitr) refines iteratively a computed least X%squares solution of Ax = b. numitr is the user supplied X%number of iterations. X%This is the linear system analog iterative refinement X%for a least squares solution. For another iterative X%refinement algorithm see LSITRN2. X%This program calls the MATCOM program LSFRQRH. X%This program implements Algorithm 7.10.1 X%of the book. X%input : Matrix A, vector b, and integer numitr X%output : vector x X X [m,n] = size(A); X xold = zeros(n,1); X for e1 = 1:numitr X r = b - A * xold; X c = lsfrqrh(A,r); X xold = xold + c; X end; X x = xold; X end; END_OF_FILE if test 774 -ne `wc -c <'chap7/lsitrn1.m'`; then echo shar: \"'chap7/lsitrn1.m'\" unpacked with wrong size! fi # end of 'chap7/lsitrn1.m' fi if test -f 'chap7/lsitrn2.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap7/lsitrn2.m'\" else echo shar: Extracting \"'chap7/lsitrn2.m'\" \(991 characters\) sed "s/^X//" >'chap7/lsitrn2.m' <<'END_OF_FILE' Xfunction x = lsitrn2(A,b,numitr) X%LSITRN2 Iterative refinement for least-squares solution. X%x = lsitrn2(A,b,numitr) refines iteratively a computed least X%squares solution x of Ax = b. numitr is the user supplied X%number of iterations. This program implements Algorithm 7.10.2 X%of the book. X%input : Matrix A, vector b, and integer numitr X%output : vector x X X [m,n] = size(A); X xold = zeros(n,1); X rold = zeros(m,1); X bigm = [eye(m,m) A;A' zeros(n,n)]; X for e1 = 1:numitr X r1 = b - ( eye(m,m) * rold + A * xold ); X r2 = A' * rold ; X [q,dp1] = qr(A); X ran = rank(dp1); X R1 = dp1(1:ran,:); X rprime = q' * r1; X r1prime = rprime(1:n); X r2prime = rprime(n+1:m); X c2prime = R1'\r2; X c2 = R1\(r1prime - c2prime); X c1 = q * [c2prime; r2prime]; X rold = rold + c1; X xold = xold + c2; X end; X x = xold; X end; END_OF_FILE if test 991 -ne `wc -c <'chap7/lsitrn2.m'`; then echo shar: \"'chap7/lsitrn2.m'\" unpacked with wrong size! fi # end of 'chap7/lsitrn2.m' fi if test -f 'chap7/lsrdqrh.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap7/lsrdqrh.m'\" else echo shar: Extracting \"'chap7/lsrdqrh.m'\" \(802 characters\) sed "s/^X//" >'chap7/lsrdqrh.m' <<'END_OF_FILE' X function [x] = lsrdqrh(A,b) X%LSRDQRH Least-squares solutions for the rank deficient problem using Householder QR X%x = lsrdqrh(A,b) computes a least-squares solution x X%to the rank - deficient overdetermined system Ax = b using X%Householder QR Factorization A. X%This program calls the MATCOM program BACKSUB. X%This program implements Algorthim 7.8.6 of the book. X%input : Matrix A and vector b X%output : vector x X X [m,n] = size(A); X ran = rank(A); X y = zeros(n,1); X [Q,R,p] = qr(A); X z1 = Q'*b; X c=z1(1:ran); X R11=R(1:ran,1:ran); X if ((ran+1) <= n) X r12=R(1:ran,ran+1:n); X y2=rand(n-ran,1); X y(ran+1:n)=y2; X z2 = c-r12*y2; X else X z2 = c; X end X y1 = backsub(R11,z2); X y(1:ran)=y1; X x=p*y; X end; X END_OF_FILE if test 802 -ne `wc -c <'chap7/lsrdqrh.m'`; then echo shar: \"'chap7/lsrdqrh.m'\" unpacked with wrong size! fi # end of 'chap7/lsrdqrh.m' fi if test -f 'chap7/mdgrsch.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap7/mdgrsch.m'\" else echo shar: Extracting \"'chap7/mdgrsch.m'\" \(655 characters\) sed "s/^X//" >'chap7/mdgrsch.m' <<'END_OF_FILE' Xfunction [Q,R] = mdgrsch(A); X%MDGRSCH Modified Gram-Schmidt for QR factorization X%[Q,R] = mdgrsch(A) computes the QR factorization of an m x n matrix A X%using the modified Gram-Schmidt method : A = QR, R is n x n upper X%triangular and Q is m x n and has orthonormal columns. X%This program implements Algorithm 7.8.4 of the book. X%input : Matrix A X%output : Matrices Q and R X X [m,n] = size(A); X Q(:,1:n) = A(:,1:n); X for k = 1 :n X R(k,k) = norm(Q(:,k)); X Q(:,k) = Q(:,k) /R(k,k); X for j = k+1:n X R(k,j) = Q(:,k)'*Q(:,j); X Q(:,j) = Q(:,j) - R(k,j) * Q(:,k); X end; X end; X end; END_OF_FILE if test 655 -ne `wc -c <'chap7/mdgrsch.m'`; then echo shar: \"'chap7/mdgrsch.m'\" unpacked with wrong size! fi # end of 'chap7/mdgrsch.m' fi if test -f 'chap7/mnudnme.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap7/mnudnme.m'\" else echo shar: Extracting \"'chap7/mnudnme.m'\" \(395 characters\) sed "s/^X//" >'chap7/mnudnme.m' <<'END_OF_FILE' Xfunction [x] = mnudnme(A,b) X%MNUDNME Minimum-norm solution for the underdetermined system using normal X%equations x = mnudnme(A,b) computes the minimum norm solution x to the full X%rank underdetermined system Ax = b using normal equations. This X%program implements Algorthim 7.9.1 of the book. X%input : Matrix A and vector b X%output : vector x X X y = (A*A') \ b; X x = A'*y; END_OF_FILE if test 395 -ne `wc -c <'chap7/mnudnme.m'`; then echo shar: \"'chap7/mnudnme.m'\" unpacked with wrong size! fi # end of 'chap7/mnudnme.m' fi if test -f 'chap7/mnudqrh.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap7/mnudqrh.m'\" else echo shar: Extracting \"'chap7/mnudqrh.m'\" \(486 characters\) sed "s/^X//" >'chap7/mnudqrh.m' <<'END_OF_FILE' Xfunction [x] = mnudqrh(A,b) X%MNUDQRH Minimum-norm solution for the underdetermined problem using QR X%x = mnudqrh(A,b) computes the minimum norm solution x to the full rank X%underdetermined system Ax = b using the Householder QR factorization of A X%This program implements Algortihm 7.9.2 of the book. X%input : Matrix A and vector b X%output : vector x X X [Q,R] = qr(A'); X ran = rank(R); X Q1 = Q(:,1:ran); X R1 = R(1:ran,:); X yr = R1'\b; X x = Q1*yr; END_OF_FILE if test 486 -ne `wc -c <'chap7/mnudqrh.m'`; then echo shar: \"'chap7/mnudqrh.m'\" unpacked with wrong size! fi # end of 'chap7/mnudqrh.m' fi if test -f 'chap7/varcovar.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap7/varcovar.m'\" else echo shar: Extracting \"'chap7/varcovar.m'\" \(889 characters\) sed "s/^X//" >'chap7/varcovar.m' <<'END_OF_FILE' Xfunction X = varcovar(A); X%VARCOVAR Variance-Covariance matrix X%X = varcovar(A) computes the variance-covariance matrix X from A : X%X = inv(A' * A) explicitly without computing A'A. X%This program implements Algorithm 7.11.1 of the book. X%input : Matrix A X%output : matrix X X X [m,n] = size(A); X [Q,R1] = qr(A); X ran = rank(R1); X R = R1(1:ran,:); X X=zeros(n,n); X ey = eye(n,n); X e = ey(:,n) ; X y = R\(e/R(n,n)) ; X X(:,n) = y; X X(n,:) = y'; X for k = n-1 :-1:1 X sum = R(k,k+1:n)*X(k+1:n,k); X X(k,k) = ( 1/R(k,k) - sum )/(R(k,k)); X for i = k-1 :-1:1 X i; X k; X sum1 = R(i,i+1:k) * X(i+1:k,k); X sum2 = R(i,k+1:n) * X(k,k+1:n)'; X X(i,k) = (sum1 + sum2)/(-R(i,i)); X X(k,i) = X(i,k); X end; X end; X end; END_OF_FILE if test 889 -ne `wc -c <'chap7/varcovar.m'`; then echo shar: \"'chap7/varcovar.m'\" unpacked with wrong size! fi # end of 'chap7/varcovar.m' fi if test ! -d 'chap8' ; then echo shar: Creating directory \"'chap8'\" mkdir 'chap8' fi if test -f 'chap8/invitr.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap8/invitr.m'\" else echo shar: Extracting \"'chap8/invitr.m'\" \(1047 characters\) sed "s/^X//" >'chap8/invitr.m' <<'END_OF_FILE' Xfunction [x,iter] = invitr(A,x0, sigma, ep, numitr) X%INVITR Inverse iteration X%[x,iter] = invitr(A,x0, sigma, ep, numitr) computes an approximation x, of the X%eigenvector corresponding to a given approximation sigma of an eigenvalue, X%using inverse iteration. x0 is the initial approximation, X%ep is the tolerance and numitr is the maximum number of iterations. X%If the iteration converged, iter is the number of iterations X%needed to converge. If the iteration did not converge, X%iter contains numitr. X%This program implements Algorithm 8.5.2 of the book. X%input : Matrix A, vector x0, scalars sigma, ep and integer numitr X%output : vector x and integer iter X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X x = x0; X for k = 1 : numitr X iter = k; X prevx = x; X xhat = (A - sigma * eye(n,n))\ x; X x = xhat/norm(xhat); X if norm( (A - sigma*eye(n,n)) * x , inf) < ep X return; X end; X end; X end; X END_OF_FILE if test 1047 -ne `wc -c <'chap8/invitr.m'`; then echo shar: \"'chap8/invitr.m'\" unpacked with wrong size! fi # end of 'chap8/invitr.m' fi if test -f 'chap8/lansym.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap8/lansym.m'\" else echo shar: Extracting \"'chap8/lansym.m'\" \(901 characters\) sed "s/^X//" >'chap8/lansym.m' <<'END_OF_FILE' Xfunction [V,T] = lansym(A,v1,k) X%LANSYM Symmetric Lanczos algorithm X% [V,T] = lansym(A,v1,k) computes a k x k symmteric tridiagonal matrix X%T and an orthonormal matrix V using the Lanczos X%algorithm. Matrix A is symmetric and v1 is a unit vector. X%Integer k is the number of basis vectors v1,v2, ... X%to be computed. This program implements Algorithm 8.12.2 X%of the book. X%Input : Matrix A, vector v1 and integer k X%Output : matrices V and T X X [m,n] = size(A); X v0 = zeros(m,1); X bet = 1; X r = v1; X vol = v0; X for j = 1:k X v = r/bet; X V(:,j) = v; X u = A * v; X rnew = u - bet * vol; X alph = v'*u; X alpha(j) = alph; X rnew = rnew - alph * v; X bet = norm(rnew); X beta(j) = bet; X r = rnew; X vol = v; X end; X T = diag(alpha(1:k))+ diag(beta(1:k-1),1)+diag(beta(1:k-1),-1); X end; END_OF_FILE if test 901 -ne `wc -c <'chap8/lansym.m'`; then echo shar: \"'chap8/lansym.m'\" unpacked with wrong size! fi # end of 'chap8/lansym.m' fi if test -f 'chap8/power.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap8/power.m'\" else echo shar: Extracting \"'chap8/power.m'\" \(1297 characters\) sed "s/^X//" >'chap8/power.m' <<'END_OF_FILE' Xfunction [lambda1,x,iter] = power(A,x0,ep,numitr) X%POWER Power method X%[lambda1,x,iter] = power(A,x0,ep,numitr) computes the dominant X%eigenvalue lambda1 and the corresponding eigenvector using the X%power method. x0 is the initial vector, ep is the X%tolerance and numitr is the X%maximum number of iterations. On output, if the power method converged, X%iter contains the iteration number needed to converge. X%If the power method did not converge, iter contains the value of numitr. X%This program implements Algorithm 8.5.1 of the book. X%input : Matrix A, vector x0, scalar ep and integer numitr X%output : Scalar lambda1, vector x, and integer iter X X t1 = cputime; X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X x = x0; X xhat = x0; X for k = 1 : numitr X iter = k; X prevxhat = xhat; X xhat = A* x; X x = xhat/max(xhat); X X if norm( (A*x - max(xhat)*x) ) < ep X lambda1 = max(xhat); X return; X end; X end; X lambda1 = max(xhat); X t2 = cputime - t1; X disp('the cpu time in seconds taken for this algorithm') X t2 X disp('the flop count for this algorithm is') X flps = flops X end; END_OF_FILE if test 1297 -ne `wc -c <'chap8/power.m'`; then echo shar: \"'chap8/power.m'\" unpacked with wrong size! fi # end of 'chap8/power.m' fi if test -f 'chap8/qritrb.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap8/qritrb.m'\" else echo shar: Extracting \"'chap8/qritrb.m'\" \(530 characters\) sed "s/^X//" >'chap8/qritrb.m' <<'END_OF_FILE' Xfunction Ak = qritrb(A,numitr) X%QRITRB Basic QR iteration X%Ak = qritrb(A,numitr) returns the matrix Ak after 'numitr' number X%of basic QR iterations. numitr is a user specified integer. X%See Section 8.9.1 of the book. X%input : Matrix A and integer numitr X%output : Matrix Ak X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X [Q,R] = qr(A); X for k = 1 : numitr X Anew = R*Q; X [Q,R] = qr(Anew); X end; X Ak = Anew; X end; END_OF_FILE if test 530 -ne `wc -c <'chap8/qritrb.m'`; then echo shar: \"'chap8/qritrb.m'\" unpacked with wrong size! fi # end of 'chap8/qritrb.m' fi if test -f 'chap8/qritrdse.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap8/qritrdse.m'\" else echo shar: Extracting \"'chap8/qritrdse.m'\" \(820 characters\) sed "s/^X//" >'chap8/qritrdse.m' <<'END_OF_FILE' Xfunction NH = qritrdse(H,numitr) X%QRITRDSE Explicit double shift QR iteration X%NH = qritrdse(H,numitr) returns an upper Hessenberg matrix NH, starting from X%a given upper Hessenberg matrix H, using explicit double-shift QR X%iteration. numitr is the number of iterations specified by the user. X%see section 8.9.5 of the book. The matrix H is overwritten by X%itself after each iteration. X%input : Matrix H and integer numitr X%output : matrix NH X X [m,n] = size(H); X if m~=n X disp('matrix H is not square') ; X return; X end; X for k = 1 : numitr X t = H(n-1,n-1) + H(n,n); X d = H(n-1,n-1) * H(n,n) - H(n,n-1) * H(n-1,n); X N = H*H - t*H + d*eye(n,n) ; X [q,r] = qr(N); X Hnew = q'*H*q; X H = Hnew; X end; X NH = H; X end; END_OF_FILE if test 820 -ne `wc -c <'chap8/qritrdse.m'`; then echo shar: \"'chap8/qritrdse.m'\" unpacked with wrong size! fi # end of 'chap8/qritrdse.m' fi if test -f 'chap8/qritrdsi.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap8/qritrdsi.m'\" else echo shar: Extracting \"'chap8/qritrdsi.m'\" \(1408 characters\) sed "s/^X//" >'chap8/qritrdsi.m' <<'END_OF_FILE' Xfunction [NH,Q] = qritrdsi(H) X%QRITRDSI One iteration step of the Implicit Double Shift QR Iteration X%[NH,Q] = qritrdsi(H) returns an upper Hessenberg matrix NH, X%starting from a given upper Hessenberg matrix H, using X%one step of implicit double shift QR iteration. Q is the X%transforming orthogonal matrix : Q'HQ = NH. X%This program calls the MATCOM program HOUSZERO. X%The program implements Algorithm 8.9.1 of the book. X%input : Matrix H X%output : Matrix NH and Q X X [m,n] = size(H); X if m~=n X disp('matrix H is not square') ; X return; X end; X Q = eye(n,n) ; X t = H(n-1,n-1) + H(n,n); X d = H(n-1,n-1) * H(n,n) - H(n,n-1) * H(n-1,n); X x = H(1,1)*H(1,1) - t* H(1,1) + d + H(1,2) * H(2,1); X y = H(2,1) * (H(1,1) + H(2,2) -t); X z = H(2,1) * H(3,2); X for k = 0: n -3 X [u,sigma] = houszero([x y z]'); X c = eye(size(u), size(u)) - 2 * u * u' / (u' * u); X Pk = eye(n,n); X Pk(k+1:k+3,k+1:k+3) = c; X Q = Q*Pk ; X H = Pk'*H*Pk; X x = H(k+2,k+1); X y = H(k+3,k+1); X if k < n-3 X z = H(k+4,k+1); X end; X end; X [u,sigma] = houszero([x y ]'); X c = eye(size(u), size(u)) - 2 * u * u' / (u' * u); X Pk = eye(n,n); X Pk(n-1:n,n-1:n) = c; X Q = Q*Pk; X H = Pk'*H*Pk; X NH = H; X end; X END_OF_FILE if test 1408 -ne `wc -c <'chap8/qritrdsi.m'`; then echo shar: \"'chap8/qritrdsi.m'\" unpacked with wrong size! fi # end of 'chap8/qritrdsi.m' fi if test -f 'chap8/qritrh.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap8/qritrh.m'\" else echo shar: Extracting \"'chap8/qritrh.m'\" \(621 characters\) sed "s/^X//" >'chap8/qritrh.m' <<'END_OF_FILE' Xfunction A = qritrh(A,numitr) X%QRITRH Hessenberg-QR iteration X%A = qritrh(A,numitr) returns a reduced Hessenberg matrix starting X%from a given upper Hessenberg matrix A using QR iterations. X%Givens rotations are used to factor A(k) into Q(k)R(k) . X%numitr is the user supplied number of iterations. X%See Section 8.9.2 of the book. X%input : Matrix A and integer numitr X%output : Matrix A X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X [Q,R] = givqr(A); X for k = 1 : numitr X A = R*Q; X [Q,R] = qr(A); X end; X end; END_OF_FILE if test 621 -ne `wc -c <'chap8/qritrh.m'`; then echo shar: \"'chap8/qritrh.m'\" unpacked with wrong size! fi # end of 'chap8/qritrh.m' fi if test -f 'chap8/qritrsse.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap8/qritrsse.m'\" else echo shar: Extracting \"'chap8/qritrsse.m'\" \(648 characters\) sed "s/^X//" >'chap8/qritrsse.m' <<'END_OF_FILE' Xfunction H = qritrsse(A,numitr) X%QRITRSSE Explicit single-shift QR iteration X%H = qritrsse(A,numitr) returns an upper Hessenberg matrix H, starting X%from an arbitrary matrix A, using single-shift explicit QR iteration. X%numitr is the user-supplied number of iterations. See Section 8.9.4 of X%the book. X%input : Matrix A and integer numitr X%output : Matrix H X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X [H,PP] = givhess(A); X for k = 0 : numitr X Hnn = H(n,n); X [Q,R] = qr(H - Hnn*eye(n,n)); X H = R * Q + Hnn*eye(n,n); X end; X end; END_OF_FILE if test 648 -ne `wc -c <'chap8/qritrsse.m'`; then echo shar: \"'chap8/qritrsse.m'\" unpacked with wrong size! fi # end of 'chap8/qritrsse.m' fi if test -f 'chap8/rayqot.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap8/rayqot.m'\" else echo shar: Extracting \"'chap8/rayqot.m'\" \(1141 characters\) sed "s/^X//" >'chap8/rayqot.m' <<'END_OF_FILE' Xfunction [sigma,x,iter] = rayqot(A,x0,ep,numitr) X%RAYQOT Rayleigh quotient iteration X%[sigma,x,iter] = rayqot(A,x0,ep,numitr) computes an approximate X%eigenpair(sigma,x) of a matrix A using Rayleigh-Quotient iteration. X%x0 is the initial approximation to the eigen vector x. ep is the tolerance X%numitr is the user-supplied number of iteration. On output, if the X%Rayleigh quotient iteration converged, iter contains the iteration number X%needed to converge. If the iteration did not converge, iter contains X%the value of numitr. X%This program implements Algorithm 8.5.3 of the book. X%input : Matrix A, vector x0, scalar ep and integer numitr X%output : Scalar sigma, vector x and integer iter X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X prevsig = 0; X x = x0; X for k = 0 : numitr X iter = k; X sigma = (x'*A*x)/(x'*x); X xhat = (A-sigma * eye(n,n))\x; X x = xhat/max(xhat); X if norm( (A-sigma * eye(n,n))*x ) < ep X return; X end X prevsig = sigma; X end; X end; END_OF_FILE if test 1141 -ne `wc -c <'chap8/rayqot.m'`; then echo shar: \"'chap8/rayqot.m'\" unpacked with wrong size! fi # end of 'chap8/rayqot.m' fi if test -f 'chap8/senseig.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap8/senseig.m'\" else echo shar: Extracting \"'chap8/senseig.m'\" \(652 characters\) sed "s/^X//" >'chap8/senseig.m' <<'END_OF_FILE' Xfunction s = senseig(A) X%SENSEIG Reciprocal of the condition numbers of eigenvalues X%s = senseig(A) computes a vector s containing the reciprocals X%of the condition numbers of the eigenvalues of a diagonalizable matrix A. X%See Section 8.7.2 of the book. X%input : Matrix A X%output : vector s X X [m,n] = size(A); X if m~=n X disp('matrix A is not square') ; X return; X end; X e = eye(n,n) ; X [X,d] = eig(A); X invX = inv(X); X for i = 1 : n X x = (X(:,i))/norm(X(:,i)) ; X invXtr=invX'; y = invXtr(:,i)/norm(invXtr(:,i)); X s(i) = abs(y' * x); X end; X end; END_OF_FILE if test 652 -ne `wc -c <'chap8/senseig.m'`; then echo shar: \"'chap8/senseig.m'\" unpacked with wrong size! fi # end of 'chap8/senseig.m' fi if test ! -d 'chap9' ; then echo shar: Creating directory \"'chap9'\" mkdir 'chap9' fi if test -f 'chap9/cholqr.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap9/cholqr.m'\" else echo shar: Extracting \"'chap9/cholqr.m'\" \(770 characters\) sed "s/^X//" >'chap9/cholqr.m' <<'END_OF_FILE' Xfunction [X,lam] = cholqr(A,B) X%CHOLQR Cholesky-QR algorithm for the symmetric definite pencil X%[X,lam] = cholqr(A,B) computes the eigenvalues and eigenvectors for the X%symmetric definite pencil (A - lambda B). Matrix X contains the eigenvectors X% and the vector lam contains the eigenvalues. X%A is symmetric and B is symmetric positive definite. X%This program calls the MATCOM program BACKSUB. X%This program implements Algorithm 9.5.1 of the book. X%input : Matrices A and B X%output : Matrix X, and vector lam. X X [m,n] = size(A); X L = chol(B); X L = L'; X D = inv(L); X C = D * A * D'; X [Y,D] = eig(C); X for i = 1:n X X(:,i) = backsub(L',Y(:,i)); X end; X lam = diag(D); X end; END_OF_FILE if test 770 -ne `wc -c <'chap9/cholqr.m'`; then echo shar: \"'chap9/cholqr.m'\" unpacked with wrong size! fi # end of 'chap9/cholqr.m' fi if test -f 'chap9/genrayqt.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap9/genrayqt.m'\" else echo shar: Extracting \"'chap9/genrayqt.m'\" \(1085 characters\) sed "s/^X//" >'chap9/genrayqt.m' <<'END_OF_FILE' Xfunction [lamda,x,iter] = genrayqt(A,B,x0,ep,numitr) X%GENRAYQT Generalized Rayleigh quotient iteration X%[lamda,x,iter] = genrayqt(A,B,x0,ep,numitr) computes approximations to X%the eigenvectors and the eigenvalues of the symmetric definite X%pencil A - lambda B. A is symmetric and B is positive definite. X%x0 is chosen such that norm(x0) = 1. numitr is the user X%supplied number of iterations. ep is the tolerance. X%If the method converged, iter contains the number of iterations X%needed to converge. If the method did not converge iter X%contains the value of numitr. X%This program implements Algorithm 9.5.4 of the book. X%input : Matrices A and B, vector x0, scalar ep and integer numitr X%output : Scalar lamda, vector x and integer iter. X X [m,n] = size(A); X x = x0; X for k = 1:numitr X iter = k; X lamda = (x'*A*x) / (x'*B*x); X xkhat = (A - lamda * B) \ (B * x); X newxk = xkhat/norm(xkhat); X if norm( (A - lamda * eye(n,n)) * newxk ) < ep X return; X end X x = newxk; X end X end END_OF_FILE if test 1085 -ne `wc -c <'chap9/genrayqt.m'`; then echo shar: \"'chap9/genrayqt.m'\" unpacked with wrong size! fi # end of 'chap9/genrayqt.m' fi if test -f 'chap9/gensturm.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap9/gensturm.m'\" else echo shar: Extracting \"'chap9/gensturm.m'\" \(814 characters\) sed "s/^X//" >'chap9/gensturm.m' <<'END_OF_FILE' Xfunction r = gensturm(A,B) X%GENSTURM Sturm sequence method X% r = gensturm(A,B) X%Given A and B n x n matrices, the algorithm computes X%the eigenvalues of the symmetric definite tridiagonal X%pencil (A - lambda B) by computing the zeros of Pn(lamda). X%Matrices A and B are both symmetric tridiagonal X%and B is positive definite. See Section 9.9.1 of the book. X%input : Matrices A and B X%output : Vector r X X X [m,n] = size(A); X p0 = 1; X p1 = [-B(1,1); A(1,1)]; X for r = 2:n X pnew1 = A(r,r) * [0;p1] - B(r,r) * [p1;0] ; X pnew2 = -(A(r,r-1)^2*[0;0;p0] -2*A(r,r-1)*B(r,r-1)*[0;p0;0]); X pnew3 = -(+B(r,r-1)^2 *[p0;0;0]); X pnew = pnew1 + pnew2 + pnew3; X p0 = p1; X p1 = pnew; X end; X r = roots(pnew); X end; END_OF_FILE if test 814 -ne `wc -c <'chap9/gensturm.m'`; then echo shar: \"'chap9/gensturm.m'\" unpacked with wrong size! fi # end of 'chap9/gensturm.m' fi if test -f 'chap9/hesstri.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap9/hesstri.m'\" else echo shar: Extracting \"'chap9/hesstri.m'\" \(1313 characters\) sed "s/^X//" >'chap9/hesstri.m' <<'END_OF_FILE' Xfunction [A,B] = hesstri(A,B) X%HESSTRI Hessenberg-Triangular Reduction X% [A,B] = hesstri(A,B) overwrites X%A with and upper hessenberg matrix by orthogonal transformations X%and B with an upper triangular matrix by orthogonal transformations. X%see Section 9.3.1 of the book X%input : Matrices A and B X%output : Matrices A and B X X X [m,n] = size(A); X Q = eye(n,n); X Z = eye(n,n); X [U,R] = qr(B); X B = R; X A = U'*A; X for j = 1:n-2 X for i = n:-1:j+2 X x = zeros(2,1); X x(1) = A(i-1,j); X x(2) = A(i,j); X [c,s] = givzero(x) ; X A1 = c * A(i-1,:) + s*A(i,:); X A2 = -s * A(i-1,:) + c*A(i,:); X A(i-1,:) = A1; X A(i,:) = A2; X B1 = c * B(i-1,:) + s*B(i,:); X B2 = -s * B(i-1,:) + c*B(i,:); X B(i-1,:) = B1; X B(i,:) = B2; X x = zeros(2,1); X x(1) = B(i,i-1); X x(2) = B(i,i); X [c,s] = givzero(x); X B4 = c * B(:,i-1) + s*B(:,i); X B3 = -s * B(:,i-1) + c *B(:,i); X B(:,i-1) = B3; X B(:,i) = B4; X A4 = c * A(:,i-1) + s*A(:,i); X A3 = -s * A(:,i-1) + c *A(:,i); X A(:,i-1) = A3; X A(:,i) = A4; X end; X end; END_OF_FILE if test 1313 -ne `wc -c <'chap9/hesstri.m'`; then echo shar: \"'chap9/hesstri.m'\" unpacked with wrong size! fi # end of 'chap9/hesstri.m' fi if test -f 'chap9/invitrgn.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap9/invitrgn.m'\" else echo shar: Extracting \"'chap9/invitrgn.m'\" \(793 characters\) sed "s/^X//" >'chap9/invitrgn.m' <<'END_OF_FILE' Xfunction [v,iter] = invitrgn(A,B,v0,lambda,numitr) X%INVITRGN Eigenvectors. X%[v,iter] = invitrgn(A,B,v0,lambda,numitr) computes the X%eigenvector v corresponding to an eigenvalue lambda of X%the generalized eigenvalue problem A - lambda B, by X%inverse iteeration. X%v0 is the initial eigenvector. X%numitr is the user supplied number of iterations. X%If the method did not converge, iter contains the value of numitr. X%This program implements Algorithm 9.4.1 of the book. X%input : Matrices A and B, vector v0, scalars lambda and integer numitr X%output : vector v, integer iter X X [m,n] = size(A); X v = v0; X for k = 1:numitr X iter = k; X vkhat = (A-lambda*B) \ (B * v); X vnew = vkhat/norm(vkhat,inf); X v = vnew; X end; X end; END_OF_FILE if test 793 -ne `wc -c <'chap9/invitrgn.m'`; then echo shar: \"'chap9/invitrgn.m'\" unpacked with wrong size! fi # end of 'chap9/invitrgn.m' fi if test -f 'chap9/lansymgn.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap9/lansymgn.m'\" else echo shar: Extracting \"'chap9/lansymgn.m'\" \(1003 characters\) sed "s/^X//" >'chap9/lansymgn.m' <<'END_OF_FILE' Xfunction [V,T] = lansymgn(A,B,v1,j) X%LANSYMGN Lanczos Algorithm for Symmetric Definite Pencil X%[V,T] = lansymgn(A,B,v1,j) constructs a symmetric tridiagonal matrix T X%and an orthonormal matrix Vj = [v1,v2, ..., vj] such that X%Vj'AVj = Tj and Vj'BVj = I(jxj). X%Matrices A and B are symmetric and B is positive definite. X%This program implements Algorithm 9.9.1 of the book. X%input : Matrices A and B, vector v1, and integer j. X%output : Matrices V and T X X [m,n] = size(A); X L = chol(B); X L=L'; X bet = 1; X r = v1; X v0 = zeros(m,1); X for i = 1:j X if i == 1 X v = v1; X else X v = B\( r /bet); X end; X V(:,i) = v; X alph = v'*(A*v - B * v0); X alpha(i) = alph; X r = A * v - alph * B * v - bet*B*v0; X bet = norm(inv(L)*r); X beta(i) = bet; X v0 = v; X end; X T = diag(alpha(1:j)) + diag(beta(1:j-1),1) + diag(beta(1:j-1),-1); END_OF_FILE if test 1003 -ne `wc -c <'chap9/lansymgn.m'`; then echo shar: \"'chap9/lansymgn.m'\" unpacked with wrong size! fi # end of 'chap9/lansymgn.m' fi if test -f 'chap9/quadeig1.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap9/quadeig1.m'\" else echo shar: Extracting \"'chap9/quadeig1.m'\" \(692 characters\) sed "s/^X//" >'chap9/quadeig1.m' <<'END_OF_FILE' Xfunction [V,lam] = quadeig1(M,D,K) X%QUADEIG1 Quadratic Eigenvalue Problem via Standard Reduction X%[V,lam] = quadeig1(M,D,K) computes the eigenvalues and X%eigenvectors of the quadratic eigenvalue problem via X%the reduction to a standard eigenvalue problem. X%The vector lam contains the eigenvalues and X%matrix V the corresponding eigenvectors. X%See Section 9.8 of the book. X%input : Matrices M,D and K X%output : Matrix V and vector lam X X [m,n] = size(M); X if m~=n X disp('matrix M is not square') ; X return; X end; X A = [zeros(n,n) eye(n,n) ; -inv(M)*K -inv(M) * D]; X [U,D] = eig(A); X V = U(1:n,:); X lam = diag(D); END_OF_FILE if test 692 -ne `wc -c <'chap9/quadeig1.m'`; then echo shar: \"'chap9/quadeig1.m'\" unpacked with wrong size! fi # end of 'chap9/quadeig1.m' fi if test -f 'chap9/quadeig2.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap9/quadeig2.m'\" else echo shar: Extracting \"'chap9/quadeig2.m'\" \(745 characters\) sed "s/^X//" >'chap9/quadeig2.m' <<'END_OF_FILE' Xfunction [V,lam] = quadeig2(M,D,K) X%QUADEIG2 Quadratic Eigenvalue Problem via Reduction to a Generalized Eigenvalue Problem X%[V,lam] = quadeig2(M,D,K) computes the eigenvalues and X%eigenvectors of the quadratic eigenvalue problem via X%reduction to a generalized eigenvalue problem. X%The vector lam contains the eigenvalues and the X%matrix V contains the corresponding eigenvectors. X%See Section 9.8 of the book. X%input : Matrices M,D and K X%output : Matrix V and vector lam X X [m,n] = size(M); X if m~=n X disp('matrix M is not square') ; X return; X end; X B = [D K; K zeros(n,n)]; X C = [-M zeros(n,n); zeros(n,n) K]; X [U,D1] = eig(B,C); X lam = diag(D1); X V = U(1:n,:); END_OF_FILE if test 745 -ne `wc -c <'chap9/quadeig2.m'`; then echo shar: \"'chap9/quadeig2.m'\" unpacked with wrong size! fi # end of 'chap9/quadeig2.m' fi if test -f 'chap9/simdiag.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chap9/simdiag.m'\" else echo shar: Extracting \"'chap9/simdiag.m'\" \(517 characters\) sed "s/^X//" >'chap9/simdiag.m' <<'END_OF_FILE' Xfunction X = simdiag(A,B) X%SIMDIAG Simultaneous diagonalization of a symmetric definite pencil X%X = simdiag(A,B) computes a nonsingular matrix X such that X%X'BX is the identity matrix and X'AX is a diagonal matrix. X%Matrices A and B are symmetric and B is positive definite. X%This program implements Algorithm 9.5.3 of the book. X%input : Matrices A and B X%output : Matrix X X X [m,n] = size(A); X [C,L] = choles(B); X Linv = inv(L); X C = Linv * A * Linv'; X [Y,d] = eig(C); X X = Linv' * Y; END_OF_FILE if test 517 -ne `wc -c <'chap9/simdiag.m'`; then echo shar: \"'chap9/simdiag.m'\" unpacked with wrong size! fi # end of 'chap9/simdiag.m' fi if test ! -d 'chapt10' ; then echo shar: Creating directory \"'chapt10'\" mkdir 'chapt10' fi if test -f 'chapt10/bidiag.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chapt10/bidiag.m'\" else echo shar: Extracting \"'chapt10/bidiag.m'\" \(726 characters\) sed "s/^X//" >'chapt10/bidiag.m' <<'END_OF_FILE' Xfunction A = bidiag(A) X%BIDIAG Reduction to bidiagonal form. X%B = bidiag(A) reduces the matrix A to bidiagonal form. X%This program calls the MATCOM programs HOUSZERO and PHOUSMUL. X%See Section 10.9 of the book. X%input : Matrix A X%output : Matrix A X X X [m,n] = size(A); X B = A; X for i = 1:n X B = A(i:m,i:n); X [m1,n1] = size(B); X x = B(1:m1,1); X [u,sigma] = houszero(x); X B = phousmul(B,u); X if i <= (n-2) X C = B(:,2:n1); X [m2,n2] = size(C); X x = C(1,1:n2); X x = x'; X [u,sigma] = houszero(x); X C = housmulp(C,u); X B(:,2:n1) = C; X end; X A(i:m,i:n) = B; X end; END_OF_FILE if test 726 -ne `wc -c <'chapt10/bidiag.m'`; then echo shar: \"'chapt10/bidiag.m'\" unpacked with wrong size! fi # end of 'chapt10/bidiag.m' fi if test -f 'chapt10/covsvd.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chapt10/covsvd.m'\" else echo shar: Extracting \"'chapt10/covsvd.m'\" \(636 characters\) sed "s/^X//" >'chapt10/covsvd.m' <<'END_OF_FILE' Xfunction C = covsvd(A) X%COVSVD Variance-Covariance matrix X%C = covsvd(A) computes the variance covariance matrix X%inv(A'*A) using the SVD of A. X% see section 10.7 of the book. X%input : Matrix A X%output : Matrix C X X [U,S,V] = svd(A); X [m,n] = size(A); X [U,S,V] = svd(A); X r = 0; X for i = 1:n X if S(i,i) > 10^(-15) X r = r+1; X else X i = n; X end; X end; X for i = 1:n X for j = 1:n X sum = 0; X for k = 1:r X sum = sum + V(i,k) * V(j,k) / (S(k,k) * S(k,k)); X end; X C(i,j) = sum; X end; X end; X end END_OF_FILE if test 636 -ne `wc -c <'chapt10/covsvd.m'`; then echo shar: \"'chapt10/covsvd.m'\" unpacked with wrong size! fi # end of 'chapt10/covsvd.m' fi if test -f 'chapt10/lsqrsvd.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chapt10/lsqrsvd.m'\" else echo shar: Extracting \"'chapt10/lsqrsvd.m'\" \(548 characters\) sed "s/^X//" >'chapt10/lsqrsvd.m' <<'END_OF_FILE' Xfunction x = lsqrsvd(A,b) X%LSQRSVD Least-squares solutions using the SVD X%x = lsqrsvd(A,b) computes the least squares solution x using the X%the SVD of A. X%This program implements Algorithm 10.8.1 of the book. X%input : Matrix A and vector b X%output : vector x X X [m,n] = size(A); X y = zeros(n,1); X [U,S,V] = svd(A); X bprime = U' * b; X for i = 1:n X if S(i,i) > 10^(-15) X y(i) = bprime(i) / S(i,i); X else X y(i) = rand(1); X end; X end; X x = V * y; X end END_OF_FILE if test 548 -ne `wc -c <'chapt10/lsqrsvd.m'`; then echo shar: \"'chapt10/lsqrsvd.m'\" unpacked with wrong size! fi # end of 'chapt10/lsqrsvd.m' fi if test -f 'chapt10/minnmsvd.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chapt10/minnmsvd.m'\" else echo shar: Extracting \"'chapt10/minnmsvd.m'\" \(617 characters\) sed "s/^X//" >'chapt10/minnmsvd.m' <<'END_OF_FILE' Xfunction x = minnmsvd(A,b) X%MINNMSVD Minimum norm least-squares solution using the SVD X%x = minnmsvd(A,b) computes the minimum norm least squares X%solution x using the SVD of A. X%This program implements Algorithm 10.8.1 of the book to X%compute the minimum-norm least-squares solution. X%input : Matrix A and vector b X%output : vector x X X [m,n] = size(A); X y = zeros(n,1); X [U,S,V] = svd(A); X bprime = U' * b; X r = 0; X for i = 1:n X if S(i,i) > 10^(-15) X y(i) = bprime(i) / S(i,i); X else X y(i) = 0; X end; X end; X x = V * y; END_OF_FILE if test 617 -ne `wc -c <'chapt10/minnmsvd.m'`; then echo shar: \"'chapt10/minnmsvd.m'\" unpacked with wrong size! fi # end of 'chapt10/minnmsvd.m' fi if test -f 'chapt10/orthproj.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'chapt10/orthproj.m'\" else echo shar: Extracting \"'chapt10/orthproj.m'\" \(971 characters\) sed "s/^X//" >'chapt10/orthproj.m' <<'END_OF_FILE' Xfunction [R,N,OR,ON] = orthproj(A) X%ORTHPROJ Orthogonal projections using the SVD. X%[R,N,OR,ON] = orthproj(A) computes the orthogonal projections X%once the orthonormal bases for the range R(A) and the null-space X%N(A) of A are obtained, using the SVD of a : [ U, S, V] = SVD(A). X%R is the projection onto R(A) = U1 * U1' X%N is the projection onto N(A) = V2 * V2' X%OR is the projection onto the orthogonal complement of R(A) = U2 * U2' X%ON is the projection onto the orthogonal complement of N(A) = V1 * V1' X%see section 10.6.2 of the book X%input : Matrix A X%output : Matrices R, N, OR and ON X X [U,S,V] = svd(A); X [m,n] = size(A); X r = 0; X for i = 1:n X if S(i,i) > 10^(-15) X r = r + 1; X else X i = n; X end; X end; X U1 = U(:,1:r); X U2 = U(:,r+1:n); X V1 = V(:,1:r); X V2 = V(:,r+1:n); X R = U1 * U1'; X N = V2 * V2'; X OR = U2 * U2'; X ON = V1 * V1'; END_OF_FILE if test 971 -ne `wc -c <'chapt10/orthproj.m'`; then echo shar: \"'chapt10/orthproj.m'\" unpacked with wrong size! fi # end of 'chapt10/orthproj.m' fi if test ! -d 'nochapt' ; then echo shar: Creating directory \"'nochapt'\" mkdir 'nochapt' fi if test -f 'nochapt/aanochp.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'nochapt/aanochp.m'\" else echo shar: Extracting \"'nochapt/aanochp.m'\" \(34 characters\) sed "s/^X//" >'nochapt/aanochp.m' <<'END_OF_FILE' X% N O C H A P T E R END_OF_FILE if test 34 -ne `wc -c <'nochapt/aanochp.m'`; then echo shar: \"'nochapt/aanochp.m'\" unpacked with wrong size! fi # end of 'nochapt/aanochp.m' fi if test -f 'nochapt/absmax.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'nochapt/absmax.m'\" else echo shar: Extracting \"'nochapt/absmax.m'\" \(405 characters\) sed "s/^X//" >'nochapt/absmax.m' <<'END_OF_FILE' X function [p] = absmax(y); X%ABSMAX The position p in a vector where the absolute maximum occurs. X%p = absmax(y) returns the postion p in the vector y X%where the absolute maximum of the vector y occurs. X%input : Vector y X%output : Integer p X X [m,n] = size(y); X s = y(1); X n1 = 1; X for k = 1:m X if abs(y(k)) > abs(s) X p = k; X s = y(k); X end; X end; X end; END_OF_FILE if test 405 -ne `wc -c <'nochapt/absmax.m'`; then echo shar: \"'nochapt/absmax.m'\" unpacked with wrong size! fi # end of 'nochapt/absmax.m' fi if test -f 'nochapt/cputime.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'nochapt/cputime.m'\" else echo shar: Extracting \"'nochapt/cputime.m'\" \(576 characters\) sed "s/^X//" >'nochapt/cputime.m' <<'END_OF_FILE' X function cputim1 = cputime X%CPUTIME Computation of the CPU time. X%cputim1 = cputime returns the cputime in seconds that X%has been used by the program since the program started. X%example : X% t1 = cputime X% your-operation X% t1 = cputime - t1 X%returns the cputime used to run your-operation. X%Since this PC version of MATLAB does not have a cputime X%I have written a program called cputime X%that gives the time. X X t1 = clock; X sum = t1(6); X sum = sum + t1(5) * 60; X sum = sum + t1(4) * 3600; X cputim1 = sum; X end; X END_OF_FILE if test 576 -ne `wc -c <'nochapt/cputime.m'`; then echo shar: \"'nochapt/cputime.m'\" unpacked with wrong size! fi # end of 'nochapt/cputime.m' fi if test -f 'nochapt/inter.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'nochapt/inter.m'\" else echo shar: Extracting \"'nochapt/inter.m'\" \(298 characters\) sed "s/^X//" >'nochapt/inter.m' <<'END_OF_FILE' X function [y,z] = inter(y,z); X%INTER Interchange two vectors X%[y,z] = inter(y,z) interchanges two vectors y and z X%so that on output y contains z and z contains y. X%input : Vectors y and z X%output : Vectors y and z X X c = y ; X y = z; X z = c; X X X X end; X END_OF_FILE if test 298 -ne `wc -c <'nochapt/inter.m'`; then echo shar: \"'nochapt/inter.m'\" unpacked with wrong size! fi # end of 'nochapt/inter.m' fi echo shar: End of shell archive. exit 0