\documentclass{article} \usepackage{axiom} \begin{document} \title{\$SPAD/src/algebra nepip.as} \author{Michael Richardson} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{NagEigenInterfacePackage} <>= +++ Author: M.G. Richardson +++ Date Created: 1996 January 12 +++ Date Last Updated: +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ This package provides Axiom-like interfaces to the NAG generalised +++ eigenvalue and eigenvector routines in the NAGlink. DF ==> DoubleFloat ; CDF ==> Complex DoubleFloat ; FFCDF ==> FormalFraction Complex DoubleFloat ; LFFCDF ==> List FormalFraction Complex DoubleFloat ; LDF ==> List DoubleFloat ; LCDF ==> List Complex DoubleFloat ; LLDF ==> List LDF ; VDF ==> Vector DoubleFloat ; LVDF ==> List VDF ; VCDF ==> Vector Complex DoubleFloat ; LVCDF ==> List VCDF ; MDF ==> Matrix DoubleFloat ; MCDF ==> Matrix Complex DoubleFloat ; INT ==> Integer ; NNI ==> NonNegativeInteger ; RCD ==> Record ; RSLT ==> Result ; STRG ==> String ; UNNRES ==> Union(a:LDF,b:LFFCDF) ; -- a & b are dummy tags RURSLV ==> RCD(eigenvalues : UNNRES, eigenvectors : LVCDF) ; NagEigenInterfacePackage: with { nagEigenvalues : (MDF,MDF,DF) -> UNNRES ; ++ nagEigenvalues(A,B,eps) returns a list of the eigenvalues #if saturn ++ $ \lambda $ #else ++ \spad{l} #endif ++ of the system #if saturn ++ $ A x = \lambda B x $ #else ++ \spad{A*x = l*B*x} #endif ++ ++ The numerical calculation is performed by one of the NAG routines ++ F02ADF and F02BJF, depending on the the form of \spad{A} and B. ++ The result is of type Union(List DoubleFloat, List FormalFraction ++ Complex DoubleFloat), the first branch resulting from F02ADF and ++ the second from F02BJF. Note that in the latter case values should ++ be inspected for numerically small numerators and denominators, ++ ratios of which may be in effect indeterminate, before the result is ++ converted to List Complex DoubleFloat. ++ ++ The parameter eps, if positive, defines a tolerance to be used in ++ recognising negligable matrix elements when F02BJF is called; setting ++ this may result in faster execution with less accuracy. ++ ++ For more detailed information, please consult the NAG manual ++ via the Browser pages for the operations f02adf and f02bjf. nagEigenvalues : (MDF,MDF) -> UNNRES ; ++ nagEigenvalues(A,B) returns a list of the eigenvalues #if saturn ++ $ \lambda $ #else ++ \spad{l} #endif ++ of the system #if saturn ++ $ A x = \lambda B x $ #else ++ \spad{A*x = l*B*x} #endif ++ ++ The numerical calculation is performed by one of the NAG routines ++ F02ADF and F02BJF, depending on the the form of \spad{A} and B. ++ The result is of type Union(List DoubleFloat, List FormalFraction ++ Complex DoubleFloat), the first branch resulting from F02ADF and ++ the second from F02BJF. Note that in the latter case values should ++ be inspected for numerically small numerators and denominators, ++ ratios of which may be in effect indeterminate, before the result is ++ converted to List Complex DoubleFloat. ++ ++ For more detailed information, please consult the NAG manual ++ via the Browser pages for the operations f02adf and f02bjf. nagEigenvectors : (MDF,MDF,DF) -> RURSLV ; ++ nagEigenvectors(A,B,eps) returns a record consisting of a list of the ++ eigenvalues #if saturn ++ $ \lambda $ #else ++ \spad{l} #endif ++ and a list of the corresponding eigenvectors of the system #if saturn ++ $ A x = \lambda B x $ #else ++ \spad{A*x = l*B*x} #endif ++ where #if saturn ++ $A$ and $B$ #else ++ \spad{A} and B #endif #if saturn ++ $B$ #else ++ B #endif ++ is positive-definite. ++ ++ The numerical calculation is performed by one of the NAG routines ++ F02AEF and F02BJF, depending on the the form of \spad{A} and B. ++ The first component of the result, \spad{eigenvalues}, ++ is of type Union(List DoubleFloat, List FormalFraction ++ Complex DoubleFloat), the first branch resulting from F02AEF and ++ the second from F02BJF. Note that in the latter case values should ++ be inspected for numerically small numerators and denominators, ++ ratios of which may be in effect indeterminate, before the result is ++ converted to List Complex DoubleFloat. ++ ++ The parameter eps, if positive, defines a tolerance to be used in ++ recognising negligable matrix elements when F02BJF is called; setting ++ this may result in faster execution with less accuracy. ++ ++ For more detailed information, please consult the NAG manual ++ via the Browser pages for the operations f02aef and f02bjf. nagEigenvectors : (MDF,MDF) -> RURSLV ; ++ nagEigenvectors(A,B) returns a record consisting of a list of the ++ eigenvalues #if saturn ++ $ \lambda $ #else ++ \spad{l} #endif ++ and a list of the corresponding eigenvectors of the system #if saturn ++ $ A x = \lambda B x $ #else ++ \spad{A*x = l*B*x} #endif ++ where #if saturn ++ $A$ and $B$ #else ++ \spad{A} and B #endif #if saturn ++ $B$ #else ++ B #endif ++ is positive-definite. ++ ++ The numerical calculation is performed by one of the NAG routines ++ F02AEF and F02BJF, depending on the the form of \spad{A} and B. ++ The first component of the result, \spad{eigenvalues}, ++ is of type Union(List DoubleFloat, List FormalFraction ++ Complex DoubleFloat), the first branch resulting from F02AEF and ++ the second from F02BJF. Note that in the latter case values should ++ be inspected for numerically small numerators and denominators, ++ ratios of which may be in effect indeterminate, before the result is ++ converted to List Complex DoubleFloat. ++ ++ For more detailed information, please consult the NAG manual ++ via the Browser pages for the operations f02aef and f02bjf. } == add { import from AnyFunctions1 INT ; import from AnyFunctions1 MDF ; import from CDF; import from ErrorFunctions ; import from MDF ; import from NagResultChecks ; import from NagEigenPackage ; import from List STRG ; import from Symbol ; import from VDF ; import from Boolean ; import from Result ; local (..)(a:INT,b:INT):Generator INT == { generate { t := a ; while (t <= b) repeat { yield t ; t := t + 1 ; } } } local ipIfail : INT := -1 ; -- First, some local functions: f02bjfEigVals(A : MDF, B : MDF, orderAB : INT, eps : DF) : LFFCDF == { -- orderAB is the common order of the square matrices A and B. local f02bjfResult : RSLT ; local numR, numI, den : LDF ; f02bjfResult := f02bjf(orderAB,orderAB,orderAB,eps, false,orderAB,A,B,ipIfail) ; den := entries(row(checkMxDF(f02bjfResult, "beta", "F02BJF"), 1)) ; numR := entries(row(retract(f02bjfResult."alfr") @ MDF, 1)) ; numI := entries(row(retract(f02bjfResult."alfi") @ MDF, 1)) ; [ (complex(r,i)/complex(d,0@DF))$FFCDF for r in numR for i in numI for d in den ] } f02bjfEigVecs(A : MDF, B : MDF, orderAB : INT, eps : DF) : RURSLV == { -- orderAB is the common order of the square matrices A and B. local size : NNI ; local f02bjfResult : RSLT ; local numR, numI, den : LDF ; local eVals : UNNRES ; local vecMat : MDF ; local eVecs : LVCDF ; local j : INT ; local thisVec, leftVec : VCDF ; size := orderAB pretend NNI ; f02bjfResult := f02bjf(orderAB,orderAB,orderAB,eps, true,orderAB,A,B,ipIfail) ; den := entries(row(checkMxDF(f02bjfResult, "beta", "F02BJF"), 1)) ; numR := entries(row(retract(f02bjfResult."alfr") @ MDF, 1)) ; numI := entries(row(retract(f02bjfResult."alfi") @ MDF, 1)) ; vecMat := retract(f02bjfResult."v") @ MDF ; -- outer [] for union type: eVals := [[(complex(r,i)/complex(d,0@DF))$FFCDF for r in numR for i in numI for d in den]] ; eVecs := [] ; j := orderAB ; while j > 0 repeat { if numI.j ~= 0$DF then { if j = 1 or numI.(j-1) = 0$DF then error("nagEigenvectors", "Inconsistent results returned from NAG routine F02BJF") ; thisVec := new(size,0) ; leftVec := new(size,0) ; for i in 1 .. orderAB repeat { thisVec.i := complex(vecMat(i,j-1),-vecMat(i,j)) ; leftVec.i := complex(vecMat(i,j-1),vecMat(i,j)) ; } eVecs := cons(leftVec,cons(thisVec,eVecs)) ; j := j - 2; } else { thisVec := new(size,0) ; for i in 1 .. orderAB repeat thisVec.i := complex(vecMat(i,j),0@DF) ; eVecs := cons(thisVec,eVecs) ; j := j - 1 ; } } [eVals,eVecs] } nagError(routine : STRG, opIfail : INT) : Exit == error ["An error was detected when calling the NAG Library routine ", routine, ". The error number (IFAIL value) is ", string(opIfail), ", please consult the NAG manual via the Browser for", " diagnostic information."] ; -- Then the exported functions: nagEigenvalues(A : MDF, B : MDF, eps : DF) : UNNRES == { -- Strategy: if either matrix is asymmetric, use F02BJF, otherwise -- try F02ADF in case B is positive-definite. -- If F02ADF gives IFAIL=1 (should happen quickly if at all), -- not positive-definite, use less efficient F02BJF. local rA, rB, cA, cB : NNI ; local orderAB, opIfail: INT ; local vals : UNNRES ; rA := nrows A ; rB := nrows B ; if rA ~= rB then error("nagEigenvalues", "the two matrices supplied are of different sizes.") ; orderAB := rA pretend INT ; if symmetric?(A) and symmetric?(B) then { f02adfResult := f02adf(orderAB,orderAB,orderAB,A,B,ipIfail) ; opIfail := retract(f02adfResult."ifail") @ INT ; if zero? opIfail then -- using [] to give union type: vals := [entries(row(retract(f02adfResult."r") @ MDF,1))] ; else if opIfail = 1 then vals := [f02bjfEigVals(A,B,orderAB,eps)] else nagError("F02BJF",opIfail) ; } else { cA := ncols A ; if cA ~= rA then error("nagEigenvalues", "the first matrix supplied is not square") ; cB := ncols B ; if cB ~= rB then error("nagEigenvalues", "the second matrix supplied is not square") ; vals := [f02bjfEigVals(A,B,orderAB,eps)] ; } vals } nagEigenvalues(A : MDF, B : MDF) : UNNRES == nagEigenvalues(A,B,0@DF) ; nagEigenvectors(A : MDF, B : MDF, eps : DF) : RURSLV == { -- Strategy: if either matrix is asymmetric, use F02BJF, otherwise -- try F02AEF in case B is positive-definite. -- If F02AEF gives IFAIL=1 (should happen quickly if at all), -- not positive-definite, use less efficient F02BJF. local rA, rB, cA, cB : NNI ; local orderAB, opIfail, j : INT ; local eVals : UNNRES ; local eVecs : LVCDF ; local vecMat : MDF ; local thisVec : VCDF ; local f02aefResult : RSLT ; local result : RURSLV ; rA := nrows A ; rB := nrows B ; if rA ~= rB then error("nagEigenvectors", "the two matrices supplied are of different sizes.") ; orderAB := rA pretend INT ; if symmetric?(A) and symmetric?(B) then { f02aefResult := f02aef(orderAB,orderAB,orderAB, orderAB,A,B,ipIfail) ; opIfail := retract(f02aefResult."ifail") @ INT ; if zero? opIfail then { -- using [] to give union type: eVals := [entries(row(retract(f02aefResult."r") @ MDF,1))] ; vecMat := retract(f02aefResult."v") @ MDF ; eVecs := [] ; j := orderAB ; while j > 0 repeat { thisVec := new(rA,0) ; for i in 1 .. orderAB repeat thisVec.i := complex(vecMat(i,j),0@DF) ; eVecs := cons(thisVec,eVecs) ; j := j - 1 ; } result := [eVals,eVecs] } else if opIfail = 1 then result := f02bjfEigVecs(A,B,orderAB,eps) else nagError("F02BJF",opIfail) ; } else { cA := ncols A ; if cA ~= rA then error("nagEigenvectors", "the first matrix supplied is not square") ; cB := ncols B ; if cB ~= rB then error("nagEigenvectors", "the second matrix supplied is not square") ; result := f02bjfEigVecs(A,B,orderAB,eps) ; } result } nagEigenvectors(A : MDF, B : MDF) : RURSLV == nagEigenvectors(A,B,0@DF) ; } #if NeverAssertThis -- Note that the conversions of results from DoubleFloat to Float -- will become unnecessary if outputGeneral is extended to apply to -- DoubleFloat quantities. )lib nrc )lib ffrac )lib nepip outputGeneral 5 mA1 := matrix [[ 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]]; mB1 := matrix[[ 1 , 3 , 4 , 1], _ [ 3 , 13 , 16 , 11], _ [ 4 , 16 , 24 , 18], _ [ 1 , 11 , 18 , 27]]; mA2 := matrix [[ 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]]; mB2 := matrix[[ 1 , 2 , -3 , 1], _ [ 1 , 3 , -5 , 4], _ [ 1 , 3 , -4 , 3], _ [ 1 , 3 , -4 , 4]]; nagEigenvalues(mA1,mB1) :: List Float -- [- 3.0,- 1.0,2.0,4.0] vv1 := nagEigenvectors(mA1,mB1); (vv1.eigenvalues) :: List Float -- [- 3.0,- 1.0,2.0,4.0] (vv1.eigenvectors) :: List Vector Complex Float -- [[- 4.35,0.05,1.0,- 0.5], [- 2.05,0.15,0.5,- 0.5], [- 3.95,0.85,0.5,- 0.5], -- [2.65,0.05,- 1.0,0.5]] nagEigenvalues(mA2,mB2) -- all components are O(1) or more so: % :: List Complex Float -- [2.0,3.0 + 4.0 %i,3.0 - 4.0 %i,4.0] vv2 := nagEigenvectors(mA2,mB2); vv2.eigenvalues -- all components are O(1) or more so: % :: List Complex Float -- [2.0,3.0 + 4.0 %i,3.0 - 4.0 %i,4.0] vv2.eigenvectors :: List Vector Complex Float -- [[0.99606,0.0056917,0.062609,0.062609], -- -- [0.94491, 0.18898 + 0.26077 E -14 %i, 0.11339 - 0.15119 %i, -- 0.11339 - 0.15119 %i] -- , -- -- [0.94491, 0.18898 - 0.26077 E -14 %i, 0.11339 + 0.15119 %i, -- 0.11339 + 0.15119 %i] -- , -- [0.98752,0.010972,- 0.032917,0.15361]] -- The same call with eps=0.0001: vv2a := nagEigenvectors(mA2,mB2,0.0001); vv2a.eigenvalues :: List Complex Float -- [1.9989,3.0003 + 3.9994 %i,3.0003 - 3.9994 %i,4.0] vv2a.eigenvectors :: List Vector Complex Float -- [[0.99605,0.0057355,0.062656,0.062656], -- -- [0.94491, 0.18899 - 0.000048882 %i, 0.11336 - 0.15119 %i, -- 0.11336 - 0.15119 %i] -- , -- -- [0.94491, 0.18899 + 0.000048882 %i, 0.11336 + 0.15119 %i, -- 0.11336 + 0.15119 %i] -- , -- [0.98751,0.011031,- 0.032912,0.15367]] mB1(1,1) := -1; -- The next test should fail on F02ADF then call F02BJF: nagEigenvalues(mA1,mB1) -- all components are O(1) or more so: % :: List Complex Float -- [3.5016,- 1.5471,0.041212 + 0.21738 %i,0.041212 - 0.21738 %i] -- Similarly, this should fail on F02AEF then call F02BJF: vv3 := nagEigenvectors(mA1,mB1); vv3.eigenvalues -- all components are O(1) or more so: % :: List Complex Float -- [3.5016,- 1.5471,0.041212 + 0.21738 %i,0.041212 - 0.21738 %i] vv3.eigenvectors :: List Vector Complex Float -- [[- 0.034577,0.63045,- 0.75202,0.1892], -- [0.17876,- 0.73845,0.047413,0.64845], -- -- [0.80838, - 0.00095133 + 0.47557 %i, - 0.20354 - 0.21737 %i, -- 0.15404 + 0.089179 %i] -- , -- -- [0.80838, - 0.00095133 - 0.47557 %i, - 0.20354 + 0.21737 %i, -- 0.15404 - 0.089179 %i] -- ] outputGeneral() output "End of tests" #endif @ \section{License} <>= --Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- - Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- -- - Redistributions in binary form must reproduce the above copyright -- notice, this list of conditions and the following disclaimer in -- the documentation and/or other materials provided with the -- distribution. -- -- - Neither the name of The Numerical ALgorithms Group Ltd. nor the -- names of its contributors may be used to endorse or promote products -- derived from this software without specific prior written permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS --IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED --TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A --PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER --OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, --EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, --PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR --PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF --LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING --NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS --SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @ <<*>>= <> -- NagEigenProblemInterfacePackage -- To test: -- sed '1,/^#if NeverAssertThis/d;/#endif/d' < nepip.as > nepip.input -- axiom -- )set nag host -- )r nepip.input #unassert saturn #include "axiom.as" <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}