unit SparseTest;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms,
  Dialogs, Basic2, StdCtrls, ComCtrls, ExtCtrls, Sparse, MtxVec, TeEngine,
  Series, TeeProcs, Chart;

type
  TfrmSparseTest = class(TBasicForm2)
    Chart1: TChart;
    Series1: TBarSeries;
    Panel4: TPanel;
    Button1: TButton;
    Label1: TLabel;
    ComboBox1: TComboBox;
    Memo1: TMemo;
    procedure FormCreate(Sender: TObject);
    procedure FormDestroy(Sender: TObject);
    procedure Button1Click(Sender: TObject);
    procedure ComboBox1Change(Sender: TObject);
  private
    { Private declarations }
    SparseA: TSparseMtx;
    A: TMtx;
    B: TVec;
    tmpTime : DWord;
    x1,x2,x3: TVec; // compare solutions
  public
    { Public declarations }
  end;

var
  frmSparseTest: TfrmSparseTest;

implementation

Uses MtxVecEdit, ansicrtl;

{$R *.dfm}

procedure TfrmSparseTest.FormCreate(Sender: TObject);
begin
  inherited;
  With RichEdit1.Lines, RichEdit1 do
  begin
    Clear;
    Add('As of MtxVec 1.5, sparse matrices are supported via the '
      + 'TmtxSparse object. Several methods known from TMtx object are '
      + 'supported.  Some features in Sparse.pas:');

    Paragraph.FirstIndent := 10;
    Paragraph.Numbering := nsBullet;
    Add('Introduces TSparseMtx class capable of handling real and  complex sparse matrices.');
    Add('Encapsulates UMF Pack v 3.2, one of the fastest direct sparse matrix solvers and interfacing MKL 5.1 BLAS III');
    Add('Features iterative sparse solvers (cg, bcg, gmres, ..)');
    Add('Compute inverse of the sparse matrix');
    Add('Many conversion routines to convert to and form different matrix formats '
      + '(triplet, banded, dense, Harwell Boeing)');
    Add('Load and save MatrixMarket file format');
    Add('Very fast banded matrix solvers already available in TMtx');
    Add('Add, subtract, multiply routines for sparse/sparse and dense/sparse matrices');
    Add('PixelDownSample routine to allow fast display of the non-zero '
      + 'pattern of the matrices of huge sizes (1000 000 x 1000 000 elements)');
    Add('Complete set of standard math routines known from TVec and TMtx. (Log, Power, Sin, Cos, etc...)'+#10+#13);

    Paragraph.FirstIndent := 0;
    Paragraph.Numbering := nsNone;
    Add('');
    Add('In this demo the the system of 256 linear equations (sparse) stystem is solved with several methods :');
    Add('');
    Add('1) Sparse matrix is converted to dense matrix and TMtx.LUSolve method is used to solve the system of equations');
    Add('2) TmtxSparse.LUSolve direct method is used to solve system of equations');
    Add('3) TmtxSparse iterative method  is used to solve system of equations');
  end;


  B := TVec.Create;
  A := TMtx.Create;
  x1 := TVec.Create;
  x2 := TVec.Create;
  x3 := TVec.Create;
  SparseA := TSparseMtx.Create;
  SparseA.ReportLevel := prlNone; { no report - fastest }
  SparseA.AutoClearReport := True; { force clear }
  try
    SparseA.LoadFromMatrixMarketFile('pde225.mtx');
  except
  end;
  { random data for Ax=b i.e. b constants }
  B.Size(SparseA.Rows, False);
  B.RandGauss;
end;

procedure TfrmSparseTest.FormDestroy(Sender: TObject);
begin
  A.Destroy;
  B.Destroy;
  x1.Destroy;
  x2.Destroy;
  x3.Destroy;
  SparseA.Destroy;
  inherited;
end;

procedure TfrmSparseTest.Button1Click(Sender: TObject);
var i : Integer;
  tmpStr: String;
begin
  Screen.Cursor := crHourGlass;
  Series1.Clear;
  { Ordinary Mtx.LuSolve }
  x1.Size(B);
  SparseA.SparseToDense(A);
  tmpTime := GetTickCount;
  for i := 1 to 100 do A.LUSolve(B,x1);
  tmpTime := GetTickCount - tmpTime;
  Series1.Add(tmpTime,'LUSolve',clRed);

  { "regular" sparse solve }
  x2.Size(B);
  tmpTime := GetTickCount;
  for i := 1 to 100 do SparseA.Solve(B,x2);
  { load report from ansicrtl RtlStrings global variable }
  Memo1.Lines.Clear();
  MtxVec.Report.Position := 0;
  Memo1.Lines.LoadFromStream(MtxVec.Report);
  tmpTime := GetTickCount - tmpTime;
  Series1.Add(tmpTime,'Sparse'+chr(13)+'Solve',clBlue);

  { iterative sparse solve }
  x3.Size(B);
  tmpTime := GetTickCount;
  SparseA.IterativeMethod := itmLUGMRES;
  for i := 1 to 100 do
  begin
       x3.Ramp; {Approximate solution}
       SparseA.SolveIterative(B,x3);
  end;
  tmpTime := GetTickCount - tmpTime;
  Series1.Add(tmpTime,'Sparse'+chr(13)+'IterativeSolve',clYellow);
  Screen.Cursor := crDefault;
  if x1.Equal(x2,1.0e-3) then tmpStr := 'LU Solve solution equal to Sparse solve solution'
  else tmpStr := 'LU Solve solution NOT equal to Sparse solve solution';

  if x1.Equal(x3,1.0e-3) then tmpStr := tmpStr + chr(13)+'LU Solve solution equal to iterative solve solution'
  else tmpStr := tmpStr + chr(13) + 'LU Solve solution NOT equal to iterative solve solution';

  if x2.Equal(x3,1.0e-3) then tmpStr := tmpStr + chr(13)+'sparse solve solution equal to iterative solve solution'
  else tmpStr := tmpStr + chr(13) + 'sparse solve solution NOT equal to iterative solve solution';
  ShowMessage(tmpStr);
end;

procedure TfrmSparseTest.ComboBox1Change(Sender: TObject);
begin
  SparseA.ReportLevel := TReportLevel(ComboBox1.ItemIndex);
end;

initialization
  RegisterClass(TfrmSparseTest);

end.
