(* :Title: Orthogonalization *)

(* :Author: John M. Novak *)

(* :Summary:
This package provides the Gram-Schmidt routine for orthogonalizing a list
of vectors or a list of functions.
*)

(* :Context: LinearAlgebra`Orthogonalization` *)

(* :Package Version: 1.2 *)

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

(* :History: 
	Version 1.0 by John M. Novak (Wolfram Research), January 1991.
	Version 1.1 by John M. Novak, February 1991. Adds ability to handle
		function spaces and ability to generate normalized basis.
    Version 1.2 by John M. Novak, February 1999. In a small cleanup for
        handling linearly dependent case, uses a Developer` context function
*)

(* :Keywords: orthogonalization, linear algebra, matrices, vector spaces *)

(* :Sources: Standard linear algebra and calculus texts *)

(* :Mathematica Version: 4.0 *)

(* :Example:
	Here is an example of using GramSchmidt with an arbitrary inner product
	in a function space:
	GramSchmidt[{1,x,x^2,x^3,x^4},
		InnerProduct->(Integrate[#1 #2,{x,-1,1}]&)]
	This generates the first five Legendre Polynomials.
*)

BeginPackage["LinearAlgebra`Orthogonalization`"]

(* Usage messages *)

Normalize::usage =
	"Normalize[vec] normalizes a vector passed in; note that an inner
	product can be specified with the option InnerProduct to allow, for
	instance, function spaces to be used."

Projection::usage =
	"Projection[v1,v2] projects a vector on another vector. Note that an
	inner product can be specified with the option InnerProduct to allow,
	for instance, function spaces."

GramSchmidt::usage =
	"GramSchmidt[vectors] performs the GramSchmidt orthogonalization process
	on a list of vectors. Note that an inner product can be specified,
	allowing, for instance, a function space to used.  Also, the option
	Normalized can be set to determine whether or not the basis is
	orthonormal."

InnerProduct::usage =
	"InnerProduct is an option to functions in the Orthogonalization
	package specifying an inner product as a pure function.  The
	default is Dot.  For example: InnerProduct->(Integrate[#1 #2,{x,-1,1}]&)."

Normalized::usage =
	"Normalized is an option for GramSchmidt that determines
	whether the basis generated is orthonormal or not.  If
	True, then it is orthonormal, False is returned otherwise."

Begin["`Private`"]

Options[Normalize] = {InnerProduct -> Dot};

Options[Projection] = {InnerProduct -> Dot};

Options[GramSchmidt] = {InnerProduct -> Dot, Normalized -> True};

Normalize[vec_,opts___] :=
	Module[{ip, ipvec},
		{ip} = {InnerProduct}/.{opts}/.Options[Normalize];
		If[Developer`ZeroQ[ipvec = ip[vec, vec]],
			vec,
			vec/Sqrt[ipvec]
		]
	]

Projection[v1_,v2_,opts___] := 
	Module[{ip, ipv2},
		{ip} = {InnerProduct}/.{opts}/.Options[Projection];
		If[Developer`ZeroQ[ipv2 = ip[v2, v2]],
			v2,
			ip[v1, v2] v2/ipv2
		]
	]

(* auxiliary internal function; projection of the vector on a list of vectors.
*)

multipleprojection[v1_, vecs_, opts___] := 
	Plus @@ Map[Projection[v1, #, opts]&, vecs]

(*
The Gram-Schmidt method works by taking a list of vectors; starting with the
first vector, it finds a new basis vector in the orthogonal set by subtracting
from each vector the projection of that vector on the basis vectors determined
so far. Note that a different inner product can be specified, allowing the
use of this function in vector spaces.
*)

GramSchmidt[vecs_List, opts___] :=
	Module[{norm, ip, ans},
		{norm, ip} = {Normalized, InnerProduct}/.{opts}/.
			Options[GramSchmidt];
		ans = Fold[Join[#1,
					{#2 - multipleprojection[#2, #1, InnerProduct->ip]}]&,
				{}, vecs];
		If[TrueQ[norm],
			Map[Normalize[#, InnerProduct -> ip]&,ans],
			ans
		]
	]

End[]

EndPackage[]
