diff options
Diffstat (limited to 'src/algebra/nepip.as.pamphlet')
| -rw-r--r-- | src/algebra/nepip.as.pamphlet | 626 |
1 files changed, 0 insertions, 626 deletions
diff --git a/src/algebra/nepip.as.pamphlet b/src/algebra/nepip.as.pamphlet deleted file mode 100644 index 4bfb2b5a..00000000 --- a/src/algebra/nepip.as.pamphlet +++ /dev/null @@ -1,626 +0,0 @@ -\documentclass{article} -\usepackage{open-axiom} -\begin{document} -\title{\$SPAD/src/algebra nepip.as} -\author{Michael Richardson} -\maketitle -\begin{abstract} -\end{abstract} -\eject -\tableofcontents -\eject -\section{NagEigenInterfacePackage} -<<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} -<<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. -@ -<<*>>= -<<license>> - --- NagEigenProblemInterfacePackage - --- To test: --- sed '1,/^#if NeverAssertThis/d;/#endif/d' < nepip.as > nepip.input --- axiom --- )set nag host <some machine running nagd> --- )r nepip.input - -#unassert saturn - -#include "axiom.as" - -<<NagEigenInterfacePackage>> -@ -\eject -\begin{thebibliography}{99} -\bibitem{1} nothing -\end{thebibliography} -\end{document} |