diff git a/books/bookvol10.4.pamphlet b/books/bookvol10.4.pamphlet
index f5526b9..addd4fa 100644
 a/books/bookvol10.4.pamphlet
+++ b/books/bookvol10.4.pamphlet
@@ 77519,929 +77519,1877 @@ MultivariateSquareFree (E,OV,R,P) : C == T where
)set message auto off
)clear all
S 1 of 1
+S 1 of 164
)show NagEigenPackage
+R
+R NagEigenPackage is a package constructor
+R Abbreviation for NagEigenPackage is NAGF02
+R This constructor is exposed in this frame.
+R Issue )edit bookvol10.4.pamphlet to see algebra source code for NAGF02
+R
+R Operations 
+R f02aaf : (Integer,Integer,Matrix(DoubleFloat),Integer) > Result
+R f02abf : (Matrix(DoubleFloat),Integer,Integer,Integer,Integer) > Result
+R f02adf : (Integer,Integer,Integer,Matrix(DoubleFloat),Matrix(DoubleFloat),Integer) > Result
+R f02aef : (Integer,Integer,Integer,Integer,Matrix(DoubleFloat),Matrix(DoubleFloat),Integer) > Result
+R f02aff : (Integer,Integer,Matrix(DoubleFloat),Integer) > Result
+R f02agf : (Integer,Integer,Integer,Integer,Matrix(DoubleFloat),Integer) > Result
+R f02ajf : (Integer,Integer,Integer,Matrix(DoubleFloat),Matrix(DoubleFloat),Integer) > Result
+R f02akf : (Integer,Integer,Integer,Integer,Integer,Matrix(DoubleFloat),Matrix(DoubleFloat),Integer) > Result
+R f02awf : (Integer,Integer,Integer,Matrix(DoubleFloat),Matrix(DoubleFloat),Integer) > Result
+R f02axf : (Matrix(DoubleFloat),Integer,Matrix(DoubleFloat),Integer,Integer,Integer,Integer,Integer) > Result
+R f02bbf : (Integer,Integer,DoubleFloat,DoubleFloat,Integer,Integer,Matrix(DoubleFloat),Integer) > Result
+R f02bjf : (Integer,Integer,Integer,DoubleFloat,Boolean,Integer,Matrix(DoubleFloat),Matrix(DoubleFloat),Integer) > Result
+R f02fjf : (Integer,Integer,DoubleFloat,Integer,Integer,Integer,Integer,Integer,Integer,Integer,Matrix(DoubleFloat),Integer,Union(fn: FileName,fp: Asp27(DOT)),Union(fn: FileName,fp: Asp28(IMAGE))) > Result
+R f02fjf : (Integer,Integer,DoubleFloat,Integer,Integer,Integer,Integer,Integer,Integer,Integer,Matrix(DoubleFloat),Integer,Union(fn: FileName,fp: Asp27(DOT)),Union(fn: FileName,fp: Asp28(IMAGE)),FileName) > Result
+R f02wef : (Integer,Integer,Integer,Integer,Integer,Boolean,Integer,Boolean,Integer,Matrix(DoubleFloat),Matrix(DoubleFloat),Integer) > Result
+R f02xef : (Integer,Integer,Integer,Integer,Integer,Boolean,Integer,Boolean,Integer,Matrix(Complex(DoubleFloat)),Matrix(Complex(DoubleFloat)),Integer) > Result
+R
E 1
)spool
)lisp (bye)
\end{chunk}
\begin{chunk}{NagEigenPackage.help}
+)clear all
This package uses the NAG Library to compute
 * eigenvalues and eigenvectors of a matrix\
 * eigenvalues and eigenvectors of generalized matrix
 * eigenvalue problems
 * singular values and singular vectors of a matrix.
+S 2 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 2
 F02  Eigenvalues and Eigenvectors Introduction  F02
 Chapter F02
 Eigenvalues and Eigenvectors
+S 3 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 3
 1. Scope of the Chapter
+S 4 of 164
+ia:=4
+R
+R
+R (3) 4
+R Type: PositiveInteger
+E 4
 This chapter is concerned with computing
+S 5 of 164
+n:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 5
  eigenvalues and eigenvectors of a matrix
+S 6 of 164
+a:Matrix SF:=
+ [[ 0.5 , 0.0 , 2.3 ,2.6 ],_
+ [ 0.0 , 0.5 ,1.4 ,0.7 ],_
+ [ 2.3 ,1.4 , 0.5 , 0.0 ],_
+ [2.6 ,0.7 , 0.0 , 0.5 ]]
+R
+R
+R (5)
+R [[0.5,0.,2.2999999999999998, 2.5999999999999996],
+R [0.,0.5, 1.3999999999999999, 0.69999999999999996],
+R [2.2999999999999998, 1.3999999999999999,0.5,0.],
+R [ 2.5999999999999996, 0.69999999999999996,0.,0.5]]
+R Type: Matrix(DoubleFloat)
+E 6
  eigenvalues and eigenvectors of generalized matrix
 eigenvalue problems
+S 7 of 164
+ result:=f02aaf(ia,n,a,1)
+E 7
  singular values and singular vectors of a matrix.
+)clear all
 2. Background to the Problems
+S 8 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 8
 2.1. Eigenvalue Problems
+S 9 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 9
 In the most usual form of eigenvalue problem we are given a
 square n by n matrix A and wish to compute (lambda) (an
 eigenvalue) and x/=0 (an eigenvector) which satisfy the equation
+S 10 of 164
+a:Matrix SF:=
+ [[ 0.5 , 0.0 , 2.3 ,2.6 ],_
+ [ 0.0 , 0.5 ,1.4 ,0.7 ],_
+ [ 2.3 ,1.4 , 0.5 , 0.0 ],_
+ [2.6 ,0.7 , 0.0 , 0.5 ]]
+R
+R
+R (3)
+R [[0.5,0.,2.2999999999999998, 2.5999999999999996],
+R [0.,0.5, 1.3999999999999999, 0.69999999999999996],
+R [2.2999999999999998, 1.3999999999999999,0.5,0.],
+R [ 2.5999999999999996, 0.69999999999999996,0.,0.5]]
+R Type: Matrix(DoubleFloat)
+E 10
 Ax=(lambda)x
+S 11 of 164
+ia:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 11
 Such problems are called 'standard' eigenvalue problems in
 contrast to 'generalized' eigenvalue problems where we wish to
 satisfy the equation
+S 12 of 164
+n:=4
+R
+R
+R (5) 4
+R Type: PositiveInteger
+E 12
 Ax=(lambda)Bx
+S 13 of 164
+iv:=4
+R
+R
+R (6) 4
+R Type: PositiveInteger
+E 13
 B also being a square n by n matrix.
+S 14 of 164
+ result:=f02abf(a,ia,n,iv,1)
+E 14
 Section 2.1.1 and Section 2.1.2 discuss, respectively, standard
 and generalized eigenvalue problems where the matrices involved
 are dense; Section 2.1.3 discusses both types of problem in the
 case where A and B are sparse (and symmetric).
+)clear all
 2.1.1. Standard eigenvalue problems
+S 15 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 15
 Some of the routines in this chapter find all the n eigenvalues,
 some find all the n eigensolutions (eigenvalues and corresponding
 eigenvectors), and some find a selected group of eigenvalues
 and/or eigenvectors. The matrix A may be:
+S 16 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 16
 (i) general (real or complex)
+S 17 of 164
+ia:=4
+R
+R
+R (3) 4
+R Type: PositiveInteger
+E 17
 (ii) real symmetric, or
+S 18 of 164
+ib:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 18
 (iii) complex Hermitian (so that if a =(alpha)+i(beta) then
 ij
 a =(alpha)i(beta)).
 ji
+S 19 of 164
+n:=4
+R
+R
+R (5) 4
+R Type: PositiveInteger
+E 19
 In all cases the computation starts with a similarity
 1
 transformation S AS=T, where S is nonsingular and is the
 product of fairly simple matrices, and T has an 'easier form'
 than A so that its eigensolutions are easily determined. The
 matrices A and T, of course, have the same eigenvalues, and if y
 is an eigenvector of T then Sy is the corresponding eigenvector
 of A.
+S 20 of 164
+a:Matrix SF:=
+ [[0.5 , 1.5 , 6.6 , 4.8 ],_
+ [1.5 , 6.5 ,16.2 , 8.6 ],_
+ [6.6 ,16.2 ,37.6 , 9.8 ],_
+ [4.8 , 8.6 , 9.8 ,17.1 ]]
+R
+R
+R (6)
+R [[0.5,1.5,6.5999999999999996,4.7999999999999998],
+R [1.5,6.5,16.199999999999999,8.5999999999999996],
+R
+R [6.5999999999999996, 16.199999999999999, 37.599999999999994,
+R 9.8000000000000007]
+R ,
+R
+R [4.7999999999999998, 8.5999999999999996, 9.8000000000000007,
+R  17.099999999999998]
+R ]
+R Type: Matrix(DoubleFloat)
+E 20
 In case (i) (general real or complex A), the selected form of T
 is an upper Hessenberg matrix (t =0 if ij>1) and S is the
 ij
 product of n2 stabilised elementary transformation matrices.
 There is no easy method of computing selected eigenvalues of a
 Hessenberg matrix, so that all eigenvalues are always calculated.
 In the real case this computation is performed via the Francis QR
 algorithm with double shifts, and in the complex case by means of
 the LR algorithm. If the eigenvectors are required they are
 computed by backsubstitution following the QR and LR algorithm.
+S 21 of 164
+b:Matrix SF:=
+ [[1 , 3 , 4 , 1 ],_
+ [3 ,13 ,16 ,11 ],_
+ [4 ,16 ,24 ,18 ],_
+ [1 ,11 ,18 ,27 ]]
+R
+R
+R +1. 3. 4. 1. +
+R  
+R 3. 13. 16. 11.
+R (7)  
+R 4. 16. 24. 18.
+R  
+R +1. 11. 18. 27.+
+R Type: Matrix(DoubleFloat)
+E 21
 In case (ii) (real and symmetric A) the selected simple form of T
 is a tridiagonal matrix (t =0 if ij>1), and S is the product
 ij
 of n2 orthogonal Householder transformation matrices. If only
 selected eigenvalues are required, they are obtained by the
 method of bisection using the Sturm sequence property, and the
 corresponding eigenvectors of T are computed by inverse
 iteration. If all eigenvalues are required, they are computed
 from T via the QL algorithm (an adaptation of the QR algorithm),
 and the corresponding eigenvectors of T are the product of the
 transformations for the QL reduction. In all cases the
 corresponding eigenvectors of A are recovered from the
 computation of x=Sy.
+S 22 of 164
+ result:=f02adf(ia,ib,n,a,b,1)
+E 22
 In case (iii) (complex Hermitian A) analogous transformations as
 in case (ii) are used. T has complex elements in offdiagonal
 positions, but a simple diagonal similarity transformation is
 then used to produce a real tridiagonal form, after which the QL
 algorithm and succeeding methods described in the previous
 paragraph are used to complete the solution.
+)clear all
 2.1.2. Generalized eigenvalue problems
+S 23 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 23
 Here we distinguish as a special case those problems in which
 both A and B are symmetric and B is positivedefinite and well
 conditioned with respect to inversion (i.e., all the eigenvalues
 of B are significantly greater than zero). Such problems can be
 satisfactorily treated by first reducing them to case (ii) of
 Section 2.1.1 and then using the methods described there to
 T
 compute the eigensolutions. If B is factorized as LL (L lower
 triangular), then Ax=(lambda)Bx is equivalent to the standard
 1 T 1 T
 symmetric problem Ry=(lambda)y, where R=L A(L ) and y=L x.
 After finding an eigenvector y of R, the required x is computed
 T
 by backsubstitution in y=L x.
+S 24 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 24
 For generalized problems of the form Ax=(lambda)Bx which do not
 fall into the special case, the QZ algorithm is provided.
+S 25 of 164
+ia:=4
+R
+R
+R (3) 4
+R Type: PositiveInteger
+E 25
 In order to appreciate the domain in which this algorithm is
 appropriate we remark first that when B is nonsingular the
 problem Ax=(lambda)Bx is fully equivalent to the problem
 1
 (B A)x=(lambda)x; both the eigenvalues and eigenvectors being
 the same. When A is nonsingular Ax=(lambda)Bx is equivalent to
 1
 the problem (A B)x=(mu)x; the eigenvalues (mu) being the
 reciprocals of the required eigenvalues and the eigenvectors
 remaining the same. In theory then, provided at least one of the
 matrices A and B is nonsingular, the generalized problem
 Ax=(lambda)Bx could be solved via the standard problem
 Cx=(lambda)x with an appropriate matrix C, and as far as economy
 of effort is concerned this is quite satisfactory. However, in
 practice, for this reduction to be satisfactory from the
 standpoint of numerical stability, one requires more than the
 1
 mere nonsingularity of A or B. It is necessary that B A (or
 1
 A B) should not only exist but that B (or A) should be well
 conditioned with respect to inversion. The nearer B (or A) is to
 1 1
 singularity the more unsatisfactory B A (or A B) will be as a
 vehicle for determining the required eigenvalues. Unfortunately
 1
 one cannot counter illconditioning in B (or A) by computing B A
 1
 (or A B) accurately to single precision using iterative
 refinement. Welldetermined eigenvalues of the original
 Ax=(lambda)Bx may be poorly determined even by the correctly
 1 1
 rounded version of B A (or A B). The situation may in some
 instances be saved by the observation that if Ax=(lambda)Bx then
 (AkB)x=((lambda)k)Bx. Hence if AkB is nonsingular we may
 1
 solve the standard problem [(AkB) B]x=(mu)x and for numerical
 stability we require only that (AkB) be wellconditioned with
 respect to inversion.

 In practice one may well be in a situation where no k is known
 for which (AkB) is wellconditioned with respect to inversion
 and indeed (AkB) may be singular for all k. The QZ algorithm is
 designed to deal directly with the problem Ax=(lambda)Bx itself
 and its performance is unaffected by singularity or near
 singularity of A, B or AkB.

 2.1.3. Sparse symmetric problems
+S 26 of 164
+ib:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 26
 If the matrices A and B are large and sparse (i.e., only a small
 proportion of the elements are nonzero), then the methods
 described in the previous Section are unsuitable, because in
 reducing the problem to a simpler form, much of the sparsity of
 the problem would be lost; hence the computing time and the
 storage required would be very large. Instead, for symmetric
 problems, the method of simultaneous iteration may be used to
 determine selected eigenvalues and the corresponding
 eigenvectors. The routine provided has been designed to handle
 both symmetric and generalized symmetric problems.
+S 27 of 164
+n:=4
+R
+R
+R (5) 4
+R Type: PositiveInteger
+E 27
 2.2. Singular Value Problems
+S 28 of 164
+iv:=4
+R
+R
+R (6) 4
+R Type: PositiveInteger
+E 28
 The singular value decomposition of an m by n real matrix A is
 given by
+S 29 of 164
+a:Matrix SF:=
+ [[0.5 , 1.5 , 6.6 , 4.8 ],_
+ [1.5 , 6.5 ,16.2 , 8.6 ],_
+ [6.6 ,16.2 ,37.6 , 9.8 ],_
+ [4.8 , 8.6 , 9.8 ,17.1 ]]
+R
+R
+R (7)
+R [[0.5,1.5,6.5999999999999996,4.7999999999999998],
+R [1.5,6.5,16.199999999999999,8.5999999999999996],
+R
+R [6.5999999999999996, 16.199999999999999, 37.599999999999994,
+R 9.8000000000000007]
+R ,
+R
+R [4.7999999999999998, 8.5999999999999996, 9.8000000000000007,
+R  17.099999999999998]
+R ]
+R Type: Matrix(DoubleFloat)
+E 29
 T
 A=QDP ,
+S 30 of 164
+b:Matrix SF:=
+ [[1 , 3 , 4 , 1 ],_
+ [3 ,13 ,16 ,11 ],_
+ [4 ,16 ,24 ,18 ],_
+ [1 ,11 ,18 ,27 ]]
+R
+R
+R +1. 3. 4. 1. +
+R  
+R 3. 13. 16. 11.
+R (8)  
+R 4. 16. 24. 18.
+R  
+R +1. 11. 18. 27.+
+R Type: Matrix(DoubleFloat)
+E 30
 where Q is an m by m orthogonal matrix, P is an n by n orthogonal
 matrix and D is an m by n diagonal matrix with nonnegative
 diagonal elements. The first k==min(m,n) columns of Q and P are
 the left and righthand singular vectors of A and the k diagonal
 elements of D are the singular values.
+S 31 of 164
+ result:=f02aef(ia,ib,n,iv,a,b,1)
+E 31
 When A is complex then the singular value decomposition is given
 by
+)clear all
 H
 A=QDP ,
+S 32 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 32
 H T
 where Q and P are unitary, P denotes the complex conjugate of P
 and D is as above for the real case.
+S 33 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 33
 If the matrix A has column means of zero, then AP is the matrix
 of principal components of A and the singular values are the
 square roots of the sample variances of the observations with
 respect to the principal components. (See also Chapter G03.)
+S 34 of 164
+ia:=4
+R
+R
+R (3) 4
+R Type: PositiveInteger
+E 34
 Routines are provided to return the singular values and vectors
 of a general real or complex matrix.
+S 35 of 164
+n:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 35
 3. Recommendations on Choice and Use of Routines
+S 36 of 164
+a:Matrix SF:=
+ [[ 1.5 ,0.1 , 4.5 ,1.5 ],_
+ [22.5 ,3.5 ,12.5 ,2.5 ],_
+ [ 2.5 ,0.3 , 4.5 ,2.5 ],_
+ [ 2.5 ,0.1 , 4.5 , 2.5 ]]
+R
+R
+R + 1.5 0.10000000000000001 4.5  1.5+
+R  
+R  22.5 3.5 12.5  2.5
+R (5)  
+R  2.5 0.29999999999999999 4.5  2.5
+R  
+R + 2.5 0.10000000000000001 4.5 2.5 +
+R Type: Matrix(DoubleFloat)
+E 36
 3.1. General Discussion
+S 37 of 164
+ result:=f02aff(ia,n,a,1)
+E 37
 There is one routine, F02FJF, which is designed for sparse
 symmetric eigenvalue problems, either standard or generalized.
 The remainder of the routines are designed for dense matrices.
+)clear all
 3.2. Eigenvalue and Eigenvector Routines
+S 38 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 38
 These reduce the matrix A to a simpler form by a similarity
 1
 transformation S AS=T where T is an upper Hessenberg or
 tridiagonal matrix, compute the eigensolutions of T, and then
 recover the eigenvectors of A via the matrix S. The eigenvectors
 are normalised so that
+S 39 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 39
 n
  2
 > x  =1
  r
 r=1
+S 40 of 164
+ia:=4
+R
+R
+R (3) 4
+R Type: PositiveInteger
+E 40
 x being the rth component of the eigenvector x, and so that the
 r
 element of largest modulus is real if x is complex. For problems
 of the type Ax=(lambda)Bx with A and B symmetric and B positive
 T
 definite, the eigenvectors are normalised so that x Bx=1, x
 always being real for such problems.
+S 41 of 164
+n:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 41
 3.3. Singular Value and Singular Vector Routines
+S 42 of 164
+ivr:=4
+R
+R
+R (5) 4
+R Type: PositiveInteger
+E 42
 These reduce the matrix A to real bidiagonal form, B say, by
 T
 orthogonal transformations Q AP=B in the real case, and by
 H
 unitary transformations Q AP=B in the complex case, and the
 singular values and vectors are computed via this bidiagonal
 form. The singular values are returned in descending order.
+S 43 of 164
+ivi:=4
+R
+R
+R (6) 4
+R Type: PositiveInteger
+E 43
 3.4. Decision Trees
+S 44 of 164
+a:Matrix SF:=
+ [[ 1.5 ,0.1 , 4.5 ,1.5 ],_
+ [22.5 ,3.5 ,12.5 ,2.5 ],_
+ [ 2.5 ,0.3 , 4.5 ,2.5 ],_
+ [ 2.5 ,0.1 , 4.5 , 2.5 ]]
+R
+R
+R + 1.5 0.10000000000000001 4.5  1.5+
+R  
+R  22.5 3.5 12.5  2.5
+R (7)  
+R  2.5 0.29999999999999999 4.5  2.5
+R  
+R + 2.5 0.10000000000000001 4.5 2.5 +
+R Type: Matrix(DoubleFloat)
+E 44
 (i) Eigenvalues and Eigenvectors
+S 45 of 164
+ result:=f02agf(ia,n,ivr,ivi,a,1)
+E 45
+)clear all
 Please see figure in printed Reference Manual
+S 46 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 46
+S 47 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 47
 (ii) Singular Values and Singular Vectors
+S 48 of 164
+iar:=4
+R
+R
+R (3) 4
+R Type: PositiveInteger
+E 48
+S 49 of 164
+iai:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 49
 Please see figure in printed Reference Manual
+S 50 of 164
+n:=4
+R
+R
+R (5) 4
+R Type: PositiveInteger
+E 50
 F02  Eigenvalues and Eigenvectors Contents  F02
 Chapter F02
+S 51 of 164
+ar:Matrix SF:=
+ [[21.0 , 0.0 ,13.6 ,0.0 ],_
+ [ 0.0 ,26.0 , 7.5 ,2.5 ],_
+ [ 2.0 , 1.68, 4.5 ,1.5 ],_
+ [ 0.0 ,2.6 ,2.7 ,2.5 ]]
+R
+R
+R + 21. 0. 13.6 0. +
+R  
+R  0. 26. 7.5 2.5
+R (6)  
+R  2. 1.6799999999999999 4.5 1.5
+R  
+R + 0.  2.5999999999999996  2.6999999999999997 2.5+
+R Type: Matrix(DoubleFloat)
+E 51
 Eigenvalues and Eigenvectors
+S 52 of 164
+ai:Matrix SF:=
+ [[5.0 ,24.6 , 10.2 ,4.0 ],_
+ [22.5 ,5.0 , 10.0 ,0.0 ],_
+ [ 1.5 , 2.24 ,5.0 , 2.0 ],_
+ [2.5 , 0.0 , 3.6 ,5.0 ]]
+R
+R
+R + 5. 24.600000000000001 10.199999999999999 4. +
+R  
+R 22.5  5.  10. 0. 
+R (7)  
+R  1.5 2.2400000000000002  5. 2. 
+R  
+R + 2.5 0. 3.5999999999999996  5.+
+R Type: Matrix(DoubleFloat)
+E 52
 F02AAF All eigenvalues of real symmetric matrix
+S 53 of 164
+ result:=f02ajf(iar,iai,n,ar,ai,1)
+E 53
 F02ABF All eigenvalues and eigenvectors of real symmetric matrix
+)clear all
 F02ADF All eigenvalues of generalized real symmetricdefinite
 eigenproblem
+S 54 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 54
 F02AEF All eigenvalues and eigenvectors of generalized real
 symmetricdefinite eigenproblem
+S 55 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 55
 F02AFF All eigenvalues of real matrix
+S 56 of 164
+iar:=4
+R
+R
+R (3) 4
+R Type: PositiveInteger
+E 56
 F02AGF All eigenvalues and eigenvectors of real matrix
+S 57 of 164
+iai:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 57
 F02AJF All eigenvalues of complex matrix
+S 58 of 164
+n:=4
+R
+R
+R (5) 4
+R Type: PositiveInteger
+E 58
 F02AKF All eigenvalues and eigenvectors of complex matrix
+S 59 of 164
+ivr:=4
+R
+R
+R (6) 4
+R Type: PositiveInteger
+E 59
 F02AWF All eigenvalues of complex Hermitian matrix
+S 60 of 164
+ivi:=4
+R
+R
+R (7) 4
+R Type: PositiveInteger
+E 60
 F02AXF All eigenvalues and eigenvectors of complex Hermitian
 matrix
+S 61 of 164
+ar:Matrix SF:=
+ [[21.0 , 0.0 ,13.6 ,0.0 ],_
+ [ 0.0 ,26.0 , 7.5 ,2.5 ],_
+ [ 2.0 , 1.68 ,4.5 ,1.5 ],_
+ [ 0.0 ,2.6 ,2.7 ,2.5 ]]
+R
+R
+R + 21. 0. 13.6 0. +
+R  
+R  0. 26. 7.5 2.5
+R (8)  
+R  2. 1.6799999999999999 4.5 1.5
+R  
+R + 0.  2.5999999999999996  2.6999999999999997 2.5+
+R Type: Matrix(DoubleFloat)
+E 61
 F02BBF Selected eigenvalues and eigenvectors of real symmetric
 matrix
+S 62 of 164
+ai:Matrix SF:=
+ [[5.0 ,24.6 ,10.2 , 4.0 ],_
+ [22.5 ,5.0 ,10.0 , 0.0 ],_
+ [ 1.5 , 2.24 ,5.0 , 2.0 ],_
+ [2.5 , 0.0 , 3.6 ,5.0 ]]
+R
+R
+R + 5. 24.600000000000001 10.199999999999999 4. +
+R  
+R 22.5  5.  10. 0. 
+R (9)  
+R  1.5 2.2400000000000002  5. 2. 
+R  
+R + 2.5 0. 3.5999999999999996  5.+
+R Type: Matrix(DoubleFloat)
+E 62
 F02BJF All eigenvalues and optionally eigenvectors of
 generalized eigenproblem by QZ algorithm, real matrices
+S 63 of 164
+ result:=f02akf(iar,iai,n,ivr,ivi,ar,ai,1)
+E 63
 F02FJF Selected eigenvalues and eigenvectors of sparse symmetric
 eigenproblem
+)clear all
 F02WEF SVD of real matrix
+S 64 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 64
 F02XEF SVD of complex matrix
+S 65 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 65
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+S 66 of 164
+iar:=4
+R
+R
+R (3) 4
+R Type: PositiveInteger
+E 66
 F02  Eigenvalue and Eigenvectors F02AAF
 F02AAF  NAG Foundation Library Routine Document
+S 67 of 164
+iai:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 67
 Note: Before using this routine, please read the Users' Note for
 your implementation to check implementationdependent details.
 The symbol (*) after a NAG routine name denotes a routine that is
 not included in the Foundation Library.
+S 68 of 164
+n:=4
+R
+R
+R (5) 4
+R Type: PositiveInteger
+E 68
 1. Purpose
+S 69 of 164
+ar:Matrix SF:=
+ [[0.5 , 0.0 , 1.84 , 2.08 ],_
+ [0.0 , 0.5 , 1.12 ,0.56 ],_
+ [1.84 , 1.12 ,0.5 , 0.0 ],_
+ [2.08 ,0.56 ,0.0 , 0.5 ]]
+R
+R
+R (6)
+R [[0.5,0.,1.8399999999999999,2.0800000000000001],
+R [0.,0.5,1.1200000000000001, 0.55999999999999994],
+R [1.8399999999999999,1.1200000000000001,0.5,0.],
+R [2.0800000000000001, 0.55999999999999994,0.,0.5]]
+R Type: Matrix(DoubleFloat)
+E 69
 F02AAF calculates all the eigenvalues of a real symmetric matrix.
+S 70 of 164
+ai:Matrix SF:=
+ [[ 0.0 , 0.0 , 1.38 ,1.56 ],_
+ [ 0.0 , 0.0 , 0.84 , 0.42 ],_
+ [1.38 ,0.84 ,0.0 , 0.0 ],_
+ [ 1.56 ,0.42 ,0.0 , 0.0 ]]
+R
+R
+R (7)
+R [[0.,0.,1.3799999999999999, 1.5599999999999998],
+R [0.,0.,0.83999999999999997,0.41999999999999998],
+R [ 1.3799999999999999, 0.84000000000000008,0.,0.],
+R [1.5600000000000001, 0.42000000000000004,0.,0.]]
+R Type: Matrix(DoubleFloat)
+E 70
 2. Specification
+S 71 of 164
+ result:=f02awf(iar,iai,n,ar,ai,1)
+E 71
 SUBROUTINE F02AAF (A, IA, N, R, E, IFAIL)
 INTEGER IA, N, IFAIL
 DOUBLE PRECISION A(IA,N), R(N), E(N)
+)clear all
 3. Description
+S 72 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 72
 This routine reduces the real symmetric matrix A to a real
 symmetric tridiagonal matrix using Householder's method. The
 eigenvalues of the tridiagonal matrix are then determined using
 the QL algorithm.
+S 73 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 73
 4. References
+S 74 of 164
+ar:Matrix SF:=
+ [[0.5 ,0.0 ,1.84 ,2.08 ],_
+ [0.0 ,0.5 ,1.12 ,0.56 ],_
+ [1.84 ,1.12 ,0.5 ,0.0 ],_
+ [2.08 ,0.56 ,0.0 ,0.5 ]]
+R
+R
+R (3)
+R [[0.5,0.,1.8399999999999999,2.0800000000000001],
+R [0.,0.5,1.1200000000000001, 0.55999999999999994],
+R [1.8399999999999999,1.1200000000000001,0.5,0.],
+R [2.0800000000000001, 0.55999999999999994,0.,0.5]]
+R Type: Matrix(DoubleFloat)
+E 74
 [1] Wilkinson J H and Reinsch C (1971) Handbook for Automatic
 Computation II, Linear Algebra. SpringerVerlag.
+S 75 of 164
+iar:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 75
 5. Parameters
+S 76 of 164
+ai:Matrix SF:=
+ [[0.0 ,0.0 ,1.38 ,1.56 ],_
+ [0.0 ,0.0 ,0.84 ,0.42 ],_
+ [1.38 ,0.84 ,0.0 ,0.0 ],_
+ [1.56 ,0.42 ,0.0 ,0.0 ]]
+R
+R
+R (5)
+R [[0.,0.,1.3799999999999999, 1.5599999999999998],
+R [0.,0.,0.83999999999999997,0.41999999999999998],
+R [ 1.3799999999999999, 0.84000000000000008,0.,0.],
+R [1.5600000000000001, 0.42000000000000004,0.,0.]]
+R Type: Matrix(DoubleFloat)
+E 76
 1: A(IA,N)  DOUBLE PRECISION array Input/Output
 On entry: the lower triangle of the n by n symmetric matrix
 A. The elements of the array above the diagonal need not be
 set. On exit: the elements of A below the diagonal are
 overwritten, and the rest of the array is unchanged.
+S 77 of 164
+iai:=4
+R
+R
+R (6) 4
+R Type: PositiveInteger
+E 77
 2: IA  INTEGER Input
 On entry:
 the first dimension of the array A as declared in the
 (sub)program from which F02AAF is called.
 Constraint: IA >= N.
+S 78 of 164
+n:=4
+R
+R
+R (7) 4
+R Type: PositiveInteger
+E 78
 3: N  INTEGER Input
 On entry: n, the order of the matrix A.
+S 79 of 164
+ivr:=4
+R
+R
+R (8) 4
+R Type: PositiveInteger
+E 79
 4: R(N)  DOUBLE PRECISION array Output
 On exit: the eigenvalues in ascending order.
+S 80 of 164
+ivi:=4
+R
+R
+R (9) 4
+R Type: PositiveInteger
+E 80
 5: E(N)  DOUBLE PRECISION array Workspace
+S 81 of 164
+ result:=f02axf(ar,iar,ai,iai,n,ivr,ivi,1)
+E 81
 6: IFAIL  INTEGER Input/Output
 On entry: IFAIL must be set to 0, 1 or 1. For users not
 familiar with this parameter (described in the Essential
 Introduction) the recommended value is 0.
+)clear all
 On exit: IFAIL = 0 unless the routine detects an error (see
 Section 6).
+S 82 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 82
 6. Error Indicators and Warnings
+S 83 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 83
 Errors detected by the routine:
+S 84 of 164
+ia:=4
+R
+R
+R (3) 4
+R Type: PositiveInteger
+E 84
 IFAIL= 1
 Failure in F02AVF(*) indicating that more than 30*N
 iterations are required to isolate all the eigenvalues.
+S 85 of 164
+n:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 85
 7. Accuracy
+S 86 of 164
+alb:=2.0
+R
+R
+R (5)  2.0
+R Type: Float
+E 86
 The accuracy of the eigenvalues depends on the sensitivity of the
 matrix to rounding errors produced in tridiagonalisation. For a
 detailed error analysis see Wilkinson and Reinsch [1] pp 222 and
 235.
+S 87 of 164
+ub:=3.0
+R
+R
+R (6) 3.0
+R Type: Float
+E 87
 8. Further Comments
+S 88 of 164
+m:=3
+R
+R
+R (7) 3
+R Type: PositiveInteger
+E 88
 3
 The time taken by the routine is approximately proportional to n
+S 89 of 164
+iv:=4
+R
+R
+R (8) 4
+R Type: PositiveInteger
+E 89
 9. Example
+S 90 of 164
+a:Matrix SF:=
+ [[0.5 ,0.0 ,2.3 ,2.6 ],_
+ [0.0 ,0.5 ,1.4 ,0.7 ],_
+ [2.3 ,1.4 ,0.5 ,0.0 ],_
+ [2.6 ,0.7 ,0.0 ,0.5 ]]
+R
+R
+R (9)
+R [[0.5,0.,2.2999999999999998, 2.5999999999999996],
+R [0.,0.5, 1.3999999999999999, 0.69999999999999996],
+R [2.2999999999999998, 1.3999999999999999,0.5,0.],
+R [ 2.5999999999999996, 0.69999999999999996,0.,0.5]]
+R Type: Matrix(DoubleFloat)
+E 90
 To calculate all the eigenvalues of the real symmetric matrix:
+S 91 of 164
+ result:=f02bbf(ia,n,alb,ub,m,iv,a,1)
+E 91
 ( 0.5 0.0 2.3 2.6)
 ( 0.0 0.5 1.4 0.7)
 ( 2.3 1.4 0.5 0.0).
 (2.6 0.7 0.0 0.5)
+)clear all
 The example program is not reproduced here. The source code for
 all example programs is distributed with the NAG Foundation
 Library software and should be available online.
+S 92 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 92
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+S 93 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 93
 F02  Eigenvalue and Eigenvectors F02ABF
 F02ABF  NAG Foundation Library Routine Document

 Note: Before using this routine, please read the Users' Note for
 your implementation to check implementationdependent details.
 The symbol (*) after a NAG routine name denotes a routine that is
 not included in the Foundation Library.
+S 94 of 164
+n:=4
+R
+R
+R (3) 4
+R Type: PositiveInteger
+E 94
 1. Purpose
+S 95 of 164
+ia:=4
+R
+R
+R (4) 4
+R Type: PositiveInteger
+E 95
 F02ABF calculates all the eigenvalues and eigenvectors of a real
 symmetric matrix.
+S 96 of 164
+ib:=4
+R
+R
+R (5) 4
+R Type: PositiveInteger
+E 96
 2. Specification
+S 97 of 164
+eps1:SF:=1.0e4
+R
+R
+R (6) 9.9999999999999991E5
+R Type: DoubleFloat
+E 97
 SUBROUTINE F02ABF (A, IA, N, R, V, IV, E, IFAIL)
 INTEGER IA, N, IV, IFAIL
 DOUBLE PRECISION A(IA,N), R(N), V(IV,N), E(N)
+S 98 of 164
+matv:=true
+R
+R
+R (7) true
+R Type: Boolean
+E 98
 3. Description
+S 99 of 164
+iv:=4
+R
+R
+R (8) 4
+R Type: PositiveInteger
+E 99
 This routine reduces the real symmetric matrix A to a real
 symmetric tridiagonal matrix by Householder's method. The
 eigenvalues and eigenvectors are calculated using the QL
 algorithm.
+S 100 of 164
+a:Matrix SF:=
+ [[3.9 ,12.5 ,34.5 ,0.5 ],_
+ [4.3 ,21.5 ,47.5 ,7.5 ],_
+ [4.3 ,21.5 ,43.5 ,3.5 ],_
+ [4.4 ,26.0 ,46.0 ,6.0 ]]
+R
+R
+R +3.8999999999999999 12.5  34.5  0.5+
+R  
+R 4.2999999999999998 21.5  47.5 7.5 
+R (9)  
+R 4.2999999999999998 21.5  43.5 3.5 
+R  
+R +4.4000000000000004 26.  46. 6. +
+R Type: Matrix(DoubleFloat)
+E 100
 4. References
+S 101 of 164
+b:Matrix SF:=
+ [[1 ,2 ,3 ,1 ],_
+ [1 ,3 ,5 ,4 ],_
+ [1 ,3 ,4 ,3 ],_
+ [1 ,3 ,4 ,4 ]]
+R
+R
+R +1. 2.  3. 1.+
+R  
+R 1. 3.  5. 4.
+R (10)  
+R 1. 3.  4. 3.
+R  
+R +1. 3.  4. 4.+
+R Type: Matrix(DoubleFloat)
+E 101
 [1] Wilkinson J H and Reinsch C (1971) Handbook for Automatic
 Computation II, Linear Algebra. SpringerVerlag.
+S 102 of 164
+ result:=f02bjf(n,ia,ib,eps1,matv,iv,a,b,1)
+E 102
 5. Parameters
+)clear all
 1: A(IA,N)  DOUBLE PRECISION array Input
 On entry: the lower triangle of the n by n symmetric matrix
 A. The elements of the array above the diagonal need not be
 set. See also Section 8.
+S 103 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 103
 2: IA  INTEGER Input
 On entry:
 the first dimension of the array A as declared in the
 (sub)program from which F02ABF is called.
 Constraint: IA >= N.
+S 104 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 104
 3: N  INTEGER Input
 On entry: n, the order of the matrix A.
+S 105 of 164
+n : Integer := 16;
+R
+R
+R Type: Integer
+E 105
 4: R(N)  DOUBLE PRECISION array Output
 On exit: the eigenvalues in ascending order.
+S 106 of 164
+k : Integer := 6
+R
+R
+R (4) 6
+R Type: Integer
+E 106
 5: V(IV,N)  DOUBLE PRECISION array Output
 On exit: the normalised eigenvectors, stored by columns;
 the ith column corresponds to the ith eigenvalue. The
 eigenvectors are normalised so that the sum of squares of
 the elements is equal to 1.
+S 107 of 164
+tol : DoubleFloat := 0.0001
+R
+R
+R (5) 9.9999999999999991E5
+R Type: DoubleFloat
+E 107
 6: IV  INTEGER Input
 On entry:
 the first dimension of the array V as declared in the
 (sub)program from which F02ABF is called.
 Constraint: IV >= N.
+S 108 of 164
+novecs : Integer := 0
+R
+R
+R (6) 0
+R Type: Integer
+E 108
 7: E(N)  DOUBLE PRECISION array Workspace
+S 109 of 164
+nrx : Integer := n
+R
+R
+R (7) 16
+R Type: Integer
+E 109
 8: IFAIL  INTEGER Input/Output
 On entry: IFAIL must be set to 0, 1 or 1. For users not
 familiar with this parameter (described in the Essential
 Introduction) the recommended value is 0.
+S 110 of 164
+lwork : Integer := 86
+R
+R
+R (8) 86
+R Type: Integer
+E 110
 On exit: IFAIL = 0 unless the routine detects an error (see
 Section 6).
+S 111 of 164
+lrwork : Integer := 1;
+R
+R
+R Type: Integer
+E 111
 6. Error Indicators and Warnings
+S 112 of 164
+liwork : Integer := 1;
+R
+R
+R Type: Integer
+E 112
 Errors detected by the routine:
+S 113 of 164
+noits : Integer := 1000
+R
+R
+R (11) 1000
+R Type: Integer
+E 113
 IFAIL= 1
 Failure in F02AMF(*) indicating that more than 30*N
 iterations are required to isolate all the eigenvalues.
+S 114 of 164
+m : Integer := 4;
+R
+R
+R Type: Integer
+E 114
 7. Accuracy
+S 115 of 164
+x :Matrix SF:=new(nrx,k,0.0);
+R
+R Compiling function G1781 with type Integer > Boolean
+R
+R Type: Matrix(DoubleFloat)
+E 115
 The eigenvectors are always accurately orthogonal but the
 accuracy of the individual eigenvectors is dependent on their
 inherent sensitivity to changes in the original matrix. For a
 detailed error analysis see Wilkinson and Reinsch [1] pp 222 and
 235.
+S 116 of 164
+ifail : Integer := 1
+R
+R
+R (14)  1
+R Type: Integer
+E 116
 8. Further Comments
+S 117 of 164
+a :Matrix FRAC INT:= new(n,n,0);
+R
+R
+R Type: Matrix(Fraction(Integer))
+E 117
 3
 The time taken by the routine is approximately proportional to n
+S 118 of 164
+a(1,1) := 1;
+R
+R
+R Type: Fraction(Integer)
+E 118
 Unless otherwise stated in the Users' Note for your
 implementation, the routine may be called with the same actual
 array supplied for parameters A and V, in which case the
 eigenvectors will overwrite the original matrix. However this is
 not standard Fortran 77, and may not work on all systems.
+S 119 of 164
+a(1,2) := 1/4;
+R
+R
+R Type: Fraction(Integer)
+E 119
 9. Example
+S 120 of 164
+a(1,5) := 1/4;
+R
+R
+R Type: Fraction(Integer)
+E 120
 To calculate all the eigenvalues and eigenvectors of the real
 symmetric matrix:
+S 121 of 164
+for i in 2..4 repeat
+ a(i,i1) := 1/4
+ a(i,i) := 1
+ a(i,i+1) := 1/4
+ a(i,i+4) := 1/4
+R
+R Type: Void
+E 121
 ( 0.5 0.0 2.3 2.6)
 ( 0.0 0.5 1.4 0.7)
 ( 2.3 1.4 0.5 0.0).
 (2.6 0.7 0.0 0.5)
+S 122 of 164
+for i in 5..n4 repeat
+ a(i,i4) := 1/4
+ a(i,i1) := 1/4
+ a(i,i) := 1
+ a(i,i+1) := 1/4
+ a(i,i+4) := 1/4
+R
+R Type: Void
+E 122
 The example program is not reproduced here. The source code for
 all example programs is distributed with the NAG Foundation
 Library software and should be available online.
+S 123 of 164
+for i in n3..n1 repeat
+ a(i,i4) := 1/4
+ a(i,i1) := 1/4
+ a(i,i) := 1
+ a(i,i+1) := 1/4
+R
+R Type: Void
+E 123
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+S 124 of 164
+a(16,16) := 1;
+R
+R
+R Type: Fraction(Integer)
+E 124
 F02  Eigenvalue and Eigenvectors F02ADF
 F02ADF  NAG Foundation Library Routine Document
+S 125 of 164
+a(16,15) := 1/4;
+R
+R
+R Type: Fraction(Integer)
+E 125
 Note: Before using this routine, please read the Users' Note for
 your implementation to check implementationdependent details.
 The symbol (*) after a NAG routine name denotes a routine that is
 not included in the Foundation Library.
+S 126 of 164
+a(16,12) := 1/4;
+R
+R
+R Type: Fraction(Integer)
+E 126
 1. Purpose
+S 127 of 164
+b:Matrix FRAC INT:= new(n,n,0);
+R
+R
+R Type: Matrix(Fraction(Integer))
+E 127
 F02ADF calculates all the eigenvalues of Ax=(lambda)Bx, where A
 is a real symmetric matrix and B is a real symmetric positive
 definite matrix.
+S 128 of 164
+b(1,1) := 1
+R
+R
+R (26) 1
+R Type: Fraction(Integer)
+E 128
 2. Specification
+S 129 of 164
+b(1,2) := 1/2
+R
+R
+R 1
+R (27)  
+R 2
+R Type: Fraction(Integer)
+E 129
 SUBROUTINE F02ADF (A, IA, B, IB, N, R, DE, IFAIL)
 INTEGER IA, IB, N, IFAIL
 DOUBLE PRECISION A(IA,N), B(IB,N), R(N), DE(N)
+S 130 of 164
+for i in 2..n1 repeat
+ b(i,i1) := 1/2
+ b(i,i) := 1
+ b(i,i+1) := 1/2
+R
+R Type: Void
+E 130
 3. Description
+S 131 of 164
+b(16,15) := 1/2
+R
+R
+R 1
+R (29)  
+R 2
+R Type: Fraction(Integer)
+E 131
 The problem is reduced to the standard symmetric eigenproblem
 using Cholesky's method to decompose B into triangular matrices,
 T
 B=LL , where L is lower triangular. Then Ax=(lambda)Bx implies
 1 T T T
 (L AL )(L x)=(lambda)(L x); hence the eigenvalues of
 Ax=(lambda)Bx are those of Py=(lambda)y where P is the symmetric
 1 T
 matrix L AL . Householder's method is used to tridiagonalise
 the matrix P and the eigenvalues are then found using the QL
 algorithm.

 4. References

 [1] Wilkinson J H and Reinsch C (1971) Handbook for Automatic
 Computation II, Linear Algebra. SpringerVerlag.

 5. Parameters

 1: A(IA,N)  DOUBLE PRECISION array Input/Output
 On entry: the upper triangle of the n by n symmetric matrix
 A. The elements of the array below the diagonal need not be
 set. On exit: the lower triangle of the array is
 overwritten. The rest of the array is unchanged.
+S 132 of 164
+b(16,16) := 1
+R
+R
+R (30) 1
+R Type: Fraction(Integer)
+E 132
 2: IA  INTEGER Input
 On entry:
 the first dimension of the array A as declared in the
 (sub)program from which F02ADF is called.
 Constraint: IA >= N.
+S 133 of 164
+c : Matrix MachineFloat := (inverse (a))*b;
+R
+R
+R Type: Matrix(MachineFloat)
+E 133
 3: B(IB,N)  DOUBLE PRECISION array Input/Output
 On entry: the upper triangle of the n by n symmetric
 positivedefinite matrix B. The elements of the array below
 the diagonal need not be set. On exit: the elements below
 the diagonal are overwritten. The rest of the array is
 unchanged.
+S 134 of 164
+bb := b :: Matrix MachineFloat
+R
+R
+R (32)
+R [[1.0, 0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
+R [ 0.5,1.0, 0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
+R [0.0, 0.5,1.0, 0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
+R [0.0,0.0, 0.5,1.0, 0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
+R [0.0,0.0,0.0, 0.5,1.0, 0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
+R [0.0,0.0,0.0,0.0, 0.5,1.0, 0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
+R [0.0,0.0,0.0,0.0,0.0, 0.5,1.0, 0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
+R [0.0,0.0,0.0,0.0,0.0,0.0, 0.5,1.0, 0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
+R [0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.5,1.0, 0.5,0.0,0.0,0.0,0.0,0.0,0.0],
+R [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.5,1.0, 0.5,0.0,0.0,0.0,0.0,0.0],
+R [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.5,1.0, 0.5,0.0,0.0,0.0,0.0],
+R [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.5,1.0, 0.5,0.0,0.0,0.0],
+R [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.5,1.0, 0.5,0.0,0.0],
+R [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.5,1.0, 0.5,0.0],
+R [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.5,1.0, 0.5],
+R [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.5,1.0]]
+R Type: Matrix(MachineFloat)
+E 134
 4: IB  INTEGER Input
 On entry:
 the first dimension of the array B as declared in the
 (sub)program from which F02ADF is called.
 Constraint: IB >= N.
+S 135 of 164
+ result:=f02fjf(n,k,tol,novecs,nrx,lwork,lrwork,liwork,m,noits,_
+ x,ifail,bb :: ASP27('DOT),c :: ASP28('IMAGE))
+E 135
 5: N  INTEGER Input
 On entry: n, the order of the matrices A and B.
+)clear all
 6: R(N)  DOUBLE PRECISION array Output
 On exit: the eigenvalues in ascending order.
+S 136 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 136
 7: DE(N)  DOUBLE PRECISION array Workspace
+S 137 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 137
 8: IFAIL  INTEGER Input/Output
 On entry: IFAIL must be set to 0, 1 or 1. For users not
 familiar with this parameter (described in the Essential
 Introduction) the recommended value is 0.
+S 138 of 164
+m := 5
+R
+R
+R (3) 5
+R Type: PositiveInteger
+E 138
 On exit: IFAIL = 0 unless the routine detects an error (see
 Section 6).
+S 139 of 164
+n := 3
+R
+R
+R (4) 3
+R Type: PositiveInteger
+E 139
 6. Error Indicators and Warnings
+S 140 of 164
+a : Matrix SF:=
+ [[ 2.0, 2.5, 2.5],_
+ [ 2.0, 2.5, 2.5],_
+ [ 1.6,0.4, 2.8],_
+ [ 2.0,0.5, 0.5],_
+ [ 1.2,0.3,2.9] ]
+R
+R
+R + 2. 2.5 2.5 +
+R  
+R  2. 2.5 2.5 
+R  
+R (5) 1.6000000000000001  0.39999999999999997 2.7999999999999998 
+R  
+R  2.  0.5 0.5 
+R  
+R + 1.2  0.30000000000000004  2.9000000000000004+
+R Type: Matrix(DoubleFloat)
+E 140
 Errors detected by the routine:
+S 141 of 164
+lda := m
+R
+R
+R (6) 5
+R Type: PositiveInteger
+E 141
 IFAIL= 1
 Failure in F01AEF(*); the matrix B is not positivedefinite
 possibly due to rounding errors.
+S 142 of 164
+ncolb := 1
+R
+R
+R (7) 1
+R Type: PositiveInteger
+E 142
 IFAIL= 2
 Failure in F02AVF(*), more than 30*N iterations are required
 to isolate all the eigenvalues.
+S 143 of 164
+b : Matrix SF:= [[ 1.1, 0.9, 0.6, 0.0, 0.8 ]]
+R
+R
+R (8)
+R [
+R [1.1000000000000001, 0.89999999999999991, 0.59999999999999998, 0.,
+R  0.79999999999999993]
+R ]
+R Type: Matrix(DoubleFloat)
+E 143
 7. Accuracy
+S 144 of 164
+ldb := 5
+R
+R
+R (9) 5
+R Type: PositiveInteger
+E 144
 In general this routine is very accurate. However, if B is ill
 conditioned with respect to inversion, the eigenvalues could be
 inaccurately determined. For a detailed error analysis see
 Wilkinson and Reinsch [1] pp 310, 222 and 235.
+S 145 of 164
+wantq := true
+R
+R
+R (10) true
+R Type: Boolean
+E 145
 8. Further Comments
+S 146 of 164
+wantp := true
+R
+R
+R (11) true
+R Type: Boolean
+E 146
 3
 The time taken by the routine is approximately proportional to n
+S 147 of 164
+ldq := 1
+R
+R
+R (12) 1
+R Type: PositiveInteger
+E 147
 9. Example
+S 148 of 164
+ldpt := n
+R
+R
+R (13) 3
+R Type: PositiveInteger
+E 148
 To calculate all the eigenvalues of the general symmetric
 eigenproblem Ax=(lambda) Bx where A is the symmetric matrix:
+S 149 of 164
+ifail := 1
+R
+R
+R (14)  1
+R Type: Integer
+E 149
 (0.5 1.5 6.6 4.8)
 (1.5 6.5 16.2 8.6)
 (6.6 16.2 37.6 9.8)
 (4.8 8.6 9.8 17.1)
+S 150 of 164
+ result:=f02wef(m,n,lda,ncolb,ldb,wantq,ldq,wantp,ldpt,a,b,ifail)
+E 150
 and B is the symmetric positivedefinite matrix:
+)clear all
+
+S 151 of 164
+showArrayValues true
+R
+R
+R (1) true
+R Type: Boolean
+E 151
 (1 3 4 1)
 (3 13 16 11)
 (4 16 24 18).
 (1 11 18 27)
+S 152 of 164
+showScalarValues true
+R
+R
+R (2) true
+R Type: Boolean
+E 152
 The example program is not reproduced here. The source code for
 all example programs is distributed with the NAG Foundation
 Library software and should be available online.
+S 153 of 164
+m:=5
+R
+R
+R (3) 5
+R Type: PositiveInteger
+E 153
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+S 154 of 164
+n:=3
+R
+R
+R (4) 3
+R Type: PositiveInteger
+E 154
 F02  Eigenvalue and Eigenvectors F02AEF
 F02AEF  NAG Foundation Library Routine Document
+S 155 of 164
+lda:=5
+R
+R
+R (5) 5
+R Type: PositiveInteger
+E 155
 Note: Before using this routine, please read the Users' Note for
 your implementation to check implementationdependent details.
 The symbol (*) after a NAG routine name denotes a routine that is
 not included in the Foundation Library.
+S 156 of 164
+ncolb:=1
+R
+R
+R (6) 1
+R Type: PositiveInteger
+E 156
 1. Purpose
+S 157 of 164
+ldb:=5
+R
+R
+R (7) 5
+R Type: PositiveInteger
+E 157
 F02AEF calculates all the eigenvalues and eigenvectors of
 Ax=(lambda)Bx, where A is a real symmetric matrix and B is a
 real symmetric positivedefinite matrix.
+S 158 of 164
+wantq:=true
+R
+R
+R (8) true
+R Type: Boolean
+E 158
 2. Specification
+S 159 of 164
+ldq:=5
+R
+R
+R (9) 5
+R Type: PositiveInteger
+E 159
 SUBROUTINE F02AEF (A, IA, B, IB, N, R, V, IV, DL, E, IFAIL)
 INTEGER IA, IB, N, IV, IFAIL
 DOUBLE PRECISION A(IA,N), B(IB,N), R(N), V(IV,N), DL(N), E
 1 (N)
+S 160 of 164
+wantp:=true
+R
+R
+R (10) true
+R Type: Boolean
+E 160
 3. Description
+S 161 of 164
+ldph:=3
+R
+R
+R (11) 3
+R Type: PositiveInteger
+E 161
 The problem is reduced to the standard symmetric eigenproblem
 using Cholesky's method to decompose B into triangular matrices
 T
 B=LL , where L is lower triangular. Then Ax=(lambda)Bx implies
 1 T T T
 (L AL )(L x)=(lambda)(L x); hence the eigenvalues of
 Ax=(lambda)Bx are those of Py=(lambda)y, where P is the symmetric
 1 T
 matrix L AL . Householder's method is used to tridiagonalise
 the matrix P and the eigenvalues are found using the QL
 algorithm. An eigenvector z of the derived problem is related to
 T
 an eigenvector x of the original problem by z=L x. The
 eigenvectors z are determined using the QL algorithm and are
 T
 normalised so that z z=1; the eigenvectors of the original
 T
 problem are then determined by solving L x=z, and are normalised
 T
 so that x Bx=1.

 4. References

 [1] Wilkinson J H and Reinsch C (1971) Handbook for Automatic
 Computation II, Linear Algebra. SpringerVerlag.

 5. Parameters

 1: A(IA,N)  DOUBLE PRECISION array Input/Output

 On entry: the upper triangle of the n by n symmetric matrix
 A. The elements of the array below the diagonal need not be
 set. On exit: the lower triangle of the array is
 overwritten. The rest of the array is unchanged. See also
 Section 8.

 2: IA  INTEGER Input
 On entry:
 the first dimension of the array A as declared in the
 (sub)program from which F02AEF is called.
 Constraint: IA >= N.

 3: B(IB,N)  DOUBLE PRECISION array Input/Output
 On entry: the upper triangle of the n by n symmetric
 positivedefinite matrix B. The elements of the array below
 the diagonal need not be set. On exit: the elements below
 the diagonal are overwritten. The rest of the array is
 unchanged.

 4: IB  INTEGER Input
 On entry:
 the first dimension of the array B as declared in the
 (sub)program from which F02AEF is called.
 Constraint: IB >= N.

 5: N  INTEGER Input
 On entry: n, the order of the matrices A and B.

 6: R(N)  DOUBLE PRECISION array Output
 On exit: the eigenvalues in ascending order.

 7: V(IV,N)  DOUBLE PRECISION array Output
 On exit: the normalised eigenvectors, stored by columns;
 the ith column corresponds to the ith eigenvalue. The
 T
 eigenvectors x are normalised so that x Bx=1. See also
 Section 8.
+S 162 of 164
+a:Matrix Complex SF:=
+ [[0.5*%i ,0.5 + 1.5*%i ,1 + 1*%i ],_
+ [0.4 + 0.3*%i ,0.9 + 1.3*%i ,0.2 + 1.4*%i ],_
+ [0.4 ,0.4 + 0.4*%i ,1.8 ],_
+ [0.3  0.4*%i ,0.1 + 0.7*%i ,0.0 ],_
+ [0.3*%i ,0.3 + 0.3*%i ,2.4*%i ]]
+R
+R
+R (12)
+R [[0.5 %i, 0.5 + 1.5 %i, 1. + %i],
+R
+R [0.40000000000000002 + 0.29999999999999999 %i,
+R 0.89999999999999991 + 1.2999999999999998 %i,
+R 0.20000000000000001 + 1.3999999999999999 %i]
+R ,
+R
+R [0.40000000000000002,  0.39999999999999997 + 0.40000000000000002 %i,
+R 1.7999999999999998]
+R ,
+R
+R [0.29999999999999999  0.40000000000000002 %i,
+R 0.10000000000000001 + 0.69999999999999996 %i, 0.]
+R ,
+R
+R [ 0.29999999999999999 %i, 0.29999999999999999 + 0.29999999999999999 %i,
+R 2.3999999999999999 %i]
+R ]
+R Type: Matrix(Complex(DoubleFloat))
+E 162
 8: IV  INTEGER Input
 On entry:
 the first dimension of the array V as declared in the
 (sub)program from which F02AEF is called.
 Constraint: IV >= N.
+S 163 of 164
+b:Matrix Complex SF:=
+ [[0.55+1.05*%i ],_
+ [0.49+0.93*%i ],_
+ [0.560.16*%i ],_
+ [0.39+0.23*%i ],_
+ [1.13+0.83*%i ]]
+R
+R
+R + 0.54999999999999993 + 1.0499999999999998 %i+
+R  
+R 0.48999999999999999 + 0.92999999999999994 %i 
+R  
+R (13) 0.56000000000000005  0.15999999999999998 %i 
+R  
+R 0.39000000000000001 + 0.22999999999999998 %i 
+R  
+R + 1.1299999999999999 + 0.82999999999999996 %i +
+R Type: Matrix(Complex(DoubleFloat))
+E 163
 9: DL(N)  DOUBLE PRECISION array Workspace
+S 164 of 164
+ result:=f02xef(m,n,lda,ncolb,ldb,wantq,ldq,wantp,ldph,a,b,1)
+E 164
 10: E(N)  DOUBLE PRECISION array Workspace
+)spool
+)lisp (bye)
+\end{chunk}
+\begin{chunk}{NagEigenPackage.help}
 11: IFAIL  INTEGER Input/Output
 On entry: IFAIL must be set to 0, 1 or 1. For users not
 familiar with this parameter (described in the Essential
 Introduction) the recommended value is 0.
 On exit: IFAIL = 0 unless the routine detects an error (see
 Section 6).
+This package uses the NAG Library to compute
+ * eigenvalues and eigenvectors of a matrix\
+ * eigenvalues and eigenvectors of generalized matrix
+ * eigenvalue problems
+ * singular values and singular vectors of a matrix.
 6. Error Indicators and Warnings
+ F02  Eigenvalues and Eigenvectors Introduction  F02
+ Chapter F02
+ Eigenvalues and Eigenvectors
 Errors detected by the routine:
+ 1. Scope of the Chapter
 IFAIL= 1
 Failure in F01AEF(*); the matrix B is not positivedefinite,
 possibly due to rounding errors.
+ This chapter is concerned with computing
 IFAIL= 2
 Failure in F02AMF(*); more than 30*N iterations are required
 to isolate all the eigenvalues.
+  eigenvalues and eigenvectors of a matrix
 7. Accuracy
+  eigenvalues and eigenvectors of generalized matrix
+ eigenvalue problems
 In general this routine is very accurate. However, if B is ill
 conditioned with respect to inversion, the eigenvectors could be
 inaccurately determined. For a detailed error analysis see
 Wilkinson and Reinsch [1] pp 310, 222 and 235.
+  singular values and singular vectors of a matrix.
 8. Further Comments
+ 2. Background to the Problems
 3
 The time taken by the routine is approximately proportional to n
+ 2.1. Eigenvalue Problems
 Unless otherwise stated in the Users' Note for your
 implementation, the routine may be called with the same actual
 array supplied for parameters A and V, in which case the
 eigenvectors will overwrite the original matrix A. However this
 is not standard Fortran 77, and may not work on all systems.
+ In the most usual form of eigenvalue problem we are given a
+ square n by n matrix A and wish to compute (lambda) (an
+ eigenvalue) and x/=0 (an eigenvector) which satisfy the equation
 9. Example
+ Ax=(lambda)x
 To calculate all the eigenvalues and eigenvectors of the general
 symmetric eigenproblem Ax=(lambda) Bx where A is the symmetric
 matrix:
+ Such problems are called 'standard' eigenvalue problems in
+ contrast to 'generalized' eigenvalue problems where we wish to
+ satisfy the equation
 (0.5 1.5 6.6 4.8)
 (1.5 6.5 16.2 8.6)
 (6.6 16.2 37.6 9.8)
 (4.8 8.6 9.8 17.1)
+ Ax=(lambda)Bx
 and B is the symmetric positivedefinite matrix:
+ B also being a square n by n matrix.
 (1 3 4 1)
 (3 13 16 11)
 (4 16 24 18).
 (1 11 18 27)
+ Section 2.1.1 and Section 2.1.2 discuss, respectively, standard
+ and generalized eigenvalue problems where the matrices involved
+ are dense; Section 2.1.3 discusses both types of problem in the
+ case where A and B are sparse (and symmetric).
 The example program is not reproduced here. The source code for
 all example programs is distributed with the NAG Foundation
 Library software and should be available online.
+ 2.1.1. Standard eigenvalue problems
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ Some of the routines in this chapter find all the n eigenvalues,
+ some find all the n eigensolutions (eigenvalues and corresponding
+ eigenvectors), and some find a selected group of eigenvalues
+ and/or eigenvectors. The matrix A may be:
 F02  Eigenvalue and Eigenvectors F02AFF
 F02AFF  NAG Foundation Library Routine Document
+ (i) general (real or complex)
 Note: Before using this routine, please read the Users' Note for
 your implementation to check implementationdependent details.
 The symbol (*) after a NAG routine name denotes a routine that is
 not included in the Foundation Library.
+ (ii) real symmetric, or
 1. Purpose
+ (iii) complex Hermitian (so that if a =(alpha)+i(beta) then
+ ij
+ a =(alpha)i(beta)).
+ ji
 F02AFF calculates all the eigenvalues of a real unsymmetric
 matrix.
+ In all cases the computation starts with a similarity
+ 1
+ transformation S AS=T, where S is nonsingular and is the
+ product of fairly simple matrices, and T has an 'easier form'
+ than A so that its eigensolutions are easily determined. The
+ matrices A and T, of course, have the same eigenvalues, and if y
+ is an eigenvector of T then Sy is the corresponding eigenvector
+ of A.
 2. Specification
+ In case (i) (general real or complex A), the selected form of T
+ is an upper Hessenberg matrix (t =0 if ij>1) and S is the
+ ij
+ product of n2 stabilised elementary transformation matrices.
+ There is no easy method of computing selected eigenvalues of a
+ Hessenberg matrix, so that all eigenvalues are always calculated.
+ In the real case this computation is performed via the Francis QR
+ algorithm with double shifts, and in the complex case by means of
+ the LR algorithm. If the eigenvectors are required they are
+ computed by backsubstitution following the QR and LR algorithm.
 SUBROUTINE F02AFF (A, IA, N, RR, RI, INTGER, IFAIL)
 INTEGER IA, N, INTGER(N), IFAIL
 DOUBLE PRECISION A(IA,N), RR(N), RI(N)
+ In case (ii) (real and symmetric A) the selected simple form of T
+ is a tridiagonal matrix (t =0 if ij>1), and S is the product
+ ij
+ of n2 orthogonal Householder transformation matrices. If only
+ selected eigenvalues are required, they are obtained by the
+ method of bisection using the Sturm sequence property, and the
+ corresponding eigenvectors of T are computed by inverse
+ iteration. If all eigenvalues are required, they are computed
+ from T via the QL algorithm (an adaptation of the QR algorithm),
+ and the corresponding eigenvectors of T are the product of the
+ transformations for the QL reduction. In all cases the
+ corresponding eigenvectors of A are recovered from the
+ computation of x=Sy.
 3. Description
+ In case (iii) (complex Hermitian A) analogous transformations as
+ in case (ii) are used. T has complex elements in offdiagonal
+ positions, but a simple diagonal similarity transformation is
+ then used to produce a real tridiagonal form, after which the QL
+ algorithm and succeeding methods described in the previous
+ paragraph are used to complete the solution.
 The matrix A is first balanced and then reduced to upper
 Hessenberg form using stabilised elementary similarity
 transformations. The eigenvalues are then found using the QR
 algorithm for real Hessenberg matrices.
+ 2.1.2. Generalized eigenvalue problems
 4. References
+ Here we distinguish as a special case those problems in which
+ both A and B are symmetric and B is positivedefinite and well
+ conditioned with respect to inversion (i.e., all the eigenvalues
+ of B are significantly greater than zero). Such problems can be
+ satisfactorily treated by first reducing them to case (ii) of
+ Section 2.1.1 and then using the methods described there to
+ T
+ compute the eigensolutions. If B is factorized as LL (L lower
+ triangular), then Ax=(lambda)Bx is equivalent to the standard
+ 1 T 1 T
+ symmetric problem Ry=(lambda)y, where R=L A(L ) and y=L x.
+ After finding an eigenvector y of R, the required x is computed
+ T
+ by backsubstitution in y=L x.
 [1] Wilkinson J H and Reinsch C (1971) Handbook for Automatic
 Computation II, Linear Algebra. SpringerVerlag.
+ For generalized problems of the form Ax=(lambda)Bx which do not
+ fall into the special case, the QZ algorithm is provided.
 5. Parameters
+ In order to appreciate the domain in which this algorithm is
+ appropriate we remark first that when B is nonsingular the
+ problem Ax=(lambda)Bx is fully equivalent to the problem
+ 1
+ (B A)x=(lambda)x; both the eigenvalues and eigenvectors being
+ the same. When A is nonsingular Ax=(lambda)Bx is equivalent to
+ 1
+ the problem (A B)x=(mu)x; the eigenvalues (mu) being the
+ reciprocals of the required eigenvalues and the eigenvectors
+ remaining the same. In theory then, provided at least one of the
+ matrices A and B is nonsingular, the generalized problem
+ Ax=(lambda)Bx could be solved via the standard problem
+ Cx=(lambda)x with an appropriate matrix C, and as far as economy
+ of effort is concerned this is quite satisfactory. However, in
+ practice, for this reduction to be satisfactory from the
+ standpoint of numerical stability, one requires more than the
+ 1
+ mere nonsingularity of A or B. It is necessary that B A (or
+ 1
+ A B) should not only exist but that B (or A) should be well
+ conditioned with respect to inversion. The nearer B (or A) is to
+ 1 1
+ singularity the more unsatisfactory B A (or A B) will be as a
+ vehicle for determining the required eigenvalues. Unfortunately
+ 1
+ one cannot counter illconditioning in B (or A) by computing B A
+ 1
+ (or A B) accurately to single precision using iterative
+ refinement. Welldetermined eigenvalues of the original
+ Ax=(lambda)Bx may be poorly determined even by the correctly
+ 1 1
+ rounded version of B A (or A B). The situation may in some
+ instances be saved by the observation that if Ax=(lambda)Bx then
+ (AkB)x=((lambda)k)Bx. Hence if AkB is nonsingular we may
+ 1
+ solve the standard problem [(AkB) B]x=(mu)x and for numerical
+ stability we require only that (AkB) be wellconditioned with
+ respect to inversion.
 1: A(IA,N)  DOUBLE PRECISION array Input/Output
 On entry: the n by n matrix A. On exit: the array is
 overwritten.
+ In practice one may well be in a situation where no k is known
+ for which (AkB) is wellconditioned with respect to inversion
+ and indeed (AkB) may be singular for all k. The QZ algorithm is
+ designed to deal directly with the problem Ax=(lambda)Bx itself
+ and its performance is unaffected by singularity or near
+ singularity of A, B or AkB.
 2: IA  INTEGER Input
 On entry:
 the dimension of the array A as declared in the (sub)program
 from which F02AFF is called.
 Constraint: IA >= N.
+ 2.1.3. Sparse symmetric problems
 3: N  INTEGER Input
 On entry: n, the order of the matrix A.
+ If the matrices A and B are large and sparse (i.e., only a small
+ proportion of the elements are nonzero), then the methods
+ described in the previous Section are unsuitable, because in
+ reducing the problem to a simpler form, much of the sparsity of
+ the problem would be lost; hence the computing time and the
+ storage required would be very large. Instead, for symmetric
+ problems, the method of simultaneous iteration may be used to
+ determine selected eigenvalues and the corresponding
+ eigenvectors. The routine provided has been designed to handle
+ both symmetric and generalized symmetric problems.
 4: RR(N)  DOUBLE PRECISION array Output
 On exit: the real parts of the eigenvalues.
+ 2.2. Singular Value Problems
 5: RI(N)  DOUBLE PRECISION array Output
 On exit: the imaginary parts of the eigenvalues.
+ The singular value decomposition of an m by n real matrix A is
+ given by
 6: INTGER(N)  INTEGER array Output
 On exit: INTGER(i) contains the number of iterations used
 to find the ith eigenvalue. If INTGER(i) is negative, the i
 th eigenvalue is the second of a pair found simultaneously.
+ T
+ A=QDP ,
 Note that the eigenvalues are found in reverse order,
 starting with the nth.
+ where Q is an m by m orthogonal matrix, P is an n by n orthogonal
+ matrix and D is an m by n diagonal matrix with nonnegative
+ diagonal elements. The first k==min(m,n) columns of Q and P are
+ the left and righthand singular vectors of A and the k diagonal
+ elements of D are the singular values.
 7: IFAIL  INTEGER Input/Output
 On entry: IFAIL must be set to 0, 1 or 1. For users not
 familiar with this parameter (described in the Essential
 Introduction) the recommended value is 0.
+ When A is complex then the singular value decomposition is given
+ by
 On exit: IFAIL = 0 unless the routine detects an error (see
 Section 6).
+ H
+ A=QDP ,
 6. Error Indicators and Warnings
+ H T
+ where Q and P are unitary, P denotes the complex conjugate of P
+ and D is as above for the real case.
 Errors detected by the routine:
+ If the matrix A has column means of zero, then AP is the matrix
+ of principal components of A and the singular values are the
+ square roots of the sample variances of the observations with
+ respect to the principal components. (See also Chapter G03.)
 IFAIL= 1
 More than 30*N iterations are required to isolate all the
 eigenvalues.
+ Routines are provided to return the singular values and vectors
+ of a general real or complex matrix.
 7. Accuracy
+ 3. Recommendations on Choice and Use of Routines
 The accuracy of the results depends on the original matrix and
 the multiplicity of the roots. For a detailed error analysis see
 Wilkinson and Reinsch [1] pp 352 and 367.
+ 3.1. General Discussion
 8. Further Comments
+ There is one routine, F02FJF, which is designed for sparse
+ symmetric eigenvalue problems, either standard or generalized.
+ The remainder of the routines are designed for dense matrices.
 3
 The time taken by the routine is approximately proportional to n
+ 3.2. Eigenvalue and Eigenvector Routines
 9. Example
+ These reduce the matrix A to a simpler form by a similarity
+ 1
+ transformation S AS=T where T is an upper Hessenberg or
+ tridiagonal matrix, compute the eigensolutions of T, and then
+ recover the eigenvectors of A via the matrix S. The eigenvectors
+ are normalised so that
 To calculate all the eigenvalues of the real matrix:
+ n
+  2
+ > x  =1
+  r
+ r=1
 ( 1.5 0.1 4.5 1.5)
 (22.5 3.5 12.5 2.5)
 ( 2.5 0.3 4.5 2.5).
 ( 2.5 0.1 4.5 2.5)
+ x being the rth component of the eigenvector x, and so that the
+ r
+ element of largest modulus is real if x is complex. For problems
+ of the type Ax=(lambda)Bx with A and B symmetric and B positive
+ T
+ definite, the eigenvectors are normalised so that x Bx=1, x
+ always being real for such problems.
 The example program is not reproduced here. The source code for
 all example programs is distributed with the NAG Foundation
 Library software and should be available online.
+ 3.3. Singular Value and Singular Vector Routines
+
+ These reduce the matrix A to real bidiagonal form, B say, by
+ T
+ orthogonal transformations Q AP=B in the real case, and by
+ H
+ unitary transformations Q AP=B in the complex case, and the
+ singular values and vectors are computed via this bidiagonal
+ form. The singular values are returned in descending order.
+
+ 3.4. Decision Trees
+
+ (i) Eigenvalues and Eigenvectors
+
+
+ Please see figure in printed Reference Manual
+
+
+ (ii) Singular Values and Singular Vectors
+
+
+ Please see figure in printed Reference Manual
+
+ F02  Eigenvalues and Eigenvectors Contents  F02
+ Chapter F02
+
+ Eigenvalues and Eigenvectors
+
+ F02AAF All eigenvalues of real symmetric matrix
+
+ F02ABF All eigenvalues and eigenvectors of real symmetric matrix
+
+ F02ADF All eigenvalues of generalized real symmetricdefinite
+ eigenproblem
+
+ F02AEF All eigenvalues and eigenvectors of generalized real
+ symmetricdefinite eigenproblem
+
+ F02AFF All eigenvalues of real matrix
+
+ F02AGF All eigenvalues and eigenvectors of real matrix
+
+ F02AJF All eigenvalues of complex matrix
+
+ F02AKF All eigenvalues and eigenvectors of complex matrix
+
+ F02AWF All eigenvalues of complex Hermitian matrix
+
+ F02AXF All eigenvalues and eigenvectors of complex Hermitian
+ matrix
+
+ F02BBF Selected eigenvalues and eigenvectors of real symmetric
+ matrix
+
+ F02BJF All eigenvalues and optionally eigenvectors of
+ generalized eigenproblem by QZ algorithm, real matrices
+
+ F02FJF Selected eigenvalues and eigenvectors of sparse symmetric
+ eigenproblem
+
+ F02WEF SVD of real matrix
+
+ F02XEF SVD of complex matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 F02  Eigenvalue and Eigenvectors F02AGF
 F02AGF  NAG Foundation Library Routine Document
+ F02  Eigenvalue and Eigenvectors F02AAF
+ F02AAF  NAG Foundation Library Routine Document
Note: Before using this routine, please read the Users' Note for
your implementation to check implementationdependent details.
@@ 78450,25 +79398,20 @@ This package uses the NAG Library to compute
1. Purpose
 F02AGF calculates all the eigenvalues and eigenvectors of a real
 unsymmetric matrix.
+ F02AAF calculates all the eigenvalues of a real symmetric matrix.
2. Specification
 SUBROUTINE F02AGF (A, IA, N, RR, RI, VR, IVR, VI, IVI,
 1 INTGER, IFAIL)
 INTEGER IA, N, IVR, IVI, INTGER(N), IFAIL
 DOUBLE PRECISION A(IA,N), RR(N), RI(N), VR(IVR,N), VI
 1 (IVI,N)
+ SUBROUTINE F02AAF (A, IA, N, R, E, IFAIL)
+ INTEGER IA, N, IFAIL
+ DOUBLE PRECISION A(IA,N), R(N), E(N)
3. Description
 The matrix A is first balanced and then reduced to upper
 Hessenberg form using real stabilised elementary similarity
 transformations. The eigenvalues and eigenvectors of the
 Hessenberg matrix are calculated using the QR algorithm. The
 eigenvectors of the Hessenberg matrix are backtransformed to
 give the eigenvectors of the original matrix A.
+ This routine reduces the real symmetric matrix A to a real
+ symmetric tridiagonal matrix using Householder's method. The
+ eigenvalues of the tridiagonal matrix are then determined using
+ the QL algorithm.
4. References
@@ 78478,57 +79421,26 @@ This package uses the NAG Library to compute
5. Parameters
1: A(IA,N)  DOUBLE PRECISION array Input/Output
 On entry: the n by n matrix A. On exit: the array is
 overwritten.
+ On entry: the lower triangle of the n by n symmetric matrix
+ A. The elements of the array above the diagonal need not be
+ set. On exit: the elements of A below the diagonal are
+ overwritten, and the rest of the array is unchanged.
2: IA  INTEGER Input
On entry:
the first dimension of the array A as declared in the
 (sub)program from which F02AGF is called.
+ (sub)program from which F02AAF is called.
Constraint: IA >= N.
3: N  INTEGER Input
On entry: n, the order of the matrix A.
 4: RR(N)  DOUBLE PRECISION array Output
 On exit: the real parts of the eigenvalues.
+ 4: R(N)  DOUBLE PRECISION array Output
+ On exit: the eigenvalues in ascending order.
 5: RI(N)  DOUBLE PRECISION array Output
 On exit: the imaginary parts of the eigenvalues.
+ 5: E(N)  DOUBLE PRECISION array Workspace
 6: VR(IVR,N)  DOUBLE PRECISION array Output
 On exit: the real parts of the eigenvectors, stored by
 columns. The ith column corresponds to the ith eigenvalue.
 The eigenvectors are normalised so that the sum of the
 squares of the moduli of the elements is equal to 1 and the
 element of largest modulus is real. This ensures that real
 eigenvalues have real eigenvectors.

 7: IVR  INTEGER Input
 On entry:
 the first dimension of the array VR as declared in the
 (sub)program from which F02AGF is called.
 Constraint: IVR >= N.

 8: VI(IVI,N)  DOUBLE PRECISION array Output
 On exit: the imaginary parts of the eigenvectors, stored by
 columns. The ith column corresponds to the ith eigenvalue.

 9: IVI  INTEGER Input
 On entry:
 the first dimension of the array VI as declared in the
 (sub)program from which F02AGF is called.
 Constraint: IVI >= N.

 10: INTGER(N)  INTEGER array Output
 On exit: INTGER(i) contains the number of iterations used
 to find the ith eigenvalue. If INTGER(i) is negative, the i
 th eigenvalue is the second of a pair found simultaneously.

 Note that the eigenvalues are found in reverse order,
 starting with the nth.

 11: IFAIL  INTEGER Input/Output
+ 6: IFAIL  INTEGER Input/Output
On entry: IFAIL must be set to 0, 1 or 1. For users not
familiar with this parameter (described in the Essential
Introduction) the recommended value is 0.
@@ 78541,14 +79453,15 @@ This package uses the NAG Library to compute
Errors detected by the routine:
IFAIL= 1
 More than 30*N iterations are required to isolate all the
 eigenvalues.
+ Failure in F02AVF(*) indicating that more than 30*N
+ iterations are required to isolate all the eigenvalues.
7. Accuracy
 The accuracy of the results depends on the original matrix and
 the multiplicity of the roots. For a detailed error analysis see
 Wilkinson and Reinsch [1] pp 352 and 390.
+ The accuracy of the eigenvalues depends on the sensitivity of the
+ matrix to rounding errors produced in tridiagonalisation. For a
+ detailed error analysis see Wilkinson and Reinsch [1] pp 222 and
+ 235.
8. Further Comments
@@ 78557,13 +79470,12 @@ This package uses the NAG Library to compute
9. Example
 To calculate all the eigenvalues and eigenvectors of the real
 matrix:
+ To calculate all the eigenvalues of the real symmetric matrix:
 ( 1.5 0.1 4.5 1.5)
 (22.5 3.5 12.5 2.5)
 ( 2.5 0.3 4.5 2.5).
 ( 2.5 0.1 4.5 2.5)
+ ( 0.5 0.0 2.3 2.6)
+ ( 0.0 0.5 1.4 0.7)
+ ( 2.3 1.4 0.5 0.0).
+ (2.6 0.7 0.0 0.5)
The example program is not reproduced here. The source code for
all example programs is distributed with the NAG Foundation
@@ 78571,8 +79483,8 @@ This package uses the NAG Library to compute
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 F02  Eigenvalue and Eigenvectors F02AJF
 F02AJF  NAG Foundation Library Routine Document
+ F02  Eigenvalue and Eigenvectors F02ABF
+ F02ABF  NAG Foundation Library Routine Document
Note: Before using this routine, please read the Users' Note for
your implementation to check implementationdependent details.
@@ 78581,21 +79493,21 @@ This package uses the NAG Library to compute
1. Purpose
 F02AJF calculates all the eigenvalues of a complex matrix.
+ F02ABF calculates all the eigenvalues and eigenvectors of a real
+ symmetric matrix.
2. Specification
 SUBROUTINE F02AJF (AR, IAR, AI, IAI, N, RR, RI, INTGER,
 1 IFAIL)
 INTEGER IAR, IAI, N, INTGER(N), IFAIL
 DOUBLE PRECISION AR(IAR,N), AI(IAI,N), RR(N), RI(N)
+ SUBROUTINE F02ABF (A, IA, N, R, V, IV, E, IFAIL)
+ INTEGER IA, N, IV, IFAIL
+ DOUBLE PRECISION A(IA,N), R(N), V(IV,N), E(N)
3. Description
 The complex matrix A is first balanced and then reduced to upper
 Hessenberg form using stabilised elementary similarity
 transformations. The eigenvalues are then found using the
 modified LR algorithm for complex Hessenberg matrices.
+ This routine reduces the real symmetric matrix A to a real
+ symmetric tridiagonal matrix by Householder's method. The
+ eigenvalues and eigenvectors are calculated using the QL
+ algorithm.
4. References
@@ 78604,38 +79516,38 @@ This package uses the NAG Library to compute
5. Parameters
 1: AR(IAR,N)  DOUBLE PRECISION array Input/Output
 On entry: the real parts of the elements of the n by n
 complex matrix A. On exit: the array is overwritten.

 2: IAR  INTEGER Input
 On entry:
 the first dimension of the array AR as declared in the
 (sub)program from which F02AJF is called.
 Constraint: IAR >= N.

 3: AI(IAI,N)  DOUBLE PRECISION array Input/Output
 On entry: the imaginary parts of the elements of the n by n
 complex matrix A. On exit: the array is overwritten.
+ 1: A(IA,N)  DOUBLE PRECISION array Input
+ On entry: the lower triangle of the n by n symmetric matrix
+ A. The elements of the array above the diagonal need not be
+ set. See also Section 8.
 4: IAI  INTEGER Input
+ 2: IA  INTEGER Input
On entry:
 the first dimension of the array AI as declared in the
 (sub)program from which F02AJF is called.
 Constraint: IAI >= N.
+ the first dimension of the array A as declared in the
+ (sub)program from which F02ABF is called.
+ Constraint: IA >= N.
 5: N  INTEGER Input
+ 3: N  INTEGER Input
On entry: n, the order of the matrix A.
 6: RR(N)  DOUBLE PRECISION array Output
 On exit: the real parts of the eigenvalues.
+ 4: R(N)  DOUBLE PRECISION array Output
+ On exit: the eigenvalues in ascending order.
 7: RI(N)  DOUBLE PRECISION array Output
 On exit: the imaginary parts of the eigenvalues.
+ 5: V(IV,N)  DOUBLE PRECISION array Output
+ On exit: the normalised eigenvectors, stored by columns;
+ the ith column corresponds to the ith eigenvalue. The
+ eigenvectors are normalised so that the sum of squares of
+ the elements is equal to 1.
 8: INTGER(N)  INTEGER array Workspace
+ 6: IV  INTEGER Input
+ On entry:
+ the first dimension of the array V as declared in the
+ (sub)program from which F02ABF is called.
+ Constraint: IV >= N.
 9: IFAIL  INTEGER Input/Output
+ 7: E(N)  DOUBLE PRECISION array Workspace
+
+ 8: IFAIL  INTEGER Input/Output
On entry: IFAIL must be set to 0, 1 or 1. For users not
familiar with this parameter (described in the Essential
Introduction) the recommended value is 0.
@@ 78648,28 +79560,37 @@ This package uses the NAG Library to compute
Errors detected by the routine:
IFAIL= 1
 More than 30*N iterations are required to isolate all the
 eigenvalues.
+ Failure in F02AMF(*) indicating that more than 30*N
+ iterations are required to isolate all the eigenvalues.
7. Accuracy
 The accuracy of the results depends on the original matrix and
 the multiplicity of the roots. For a detailed error analysis see
 Wilkinson and Reinsch [1] pp 352 and 401.
+ The eigenvectors are always accurately orthogonal but the
+ accuracy of the individual eigenvectors is dependent on their
+ inherent sensitivity to changes in the original matrix. For a
+ detailed error analysis see Wilkinson and Reinsch [1] pp 222 and
+ 235.
8. Further Comments
3
The time taken by the routine is approximately proportional to n
+ Unless otherwise stated in the Users' Note for your
+ implementation, the routine may be called with the same actual
+ array supplied for parameters A and V, in which case the
+ eigenvectors will overwrite the original matrix. However this is
+ not standard Fortran 77, and may not work on all systems.
+
9. Example
 To calculate all the eigenvalues of the complex matrix:
+ To calculate all the eigenvalues and eigenvectors of the real
+ symmetric matrix:
 (21.05.0i 24.60i 13.6+10.2i 4.0i)
 ( 22.5i 26.005.00i 7.510.0i 2.5 )
 ( 2.0+1.5i 1.68+2.24i 4.55.0i 1.5+2.0i).
 ( 2.5i 2.60 2.7+3.6i 2.55.0i)
+ ( 0.5 0.0 2.3 2.6)
+ ( 0.0 0.5 1.4 0.7)
+ ( 2.3 1.4 0.5 0.0).
+ (2.6 0.7 0.0 0.5)
The example program is not reproduced here. The source code for
all example programs is distributed with the NAG Foundation
@@ 78677,8 +79598,8 @@ This package uses the NAG Library to compute
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 F02  Eigenvalue and Eigenvectors F02AKF
 F02AKF  NAG Foundation Library Routine Document
+ F02  Eigenvalue and Eigenvectors F02ADF
+ F02ADF  NAG Foundation Library Routine Document
Note: Before using this routine, please read the Users' Note for
your implementation to check implementationdependent details.
@@ 78687,25 +79608,29 @@ This package uses the NAG Library to compute
1. Purpose
 F02AKF calculates all the eigenvalues and eigenvectors of a
 complex matrix.
+ F02ADF calculates all the eigenvalues of Ax=(lambda)Bx, where A
+ is a real symmetric matrix and B is a real symmetric positive
+ definite matrix.
2. Specification
 SUBROUTINE F02AKF (AR, IAR, AI, IAI, N, RR, RI, VR, IVR,
 1 VI, IVI, INTGER, IFAIL)
 INTEGER IAR, IAI, N, IVR, IVI, INTGER(N), IFAIL
 DOUBLE PRECISION AR(IAR,N), AI(IAI,N), RR(N), RI(N), VR
 1 (IVR,N), VI(IVI,N)
+ SUBROUTINE F02ADF (A, IA, B, IB, N, R, DE, IFAIL)
+ INTEGER IA, IB, N, IFAIL
+ DOUBLE PRECISION A(IA,N), B(IB,N), R(N), DE(N)
3. Description
 The complex matrix A is first balanced and then reduced to upper
 Hessenberg form by stabilised elementary similarity
 transformations. The eigenvalues and eigenvectors of the
 Hessenberg matrix are calculated using the LR algorithm. The
 eigenvectors of the Hessenberg matrix are backtransformed to
 give the eigenvectors of the original matrix.
+ The problem is reduced to the standard symmetric eigenproblem
+ using Cholesky's method to decompose B into triangular matrices,
+ T
+ B=LL , where L is lower triangular. Then Ax=(lambda)Bx implies
+ 1 T T T
+ (L AL )(L x)=(lambda)(L x); hence the eigenvalues of
+ Ax=(lambda)Bx are those of Py=(lambda)y where P is the symmetric
+ 1 T
+ matrix L AL . Householder's method is used to tridiagonalise
+ the matrix P and the eigenvalues are then found using the QL
+ algorithm.
4. References
@@ 78714,61 +79639,40 @@ This package uses the NAG Library to compute
5. Parameters
 1: AR(IAR,N)  DOUBLE PRECISION array Input/Output
 On entry: the real parts of the elements of the n by n
 complex matrix A. On exit: the array is overwritten.
+ 1: A(IA,N)  DOUBLE PRECISION array Input/Output
+ On entry: the upper triangle of the n by n symmetric matrix
+ A. The elements of the array below the diagonal need not be
+ set. On exit: the lower triangle of the array is
+ overwritten. The rest of the array is unchanged.
 2: IAR  INTEGER Input
+ 2: IA  INTEGER Input
On entry:
 the first dimension of the array AR as declared in the
 (sub)program from which F02AKF is called.
 Constraint: IAR >= N.
+ the first dimension of the array A as declared in the
+ (sub)program from which F02ADF is called.
+ Constraint: IA >= N.
 3: AI(IAI,N)  DOUBLE PRECISION array Input/Output
 On entry: the imaginary parts of the elements of the n by n
 complex matrix A. On exit: the array is overwritten.
+ 3: B(IB,N)  DOUBLE PRECISION array Input/Output
+ On entry: the upper triangle of the n by n symmetric
+ positivedefinite matrix B. The elements of the array below
+ the diagonal need not be set. On exit: the elements below
+ the diagonal are overwritten. The rest of the array is
+ unchanged.
 4: IAI  INTEGER Input
+ 4: IB  INTEGER Input
On entry:
 the first dimension of the array AI as declared in the
 (sub)program from which F02AKF is called.
 Constraint: IAI >= N.
+ the first dimension of the array B as declared in the
+ (sub)program from which F02ADF is called.
+ Constraint: IB >= N.
5: N  INTEGER Input
 On entry: n, the order of the matrix A.

 6: RR(N)  DOUBLE PRECISION array Output
 On exit: the real parts of the eigenvalues.

 7: RI(N)  DOUBLE PRECISION array Output
 On exit: the imaginary parts of the eigenvalues.

 8: VR(IVR,N)  DOUBLE PRECISION array Output
 On exit: the real parts of the eigenvectors, stored by
 columns. The ith column corresponds to the ith eigenvalue.
 The eigenvectors are normalised so that the sum of squares
 of the moduli of the elements is equal to 1 and the element
 of largest modulus is real.

 9: IVR  INTEGER Input
 On entry:
 the first dimension of the array VR as declared in the
 (sub)program from which F02AKF is called.
 Constraint: IVR >= N.

 10: VI(IVI,N)  DOUBLE PRECISION array Output
 On exit: the imaginary parts of the eigenvectors, stored by
 columns. The ith column corresponds to the ith eigenvalue.
+ On entry: n, the order of the matrices A and B.
 11: IVI  INTEGER Input
 On entry:
 the first dimension of the array VI as declared in the
 (sub)program from which F02AKF is called.
 Constraint: IVI >= N.
+ 6: R(N)  DOUBLE PRECISION array Output
+ On exit: the eigenvalues in ascending order.
 12: INTGER(N)  INTEGER array Workspace
+ 7: DE(N)  DOUBLE PRECISION array Workspace
 13: IFAIL  INTEGER Input/Output
+ 8: IFAIL  INTEGER Input/Output
On entry: IFAIL must be set to 0, 1 or 1. For users not
familiar with this parameter (described in the Essential
Introduction) the recommended value is 0.
@@ 78781,14 +79685,19 @@ This package uses the NAG Library to compute
Errors detected by the routine:
IFAIL= 1
 More than 30*N iterations are required to isolate all the
 eigenvalues.
+ Failure in F01AEF(*); the matrix B is not positivedefinite
+ possibly due to rounding errors.
+
+ IFAIL= 2
+ Failure in F02AVF(*), more than 30*N iterations are required
+ to isolate all the eigenvalues.
7. Accuracy
 The accuracy of the results depends on the conditioning of the
 original matrix and the multiplicity of the roots. For a detailed
 error analysis see Wilkinson and Reinsch [1] pp 352 and 390.
+ In general this routine is very accurate. However, if B is ill
+ conditioned with respect to inversion, the eigenvalues could be
+ inaccurately determined. For a detailed error analysis see
+ Wilkinson and Reinsch [1] pp 310, 222 and 235.
8. Further Comments
@@ 78797,13 +79706,20 @@ This package uses the NAG Library to compute
9. Example
 To calculate all the eigenvalues and eigenvectors of the complex
 matrix:
+ To calculate all the eigenvalues of the general symmetric
+ eigenproblem Ax=(lambda) Bx where A is the symmetric matrix:
 (21.05.0i 24.60i 13.6+10.2i 4.0i)
 ( 22.5i 26.005.00i 7.510.0i 2.5 )
 ( 2.0+1.5i 1.68+2.24i 4.55.0i 1.5+2.0i).
 ( 2.5i 2.60 2.7+3.6i 2.55.0i)
+ (0.5 1.5 6.6 4.8)
+ (1.5 6.5 16.2 8.6)
+ (6.6 16.2 37.6 9.8)
+ (4.8 8.6 9.8 17.1)
+
+ and B is the symmetric positivedefinite matrix:
+
+ (1 3 4 1)
+ (3 13 16 11)
+ (4 16 24 18).
+ (1 11 18 27)
The example program is not reproduced here. The source code for
all example programs is distributed with the NAG Foundation
@@ 78811,8 +79727,8 @@ This package uses the NAG Library to compute
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 F02  Eigenvalue and Eigenvectors F02AWF
 F02AWF  NAG Foundation Library Routine Document
+ F02  Eigenvalue and Eigenvectors F02AEF
+ F02AEF  NAG Foundation Library Routine Document
Note: Before using this routine, please read the Users' Note for
your implementation to check implementationdependent details.
@@ 78821,76 +79737,101 @@ This package uses the NAG Library to compute
1. Purpose
 F02AWF calculates all the eigenvalues of a complex Hermitian
 matrix.
+ F02AEF calculates all the eigenvalues and eigenvectors of
+ Ax=(lambda)Bx, where A is a real symmetric matrix and B is a
+ real symmetric positivedefinite matrix.
2. Specification
 SUBROUTINE F02AWF (AR, IAR, AI, IAI, N, R, WK1, WK2, WK3,
 1 IFAIL)
 INTEGER IAR, IAI, N, IFAIL
 DOUBLE PRECISION AR(IAR,N), AI(IAI,N), R(N), WK1(N),
 1 WK2(N), WK3(N)
+ SUBROUTINE F02AEF (A, IA, B, IB, N, R, V, IV, DL, E, IFAIL)
+ INTEGER IA, IB, N, IV, IFAIL
+ DOUBLE PRECISION A(IA,N), B(IB,N), R(N), V(IV,N), DL(N), E
+ 1 (N)
3. Description
 The complex Hermitian matrix A is first reduced to a real
 tridiagonal matrix by n2 unitary transformations, and a
 subsequent diagonal transformation. The eigenvalues are then
 derived using the QL algorithm, an adaptation of the QR
 algorithm.
+ The problem is reduced to the standard symmetric eigenproblem
+ using Cholesky's method to decompose B into triangular matrices
+ T
+ B=LL , where L is lower triangular. Then Ax=(lambda)Bx implies
+ 1 T T T
+ (L AL )(L x)=(lambda)(L x); hence the eigenvalues of
+ Ax=(lambda)Bx are those of Py=(lambda)y, where P is the symmetric
+ 1 T
+ matrix L AL . Householder's method is used to tridiagonalise
+ the matrix P and the eigenvalues are found using the QL
+ algorithm. An eigenvector z of the derived problem is related to
+ T
+ an eigenvector x of the original problem by z=L x. The
+ eigenvectors z are determined using the QL algorithm and are
+ T
+ normalised so that z z=1; the eigenvectors of the original
+ T
+ problem are then determined by solving L x=z, and are normalised
+ T
+ so that x Bx=1.
4. References
 [1] Peters G (1967) NPL Algorithms Library. Document No.
 F1/04/A.

 [2] Wilkinson J H and Reinsch C (1971) Handbook for Automatic
+ [1] Wilkinson J H and Reinsch C (1971) Handbook for Automatic
Computation II, Linear Algebra. SpringerVerlag.
5. Parameters
 1: AR(IAR,N)  DOUBLE PRECISION array Input/Output
 On entry: the real parts of the elements of the lower
 triangle of the n by n complex Hermitian matrix A. Elements
 of the array above the diagonal need not be set. On exit:
 the array is overwritten.
+ 1: A(IA,N)  DOUBLE PRECISION array Input/Output
 2: IAR  INTEGER Input
+ On entry: the upper triangle of the n by n symmetric matrix
+ A. The elements of the array below the diagonal need not be
+ set. On exit: the lower triangle of the array is
+ overwritten. The rest of the array is unchanged. See also
+ Section 8.
+
+ 2: IA  INTEGER Input
On entry:
 the first dimension of the array AR as declared in the
 (sub)program from which F02AWF is called.
 Constraint: IAR >= N.
+ the first dimension of the array A as declared in the
+ (sub)program from which F02AEF is called.
+ Constraint: IA >= N.
 3: AI(IAI,N)  DOUBLE PRECISION array Input/Output
 On entry: the imaginary parts of the elements of the lower
 triangle of the n by n complex Hermitian matrix A. Elements
 of the array above the diagonal need not be set. On exit:
 the array is overwritten.
+ 3: B(IB,N)  DOUBLE PRECISION array Input/Output
+ On entry: the upper triangle of the n by n symmetric
+ positivedefinite matrix B. The elements of the array below
+ the diagonal need not be set. On exit: the elements below
+ the diagonal are overwritten. The rest of the array is
+ unchanged.
 4: IAI  INTEGER Input
+ 4: IB  INTEGER Input
On entry:
 the first dimension of the array AI as declared in the
 (sub)program from which F02AWF is called.
 Constraint: IAI >= N.
+ the first dimension of the array B as declared in the
+ (sub)program from which F02AEF is called.
+ Constraint: IB >= N.
5: N  INTEGER Input
 On entry: n, the order of the complex Hermitian matrix, A.
+ On entry: n, the order of the matrices A and B.
6: R(N)  DOUBLE PRECISION array Output
On exit: the eigenvalues in ascending order.
 7: WK1(N)  DOUBLE PRECISION array Workspace
+ 7: V(IV,N)  DOUBLE PRECISION array Output
+ On exit: the normalised eigenvectors, stored by columns;
+ the ith column corresponds to the ith eigenvalue. The
+ T
+ eigenvectors x are normalised so that x Bx=1. See also
+ Section 8.
 8: WK2(N)  DOUBLE PRECISION array Workspace
+ 8: IV  INTEGER Input
+ On entry:
+ the first dimension of the array V as declared in the
+ (sub)program from which F02AEF is called.
+ Constraint: IV >= N.
 9: WK3(N)  DOUBLE PRECISION array Workspace
+ 9: DL(N)  DOUBLE PRECISION array Workspace
 10: IFAIL  INTEGER Input/Output
+ 10: E(N)  DOUBLE PRECISION array Workspace
+
+ 11: IFAIL  INTEGER Input/Output
On entry: IFAIL must be set to 0, 1 or 1. For users not
familiar with this parameter (described in the Essential
Introduction) the recommended value is 0.

On exit: IFAIL = 0 unless the routine detects an error (see
Section 6).
@@ 78899,27 +79840,48 @@ This package uses the NAG Library to compute
Errors detected by the routine:
IFAIL= 1
 More than 30*N iterations are required to isolate all the
 eigenvalues.
+ Failure in F01AEF(*); the matrix B is not positivedefinite,
+ possibly due to rounding errors.
+
+ IFAIL= 2
+ Failure in F02AMF(*); more than 30*N iterations are required
+ to isolate all the eigenvalues.
7. Accuracy
 For a detailed error analysis see Peters [1] page 3 and Wilkinson
 and Reinsch [2] page 235.
+ In general this routine is very accurate. However, if B is ill
+ conditioned with respect to inversion, the eigenvectors could be
+ inaccurately determined. For a detailed error analysis see
+ Wilkinson and Reinsch [1] pp 310, 222 and 235.
8. Further Comments
3
The time taken by the routine is approximately proportional to n
+ Unless otherwise stated in the Users' Note for your
+ implementation, the routine may be called with the same actual
+ array supplied for parameters A and V, in which case the
+ eigenvectors will overwrite the original matrix A. However this
+ is not standard Fortran 77, and may not work on all systems.
+
9. Example
 To calculate all the eigenvalues of the complex Hermitian matrix:
+ To calculate all the eigenvalues and eigenvectors of the general
+ symmetric eigenproblem Ax=(lambda) Bx where A is the symmetric
+ matrix:
 (0.50 0.00 1.84+1.38i 2.081.56i)
 (0.00 0.50 1.12+0.84i 0.56+0.42i)
 (1.841.38i 1.120.84i 0.50 0.00 ).
 (2.08+1.56i 0.560.42i 0.00 0.50 )
+ (0.5 1.5 6.6 4.8)
+ (1.5 6.5 16.2 8.6)
+ (6.6 16.2 37.6 9.8)
+ (4.8 8.6 9.8 17.1)
+
+ and B is the symmetric positivedefinite matrix:
+
+ (1 3 4 1)
+ (3 13 16 11)
+ (4 16 24 18).
+ (1 11 18 27)
The example program is not reproduced here. The source code for
all example programs is distributed with the NAG Foundation
@@ 78927,8 +79889,8 @@ This package uses the NAG Library to compute
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 F02  Eigenvalue and Eigenvectors F02AXF
 F02AXF  NAG Foundation Library Routine Document
+ F02  Eigenvalue and Eigenvectors F02AFF
+ F02AFF  NAG Foundation Library Routine Document
Note: Before using this routine, please read the Users' Note for
your implementation to check implementationdependent details.
@@ 78937,96 +79899,57 @@ This package uses the NAG Library to compute
1. Purpose
 F02AXF calculates all the eigenvalues and eigenvectors of a
 complex Hermitian matrix.
+ F02AFF calculates all the eigenvalues of a real unsymmetric
+ matrix.
2. Specification
 SUBROUTINE F02AXF (AR, IAR, AI, IAI, N, R, VR, IVR, VI,
 1 IVI, WK1, WK2, WK3, IFAIL)
 INTEGER IAR, IAI, N, IVR, IVI, IFAIL
 DOUBLE PRECISION AR(IAR,N), AI(IAI,N), R(N), VR(IVR,N), VI
 1 (IVI,N), WK1(N), WK2(N), WK3(N)
+ SUBROUTINE F02AFF (A, IA, N, RR, RI, INTGER, IFAIL)
+ INTEGER IA, N, INTGER(N), IFAIL
+ DOUBLE PRECISION A(IA,N), RR(N), RI(N)
3. Description
 The complex Hermitian matrix is first reduced to a real
 tridiagonal matrix by n2 unitary transformations and a
 subsequent diagonal transformation. The eigenvalues and
 eigenvectors are then derived using the QL algorithm, an
 adaptation of the QR algorithm.
+ The matrix A is first balanced and then reduced to upper
+ Hessenberg form using stabilised elementary similarity
+ transformations. The eigenvalues are then found using the QR
+ algorithm for real Hessenberg matrices.
4. References
 [1] Peters G (1967) NPL Algorithms Library. Document No.
 F2/03/A.

 [2] Peters G (1967) NPL Algorithms Library. Document No.
 F1/04/A.
+ [1] Wilkinson J H and Reinsch C (1971) Handbook for Automatic
+ Computation II, Linear Algebra. SpringerVerlag.
5. Parameters
 1: AR(IAR,N)  DOUBLE PRECISION array Input
 On entry: the real parts of the elements of the lower
 triangle of the n by n complex Hermitian matrix A. Elements
 of the array above the diagonal need not be set. See also
 Section 8.
+ 1: A(IA,N)  DOUBLE PRECISION array Input/Output
+ On entry: the n by n matrix A. On exit: the array is
+ overwritten.
 2: IAR  INTEGER Input
+ 2: IA  INTEGER Input
On entry:
 the first dimension of the array AR as declared in the
 (sub)program from which F02AXF is called.
 Constraint: IAR >= N.
+ the dimension of the array A as declared in the (sub)program
+ from which F02AFF is called.
+ Constraint: IA >= N.
 3: AI(IAI,N)  DOUBLE PRECISION array Input
 On entry: the imaginary parts of the elements of the lower
 triangle of the n by n complex Hermitian matrix A. Elements
 of the array above the diagonal need not be set. See also
 Section 8.
+ 3: N  INTEGER Input
+ On entry: n, the order of the matrix A.
 4: IAI  INTEGER Input
 On entry:
 the first dimension of the array AI as declared in the
 (sub)program from which F02AXF is called.
 Constraint: IAI >= N.
+ 4: RR(N)  DOUBLE PRECISION array Output
+ On exit: the real parts of the eigenvalues.
 5: N  INTEGER Input
 On entry: n, the order of the matrix, A.
+ 5: RI(N)  DOUBLE PRECISION array Output
+ On exit: the imaginary parts of the eigenvalues.
 6: R(N)  DOUBLE PRECISION array Output
 On exit: the eigenvalues in ascending order.
+ 6: INTGER(N)  INTEGER array Output
+ On exit: INTGER(i) contains the number of iterations used
+ to find the ith eigenvalue. If INTGER(i) is negative, the i
+ th eigenvalue is the second of a pair found simultaneously.
 7: VR(IVR,N)  DOUBLE PRECISION array Output
 On exit: the real parts of the eigenvectors, stored by
 columns. The ith column corresponds to the ith eigenvector.
 The eigenvectors are normalised so that the sum of the
 squares of the moduli of the elements is equal to 1 and the
 element of largest modulus is real. See also Section 8.

 8: IVR  INTEGER Input
 On entry:
 the first dimension of the array VR as declared in the
 (sub)program from which F02AXF is called.
 Constraint: IVR >= N.

 9: VI(IVI,N)  DOUBLE PRECISION array Output
 On exit: the imaginary parts of the eigenvectors, stored by
 columns. The ith column corresponds to the ith eigenvector.
 See also Section 8.

 10: IVI  INTEGER Input
 On entry:
 the first dimension of the array VI as declared in the
 (sub)program from which F02AXF is called.
 Constraint: IVI >= N.

 11: WK1(N)  DOUBLE PRECISION array Workspace

 12: WK2(N)  DOUBLE PRECISION array Workspace

 13: WK3(N)  DOUBLE PRECISION array Workspace
+ Note that the eigenvalues are found in reverse order,
+ starting with the nth.
 14: IFAIL  INTEGER Input/Output
+ 7: IFAIL  INTEGER Input/Output
On entry: IFAIL must be set to 0, 1 or 1. For users not
familiar with this parameter (described in the Essential
Introduction) the recommended value is 0.
@@ 79042,38 +79965,25 @@ This package uses the NAG Library to compute
More than 30*N iterations are required to isolate all the
eigenvalues.
 IFAIL= 2
 The diagonal elements of AI are not all zero, i.e., the
 complex matrix is not Hermitian.

7. Accuracy
 The eigenvectors are always accurately orthogonal but the
 accuracy of the individual eigenvalues and eigenvectors is
 dependent on their inherent sensitivity to small changes in the
 original matrix. For a detailed error analysis see Peters [1]
 page 3 and [2] page 3.
+ The accuracy of the results depends on the original matrix and
+ the multiplicity of the roots. For a detailed error analysis see
+ Wilkinson and Reinsch [1] pp 352 and 367.
8. Further Comments
3
The time taken by the routine is approximately proportional to n
 Unless otherwise stated in the implementation document, the
 routine may be called with the same actual array supplied for
 parameters AR and VR, and for AI and VI, in which case the
 eigenvectors will overwrite the original matrix A. However this
 is not standard Fortran 77, and may not work on all systems.

9. Example
 To calculate the eigenvalues and eigenvectors of the complex
 Hermitian matrix:
+ To calculate all the eigenvalues of the real matrix:
 (0.50 0.00 1.84+1.38i 2.081.56i)
 (0.00 0.50 1.12+0.84i 0.56+0.42i)
 (1.841.38i 1.120.84i 0.50 0.00 ).
 (2.08+1.56i 0.560.42i 0.00 0.50 )
+ ( 1.5 0.1 4.5 1.5)
+ (22.5 3.5 12.5 2.5)
+ ( 2.5 0.3 4.5 2.5).
+ ( 2.5 0.1 4.5 2.5)
The example program is not reproduced here. The source code for
all example programs is distributed with the NAG Foundation
@@ 79081,8 +79991,8 @@ This package uses the NAG Library to compute
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 F02  Eigenvalue and Eigenvectors F02BBF
 F02BBF  NAG Foundation Library Routine Document
+ F02  Eigenvalue and Eigenvectors F02AGF
+ F02AGF  NAG Foundation Library Routine Document
Note: Before using this routine, please read the Users' Note for
your implementation to check implementationdependent details.
@@ 79091,28 +80001,25 @@ This package uses the NAG Library to compute
1. Purpose
 F02BBF calculates selected eigenvalues and eigenvectors of a real
 symmetric matrix by reduction to tridiagonal form, bisection and
 inverse iteration, where the selected eigenvalues lie within a
 given interval.
+ F02AGF calculates all the eigenvalues and eigenvectors of a real
+ unsymmetric matrix.
2. Specification
 SUBROUTINE F02BBF (A, IA, N, ALB, UB, M, MM, R, V, IV, D,
 1 E, E2, X, G, C, ICOUNT, IFAIL)
 INTEGER IA, N, M, MM, IV, ICOUNT(M), IFAIL
 DOUBLE PRECISION A(IA,N), ALB, UB, R(M), V(IV,M), D(N), E
 1 (N), E2(N), X(N,7), G(N)
 LOGICAL C(N)
+ SUBROUTINE F02AGF (A, IA, N, RR, RI, VR, IVR, VI, IVI,
+ 1 INTGER, IFAIL)
+ INTEGER IA, N, IVR, IVI, INTGER(N), IFAIL
+ DOUBLE PRECISION A(IA,N), RR(N), RI(N), VR(IVR,N), VI
+ 1 (IVI,N)
3. Description
 The real symmetric matrix A is reduced to a symmetric tridiagonal
 matrix T by Householder's method. The eigenvalues which lie
 within a given interval [l,u], are calculated by the method of
 bisection. The corresponding eigenvectors of T are calculated by
 inverse iteration. A backtransformation is then performed to
 obtain the eigenvectors of the original matrix A.
+ The matrix A is first balanced and then reduced to upper
+ Hessenberg form using real stabilised elementary similarity
+ transformations. The eigenvalues and eigenvectors of the
+ Hessenberg matrix are calculated using the QR algorithm. The
+ eigenvectors of the Hessenberg matrix are backtransformed to
+ give the eigenvectors of the original matrix A.
4. References
@@ 79122,67 +80029,57 @@ This package uses the NAG Library to compute
5. Parameters
1: A(IA,N)  DOUBLE PRECISION array Input/Output
 On entry: the lower triangle of the n by n symmetric matrix
 A. The elements of the array above the diagonal need not be
 set. On exit: the elements of A below the diagonal are
 overwritten, and the rest of the array is unchanged.
+ On entry: the n by n matrix A. On exit: the array is
+ overwritten.
2: IA  INTEGER Input
On entry:
the first dimension of the array A as declared in the
 (sub)program from which F02BBF is called.
+ (sub)program from which F02AGF is called.
Constraint: IA >= N.
3: N  INTEGER Input
On entry: n, the order of the matrix A.
 4: ALB  DOUBLE PRECISION Input

 5: UB  DOUBLE PRECISION Input
 On entry: l and u, the lower and upper endpoints of the
 interval within which eigenvalues are to be calculated.

 6: M  INTEGER Input
 On entry: an upper bound for the number of eigenvalues
 within the interval.

 7: MM  INTEGER Output
 On exit: the actual number of eigenvalues within the
 interval.
+ 4: RR(N)  DOUBLE PRECISION array Output
+ On exit: the real parts of the eigenvalues.
 8: R(M)  DOUBLE PRECISION array Output
 On exit: the eigenvalues, not necessarily in ascending
 order.
+ 5: RI(N)  DOUBLE PRECISION array Output
+ On exit: the imaginary parts of the eigenvalues.
 9: V(IV,M)  DOUBLE PRECISION array Output
 On exit: the eigenvectors, stored by columns. The ith
 column corresponds to the ith eigenvalue. The eigenvectors
 are normalised so that the sum of the squares of the
 elements are equal to 1.
+ 6: VR(IVR,N)  DOUBLE PRECISION array Output
+ On exit: the real parts of the eigenvectors, stored by
+ columns. The ith column corresponds to the ith eigenvalue.
+ The eigenvectors are normalised so that the sum of the
+ squares of the moduli of the elements is equal to 1 and the
+ element of largest modulus is real. This ensures that real
+ eigenvalues have real eigenvectors.
 10: IV  INTEGER Input
+ 7: IVR  INTEGER Input
On entry:
 the first dimension of the array V as declared in the
 (sub)program from which F02BBF is called.
 Constraint: IV >= N.

 11: D(N)  DOUBLE PRECISION array Workspace

 12: E(N)  DOUBLE PRECISION array Workspace

 13: E2(N)  DOUBLE PRECISION array Workspace
+ the first dimension of the array VR as declared in the
+ (sub)program from which F02AGF is called.
+ Constraint: IVR >= N.
 14: X(N,7)  DOUBLE PRECISION array Workspace
+ 8: VI(IVI,N)  DOUBLE PRECISION array Output
+ On exit: the imaginary parts of the eigenvectors, stored by
+ columns. The ith column corresponds to the ith eigenvalue.
 15: G(N)  DOUBLE PRECISION array Workspace
+ 9: IVI  INTEGER Input
+ On entry:
+ the first dimension of the array VI as declared in the
+ (sub)program from which F02AGF is called.
+ Constraint: IVI >= N.
 16: C(N)  LOGICAL array Workspace
+ 10: INTGER(N)  INTEGER array Output
+ On exit: INTGER(i) contains the number of iterations used
+ to find the ith eigenvalue. If INTGER(i) is negative, the i
+ th eigenvalue is the second of a pair found simultaneously.
 17: ICOUNT(M)  INTEGER array Output
 On exit: ICOUNT(i) contains the number of iterations for
 the ith eigenvalue.
+ Note that the eigenvalues are found in reverse order,
+ starting with the nth.
 18: IFAIL  INTEGER Input/Output
+ 11: IFAIL  INTEGER Input/Output
On entry: IFAIL must be set to 0, 1 or 1. For users not
familiar with this parameter (described in the Essential
Introduction) the recommended value is 0.
@@ 79195,40 +80092,29 @@ This package uses the NAG Library to compute
Errors detected by the routine:
IFAIL= 1
 M is less than the number of eigenvalues in the given
 interval. On exit MM contains the number of eigenvalues in
 the interval. Rerun with this value for M.

 IFAIL= 2
 More than 5 iterations are required to determine any one
 eigenvector.
+ More than 30*N iterations are required to isolate all the
+ eigenvalues.
7. Accuracy
 There is no guarantee of the accuracy of the eigenvectors as the
 results depend on the original matrix and the multiplicity of the
 roots. For a detailed error analysis see Wilkinson and Reinsch
 [1] pp 222 and 436.
+ The accuracy of the results depends on the original matrix and
+ the multiplicity of the roots. For a detailed error analysis see
+ Wilkinson and Reinsch [1] pp 352 and 390.
8. Further Comments
3
The time taken by the routine is approximately proportional to n
 This subroutine should only be used when less than 25% of the
 eigenvalues and the corresponding eigenvectors are required. Also
 this subroutine is less efficient with matrices which have
 multiple eigenvalues.

9. Example
 To calculate the eigenvalues lying between 2.0 and 3.0, and the
 corresponding eigenvectors of the real symmetric matrix:
+ To calculate all the eigenvalues and eigenvectors of the real
+ matrix:
 ( 0.5 0.0 2.3 2.6)
 ( 0.0 0.5 1.4 0.7)
 ( 2.3 1.4 0.5 0.0).
 (2.6 0.7 0.0 0.5)
+ ( 1.5 0.1 4.5 1.5)
+ (22.5 3.5 12.5 2.5)
+ ( 2.5 0.3 4.5 2.5).
+ ( 2.5 0.1 4.5 2.5)
The example program is not reproduced here. The source code for
all example programs is distributed with the NAG Foundation
@@ 79236,8 +80122,8 @@ This package uses the NAG Library to compute
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 F02  Eigenvalue and Eigenvectors F02BJF
 F02BJF  NAG Foundation Library Routine Document
+ F02  Eigenvalue and Eigenvectors F02AJF
+ F02AJF  NAG Foundation Library Routine Document
Note: Before using this routine, please read the Users' Note for
your implementation to check implementationdependent details.
@@ 79246,146 +80132,61 @@ This package uses the NAG Library to compute
1. Purpose
 F02BJF calculates all the eigenvalues and, if required, all the
 eigenvectors of the generalized eigenproblem Ax=(lambda)Bx
 where A and B are real, square matrices, using the QZ algorithm.
+ F02AJF calculates all the eigenvalues of a complex matrix.
2. Specification
 SUBROUTINE F02BJF (N, A, IA, B, IB, EPS1, ALFR, ALFI,
 1 BETA, MATV, V, IV, ITER, IFAIL)
 INTEGER N, IA, IB, IV, ITER(N), IFAIL
 DOUBLE PRECISION A(IA,N), B(IB,N), EPS1, ALFR(N), ALFI(N),
 1 BETA(N), V(IV,N)
 LOGICAL MATV
+ SUBROUTINE F02AJF (AR, IAR, AI, IAI, N, RR, RI, INTGER,
+ 1 IFAIL)
+ INTEGER IAR, IAI, N, INTGER(N), IFAIL
+ DOUBLE PRECISION AR(IAR,N), AI(IAI,N), RR(N), RI(N)
3. Description
 All the eigenvalues and, if required, all the eigenvectors of the
 generalized eigenproblem Ax=(lambda)Bx where A and B are real,
 square matrices, are determined using the QZ algorithm. The QZ
 algorithm consists of 4 stages:

 (a) A is reduced to upper Hessenberg form and at the same time
 B is reduced to upper triangular form.
+ The complex matrix A is first balanced and then reduced to upper
+ Hessenberg form using stabilised elementary similarity
+ transformations. The eigenvalues are then found using the
+ modified LR algorithm for complex Hessenberg matrices.
 (b) A is further reduced to quasitriangular form while the
 triangular form of B is maintained.
+ 4. References
 (c) The quasitriangular form of A is reduced to triangular
 form and the eigenvalues extracted. This routine does not
 actually produce the eigenvalues (lambda) , but instead
 j
 returns (alpha) and (beta) such that
+ [1] Wilkinson J H and Reinsch C (1971) Handbook for Automatic
+ Computation II, Linear Algebra. SpringerVerlag.
 j j
 (lambda) =(alpha) /(beta) , j=1,2,...,.n
 j j j
 The division by (beta) becomes the responsibility of the
 j
 user's program, since (beta) may be zero indicating an
 j
 infinite eigenvalue. Pairs of complex eigenvalues occur
 with (alpha) /(beta) and (alpha) /(beta) complex
 j j j+1 j+1
+ 5. Parameters
 conjugates, even though (alpha) and (alpha) are not
 j j+1
 conjugate.
+ 1: AR(IAR,N)  DOUBLE PRECISION array Input/Output
+ On entry: the real parts of the elements of the n by n
+ complex matrix A. On exit: the array is overwritten.
 (d) If the eigenvectors are required (MATV = .TRUE.), they are
 obtained from the triangular matrices and then transformed
 back into the original coordinate system.

 4. References

 [1] Moler C B and Stewart G W (1973) An Algorithm for
 Generalized Matrix Eigenproblems. SIAM J. Numer. Anal. 10
 241256.

 [2] Ward R C (1975) The Combination Shift QZ Algorithm. SIAM J.
 Numer. Anal. 12 835853.

 [3] Wilkinson J H (1979) Kronecker's Canonical Form and the QZ
 Algorithm. Linear Algebra and Appl. 28 285303.

 5. Parameters

 1: N  INTEGER Input
 On entry: n, the order of the matrices A and B.

 2: A(IA,N)  DOUBLE PRECISION array Input/Output
 On entry: the n by n matrix A. On exit: the array is
 overwritten.

 3: IA  INTEGER Input
+ 2: IAR  INTEGER Input
On entry:
 the first dimension of the array A as declared in the
 (sub)program from which F02BJF is called.
 Constraint: IA >= N.
+ the first dimension of the array AR as declared in the
+ (sub)program from which F02AJF is called.
+ Constraint: IAR >= N.
 4: B(IB,N)  DOUBLE PRECISION array Input/Output
 On entry: the n by n matrix B. On exit: the array is
 overwritten.
+ 3: AI(IAI,N)  DOUBLE PRECISION array Input/Output
+ On entry: the imaginary parts of the elements of the n by n
+ complex matrix A. On exit: the array is overwritten.
 5: IB  INTEGER Input
+ 4: IAI  INTEGER Input
On entry:
 the first dimension of the array B as declared in the
 (sub)program from which F02BJF is called.
 Constraint: IB >= N.

 6: EPS1  DOUBLE PRECISION Input
 On entry: the tolerance used to determine negligible
 elements. If EPS1 > 0.0, an element will be considered
 negligible if it is less than EPS1 times the norm of its
 matrix. If EPS1 <= 0.0, machine precision is used in place
 of EPS1. A positive value of EPS1 may result in faster
 execution but less accurate results.

 7: ALFR(N)  DOUBLE PRECISION array Output

 8: ALFI(N)  DOUBLE PRECISION array Output
 On exit: the real and imaginary parts of (alpha) , for
 j
 j=1,2,...,n.

 9: BETA(N)  DOUBLE PRECISION array Output
 On exit: (beta) , for j=1,2,...,n.
 j

 10: MATV  LOGICAL Input
 On entry: MATV must be set .TRUE. if the eigenvectors are
 required, otherwise .FALSE..

 11: V(IV,N)  DOUBLE PRECISION array Output
 On exit: if MATV = .TRUE., then
 (i)if the jth eigenvalue is real, the jth column of V
 contains its eigenvector;
+ the first dimension of the array AI as declared in the
+ (sub)program from which F02AJF is called.
+ Constraint: IAI >= N.
 (ii) if the jth and (j+1)th eigenvalues form a complex
 pair, the jth and (j+1)th columns of V contain the
 real and imaginary parts of the eigenvector associated
 with the first eigenvalue of the pair. The conjugate
 of this vector is the eigenvector for the conjugate
 eigenvalue.
 Each eigenvector is normalised so that the component of
 largest modulus is real and the sum of squares of the moduli
 equal one.
+ 5: N  INTEGER Input
+ On entry: n, the order of the matrix A.
 If MATV = .FALSE., V is not used.
+ 6: RR(N)  DOUBLE PRECISION array Output
+ On exit: the real parts of the eigenvalues.
 12: IV  INTEGER Input
 On entry:
 the first dimension of the array V as declared in the
 (sub)program from which F02BJF is called.
 Constraint: IV >= N.
+ 7: RI(N)  DOUBLE PRECISION array Output
+ On exit: the imaginary parts of the eigenvalues.
 13: ITER(N)  INTEGER array Output
 On exit: ITER(j) contains the number of iterations needed
 to obtain the jth eigenvalue. Note that the eigenvalues are
 obtained in reverse order, starting with the nth.
+ 8: INTGER(N)  INTEGER array Workspace
 14: IFAIL  INTEGER Input/Output
+ 9: IFAIL  INTEGER Input/Output
On entry: IFAIL must be set to 0, 1 or 1. For users not
familiar with this parameter (described in the Essential
Introduction) the recommended value is 0.
@@ 79397,59 +80198,29 @@ This package uses the NAG Library to compute
Errors detected by the routine:
 IFAIL= i
 More than 30*N iterations are required to determine all the
 diagonal 1 by 1 or 2 by 2 blocks of the quasitriangular
 form in the second step of the QZ algorithm. IFAIL is set to
 the index i of the eigenvalue at which this failure occurs.
 If the soft failure option is used, (alpha) and (beta) are
 j j
 correct for j=i+1,i+2,...,n, but V does not contain any
 correct eigenvectors.
+ IFAIL= 1
+ More than 30*N iterations are required to isolate all the
+ eigenvalues.
7. Accuracy
 The computed eigenvalues are always exact for a problem
 (A+E)x=(lambda)(B+F)x where E/A and F/B
 are both of the order of max(EPS1,(epsilon)), EPS1 being defined
 as in Section 5 and (epsilon) being the machine precision.

 Note: interpretation of results obtained with the QZ algorithm
 often requires a clear understanding of the effects of small
 changes in the original data. These effects are reviewed in
 Wilkinson [3], in relation to the significance of small values of
 (alpha) and (beta) . It should be noted that if (alpha) and
 j j j
 (beta) are both small for any j, it may be that no reliance can
 j
 be placed on any of the computed eigenvalues
 (lambda) =(alpha) /(beta) . The user is recommended to study [3]
 i i i
 and, if in difficulty, to seek expert advice on determining the
 sensitivity of the eigenvalues to perturbations in the data.
+ The accuracy of the results depends on the original matrix and
+ the multiplicity of the roots. For a detailed error analysis see
+ Wilkinson and Reinsch [1] pp 352 and 401.
8. Further Comments
3
The time taken by the routine is approximately proportional to n
 and also depends on the value chosen for parameter EPS1.
9. Example
 To find all the eigenvalues and eigenvectors of Ax=(lambda) Bx
 where

 (3.9 12.5 34.5 0.5)
 (4.3 21.5 47.5 7.5)
 A=(4.3 21.5 43.5 3.5)
 (4.4 26.0 46.0 6.0)

 and
+ To calculate all the eigenvalues of the complex matrix:
 (1 2 3 1)
 (1 3 5 4)
 B=(1 3 4 3).
 (1 3 4 4)
+ (21.05.0i 24.60i 13.6+10.2i 4.0i)
+ ( 22.5i 26.005.00i 7.510.0i 2.5 )
+ ( 2.0+1.5i 1.68+2.24i 4.55.0i 1.5+2.0i).
+ ( 2.5i 2.60 2.7+3.6i 2.55.0i)
The example program is not reproduced here. The source code for
all example programs is distributed with the NAG Foundation
@@ 79457,8 +80228,8 @@ This package uses the NAG Library to compute
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 F02  Eigenvalue and Eigenvectors F02FJF
 F02FJF  NAG Foundation Library Routine Document
+ F02  Eigenvalue and Eigenvectors F02AKF
+ F02AKF  NAG Foundation Library Routine Document
Note: Before using this routine, please read the Users' Note for
your implementation to check implementationdependent details.
@@ 79467,653 +80238,393 @@ This package uses the NAG Library to compute
1. Purpose
 To find eigenvalues and eigenvectors of a real sparse symmetric
 or generalized symmetric eigenvalue problem.
+ F02AKF calculates all the eigenvalues and eigenvectors of a
+ complex matrix.
2. Specification
 SUBROUTINE F02FJF (N, M, K, NOITS, TOL, DOT, IMAGE, MONIT,
 1 NOVECS, X, NRX, D, WORK, LWORK, RWORK,
 2 LRWORK, IWORK, LIWORK, IFAIL)
 INTEGER N, M, K, NOITS, NOVECS, NRX, LWORK,
 1 LRWORK, IWORK(LIWORK), LIWORK, IFAIL
 DOUBLE PRECISION TOL, DOT, X(NRX,K), D(K), WORK(LWORK),
 1 RWORK(LRWORK)
 EXTERNAL DOT, IMAGE, MONIT
+ SUBROUTINE F02AKF (AR, IAR, AI, IAI, N, RR, RI, VR, IVR,
+ 1 VI, IVI, INTGER, IFAIL)
+ INTEGER IAR, IAI, N, IVR, IVI, INTGER(N), IFAIL
+ DOUBLE PRECISION AR(IAR,N), AI(IAI,N), RR(N), RI(N), VR
+ 1 (IVR,N), VI(IVI,N)
3. Description
 F02FJF finds the m eigenvalues of largest absolute value and the
 corresponding eigenvectors for the real eigenvalue problem
+ The complex matrix A is first balanced and then reduced to upper
+ Hessenberg form by stabilised elementary similarity
+ transformations. The eigenvalues and eigenvectors of the
+ Hessenberg matrix are calculated using the LR algorithm. The
+ eigenvectors of the Hessenberg matrix are backtransformed to
+ give the eigenvectors of the original matrix.
 Cx=(lambda)x (1)
+ 4. References
 where C is an n by n matrix such that
+ [1] Wilkinson J H and Reinsch C (1971) Handbook for Automatic
+ Computation II, Linear Algebra. SpringerVerlag.
 T
 BC=C B (2)
+ 5. Parameters
 for a given positivedefinite matrix B. C is said to be B
 symmetric. Different specifications of C allow for the solution
 of a variety of eigenvalue problems. For example, when
+ 1: AR(IAR,N)  DOUBLE PRECISION array Input/Output
+ On entry: the real parts of the elements of the n by n
+ complex matrix A. On exit: the array is overwritten.
 T
 C=A and B=I where A=A
+ 2: IAR  INTEGER Input
+ On entry:
+ the first dimension of the array AR as declared in the
+ (sub)program from which F02AKF is called.
+ Constraint: IAR >= N.
 the routine finds the m eigenvalues of largest absolute magnitude
 for the standard symmetric eigenvalue problem
+ 3: AI(IAI,N)  DOUBLE PRECISION array Input/Output
+ On entry: the imaginary parts of the elements of the n by n
+ complex matrix A. On exit: the array is overwritten.
 Ax=(lambda)x. (3)
+ 4: IAI  INTEGER Input
+ On entry:
+ the first dimension of the array AI as declared in the
+ (sub)program from which F02AKF is called.
+ Constraint: IAI >= N.
 The routine is intended for the case where A is sparse.
+ 5: N  INTEGER Input
+ On entry: n, the order of the matrix A.
 As a second example, when
+ 6: RR(N)  DOUBLE PRECISION array Output
+ On exit: the real parts of the eigenvalues.
 1
 C=B A
+ 7: RI(N)  DOUBLE PRECISION array Output
+ On exit: the imaginary parts of the eigenvalues.
 where
+ 8: VR(IVR,N)  DOUBLE PRECISION array Output
+ On exit: the real parts of the eigenvectors, stored by
+ columns. The ith column corresponds to the ith eigenvalue.
+ The eigenvectors are normalised so that the sum of squares
+ of the moduli of the elements is equal to 1 and the element
+ of largest modulus is real.
 T
 A=A
+ 9: IVR  INTEGER Input
+ On entry:
+ the first dimension of the array VR as declared in the
+ (sub)program from which F02AKF is called.
+ Constraint: IVR >= N.
 the routine finds the m eigenvalues of largest absolute magnitude
 for the generalized symmetric eigenvalue problem
+ 10: VI(IVI,N)  DOUBLE PRECISION array Output
+ On exit: the imaginary parts of the eigenvectors, stored by
+ columns. The ith column corresponds to the ith eigenvalue.
 Ax=(lambda)Bx. (4)
+ 11: IVI  INTEGER Input
+ On entry:
+ the first dimension of the array VI as declared in the
+ (sub)program from which F02AKF is called.
+ Constraint: IVI >= N.
 The routine is intended for the case where A and B are sparse.
+ 12: INTGER(N)  INTEGER array Workspace
 The routine does not require C explicitly, but C is specified via
 a usersupplied routine IMAGE which, given an n element vector z,
 computes the image w given by
+ 13: IFAIL  INTEGER Input/Output
+ On entry: IFAIL must be set to 0, 1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
 w=Cz.
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
 1
 For instance, in the above example, where C=B A, routine IMAGE
 will need to solve the positivedefinite system of equations
 Bw=Az for w.
+ 6. Error Indicators and Warnings
 To find the m eigenvalues of smallest absolute magnitude of (3)
 1
 we can choose C=A and hence find the reciprocals of the
 required eigenvalues, so that IMAGE will need to solve Aw=z for
 1
 w, and correspondingly for (4) we can choose C=A B and solve
 Aw=Bz for w.
+ Errors detected by the routine:
 A table of examples of choice of IMAGE is given in Table 3.1. It
 should be remembered that the routine also returns the
 corresponding eigenvectors and that B is positivedefinite.
 Throughout A is assumed to be symmetric and, where necessary,
 nonsingularity is also assumed.
+ IFAIL= 1
+ More than 30*N iterations are required to isolate all the
+ eigenvalues.
 Eigenvalues Problem
 Required
+ 7. Accuracy
 Ax=(lambda)x (B=I)Ax=(lambda)Bx ABx=(lambda)x
+ The accuracy of the results depends on the conditioning of the
+ original matrix and the multiplicity of the roots. For a detailed
+ error analysis see Wilkinson and Reinsch [1] pp 352 and 390.
 Largest Compute Solve Compute
 w=Az Bw=Az w=ABz
+ 8. Further Comments
 Smallest Solve Solve Solve
 (Find Aw=z Aw=Bz Av=z, Bw=(nu)
 1/(lambda))
+ 3
+ The time taken by the routine is approximately proportional to n
 Furthest Compute Solve Compute
 from w=(A(sigma)I)z Bw=(A(sigma)B)z w=(AB(sigma)I)z
 (sigma)
 (Find
 (lambda)
 (sigma))
+ 9. Example
 Closest to Solve Solve Solve
 (sigma) (A(sigma)I)w=z (A(sigma)B)w=Bz (AB(sigma)I)w=z
 (Find 1/((
 lambda)
 (sigma)))
+ To calculate all the eigenvalues and eigenvectors of the complex
+ matrix:
+ (21.05.0i 24.60i 13.6+10.2i 4.0i)
+ ( 22.5i 26.005.00i 7.510.0i 2.5 )
+ ( 2.0+1.5i 1.68+2.24i 4.55.0i 1.5+2.0i).
+ ( 2.5i 2.60 2.7+3.6i 2.55.0i)
 Table 3.1
 The Requirement of IMAGE for Various Problems
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available online.
 The matrix B also need not be supplied explicitly, but is
 specified via a usersupplied routine DOT which, given n element
 T
 vectors z and w, computes the generalized dot product w Bz.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 F02FJF is based upon routine SIMITZ (see Nikolai [1]), which is
 itself a derivative of the Algol procedure ritzit (see
 Rutishauser [4]), and uses the method of simultaneous (subspace)
 iteration. (See Parlett [2] for description, analysis and advice
 on the use of the method.)
+ F02  Eigenvalue and Eigenvectors F02AWF
+ F02AWF  NAG Foundation Library Routine Document
 The routine performs simultaneous iteration on k>m vectors.
 Initial estimates to p<=k eigenvectors, corresponding to the p
 eigenvalues of C of largest absolute value, may be supplied by
 the user to F02FJF. When possible k should be chosen so that the
 kth eigenvalue is not too close to the m required eigenvalues,
 but if k is initially chosen too small then F02FJF may be re
 entered, supplying approximations to the k eigenvectors found so
 far and with k then increased.
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementationdependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
 At each major iteration F02FJF solves an r by r (r<=k) eigenvalue
 subproblem in order to obtain an approximation to the
 eigenvalues for which convergence has not yet occurred. This
 approximation is refined by Chebyshev acceleration.
+ 1. Purpose
 4. References
+ F02AWF calculates all the eigenvalues of a complex Hermitian
+ matrix.
 [1] Nikolai P J (1979) Algorithm 538: Eigenvectors and
 eigenvalues of real generalized symmetric matrices by
 simultaneous iteration. ACM Trans. Math. Softw. 5 118125.
+ 2. Specification
 [2] Parlett B N (1980) The Symmetric Eigenvalue Problem.
 PrenticeHall.
+ SUBROUTINE F02AWF (AR, IAR, AI, IAI, N, R, WK1, WK2, WK3,
+ 1 IFAIL)
+ INTEGER IAR, IAI, N, IFAIL
+ DOUBLE PRECISION AR(IAR,N), AI(IAI,N), R(N), WK1(N),
+ 1 WK2(N), WK3(N)
 [3] Rutishauser H (1969) Computational aspects of F L Bauer's
 simultaneous iteration method. Num. Math. 13 413.
+ 3. Description
 [4] Rutishauser H (1970) Simultaneous iteration method for
 symmetric matrices. Num. Math. 16 205223.
+ The complex Hermitian matrix A is first reduced to a real
+ tridiagonal matrix by n2 unitary transformations, and a
+ subsequent diagonal transformation. The eigenvalues are then
+ derived using the QL algorithm, an adaptation of the QR
+ algorithm.
 5. Parameters
+ 4. References
 1: N  INTEGER Input
 On entry: n, the order of the matrix C. Constraint: N >= 1.
+ [1] Peters G (1967) NPL Algorithms Library. Document No.
+ F1/04/A.
 2: M  INTEGER Input/Output
 On entry: m, the number of eigenvalues required.
 '
 Constraint: M >= 1. On exit: m, the number of eigenvalues
 actually found. It is equal to m if IFAIL = 0 on exit, and
 is less than m if IFAIL = 2, 3 or 4. See Section 6 and
 Section 8 for further information.
+ [2] Wilkinson J H and Reinsch C (1971) Handbook for Automatic
+ Computation II, Linear Algebra. SpringerVerlag.
 3: K  INTEGER Input
 On entry: the number of simultaneous iteration vectors to be
 used. Too small a value of K may inhibit convergence, while
 a larger value of K incurs additional storage and additional
 work per iteration. Suggested value: K = M + 4 will often be
 a reasonable choice in the absence of better information.
 Constraint: M < K <= N.
+ 5. Parameters
 4: NOITS  INTEGER Input/Output
 On entry: the maximum number of major iterations (eigenvalue
 subproblems) to be performed. If NOITS <= 0, then the value
 100 is used in place of NOITS. On exit: the number of
 iterations actually performed.
+ 1: AR(IAR,N)  DOUBLE PRECISION array Input/Output
+ On entry: the real parts of the elements of the lower
+ triangle of the n by n complex Hermitian matrix A. Elements
+ of the array above the diagonal need not be set. On exit:
+ the array is overwritten.
 5: TOL  DOUBLE PRECISION Input
 On entry: a relative tolerance to be used in accepting
 eigenvalues and eigenvectors. If the eigenvalues are
 required to about t significant figures, then TOL should be
 t
 set to about 10 . d is accepted as an eigenvalue as soon
 i
 as two successive approximations to d differ by less than
 i
 ~ ~
 (d *TOL)/10, where d is the latest aproximation to d .
 i i i
 Once an eigenvalue has been accepted, then an eigenvector is
 accepted as soon as (d f )/(d d )= N.
 6: DOT  DOUBLE PRECISION FUNCTION, supplied by the user.
 External Procedure
 T
 DOT must return the value w Bz for given vectors w and z.
 For the standard eigenvalue problem, where B=I, DOT must
 T
 return the dot product w z.
+ 3: AI(IAI,N)  DOUBLE PRECISION array Input/Output
+ On entry: the imaginary parts of the elements of the lower
+ triangle of the n by n complex Hermitian matrix A. Elements
+ of the array above the diagonal need not be set. On exit:
+ the array is overwritten.
 Its specification is:
+ 4: IAI  INTEGER Input
+ On entry:
+ the first dimension of the array AI as declared in the
+ (sub)program from which F02AWF is called.
+ Constraint: IAI >= N.
 DOUBLE PRECISION FUNCTION DOT (IFLAG, N, Z, W,
 1 RWORK, LRWORK,
 2 IWORK, LIWORK)
 INTEGER IFLAG, N, LRWORK, IWORK(LIWORK),
 1 LIWORK
 DOUBLE PRECISION Z(N), W(N), RWORK(LRWORK)
+ 5: N  INTEGER Input
+ On entry: n, the order of the complex Hermitian matrix, A.
 1: IFLAG  INTEGER Input/Output
 On entry: IFLAG is always nonnegative. On exit: IFLAG
 may be used as a flag to indicate a failure in the
 T
 computation of w Bz. If IFLAG is negative on exit from
 DOT, then F02FJF will exit immediately with IFAIL set
 to IFLAG. Note that in this case DOT must still be
 assigned a value.
+ 6: R(N)  DOUBLE PRECISION array Output
+ On exit: the eigenvalues in ascending order.
 2: N  INTEGER Input
 On entry: the number of elements in the vectors z and w
 and the order of the matrix B.
+ 7: WK1(N)  DOUBLE PRECISION array Workspace
 3: Z(N)  DOUBLE PRECISION array Input
 T
 On entry: the vector z for which w Bz is required.
+ 8: WK2(N)  DOUBLE PRECISION array Workspace
 4: W(N)  DOUBLE PRECISION array Input
 T
 On entry: the vector w for which w Bz is required.
+ 9: WK3(N)  DOUBLE PRECISION array Workspace
 5: RWORK(LRWORK)  DOUBLE PRECISION array User Workspace
+ 10: IFAIL  INTEGER Input/Output
+ On entry: IFAIL must be set to 0, 1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
 6: LRWORK  INTEGER Input
+ 6. Error Indicators and Warnings
+ Errors detected by the routine:
 7: IWORK(LIWORK)  INTEGER array User Workspace
+ IFAIL= 1
+ More than 30*N iterations are required to isolate all the
+ eigenvalues.
+ 7. Accuracy
 8: LIWORK  INTEGER Input
 DOT is called from F02FJF with the parameters RWORK,
 LRWORK, IWORK and LIWORK as supplied to F02FJF. The
 user is free to use the arrays RWORK and IWORK to
 supply information to DOT and to IMAGE as an
 alternative to using COMMON.
 DOT must be declared as EXTERNAL in the (sub)program
 from which F02FJF is called. Parameters denoted as
 Input must not be changed by this procedure.
+ For a detailed error analysis see Peters [1] page 3 and Wilkinson
+ and Reinsch [2] page 235.
 7: IMAGE  SUBROUTINE, supplied by the user.
 External Procedure
 IMAGE must return the vector w=Cz for a given vector z.

 Its specification is:
+ 8. Further Comments
 SUBROUTINE IMAGE (IFLAG, N, Z, W, RWORK, LRWORK,
 1 IWORK, LIWORK)
 INTEGER IFLAG, N, LRWORK, IWORK(LIWORK),
 1 LIWORK
 DOUBLE PRECISION Z(N), W(N), RWORK(LRWORK)
+ 3
+ The time taken by the routine is approximately proportional to n
 1: IFLAG  INTEGER Input/Output
 On entry: IFLAG is always nonnegative. On exit: IFLAG
 may be used as a flag to indicate a failure in the
 computation of w. If IFLAG is negative on exit from
 IMAGE, then F02FJF will exit immediately with IFAIL set
 to IFLAG.
+ 9. Example
 2: N  INTEGER Input
 On entry: n, the number of elements in the vectors w
 and z, and the order of the matrix C.
+ To calculate all the eigenvalues of the complex Hermitian matrix:
 3: Z(N)  DOUBLE PRECISION array Input
 On entry: the vector z for which Cz is required.
+ (0.50 0.00 1.84+1.38i 2.081.56i)
+ (0.00 0.50 1.12+0.84i 0.56+0.42i)
+ (1.841.38i 1.120.84i 0.50 0.00 ).
+ (2.08+1.56i 0.560.42i 0.00 0.50 )
 4: W(N)  DOUBLE PRECISION array Output
 On exit: the vector w=Cz.
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available online.
 5: RWORK(LRWORK)  DOUBLE PRECISION array User Workspace
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 6: LRWORK  INTEGER Input
+ F02  Eigenvalue and Eigenvectors F02AXF
+ F02AXF  NAG Foundation Library Routine Document
 7: IWORK(LIWORK)  INTEGER array User Workspace
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementationdependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
 8: LIWORK  INTEGER Input
 IMAGE is called from F02FJF with the parameters RWORK,
 LRWORK, IWORK and LIWORK as supplied to F02FJF. The
 user is free to use the arrays RWORK and IWORK to
 supply information to IMAGE and DOT as an alternative
 to using COMMON.
 IMAGE must be declared as EXTERNAL in the (sub)program
 from which F02FJF is called. Parameters denoted as
 Input must not be changed by this procedure.
+ 1. Purpose
 8: MONIT  SUBROUTINE, supplied by the user.
 External Procedure
 MONIT is used to monitor the progress of F02FJF. MONIT may
 be the dummy subroutine F02FJZ if no monitoring is actually
 required. (F02FJZ is included in the NAG Foundation Library
 and so need not be supplied by the user. The routine name
 F02FJZ may be implementation dependent: see the Users' Note
 for your implementation for details.) MONIT is called after
 the solution of each eigenvalue subproblem and also just
 prior to return from F02FJF. The parameters ISTATE and
 NEXTIT allow selective printing by MONIT.
+ F02AXF calculates all the eigenvalues and eigenvectors of a
+ complex Hermitian matrix.
 Its specification is:
+ 2. Specification
 SUBROUTINE MONIT (ISTATE, NEXTIT, NEVALS,
 1 NEVECS, K, F, D)
 INTEGER ISTATE, NEXTIT, NEVALS, NEVECS,
 1 K
 DOUBLE PRECISION F(K), D(K)
+ SUBROUTINE F02AXF (AR, IAR, AI, IAI, N, R, VR, IVR, VI,
+ 1 IVI, WK1, WK2, WK3, IFAIL)
+ INTEGER IAR, IAI, N, IVR, IVI, IFAIL
+ DOUBLE PRECISION AR(IAR,N), AI(IAI,N), R(N), VR(IVR,N), VI
+ 1 (IVI,N), WK1(N), WK2(N), WK3(N)
 1: ISTATE  INTEGER Input
 On entry: ISTATE specifies the state of F02FJF and will
 have values as follows:
 ISTATE = 0
 No eigenvalue or eigenvector has just been
 accepted.
+ 3. Description
 ISTATE = 1
 One or more eigenvalues have been accepted since
 the last call to MONIT.
+ The complex Hermitian matrix is first reduced to a real
+ tridiagonal matrix by n2 unitary transformations and a
+ subsequent diagonal transformation. The eigenvalues and
+ eigenvectors are then derived using the QL algorithm, an
+ adaptation of the QR algorithm.
 ISTATE = 2
 One or more eigenvectors have been accepted since
 the last call to MONIT.
+ 4. References
 ISTATE = 3
 One or more eigenvalues and eigenvectors have
 been accepted since the last call to MONIT.
+ [1] Peters G (1967) NPL Algorithms Library. Document No.
+ F2/03/A.
 ISTATE = 4
 Return from F02FJF is about to occur.
+ [2] Peters G (1967) NPL Algorithms Library. Document No.
+ F1/04/A.
 2: NEXTIT  INTEGER Input
 On entry: the number of the next iteration.
+ 5. Parameters
 3: NEVALS  INTEGER Input
 On entry: the number of eigenvalues accepted so far.
+ 1: AR(IAR,N)  DOUBLE PRECISION array Input
+ On entry: the real parts of the elements of the lower
+ triangle of the n by n complex Hermitian matrix A. Elements
+ of the array above the diagonal need not be set. See also
+ Section 8.
 4: NEVECS  INTEGER Input
 On entry: the number of eigenvectors accepted so far.
+ 2: IAR  INTEGER Input
+ On entry:
+ the first dimension of the array AR as declared in the
+ (sub)program from which F02AXF is called.
+ Constraint: IAR >= N.
 5: K  INTEGER Input
 On entry: k, the number of simultaneous iteration
 vectors.
+ 3: AI(IAI,N)  DOUBLE PRECISION array Input
+ On entry: the imaginary parts of the elements of the lower
+ triangle of the n by n complex Hermitian matrix A. Elements
+ of the array above the diagonal need not be set. See also
+ Section 8.
 6: F(K)  DOUBLE PRECISION array Input
 On entry: a vector of error quantities measuring the
 state of convergence of the simultaneous iteration
 vectors. See the parameter TOL of F02FJF above and
 Section 8 for further details. Each element of F is
 initially set to the value 4.0 and an element remains
 at 4.0 until the corresponding vector is tested.
+ 4: IAI  INTEGER Input
+ On entry:
+ the first dimension of the array AI as declared in the
+ (sub)program from which F02AXF is called.
+ Constraint: IAI >= N.
 7: D(K)  DOUBLE PRECISION array Input
 On entry: D(i) contains the latest approximation to the
 absolute value of the ith eigenvalue of C.
 MONIT must be declared as EXTERNAL in the (sub)program
 from which F02FJF is called. Parameters denoted as
 Input must not be changed by this procedure.
+ 5: N  INTEGER Input
+ On entry: n, the order of the matrix, A.
 9: NOVECS  INTEGER Input
 On entry: the number of approximate vectors that are being
 supplied in X. If NOVECS is outside the range (0,K), then
 the value 0 is used in place of NOVECS.
+ 6: R(N)  DOUBLE PRECISION array Output
+ On exit: the eigenvalues in ascending order.
 10: X(NRX,K)  DOUBLE PRECISION array Input/Output
 On entry: if 0 < NOVECS <= K, the first NOVECS columns of X
 must contain approximations to the eigenvectors
 corresponding to the NOVECS eigenvalues of largest absolute
 value of C. Supplying approximate eigenvectors can be useful
 when reasonable approximations are known, or when the
 routine is being restarted with a larger value of K.
 Otherwise it is not necessary to supply approximate vectors,
 as simultaneous iteration vectors will be generated randomly
 by the routine. On exit: if IFAIL = 0, 2, 3 or 4, the first
 m' columns contain the eigenvectors corresponding to the
 eigenvalues returned in the first m' elements of D (see
 below); and the next km'1 columns contain approximations
 to the eigenvectors corresponding to the approximate
 eigenvalues returned in the next km'1 elements of D. Here
 m' is the value returned in M (see above), the number of
 eigenvalues actually found. The kth column is used as
 workspace.
+ 7: VR(IVR,N)  DOUBLE PRECISION array Output
+ On exit: the real parts of the eigenvectors, stored by
+ columns. The ith column corresponds to the ith eigenvector.
+ The eigenvectors are normalised so that the sum of the
+ squares of the moduli of the elements is equal to 1 and the
+ element of largest modulus is real. See also Section 8.
 11: NRX  INTEGER Input
+ 8: IVR  INTEGER Input
On entry:
 the first dimension of the array X as declared in the
 (sub)program from which F02FJF is called.
 Constraint: NRX >= N.

 12: D(K)  DOUBLE PRECISION array Output
 On exit: if IFAIL = 0, 2, 3 or 4, the first m' elements
 contain the first m' eigenvalues in decreasing order of
 magnitude; and the next km'1 elements contain
 approximations to the next km'1 eigenvalues. Here m' is
 the value returned in M (see above), the number of
 eigenvalues actually found. D(k) contains the value e where
 (e,e) is the latest interval over which Chebyshev
 acceleration is performed.

 13: WORK(LWORK)  DOUBLE PRECISION array Workspace

 14: LWORK  INTEGER Input
 On entry: the length of the array WORK, as declared in the
 (sub)program from which F02FJF is called. Constraint:
 LWORK>=3*K+max(K*K,2*N).
+ the first dimension of the array VR as declared in the
+ (sub)program from which F02AXF is called.
+ Constraint: IVR >= N.
 15: RWORK(LRWORK)  DOUBLE PRECISION array User Workspace
 RWORK is not used by F02FJF, but is passed directly to
 routines DOT and IMAGE and may be used to supply information
 to these routines.
+ 9: VI(IVI,N)  DOUBLE PRECISION array Output
+ On exit: the imaginary parts of the eigenvectors, stored by
+ columns. The ith column corresponds to the ith eigenvector.
+ See also Section 8.
 16: LRWORK  INTEGER Input
 On entry: the length of the array RWORK, as declared in the
 (sub)program from which F02FJF is called. Constraint: LRWORK
 >= 1.
+ 10: IVI  INTEGER Input
+ On entry:
+ the first dimension of the array VI as declared in the
+ (sub)program from which F02AXF is called.
+ Constraint: IVI >= N.
 17: IWORK(LIWORK)  INTEGER array User Workspace
 IWORK is not used by F02FJF, but is passed directly to
 routines DOT and IMAGE and may be used to supply information
 to these routines.
+ 11: WK1(N)  DOUBLE PRECISION array Workspace
 18: LIWORK  INTEGER Input
 On entry: the length of the array IWORK, as declared in the
 (sub)program from which F02FJF is called. Constraint: LIWORK
 >= 1.
+ 12: WK2(N)  DOUBLE PRECISION array Workspace
 19: IFAIL  INTEGER Input/Output
 On entry: IFAIL must be set to 0, 1 or 1. Users who are
 unfamiliar with this parameter should refer to the Essential
 Introduction for details.
+ 13: WK3(N)  DOUBLE PRECISION array Workspace
 On exit: IFAIL = 0 unless the routine detects an error or
 gives a warning (see Section 6).
+ 14: IFAIL  INTEGER Input/Output
+ On entry: IFAIL must be set to 0, 1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
 For this routine, because the values of output parameters
 may be useful even if IFAIL /=0 on exit, users are
 recommended to set IFAIL to 1 before entry. It is then
 essential to test the value of IFAIL on exit. To suppress
 the output of an error message when soft failure occurs, set
 IFAIL to 1.
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
6. Error Indicators and Warnings
 Errors or warnings specified by the routine:
+ Errors detected by the routine:

 IFAIL< 0
 A negative value of IFAIL indicates an exit from F02FJF
 because the user has set IFLAG negative in DOT or IMAGE. The
 value of IFAIL will be the same as the user's setting of
 IFLAG.

 IFAIL= 1
 On entry N < 1,

 or M < 1,

 or M >= K,

 or K > N,

 or NRX < N,

 or LWORK <3*K+max(K*K*N),

 or LRWORK < 1,

 or LIWORK < 1.
+ IFAIL= 1
+ More than 30*N iterations are required to isolate all the
+ eigenvalues.
IFAIL= 2
 Not all the requested eigenvalues and vectors have been
 obtained. Approximations to the rth eigenvalue are
 oscillating rapidly indicating that severe cancellation is
 occurring in the rth eigenvector and so M is returned as (r
 1). A restart with a larger value of K may permit
 convergence.

 IFAIL= 3
 Not all the requested eigenvalues and vectors have been
 obtained. The rate of convergence of the remaining
 eigenvectors suggests that more than NOITS iterations would
 be required and so the input value of M has been reduced. A
 restart with a larger value of K may permit convergence.

 IFAIL= 4
 Not all the requested eigenvalues and vectors have been
 obtained. NOITS iterations have been performed. A restart,
 possibly with a larger value of K, may permit convergence.

 IFAIL= 5
 This error is very unlikely to occur, but indicates that
 convergence of the eigenvalue subproblem has not taken
 place. Restarting with a different set of approximate
 vectors may allow convergence. If this error occurs the user
 should check carefully that F02FJF is being called
 correctly.
+ The diagonal elements of AI are not all zero, i.e., the
+ complex matrix is not Hermitian.
7. Accuracy
 Eigenvalues and eigenvectors will normally be computed to the
 accuracy requested by the parameter TOL, but eigenvectors
 corresponding to small or to close eigenvalues may not always be
 computed to the accuracy requested by the parameter TOL. Use of
 the routine MONIT to monitor acceptance of eigenvalues and
 eigenvectors is recommended.
+ The eigenvectors are always accurately orthogonal but the
+ accuracy of the individual eigenvalues and eigenvectors is
+ dependent on their inherent sensitivity to small changes in the
+ original matrix. For a detailed error analysis see Peters [1]
+ page 3 and [2] page 3.
8. Further Comments
 The time taken by the routine will be principally determined by
 the time taken to solve the eigenvalue subproblem and the time
 taken by the routines DOT and IMAGE. The time taken to solve an
 2
 eigenvalue subproblem is approximately proportional to nk . It
 is important to be aware that several calls to DOT and IMAGE may
 occur on each major iteration.

 As can be seen from Table 3.1, many applications of F02FJF will
 require routine IMAGE to solve a system of linear equations. For
 example, to find the smallest eigenvalues of Ax=(lambda)Bx, IMAGE
 needs to solve equations of the form Aw=Bz for w and routines
 from Chapters F01 and F04 of the NAG Foundation Library will
 frequently be useful in this context. In particular, if A is a
 positivedefinite variable band matrix, F04MCF may be used after
 A has been factorized by F01MCF. Thus factorization need be
 performed only once prior to calling F02FJF. An illustration of
 this type of use is given in the example program in Section 9.

 ~
 An approximation d , to the ith eigenvalue, is accepted as soon
 h
 ~
 as d and the previous approximation differ by less than
 h
 ~
 d *TOL/10. Eigenvectors are accepted in groups corresponding to
 h
 clusters of eigenvalues that are equal, or nearly equal, in
 absolute value and that have already been accepted. If d is the
 r
 last eigenvalue in such a group and we define the residual r as
 j

 r =Cx y
 j j r

 where y is the projection of Cx , with respect to B, onto the
 r j
 space spanned by x ,x ,...,x and x is the current approximation
 1 2 r j
 to the jth eigenvector, then the value f returned in MONIT is
 i
 given by

 2 T
 f =maxr  /Cx  x =x Bx
 i j B j B B

 and each vector in the group is accepted as an eigenvector if

 (d f )/(d e)n,
+ 5. Parameters
 D=S, m=n,
+ 1: A(IA,N)  DOUBLE PRECISION array Input/Output
+ On entry: the lower triangle of the n by n symmetric matrix
+ A. The elements of the array above the diagonal need not be
+ set. On exit: the elements of A below the diagonal are
+ overwritten, and the rest of the array is unchanged.
 D=(S 0), m= N.
 Q is an m by m orthogonal matrix, P is an n by n orthogonal
 matrix and S is a min(m,n) by min(m,n) diagonal matrix with non
 negative diagonal elements, sv ,sv ,...,sv , ordered such
 1 2 min(m,n)
 that
+ 3: N  INTEGER Input
+ On entry: n, the order of the matrix A.
 sv >=sv >=...>=sv >=0.
 1 2 min(m,n)
+ 4: ALB  DOUBLE PRECISION Input
 The first min(m,n) columns of Q are the lefthand singular
 vectors of A, the diagonal elements of S are the singular values
 of A and the first min(m,n) columns of P are the righthand
 singular vectors of A.

 Either or both of the lefthand and righthand singular vectors
 of A may be requested and the matrix C given by
+ 5: UB  DOUBLE PRECISION Input
+ On entry: l and u, the lower and upper endpoints of the
+ interval within which eigenvalues are to be calculated.
 T
 C=Q B,
+ 6: M  INTEGER Input
+ On entry: an upper bound for the number of eigenvalues
+ within the interval.
 where B is an m by ncolb given matrix, may also be requested.
+ 7: MM  INTEGER Output
+ On exit: the actual number of eigenvalues within the
+ interval.
 The routine obtains the singular value decomposition by first
 reducing A to upper triangular form by means of Householder
 transformations, from the left when m>=n and from the right when
 m= N.
 T
 KK =I,
+ 11: D(N)  DOUBLE PRECISION array Workspace
 (so that K has elements +1 or 1 on the diagonal)
+ 12: E(N)  DOUBLE PRECISION array Workspace
 then
+ 13: E2(N)  DOUBLE PRECISION array Workspace
 T
 A=(QK)D(PK)
+ 14: X(N,7)  DOUBLE PRECISION array Workspace
 is also a singular value decomposition of A.
+ 15: G(N)  DOUBLE PRECISION array Workspace
 4. References
+ 16: C(N)  LOGICAL array Workspace
 [1] Dongarra J J, Moler C B, Bunch J R and Stewart G W (1979)
 LINPACK Users' Guide. SIAM, Philadelphia.
+ 17: ICOUNT(M)  INTEGER array Output
+ On exit: ICOUNT(i) contains the number of iterations for
+ the ith eigenvalue.
 [2] Hammarling S (1985) The Singular Value Decomposition in
 Multivariate Statistics. ACM Signum Newsletter. 20, 3 225.
+ 18: IFAIL  INTEGER Input/Output
+ On entry: IFAIL must be set to 0, 1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
 [3] Wilkinson J H (1978) Singular Value Decomposition  Basic
 Aspects. Numerical Software  Needs and Availability. (ed D
 A H Jacobs) Academic Press.
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
 5. Parameters
+ 6. Error Indicators and Warnings
 1: M  INTEGER Input
 On entry: the number of rows, m, of the matrix A.
 Constraint: M >= 0.
+ Errors detected by the routine:
 When M = 0 then an immediate return is effected.
+ IFAIL= 1
+ M is less than the number of eigenvalues in the given
+ interval. On exit MM contains the number of eigenvalues in
+ the interval. Rerun with this value for M.
 2: N  INTEGER Input
 On entry: the number of columns, n, of the matrix A.
 Constraint: N >= 0.
+ IFAIL= 2
+ More than 5 iterations are required to determine any one
+ eigenvector.
 When N = 0 then an immediate return is effected.
+ 7. Accuracy
 3: A(LDA,*)  DOUBLE PRECISION array Input/Output
 Note: the second dimension of the array A must be at least
 max(1,N).
 On entry: the leading m by n part of the array A must
 contain the matrix A whose singular value decomposition is
 required. On exit: if M >= N and WANTQ = .TRUE., then the
 leading m by n part of A will contain the first n columns of
 the orthogonal matrix Q.
+ There is no guarantee of the accuracy of the eigenvectors as the
+ results depend on the original matrix and the multiplicity of the
+ roots. For a detailed error analysis see Wilkinson and Reinsch
+ [1] pp 222 and 436.
 If M < N and WANTP = .TRUE., then the leading m by n part of
 T
 A will contain the first m rows of the orthogonal matrix P .
+ 8. Further Comments
 If M >= N and WANTQ = .FALSE. and WANTP = .TRUE., then the
 leading n by n part of A will contain the first n rows of
 T
 the orthogonal matrix P .
+ 3
+ The time taken by the routine is approximately proportional to n
 Otherwise the leading m by n part of A is used as internal
 workspace.
+ This subroutine should only be used when less than 25% of the
+ eigenvalues and the corresponding eigenvectors are required. Also
+ this subroutine is less efficient with matrices which have
+ multiple eigenvalues.
 4: LDA  INTEGER Input
 On entry:
 the first dimension of the array A as declared in the
 (sub)program from which F02WEF is called.
 Constraint: LDA >= max(1,M).
+ 9. Example
 5: NCOLB  INTEGER Input
 On entry: ncolb, the number of columns of the matrix B.
 When NCOLB = 0 the array B is not referenced. Constraint:
 NCOLB >= 0.
+ To calculate the eigenvalues lying between 2.0 and 3.0, and the
+ corresponding eigenvectors of the real symmetric matrix:
 6: B(LDB,*)  DOUBLE PRECISION array Input/Output
 Note: the second dimension of the array B must be at least
 max(1,ncolb) On entry: if NCOLB > 0, the leading m by ncolb
 part of the array B must contain the matrix to be
 transformed. On exit: B is overwritten by the m by ncolb
 T
 matrix Q B.
+ ( 0.5 0.0 2.3 2.6)
+ ( 0.0 0.5 1.4 0.7)
+ ( 2.3 1.4 0.5 0.0).
+ (2.6 0.7 0.0 0.5)
 7: LDB  INTEGER Input
 On entry:
 the first dimension of the array B as declared in the
 (sub)program from which F02WEF is called.
 Constraint: if NCOLB > 0 then LDB >= max(1,M).
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available online.
 8: WANTQ  LOGICAL Input
 On entry: WANTQ must be .TRUE., if the lefthand singular
 vectors are required. If WANTQ = .FALSE., then the array Q
 is not referenced.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 9: Q(LDQ,*)  DOUBLE PRECISION array Output
 Note: the second dimension of the array Q must be at least
 max(1,M).
 On exit: if M < N and WANTQ = .TRUE., the leading m by m
 part of the array Q will contain the orthogonal matrix Q.
 Otherwise the array Q is not referenced.
+ F02  Eigenvalue and Eigenvectors F02BJF
+ F02BJF  NAG Foundation Library Routine Document
 10: LDQ  INTEGER Input
 On entry:
 the first dimension of the array Q as declared in the
 (sub)program from which F02WEF is called.
 Constraint: if M < N and WANTQ = .TRUE., LDQ >= max(1,M).
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementationdependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
 11: SV(*)  DOUBLE PRECISION array Output
 Note: the length of SV must be at least min(M,N). On exit:
 the min(M,N) diagonal elements of the matrix S.
+ 1. Purpose
 12: WANTP  LOGICAL Input
 On entry: WANTP must be .TRUE. if the righthand singular
 vectors are required. If WANTP = .FALSE., then the array PT
 is not referenced.
+ F02BJF calculates all the eigenvalues and, if required, all the
+ eigenvectors of the generalized eigenproblem Ax=(lambda)Bx
+ where A and B are real, square matrices, using the QZ algorithm.
 13: PT(LDPT,*)  DOUBLE PRECISION array Output
 Note: the second dimension of the array PT must be at least
 max(1,N).
 On exit: if M >= N and WANTQ and WANTP are .TRUE., the
 leading n by n part of the array PT will contain the
 T
 orthogonal matrix P . Otherwise the array PT is not
 referenced.
+ 2. Specification
 14: LDPT  INTEGER Input
 On entry:
 the first dimension of the array PT as declared in the
 (sub)program from which F02WEF is called.
 Constraint: if M >= N and WANTQ and WANTP are .TRUE., LDPT
 >= max(1,N).
+ SUBROUTINE F02BJF (N, A, IA, B, IB, EPS1, ALFR, ALFI,
+ 1 BETA, MATV, V, IV, ITER, IFAIL)
+ INTEGER N, IA, IB, IV, ITER(N), IFAIL
+ DOUBLE PRECISION A(IA,N), B(IB,N), EPS1, ALFR(N), ALFI(N),
+ 1 BETA(N), V(IV,N)
+ LOGICAL MATV
 15: WORK(*)  DOUBLE PRECISION array Output
 Note: the length of WORK must be at least max(1,lwork),
 where lwork must be as given in the following table:
+ 3. Description
 M >= N
 WANTQ is .TRUE. and WANTP = .TRUE.
 2
 lwork=max(N +5*(N1),N+NCOLB,4)
+ All the eigenvalues and, if required, all the eigenvectors of the
+ generalized eigenproblem Ax=(lambda)Bx where A and B are real,
+ square matrices, are determined using the QZ algorithm. The QZ
+ algorithm consists of 4 stages:
 WANTQ = .TRUE. and WANTP = .FALSE.
 2
 lwork=max(N +4*(N1),N+NCOLB,4)
+ (a) A is reduced to upper Hessenberg form and at the same time
+ B is reduced to upper triangular form.
 WANTQ = .FALSE. and WANTP = .TRUE.
 lwork=max(3*(N1),2) when NCOLB = 0
+ (b) A is further reduced to quasitriangular form while the
+ triangular form of B is maintained.
 lwork=max(5*(N1),2) when NCOLB > 0
+ (c) The quasitriangular form of A is reduced to triangular
+ form and the eigenvalues extracted. This routine does not
+ actually produce the eigenvalues (lambda) , but instead
+ j
+ returns (alpha) and (beta) such that
 WANTQ = .FALSE. and WANTP = .FALSE.
 lwork=max(2*(N1),2) when NCOLB = 0
+ j j
+ (lambda) =(alpha) /(beta) , j=1,2,...,.n
+ j j j
+ The division by (beta) becomes the responsibility of the
+ j
+ user's program, since (beta) may be zero indicating an
+ j
+ infinite eigenvalue. Pairs of complex eigenvalues occur
+ with (alpha) /(beta) and (alpha) /(beta) complex
+ j j j+1 j+1
 lwork=max(3*(N1),2) when NCOLB > 0
+ conjugates, even though (alpha) and (alpha) are not
+ j j+1
+ conjugate.
 M < N
 WANTQ = .TRUE. and WANTP = .TRUE.
 2
 lwork=max(M +5*(M1),2)
+ (d) If the eigenvectors are required (MATV = .TRUE.), they are
+ obtained from the triangular matrices and then transformed
+ back into the original coordinate system.
 WANTQ = .TRUE. and WANTP = .FALSE.
 lwork=max(3*(M1),1)
+ 4. References
 WANTQ = .FALSE. and WANTP = .TRUE.
 2
 lwork=max(M +3*(M1),2) when NCOLB = 0
+ [1] Moler C B and Stewart G W (1973) An Algorithm for
+ Generalized Matrix Eigenproblems. SIAM J. Numer. Anal. 10
+ 241256.
 2
 lwork=max(M +5*(M1),2) when NCOLB > 0
+ [2] Ward R C (1975) The Combination Shift QZ Algorithm. SIAM J.
+ Numer. Anal. 12 835853.
 WANTQ = .FALSE. and WANTP = .FALSE.
 lwork=max(2*(M1),1) when NCOLB = 0
+ [3] Wilkinson J H (1979) Kronecker's Canonical Form and the QZ
+ Algorithm. Linear Algebra and Appl. 28 285303.
 lwork=max(3*(M1),1) when NCOLB > 0
 On exit: WORK(min(M,N)) contains the total number of
 iterations taken by the R algorithm.
+ 5. Parameters
 The rest of the array is used as workspace.
+ 1: N  INTEGER Input
+ On entry: n, the order of the matrices A and B.
 16: IFAIL  INTEGER Input/Output
 On entry: IFAIL must be set to 0, 1 or 1. For users not
 familiar with this parameter (described in the Essential
 Introduction) the recommended value is 0.
+ 2: A(IA,N)  DOUBLE PRECISION array Input/Output
+ On entry: the n by n matrix A. On exit: the array is
+ overwritten.
 On exit: IFAIL = 0 unless the routine detects an error (see
 Section 6).
+ 3: IA  INTEGER Input
+ On entry:
+ the first dimension of the array A as declared in the
+ (sub)program from which F02BJF is called.
+ Constraint: IA >= N.
 6. Error Indicators and Warnings
+ 4: B(IB,N)  DOUBLE PRECISION array Input/Output
+ On entry: the n by n matrix B. On exit: the array is
+ overwritten.
 Errors detected by the routine:
+ 5: IB  INTEGER Input
+ On entry:
+ the first dimension of the array B as declared in the
+ (sub)program from which F02BJF is called.
+ Constraint: IB >= N.
 If on entry IFAIL = 0 or 1, explanatory error messages are
 output on the current error message unit (as defined by X04AAF).
+ 6: EPS1  DOUBLE PRECISION Input
+ On entry: the tolerance used to determine negligible
+ elements. If EPS1 > 0.0, an element will be considered
+ negligible if it is less than EPS1 times the norm of its
+ matrix. If EPS1 <= 0.0, machine precision is used in place
+ of EPS1. A positive value of EPS1 may result in faster
+ execution but less accurate results.
 IFAIL=1
 One or more of the following conditions holds:
 M < 0,
+ 7: ALFR(N)  DOUBLE PRECISION array Output
 N < 0,
+ 8: ALFI(N)  DOUBLE PRECISION array Output
+ On exit: the real and imaginary parts of (alpha) , for
+ j
+ j=1,2,...,n.
 LDA < M,
+ 9: BETA(N)  DOUBLE PRECISION array Output
+ On exit: (beta) , for j=1,2,...,n.
+ j
 NCOLB < 0,
+ 10: MATV  LOGICAL Input
+ On entry: MATV must be set .TRUE. if the eigenvectors are
+ required, otherwise .FALSE..
 LDB < M and NCOLB > 0,
+ 11: V(IV,N)  DOUBLE PRECISION array Output
+ On exit: if MATV = .TRUE., then
+ (i)if the jth eigenvalue is real, the jth column of V
+ contains its eigenvector;
 LDQ < M and M < N and WANTQ = .TRUE.,
+ (ii) if the jth and (j+1)th eigenvalues form a complex
+ pair, the jth and (j+1)th columns of V contain the
+ real and imaginary parts of the eigenvector associated
+ with the first eigenvalue of the pair. The conjugate
+ of this vector is the eigenvector for the conjugate
+ eigenvalue.
+ Each eigenvector is normalised so that the component of
+ largest modulus is real and the sum of squares of the moduli
+ equal one.
 LDPT < N and M >= N and WANTQ = .TRUE., and WANTP = .
 TRUE..
+ If MATV = .FALSE., V is not used.
 IFAIL> 0
 The QR algorithm has failed to converge in 50*min(m,n)
 iterations. In this case SV(1), SV(2),..., SV(IFAIL) may not
 have been found correctly and the remaining singular values
 may not be the smallest. The matrix A will nevertheless have
 T
 been factorized as A=QEP , where the leading min(m,n) by
 min(m,n) part of E is a bidiagonal matrix with SV(1), SV(2),
 ..., SV(min(m,n)) as the diagonal elements and WORK(1), WORK
 (2),..., WORK(min(m,n)1) as the superdiagonal elements.
+ 12: IV  INTEGER Input
+ On entry:
+ the first dimension of the array V as declared in the
+ (sub)program from which F02BJF is called.
+ Constraint: IV >= N.
 This failure is not likely to occur.
+ 13: ITER(N)  INTEGER array Output
+ On exit: ITER(j) contains the number of iterations needed
+ to obtain the jth eigenvalue. Note that the eigenvalues are
+ obtained in reverse order, starting with the nth.
 7. Accuracy
+ 14: IFAIL  INTEGER Input/Output
+ On entry: IFAIL must be set to 0, 1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
 The computed factors Q, D and P satisfy the relation
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
 T
 QDP =A+E,
+ 6. Error Indicators and Warnings
 where
+ Errors detected by the routine:
 E<=c(epsilon)A,
+ IFAIL= i
+ More than 30*N iterations are required to determine all the
+ diagonal 1 by 1 or 2 by 2 blocks of the quasitriangular
+ form in the second step of the QZ algorithm. IFAIL is set to
+ the index i of the eigenvalue at which this failure occurs.
+ If the soft failure option is used, (alpha) and (beta) are
+ j j
+ correct for j=i+1,i+2,...,n, but V does not contain any
+ correct eigenvectors.
 (epsilon) being the machine precision, c is a modest function of
 m and n and . denotes the spectral (two) norm. Note that
 A=sv .
 1
+ 7. Accuracy
 8. Further Comments
+ The computed eigenvalues are always exact for a problem
+ (A+E)x=(lambda)(B+F)x where E/A and F/B
+ are both of the order of max(EPS1,(epsilon)), EPS1 being defined
+ as in Section 5 and (epsilon) being the machine precision.
 Following the use of this routine the rank of A may be estimated
 by a call to the INTEGER FUNCTION F06KLF(*). The statement:
+ Note: interpretation of results obtained with the QZ algorithm
+ often requires a clear understanding of the effects of small
+ changes in the original data. These effects are reviewed in
+ Wilkinson [3], in relation to the significance of small values of
+ (alpha) and (beta) . It should be noted that if (alpha) and
+ j j j
+ (beta) are both small for any j, it may be that no reliance can
+ j
+ be placed on any of the computed eigenvalues
+ (lambda) =(alpha) /(beta) . The user is recommended to study [3]
+ i i i
+ and, if in difficulty, to seek expert advice on determining the
+ sensitivity of the eigenvalues to perturbations in the data.
 IRANK = F06KLF(MIN(M, N), SV, 1, TOL)
+ 8. Further Comments
 returns the value (k1) in IRANK, where k is the smallest integer
 for which SV(k)n,
+ T
+ BC=C B (2)
 D=S, m=n,
+ for a given positivedefinite matrix B. C is said to be B
+ symmetric. Different specifications of C allow for the solution
+ of a variety of eigenvalue problems. For example, when
 D=(S 0), m=sv >=...>=sv >=0.
 1 2 min(m,n)

 The first min(m,n) columns of Q are the lefthand singular
 vectors of A, the diagonal elements of S are the singular values
 of A and the first min(m,n) columns of P are the righthand
 singular vectors of A.

 Either or both of the lefthand and righthand singular vectors
 of A may be requested and the matrix C given by

 H
 C=Q B,

 where B is an m by ncolb given matrix, may also be requested.

 The routine obtains the singular value decomposition by first
 reducing A to upper triangular form by means of Householder
 transformations, from the left when m>=n and from the right when
 m= 0.
+ To find the m eigenvalues of smallest absolute magnitude of (3)
+ 1
+ we can choose C=A and hence find the reciprocals of the
+ required eigenvalues, so that IMAGE will need to solve Aw=z for
+ 1
+ w, and correspondingly for (4) we can choose C=A B and solve
+ Aw=Bz for w.
 When M = 0 then an immediate return is effected.
+ A table of examples of choice of IMAGE is given in Table 3.1. It
+ should be remembered that the routine also returns the
+ corresponding eigenvectors and that B is positivedefinite.
+ Throughout A is assumed to be symmetric and, where necessary,
+ nonsingularity is also assumed.
