path: root/src/algebra/nepip.as.pamphlet

\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}
<<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}