using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Diagnostics;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;

using Mol;
/* 使用例：Hirbert行列 A を用いた Ax=b を解きます（bの全要素は1）。
            // LuDecomp(A,n) で A の要素は LU 要素で上書きされます。
            // 元の A[i,j] 要素を保存し、検算に使用する場合は * の付いたコメントを参考(A.FieldIndexを使用)にしてください。
            int n = 1000;
            BigMatrix A = new BigMatrix();
            A.OpenRenew("Test.ldb");
            DateTime dt = DateTime.Now;
            for (int i = 0; i < n; ++i)
            {
                for (int j = 0; j < n; ++j)
                {
                    double Aij = 1.0 / (i + j + 1);
                    A[i, j] = Aij;
                    // * A[i,j] 要素を２つ確保。先頭要素は LU 分解で上書きされるため。
                    // A.AddRecord(new int[] { i, j }, new double[] { Aij, Aij });
                }
                
            }
            T.P("書込み時間   = " + (DateTime.Now - dt));
            dt = DateTime.Now;
            int [] ps = BigMatrix.LuDecomp(A,n);
            T.P("因子分解時間 = " + (DateTime.Now - dt));
            double [] b = new double[n];
            for(int i=0;i<n;++i) b[i] = 1.0;
            dt = DateTime.Now;
            double[] x = BigMatrix.LuSolve(A, b, ps, n);
            T.P("解法時間 　  = " + (DateTime.Now - dt));
            double d = 0;
            // * FieldIndex のデフォルト値は 0
            // A.FieldIndex = 1; //* 検算には、元（２番目のデータ配列要素）の A[i,j]要素を使用 
            for (int i = 0; i < n; ++i)
            {
                double s = 0;
                for (int j = 0; j < n; ++j)
                {
                    // A[] は LU で上書きされているので検算に使用できない。
                    double Aij = (1.0 / (i + j + 1));
                    // *データフィールド配列２番目、元の A[i,j] を使用して検算。
                    // double Aij = A[i, j];

                    s += Aij * x[j];
                }
                s = b[i] - s;
                d += s * s;
            }
            T.P("残差平方和 　= " + d);
            A.Verify(false);
            A.Close();
 */

namespace Examples
{
    /// <summary>
    /// 
    /// Copyright (C) 2014 Shigeo Kobayashi. All rights reserved.
    /// 
    /// メモリー上に作成できないような巨大行列(BigMatrix) A を用いて連立方程式 Ax=b を解く例<br/>
    /// 
    /// ※このソースコードは自由に利用してかまいません。<br/>
    /// ※本例は、あくまでも LightDB の「例」であり、現実的な使用に耐えるかどうかの考察は一切していません。<br/>
    /// 　一般には、Mol に用意されている、対称・疎・バンド行列等様々な、連立方程式解法を選択すべきです。
    /// 　メモリーサイズの関係で解くことのできない巨大行列が本例を使用する唯一の選択肢ですが、
    /// 　速度と精度の点で現実的な解が得られるかどうかは保証の限りではありません。
    /// </summary>
    public class BigMatrix : LightDB<int, double>
    {
        /// <summary>
        /// コンストラクター
        /// </summary>
        public BigMatrix()
            : base(2, 2) // データ要素数、Lu 分解では元の要素が上書きされるため、元の要素保存用に要素数２。、
        {

        }

        /// <summary>
        /// 線形連立方程式 Ax = b (x と b は一次元のベクトル)を解く前処理として、正方行列 A を LU 分解します。
        /// A はは結果(LU)で上書きされます。
        /// </summary>
        /// <param name="A">係数行列／LU分解の結果</param>
        /// <param name="n">A のサイズ</param>
        /// <returns>ピボット交換用のインデックス配列(後のLuSolve()で使用)</returns>
        public static int[] LuDecomp(BigMatrix A, int n)
        {
            int[] ps = new int[n];
            double[] scales = new double[n];

            // ピボット要素の選択（行単位）
            for (int i = 0; i < n; ++i)
            {
                ps[i] = i;
                double nrmrow = 0;
                for (int j = 0; j < n; ++j)
                {
                    double v = Math.Abs(A[i, j]);
                    if (v > nrmrow) nrmrow = v;
                }
                if (nrmrow<=double.Epsilon) throw new Exception("Singular matrix in LuDecomp(phase-1):" + i);
                scales[i] = 1.0/ nrmrow;
            }

            // 消去法による三角化
            int n1 = n - 1;
            for (int k = 0; k < n1; ++k)
            {
                double biggst = 0;
                int    pividx = 0;
                for (int i = k; i < n; ++i)
                {
                    double size = Math.Abs(A[ps[i], k] * scales[ps[i]]);
                    if (size > biggst)
                    {
                        biggst = size;
                        pividx = i;
                    }
                }
                if (biggst<=double.Epsilon) throw new Exception("Singular matrix in LuDecomp(phase-2):" + k);
                if (pividx != k)
                {
                    int j = ps[k];
                    ps[k] = ps[pividx];
                    ps[pividx] = j;
                }
                double pivot = A[ps[k], k];
                for (int i = (k + 1); i < n; ++i)
                {
                    double EM = A[ps[i], k] / pivot;
                    A[ps[i], k] = EM;
                    for (int j = (k + 1); j < n; ++j)
                    {
                        A[ps[i], j] -= EM * A[ps[k], j];
                    }
                }
            }
            if (Math.Abs(A[ps[n1], n1])<=double.Epsilon) throw new Exception("Singular matrix in LuDecomp(phase-3)" + n1);
            return ps;
        }

        /// <summary>
        /// LuDecomp()で LU 分解された A を使用して線形連立式 Ax = b を解きます。
        /// </summary>
        /// <param name="A">LU 分解された係数行列</param>
        /// <param name="b">定数ベクトル</param>
        /// <param name="ps">ピボット変換用インデックス</param>
        /// <param name="n">行列のサイズ</param>
        /// <returns>解ベクトル x</returns>
        public static double[] LuSolve(BigMatrix A, double[] b, int[] ps, int n)
        {
            double[] x = new double[n];

            for (int i = 0; i < n; ++i)
            {
                double dot = 0;
                for (int j = 0; j < i; ++j)
                {
                    dot -= A[ps[i], j] * x[j];
                }
                x[i] = b[ps[i]] + dot;
            }
            for (int i = n - 1; i >= 0; --i)
            {
                double dot = 0;
                for (int j = (i + 1); j < n; ++j)
                {
                    dot += A[ps[i], j] * x[j];
                }
                dot = x[i] - dot;
                x[i] = dot / A[ps[i], i];
            }
            return x;
        }
    }
}