(* :Name: LinearAlgebra`Vandermonde` *)

(* :Title: Solver for linear Vandermonde systems of equations *)

(* :Author: Mark Sofroniou *)

(* :Summary:
 This package solves a linear Vandermonde system of equations using
 5*n^2/2 flops. It also solves a transposed system in the same complexity.
*)

(* :Context: LinearAlgebra`VandermondeSolve` *)

(* :Package Version: 1.0 *)

(* :Copyright: Copyright 1997-1999, Wolfram Research, Inc. *)

(* :History:
 Version 1.0 by Mark Sofroniou. *)

(* :Keywords:
 Linear Algebra, Vandermonde matrix. *)

(* :Source:
 Gene H. Golub and C. F. Van Loan,
 Matrix Computations 3rd Edition,
 Johns Hopkins University Press, 1996.
*)

(* :Mathematica Version: 3.1 *)

(* :Limitation:
 Modified methods exist for confluent systems, but these have not yet
 been implemented.
 This package works for symbolic inputs but the results can become large.
 The option Simplify can be applied to reduce the size of intermediary
 computations.
 Fast Fourier Transform methods exist and are more efficient for certain
 types of Vandermonde systems.
*)

(* :Discussion:
 Vandermonde systems arise in many approximation and interpolation problems.
 However, the rows and/or columns of Vandermonde matrices can be almost
 collinear resulting in an ill conditioned system and hence numerical
 instability using conventional linear solvers.
 The key to constructing a more efficient and stable solver is to recognize
 the correspondence between solving a linear system and polynomial
 interpolation.
*)

BeginPackage["LinearAlgebra`Vandermonde`"];

VandermondeMatrix::usage =
"VandermondeMatrix[{x1,...,xn}] constructs the Vandermonde matrix
V(x1,...,xn) whose rows are {x1,...,xn}^i, for i=0,...,n-1.";

VandermondeSolve::usage =
"VandermondeSolve[{x1,...,xn}, {b1,...,bn}, opts] solves the
linear system V*x = b for unknowns x, where V(x1,...,xn)
is a Vandermonde matrix.";

VandermondeTransposeSolve::usage =
"VandermondeTransposeSolve[{x1,...,xn}, {b1,...,bn}, opts] solves the
linear system V^T*x = b for unknowns x, where V(x1,...,xn)^T
is the transpose of a Vandermonde matrix.";

Unprotect[VandermondeMatrix, VandermondeSolve, VandermondeTransposeSolve];

Begin["`Private`"];

(* Utility routines and predicates *)

(* Developer context functions *)

pack = Developer`ToPackedArray;
packedQ = Developer`PackedArrayQ;

(* Test for determining when to switching to symbolic code *)

NumberVectorQ[x_?packedQ]:= TensorRank[ x ] == 1;
NumberVectorQ[{__?NumberQ}] = True;
NumberVectorQ[_] = False;
  
(* Get the type of the input vector for compilation *)

PackedType[l_] := 
  Which[
    packedQ[ l, Integer], Integer,
    packedQ[ l, Real], Real,
    _, Complex
  ];


(**** Vandermonde matrix constructor ****)

VandermondeMatrix[x_?VectorQ]:= pack[ Table[ x^i, {i, 0, Length[x]-1}] ];

(**** Argument and option testing ****)

VandermondeSolve::conf =
"The argument `1` in VandermondeSolve contains identical elements.
Unable to solve the resulting confluent system.";

VandermondeSolve::args =
"The arguments `1` and `2` in VandermondeSolve did not evaluate to vectors
of the same length";

VandermondeSolve::opts =
"The option `1` in VandermondeSolve did not evaluate to `2`";

VandermondeTransposeSolve::conf =
"The argument `1` in VandermondeTransposeSolve contains identical elements.
Unable to solve the resulting confluent system.";

VandermondeTransposeSolve::args =
"The arguments `1` and `2` in VandermondeTransposeSolve did not evaluate
to vectors of the same length";

VandermondeTransposeSolve::opts =
"The option `1` in VandermondeTransposeSolve did not evaluate to `2`";

VStest[ vsfun_, bool_, msgvals_ ]:=
  If[ bool,
    True,
    Message[ MessageName[ vsfun, "opts" ], Apply[Sequence, msgvals ] ]; False
  ];

VSmessages = {
  {Compiled, "True, False or Automatic"},
  {Simplify, "a symbol or a function"}
  };

VSOptionsTest[ vsfun_, x_, b_, opts___]:=
  Module[{opttests, optvals},

    (* Check arguments *)

    If[ !VectorQ[x] || !VectorQ[b] || Length[x]=!=Length[b] ||
        !MatchQ[Length[x],_Integer?Positive],
      Message[ MessageName[ vsfun, "args" ], x, b ];
      Return[ $Failed ]
    ];

    (* Check confluence *)

    If[ Length[ Union[x] ] =!= Length[x],
      Message[ MessageName[ vsfun, "conf" ], x ]; Return[ $Failed ] ];

    (* Check options *)

    optvals = {Compiled, Simplify} /. {opts} /.  Options[vsfun];

    opttests = {
      MatchQ[ Part[optvals,1], True|False|Automatic ], (* Compiled *)
      MatchQ[ Part[optvals,2], _Symbol|_Function ] (* Simplify *)
    };

    If[ Apply[ And, MapThread[ VStest[ vsfun, #1, #2 ]&, {opttests, VSmessages}] ],
      optvals,
      $Failed
    ]
  ];

Options[VandermondeSolve] = Options[VandermondeTransposeSolve] =
  {Compiled -> Automatic, Simplify -> Together};

(**** User level interface functions ****)

VandermondeSolve[x_, b_, opts___?OptionQ]:=
  Module[{res},
    res /; ((res = VSOptionsTest[ VandermondeSolve, x, b, opts ]) =!= $Failed) &&
                      (res = VSmain[ VandermondeSolve, x, b, res ]; True)
  ];

VandermondeTransposeSolve[x_, b_, opts___?OptionQ]:=
  Module[{res},
    res /; ((res = VSOptionsTest[ VandermondeTransposeSolve, x, b, opts ]) =!= $Failed) &&
                      (res = VSmain[ VandermondeTransposeSolve, x, b, res ]; True)
  ];

(**** Available solvers ****)

VSsolver[ VandermondeSolve, Compiled ] = cvsol;
VSsolver[ VandermondeSolve, Uncompiled ] = gvsol;
VSsolver[ VandermondeTransposeSolve, Compiled ] = cvtsol;
VSsolver[ VandermondeTransposeSolve, Uncompiled ] = gvtsol;

(**** Dispatching function for compiled solver or general solver ****)

VandermondeSolve::nocomp =
"Warning: Unable to use compiled code, proceeding with uncompiled version.";

VandermondeTransposeSolve::nocomp =
"Warning: Unable to use compiled code, proceeding with uncompiled version.";

VSmain[vsfun_, x_, b_, {compQ_, simp_}]:=
  Module[{pb, px},
    px = pack[ {x, b}, Union ]; (* Pack together to least common type *)
    If[ compQ=!=False && packedQ[px],
      pb = Last[px];
      px = First[px];
      VSsolver[vsfun, Compiled][PackedType[px]][px, pb],
      If[ compQ===True, Message[ MessageName[ vsfun, "nocomp" ] ] ];
      VSsolver[vsfun, Uncompiled][ x, b, simp ]
    ]
  ];

(* Algorithm 4.6.2 from Golub and Van Loan. Doing multiple part
 extractions and assignments in one go using Take speeds things
 up considerably. *)

(* Make and store a compiled function according to the input type. *)

cvsol[type_]:=
  cvsol[type] =
    Compile[{{a, Blank[type], 1}, {b, Blank[type], 1}},
      Module[{x = b, n = Length[a]},
        Do[
          Part[x, Internal`TakePart[i+1, n]] -= 
            Part[a, i]*Take[x, {i, n-1}],
        {i, n-1}];
        Do[
          Part[x, Internal`TakePart[i+1, n]] /=
            (Take[a, {i+1, n}] - Take[a, {1, n-i}]);
          Part[x, Internal`TakePart[i, n-1]] = 
            (Take[x, {i, n-1}] - Take[x, {i+1, n}]),
        {i, n-1, 1, -1}];
        x
	  ]
    ];

(*
 This routine also works for exact entries by applying a
 simplification function to intermediary steps.
*)

gvsol[a_, b_, s_]:=
  Module[{x = b, n = Length[a], simp},
    simp =
      If[ NumberVectorQ[a] && NumberVectorQ[b] || SameQ[s, None], #&, s ];
    Do[
      Part[x, Internal`TakePart[i+1, n]] -= 
        Part[a, i]*Take[x, {i, n-1}],
    {i, n-1}];
    Do[
      Part[x, Internal`TakePart[i+1, n]] /=
        (Take[a, {i+1, n}] - Take[a, {1, n-i}]);
      Part[x, Internal`TakePart[i, n-1]] = 
        simp[ Take[x, {i, n-1}] - Take[x, {i+1, n}] ],
    {i, n-1, 1, -1}];
    x
  ];

(**** Vandermonde transpose solver ****)

(* Algorithm 4.6.1 from Golub and Van Loan. Doing multiple part
 extractions and assignments in one go using Take speeds things
 up considerably. *)

(* Make and store a compiled function according to the input type. *)

cvtsol[type_]:=
  cvtsol[type] =
    Compile[{{a, Blank[type], 1}, {b, Blank[type], 1}},
      Module[{x = b, n = Length[a]},
        Do[
          Part[x, Internal`TakePart[i+1, n]] =
            (Take[x, {i+1, n}] - Take[x, {i, n-1}])/
              (Take[a, {i+1, n}] - Take[a, {1, n-i}]),
        {i, n-1}];
        Do[
          Part[x, Internal`TakePart[i, n-1]] =
            (Take[x, {i, n-1}] - Take[x,{i+1, n}]*Part[a, i]),
        {i, n-1, 1, -1}];
        x
      ]
    ];

(*
 This routine also works for exact entries by applying a
 simplification function to intermediary steps.
*)

gvtsol[a_, b_, s_]:=
  Module[{x = b, n = Length[a], simp},
    simp = 
      If[ NumberVectorQ[a] && NumberVectorQ[x] || SameQ[s, None], #&, s ];
    Do[
      Part[x, Internal`TakePart[i+1, n]] =
	    simp[ (Take[x, {i+1, n}] - Take[x, {i, n-1}])/
              (Take[a, {i+1, n}] - Take[a, {1, n-i}]) ],
    {i, n-1}];
    Do[
      Part[x, Internal`TakePart[i, n-1]] =
	    simp[ Take[x, {i, n-1}] - Take[x,{i+1, n}]*Part[a, i] ],
    {i, n-1, 1, -1}];
    x
  ];

End[ ]; (* End `Private` Context. *)

SetAttributes[
{VandermondeMatrix, VandermondeSolve, VandermondeTransposeSolve},
ReadProtected];

Protect[VandermondeMatrix, VandermondeSolve, VandermondeTransposeSolve];

EndPackage[ ]; (* End package Context. *)

(* Examples:
VandermondeMatrix[ {1, 2, 3, 4} ]
VandermondeSolve[ {1, 2, 3, 4}, {0, -1, 3, 35} ]
VandermondeSolve[ {1., 2., 3., 4.}, {0, -1, 3, 35} ]
VandermondeSolve[ {1, 2, 3, 4}, {0, -1, 3, 35} , Compiled -> False ]
VandermondeSolve[ Range[1., 100., 1.], 1/Range[1., 100., 1.] ]
VandermondeSolve[ {a, b, c}, {d, e, f} ]
VandermondeTransposeSolve[ {1, 2, 3, 4},{10, 26, 58, 112} ]
VandermondeTransposeSolve[ {a, b, c}, {d, e, f} ]
*)


