(*
Copyright (C) 2004 Evgenii Rudnyi, http://Evgenii.Rudnyi.Ru

This software is a copyrighted work licensed under the terms, described
in the file "FREE_LICENSE". 
*)

(* Title: Slicot.m *)
(* Context: Imtek`Slicot*)
(* Author: Evgenii Rudnyi *)
(* Date: December 2004, Freiburg *)
(* Summary: This is a package for calling Slicot for model reduciton.*)
(* Package Copyright: GNU GPL *)
(* Package Version: 1.0 *)
(* Mathematica Version: 5.01 *)
(* History: *)
(* Keywords: *)
(* Sources: *)
(* Warnings: *)
(* Limitations: *)
(* Discussion: *)
(* Requirements: *)
(* Examples: See documentations *)
(* *)

(* Whereever the GNU GPL is not applicable, 
  the software should be used in the same spirit. *)

(* Users of this code must verify correctness for their application. *)
(* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
  MA 02111-1307 USA *)
(* Disclaimer: *)
(* <one line to give the program's name and a brief idea of what it does.> *)
(* Copyright (C) 2004 Evgenii Rudnyi*)

(* This program is free software; *)
(* you can redistribute it and or modify it under the terms of the GNU General Public License *)
(* as published by the Free Software Foundation; either version 2 of the License,*)
(* or (at your option) any later version.This program is distributed in the \
hope that *)
(* it will be useful, but WITHOUT ANY WARRANTY; *)
(* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *)
(* See the GNU General Public License for more details. 
      You should have received a copy of *)
(* the GNU General Public License along with this program; if not, 
  write to the *)
(* Free Software Foundation,Inc.,59 Temple Place,Suite 330,Boston,
  MA 02111-1307 USA *)


BeginPackage["Imtek`Slicot`", {"Imtek`Post4MOR`", "Imtek`Slicot`LowLevel`"}];

CallSlicot::usage="CallSlicot[sys, dim, ops] gives a reduced system of dimension dim. CallSlicot[sys, precision, ops] gives a reduced system with given precision. sys must be an explicit system of first order. Function returns a list {reduced system, Hankel singular values, {a used and optimal setting for TemporaryStorage}}. Options: BalancingSquareRoot, Equlibration, Tolerance, TemporaryStorage, Method (\"BTA\", \"SPA\", \"HNA\").";
BalancingSquareRoot::usage="an option for CallSlicot: True or False.";
Equlibration::usage="an option for CallSlicot: True or False.";
Tolerance::usage="an option for CallSlicot if Method is SPA or HNA: Real value.";
TemporaryStorage::usage="an option for CallSlicot: Integer.";

Begin["`Private`"];

Options[CallSlicot] = {BalancingSquareRoot->True, Equlibration->False, Tolerance->0., TemporaryStorage->0, Method->"BTA"};

CallSlicot[sys_DynamicSystem, dim_Integer, ops___] := 
	callSlicot[sys, dim, 0., "F", ops] /; ExplicitSystemQ[sys]

CallSlicot[sys_DynamicSystem, precision_Real, ops___] := 
	callSlicot[sys, 0, precision, "A", ops] /; ExplicitSystemQ[sys]


callSlicot[sys_, dim_, tol1_, ordsel_, ops___] := Module[
	{meth, res, job, equil, n, m, p, matA, matB, matC, matD, tol, dwork},
	meth = Method /. {ops} /. Options[CallSlicot];
	job = BalancingSquareRoot /. {ops} /. Options[CallSlicot] /. {True->"B", False->"N"};
	equil = Equlibration /. {ops} /. Options[CallSlicot] /. {True->"S", False->"N"};
	n = Length[MatrixK[sys]];
	m = Last[Dimensions[MatrixB[sys]]];
	p = Length[MatrixC[sys]];
	matA = Flatten[Transpose[Normal[N[MatrixA[sys]]]]];
	matB = Flatten[Transpose[Normal[N[MatrixB[sys]]]]];
	matC = Flatten[Transpose[Normal[N[MatrixC[sys]]]]];
	matD = Flatten[Transpose[Normal[N[MatrixD[sys]]]]];
	tol2 = Tolerance /. {ops} /. Options[CallSlicot];
	dwork = TemporaryStorage /. {ops} /. Options[CallSlicot];
	res = Switch[meth,
		"BTA",
			RunBTA[job, equil, ordsel, n, m, p, dim, matA, matB, matC, tol1, dwork],
		"SPA",
			RunSPA[job, equil, ordsel, n, m, p, dim, matA, matB, matC, matD, tol1, tol2, dwork],
		"HNA",
			RunHNA[equil, ordsel, n, m, p, dim, matA, matB, matC, matD, tol1, tol2, dwork],
		_,
			Throw["Unknown method "<>meth]
	];
	Switch[Length[res],
		6,
			If[res[[1]] === 0, Throw["Dimension of reduced model is 0"]];
			{MakeFirstOrderSystem[
				Partition[res[[2]], res[[1]]],
				Partition[res[[3]], m],
				Partition[res[[4]], res[[1]]],
				OutputNames[sys],
				Imtek`Post4MOR`Private`originalMatrixD[sys]
			],
			res[[5]],
			res[[6]]},
		7,
			If[res[[1]] === 0, Throw["Dimension of reduced model is 0"]];
			{MakeFirstOrderSystem[
				Partition[res[[2]], res[[1]]],
				Partition[res[[3]], m],
				Partition[res[[4]], res[[1]]],
				OutputNames[sys],
				Partition[res[[5]], m]
			],
			res[[6]],
			res[[7]]},
		_,
			Throw[res]
	]
]

End[];

EndPackage[];