Don't like ads? PRO users don't see any ads ;-)
Guest

Untitled

By: a guest on Jul 15th, 2012  |  syntax: None  |  size: 20.42 KB  |  hits: 12  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. using System;
  2. using System.Text.RegularExpressions;
  3.  
  4. namespace MatrixApplication.Common
  5. {
  6.     public class Matrix
  7.     {
  8.         public int rows;
  9.         public int cols;
  10.         public double[,] mat;
  11.  
  12.         public Matrix L;
  13.         public Matrix U;
  14.         private int[] pi;
  15.         private double detOfP = 1;
  16.  
  17.         public Matrix(int iRows, int iCols)         // Matrix Class constructor
  18.         {
  19.             rows = iRows;
  20.             cols = iCols;
  21.             mat = new double[rows, cols];
  22.         }
  23.  
  24.         public Boolean IsSquare()
  25.         {
  26.             return (rows == cols);
  27.         }
  28.  
  29.         public double this[int iRow, int iCol]      // Access this matrix as a 2D array
  30.         {
  31.             get { return mat[iRow, iCol]; }
  32.             set { mat[iRow, iCol] = value; }
  33.         }
  34.  
  35.         public Matrix GetCol(int k)
  36.         {
  37.             Matrix m = new Matrix(rows, 1);
  38.             for (int i = 0; i < rows; i++) m[i, 0] = mat[i, k];
  39.             return m;
  40.         }
  41.  
  42.         public void SetCol(Matrix v, int k)
  43.         {
  44.             for (int i = 0; i < rows; i++) mat[i, k] = v[i, 0];
  45.         }
  46.  
  47.         public void MakeLU()                        // Function for LU decomposition
  48.         {
  49.             if (!IsSquare()) throw new MException("The matrix is not square!");
  50.             L = IdentityMatrix(rows, cols);
  51.             U = Duplicate();
  52.  
  53.             pi = new int[rows];
  54.             for (int i = 0; i < rows; i++) pi[i] = i;
  55.  
  56.             double p = 0;
  57.             double pom2;
  58.             int k0 = 0;
  59.             int pom1 = 0;
  60.  
  61.             for (int k = 0; k < cols - 1; k++)
  62.             {
  63.                 p = 0;
  64.                 for (int i = k; i < rows; i++)      // find the row with the biggest pivot
  65.                 {
  66.                     if (Math.Abs(U[i, k]) > p)
  67.                     {
  68.                         p = Math.Abs(U[i, k]);
  69.                         k0 = i;
  70.                     }
  71.                 }
  72.                 if (p == 0) // Matrix is singular!
  73.                     throw new MException("The matrix is singular!");
  74.  
  75.                 pom1 = pi[k]; pi[k] = pi[k0]; pi[k0] = pom1;    // switch two rows in permutation matrix
  76.  
  77.                 for (int i = 0; i < k; i++)
  78.                 {
  79.                     pom2 = L[k, i]; L[k, i] = L[k0, i]; L[k0, i] = pom2;
  80.                 }
  81.  
  82.                 if (k != k0) detOfP *= -1;
  83.  
  84.                 for (int i = 0; i < cols; i++)                  // Switch rows in U
  85.                 {
  86.                     pom2 = U[k, i]; U[k, i] = U[k0, i]; U[k0, i] = pom2;
  87.                 }
  88.  
  89.                 for (int i = k + 1; i < rows; i++)
  90.                 {
  91.                     L[i, k] = U[i, k] / U[k, k];
  92.                     for (int j = k; j < cols; j++)
  93.                         U[i, j] = U[i, j] - L[i, k] * U[k, j];
  94.                 }
  95.             }
  96.         }
  97.  
  98.  
  99.         public Matrix SolveWith(Matrix v)                        // Function solves Ax = v in confirmity with solution vector "v"
  100.         {
  101.             if (rows != cols) throw new MException("The matrix is not square!");
  102.             if (rows != v.rows) throw new MException("Wrong number of results in solution vector!");
  103.             if (L == null) MakeLU();
  104.  
  105.             Matrix b = new Matrix(rows, cols);
  106.             for (int i = 0; i < rows; i++) b[pi[i], 0] = v[i, 0];   // switch two items in "v" due to permutation matrix
  107.  
  108.             Matrix z = SubsForth(L, b);
  109.             Matrix x = SubsBack(U, z);
  110.  
  111.             return x;
  112.         }
  113.  
  114.         public Matrix Invert()                                   // Function returns the inverted matrix
  115.         {
  116.             if (L == null) MakeLU();
  117.  
  118.             Matrix inv = new Matrix(rows, cols);
  119.  
  120.             for (int i = 0; i < rows; i++)
  121.             {
  122.                 Matrix Ei = Matrix.ZeroMatrix(rows, 1);
  123.                 Ei[pi[pi[i]], 0] = 1;
  124.                 Matrix col = SolveWith(Ei);
  125.                 inv.SetCol(col, i);
  126.             }
  127.             return inv;
  128.         }
  129.  
  130.  
  131.         public double Det()                         // Function for determinant
  132.         {
  133.             if (L == null) MakeLU();
  134.             double det = detOfP;
  135.             for (int i = 0; i < rows; i++) det *= U[i, i];
  136.             return det;
  137.         }
  138.  
  139.         public Matrix GetP()                        // Function returns permutation matrix "P" due to permutation vector "pi"
  140.         {
  141.             if (L == null) MakeLU();
  142.  
  143.             Matrix matrix = ZeroMatrix(rows, cols);
  144.             for (int i = 0; i < rows; i++) matrix[i, pi[i]] = 1;
  145.             return matrix;
  146.         }
  147.  
  148.         public Matrix Duplicate()                   // Function returns the copy of this matrix
  149.         {
  150.             Matrix matrix = new Matrix(rows, cols);
  151.             for (int i = 0; i < rows; i++)
  152.                 for (int j = 0; j < cols; j++)
  153.                     matrix[i, j] = mat[i, j];
  154.             return matrix;
  155.         }
  156.  
  157.         public static Matrix SubsForth(Matrix A, Matrix b)          // Function solves Ax = b for A as a lower triangular matrix
  158.         {
  159.             if (A.L == null) A.MakeLU();
  160.             int n = A.rows;
  161.             Matrix x = new Matrix(n, 1);
  162.  
  163.             for (int i = 0; i < n; i++)
  164.             {
  165.                 x[i, 0] = b[i, 0];
  166.                 for (int j = 0; j < i; j++) x[i, 0] -= A[i, j] * x[j, 0];
  167.                 x[i, 0] = x[i, 0] / A[i, i];
  168.             }
  169.             return x;
  170.         }
  171.  
  172.         public static Matrix SubsBack(Matrix A, Matrix b)           // Function solves Ax = b for A as an upper triangular matrix
  173.         {
  174.             if (A.L == null) A.MakeLU();
  175.             int n = A.rows;
  176.             Matrix x = new Matrix(n, 1);
  177.  
  178.             for (int i = n - 1; i > -1; i--)
  179.             {
  180.                 x[i, 0] = b[i, 0];
  181.                 for (int j = n - 1; j > i; j--) x[i, 0] -= A[i, j] * x[j, 0];
  182.                 x[i, 0] = x[i, 0] / A[i, i];
  183.             }
  184.             return x;
  185.         }
  186.  
  187.         public static Matrix ZeroMatrix(int iRows, int iCols)       // Function generates the zero matrix
  188.         {
  189.             Matrix matrix = new Matrix(iRows, iCols);
  190.             for (int i = 0; i < iRows; i++)
  191.                 for (int j = 0; j < iCols; j++)
  192.                     matrix[i, j] = 0;
  193.             return matrix;
  194.         }
  195.  
  196.         public static Matrix IdentityMatrix(int iRows, int iCols)   // Function generates the identity matrix
  197.         {
  198.             Matrix matrix = ZeroMatrix(iRows, iCols);
  199.             for (int i = 0; i < Math.Min(iRows, iCols); i++)
  200.                 matrix[i, i] = 1;
  201.             return matrix;
  202.         }
  203.  
  204.         public static Matrix RandomMatrix(int iRows, int iCols, int dispersion)       // Function generates the zero matrix
  205.         {
  206.             Random random = new Random();
  207.             Matrix matrix = new Matrix(iRows, iCols);
  208.             for (int i = 0; i < iRows; i++)
  209.                 for (int j = 0; j < iCols; j++)
  210.                     matrix[i, j] = random.Next(-dispersion, dispersion);
  211.             return matrix;
  212.         }
  213.  
  214.         public static Matrix Parse(string s)                        // Function parses the matrix from string
  215.         {
  216.             string[] rows = Regex.Split(s, "\r\n");
  217.             string[] nums = rows[0].Split(' ');
  218.             Matrix matrix = new Matrix(rows.Length, nums.Length);
  219.             try
  220.             {
  221.                 for (int i = 0; i < rows.Length; i++)
  222.                 {
  223.                     nums = rows[i].Split(' ');
  224.                     for (int j = 0; j < nums.Length; j++) matrix[i, j] = double.Parse(nums[j]);
  225.                 }
  226.             }
  227.             catch (FormatException exc) { throw new MException("Wrong input format!"); }
  228.             return matrix;
  229.         }
  230.  
  231.         public override string ToString()                           // Function returns matrix as a string
  232.         {
  233.             string s = "";
  234.             for (int i = 0; i < rows; i++)
  235.             {
  236.                 for (int j = 0; j < cols; j++) s += String.Format("{0,5:0.00}", mat[i, j]) + " ";
  237.                 s += "\r\n";
  238.             }
  239.             return s;
  240.         }
  241.  
  242.         public static Matrix Transpose(Matrix m)              // Matrix transpose
  243.         {
  244.             Matrix t = new Matrix(m.cols, m.rows);
  245.             for (int i = 0; i < m.rows; i++)
  246.                 for (int j = 0; j < m.cols; j++)
  247.                     t[i, j] = m[j, i];
  248.             return t;
  249.         }
  250.  
  251.         public static Matrix Power(Matrix m, int pow)           // Power matrix to exponent
  252.         {
  253.             if (pow == 0) return IdentityMatrix(m.rows, m.cols);
  254.             if (pow == 1) return m.Duplicate();
  255.             if (pow == -1) return m.Invert();
  256.  
  257.             Matrix x;
  258.             if (pow < 0) { x = m.Invert(); pow *= -1; }
  259.             else x = m.Duplicate();
  260.  
  261.             Matrix ret = IdentityMatrix(m.rows, m.cols);
  262.             while (pow != 0)
  263.             {
  264.                 if ((pow & 1) == 1) ret *= x;
  265.                 x *= x;
  266.                 pow >>= 1;
  267.             }
  268.             return ret;
  269.         }
  270.  
  271.         private static void SafeAplusBintoC(Matrix A, int xa, int ya, Matrix B, int xb, int yb, Matrix C, int size)
  272.         {
  273.             for (int i = 0; i < size; i++)          // rows
  274.                 for (int j = 0; j < size; j++)     // cols
  275.                 {
  276.                     C[i, j] = 0;
  277.                     if (xa + j < A.cols && ya + i < A.rows) C[i, j] += A[ya + i, xa + j];
  278.                     if (xb + j < B.cols && yb + i < B.rows) C[i, j] += B[yb + i, xb + j];
  279.                 }
  280.         }
  281.  
  282.         private static void SafeAminusBintoC(Matrix A, int xa, int ya, Matrix B, int xb, int yb, Matrix C, int size)
  283.         {
  284.             for (int i = 0; i < size; i++)          // rows
  285.                 for (int j = 0; j < size; j++)     // cols
  286.                 {
  287.                     C[i, j] = 0;
  288.                     if (xa + j < A.cols && ya + i < A.rows) C[i, j] += A[ya + i, xa + j];
  289.                     if (xb + j < B.cols && yb + i < B.rows) C[i, j] -= B[yb + i, xb + j];
  290.                 }
  291.         }
  292.  
  293.         private static void SafeACopytoC(Matrix A, int xa, int ya, Matrix C, int size)
  294.         {
  295.             for (int i = 0; i < size; i++)          // rows
  296.                 for (int j = 0; j < size; j++)     // cols
  297.                 {
  298.                     C[i, j] = 0;
  299.                     if (xa + j < A.cols && ya + i < A.rows) C[i, j] += A[ya + i, xa + j];
  300.                 }
  301.         }
  302.  
  303.         private static void AplusBintoC(Matrix A, int xa, int ya, Matrix B, int xb, int yb, Matrix C, int size)
  304.         {
  305.             for (int i = 0; i < size; i++)          // rows
  306.                 for (int j = 0; j < size; j++) C[i, j] = A[ya + i, xa + j] + B[yb + i, xb + j];
  307.         }
  308.  
  309.         private static void AminusBintoC(Matrix A, int xa, int ya, Matrix B, int xb, int yb, Matrix C, int size)
  310.         {
  311.             for (int i = 0; i < size; i++)          // rows
  312.                 for (int j = 0; j < size; j++) C[i, j] = A[ya + i, xa + j] - B[yb + i, xb + j];
  313.         }
  314.  
  315.         private static void ACopytoC(Matrix A, int xa, int ya, Matrix C, int size)
  316.         {
  317.             for (int i = 0; i < size; i++)          // rows
  318.                 for (int j = 0; j < size; j++) C[i, j] = A[ya + i, xa + j];
  319.         }
  320.  
  321.         private static Matrix StrassenMultiply(Matrix A, Matrix B)                // Smart matrix multiplication
  322.         {
  323.             if (A.cols != B.rows) throw new MException("Wrong dimension of matrix!");
  324.  
  325.             Matrix R;
  326.  
  327.             int msize = Math.Max(Math.Max(A.rows, A.cols), Math.Max(B.rows, B.cols));
  328.  
  329.             if (msize < 32)
  330.             {
  331.                 R = ZeroMatrix(A.rows, B.cols);
  332.                 for (int i = 0; i < R.rows; i++)
  333.                     for (int j = 0; j < R.cols; j++)
  334.                         for (int k = 0; k < A.cols; k++)
  335.                             R[i, j] += A[i, k] * B[k, j];
  336.                 return R;
  337.             }
  338.  
  339.             int size = 1; int n = 0;
  340.             while (msize > size) { size *= 2; n++; };
  341.             int h = size / 2;
  342.  
  343.  
  344.             Matrix[,] mField = new Matrix[n, 9];
  345.  
  346.             /*
  347.              *  8x8, 8x8, 8x8, ...
  348.              *  4x4, 4x4, 4x4, ...
  349.              *  2x2, 2x2, 2x2, ...
  350.              *  . . .
  351.              */
  352.  
  353.             int z;
  354.             for (int i = 0; i < n - 4; i++)          // rows
  355.             {
  356.                 z = (int)Math.Pow(2, n - i - 1);
  357.                 for (int j = 0; j < 9; j++) mField[i, j] = new Matrix(z, z);
  358.             }
  359.  
  360.             SafeAplusBintoC(A, 0, 0, A, h, h, mField[0, 0], h);
  361.             SafeAplusBintoC(B, 0, 0, B, h, h, mField[0, 1], h);
  362.             StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 1], 1, mField); // (A11 + A22) * (B11 + B22);
  363.  
  364.             SafeAplusBintoC(A, 0, h, A, h, h, mField[0, 0], h);
  365.             SafeACopytoC(B, 0, 0, mField[0, 1], h);
  366.             StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 2], 1, mField); // (A21 + A22) * B11;
  367.  
  368.             SafeACopytoC(A, 0, 0, mField[0, 0], h);
  369.             SafeAminusBintoC(B, h, 0, B, h, h, mField[0, 1], h);
  370.             StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 3], 1, mField); //A11 * (B12 - B22);
  371.  
  372.             SafeACopytoC(A, h, h, mField[0, 0], h);
  373.             SafeAminusBintoC(B, 0, h, B, 0, 0, mField[0, 1], h);
  374.             StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 4], 1, mField); //A22 * (B21 - B11);
  375.  
  376.             SafeAplusBintoC(A, 0, 0, A, h, 0, mField[0, 0], h);
  377.             SafeACopytoC(B, h, h, mField[0, 1], h);
  378.             StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 5], 1, mField); //(A11 + A12) * B22;
  379.  
  380.             SafeAminusBintoC(A, 0, h, A, 0, 0, mField[0, 0], h);
  381.             SafeAplusBintoC(B, 0, 0, B, h, 0, mField[0, 1], h);
  382.             StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 6], 1, mField); //(A21 - A11) * (B11 + B12);
  383.  
  384.             SafeAminusBintoC(A, h, 0, A, h, h, mField[0, 0], h);
  385.             SafeAplusBintoC(B, 0, h, B, h, h, mField[0, 1], h);
  386.             StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 7], 1, mField); // (A12 - A22) * (B21 + B22);
  387.  
  388.             R = new Matrix(A.rows, B.cols);                  // result
  389.  
  390.             /// C11
  391.             for (int i = 0; i < Math.Min(h, R.rows); i++)          // rows
  392.                 for (int j = 0; j < Math.Min(h, R.cols); j++)     // cols
  393.                     R[i, j] = mField[0, 1 + 1][i, j] + mField[0, 1 + 4][i, j] - mField[0, 1 + 5][i, j] + mField[0, 1 + 7][i, j];
  394.  
  395.             /// C12
  396.             for (int i = 0; i < Math.Min(h, R.rows); i++)          // rows
  397.                 for (int j = h; j < Math.Min(2 * h, R.cols); j++)     // cols
  398.                     R[i, j] = mField[0, 1 + 3][i, j - h] + mField[0, 1 + 5][i, j - h];
  399.  
  400.             /// C21
  401.             for (int i = h; i < Math.Min(2 * h, R.rows); i++)          // rows
  402.                 for (int j = 0; j < Math.Min(h, R.cols); j++)     // cols
  403.                     R[i, j] = mField[0, 1 + 2][i - h, j] + mField[0, 1 + 4][i - h, j];
  404.  
  405.             /// C22
  406.             for (int i = h; i < Math.Min(2 * h, R.rows); i++)          // rows
  407.                 for (int j = h; j < Math.Min(2 * h, R.cols); j++)     // cols
  408.                     R[i, j] = mField[0, 1 + 1][i - h, j - h] - mField[0, 1 + 2][i - h, j - h] + mField[0, 1 + 3][i - h, j - h] + mField[0, 1 + 6][i - h, j - h];
  409.  
  410.             return R;
  411.         }
  412.  
  413.         // function for square matrix 2^N x 2^N
  414.  
  415.         private static void StrassenMultiplyRun(Matrix A, Matrix B, Matrix C, int l, Matrix[,] f)    // A * B into C, level of recursion, matrix field
  416.         {
  417.             int size = A.rows;
  418.             int h = size / 2;
  419.  
  420.             if (size < 32)
  421.             {
  422.                 for (int i = 0; i < C.rows; i++)
  423.                     for (int j = 0; j < C.cols; j++)
  424.                     {
  425.                         C[i, j] = 0;
  426.                         for (int k = 0; k < A.cols; k++) C[i, j] += A[i, k] * B[k, j];
  427.                     }
  428.                 return;
  429.             }
  430.  
  431.             AplusBintoC(A, 0, 0, A, h, h, f[l, 0], h);
  432.             AplusBintoC(B, 0, 0, B, h, h, f[l, 1], h);
  433.             StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 1], l + 1, f); // (A11 + A22) * (B11 + B22);
  434.  
  435.             AplusBintoC(A, 0, h, A, h, h, f[l, 0], h);
  436.             ACopytoC(B, 0, 0, f[l, 1], h);
  437.             StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 2], l + 1, f); // (A21 + A22) * B11;
  438.  
  439.             ACopytoC(A, 0, 0, f[l, 0], h);
  440.             AminusBintoC(B, h, 0, B, h, h, f[l, 1], h);
  441.             StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 3], l + 1, f); //A11 * (B12 - B22);
  442.  
  443.             ACopytoC(A, h, h, f[l, 0], h);
  444.             AminusBintoC(B, 0, h, B, 0, 0, f[l, 1], h);
  445.             StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 4], l + 1, f); //A22 * (B21 - B11);
  446.  
  447.             AplusBintoC(A, 0, 0, A, h, 0, f[l, 0], h);
  448.             ACopytoC(B, h, h, f[l, 1], h);
  449.             StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 5], l + 1, f); //(A11 + A12) * B22;
  450.  
  451.             AminusBintoC(A, 0, h, A, 0, 0, f[l, 0], h);
  452.             AplusBintoC(B, 0, 0, B, h, 0, f[l, 1], h);
  453.             StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 6], l + 1, f); //(A21 - A11) * (B11 + B12);
  454.  
  455.             AminusBintoC(A, h, 0, A, h, h, f[l, 0], h);
  456.             AplusBintoC(B, 0, h, B, h, h, f[l, 1], h);
  457.             StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 7], l + 1, f); // (A12 - A22) * (B21 + B22);
  458.  
  459.             /// C11
  460.             for (int i = 0; i < h; i++)          // rows
  461.                 for (int j = 0; j < h; j++)     // cols
  462.                     C[i, j] = f[l, 1 + 1][i, j] + f[l, 1 + 4][i, j] - f[l, 1 + 5][i, j] + f[l, 1 + 7][i, j];
  463.  
  464.             /// C12
  465.             for (int i = 0; i < h; i++)          // rows
  466.                 for (int j = h; j < size; j++)     // cols
  467.                     C[i, j] = f[l, 1 + 3][i, j - h] + f[l, 1 + 5][i, j - h];
  468.  
  469.             /// C21
  470.             for (int i = h; i < size; i++)          // rows
  471.                 for (int j = 0; j < h; j++)     // cols
  472.                     C[i, j] = f[l, 1 + 2][i - h, j] + f[l, 1 + 4][i - h, j];
  473.  
  474.             /// C22
  475.             for (int i = h; i < size; i++)          // rows
  476.                 for (int j = h; j < size; j++)     // cols
  477.                     C[i, j] = f[l, 1 + 1][i - h, j - h] - f[l, 1 + 2][i - h, j - h] + f[l, 1 + 3][i - h, j - h] + f[l, 1 + 6][i - h, j - h];
  478.         }
  479.  
  480.         public static Matrix StupidMultiply(Matrix m1, Matrix m2)                  // Stupid matrix multiplication
  481.         {
  482.             if (m1.cols != m2.rows) throw new MException("Wrong dimensions of matrix!");
  483.  
  484.             Matrix result = ZeroMatrix(m1.rows, m2.cols);
  485.             for (int i = 0; i < result.rows; i++)
  486.                 for (int j = 0; j < result.cols; j++)
  487.                     for (int k = 0; k < m1.cols; k++)
  488.                         result[i, j] += m1[i, k] * m2[k, j];
  489.             return result;
  490.         }
  491.         private static Matrix Multiply(double n, Matrix m)                          // Multiplication by constant n
  492.         {
  493.             Matrix r = new Matrix(m.rows, m.cols);
  494.             for (int i = 0; i < m.rows; i++)
  495.                 for (int j = 0; j < m.cols; j++)
  496.                     r[i, j] = m[i, j] * n;
  497.             return r;
  498.         }
  499.         private static Matrix Add(Matrix m1, Matrix m2)
  500.         {
  501.             if (m1.rows != m2.rows || m1.cols != m2.cols) throw new MException("Matrices must have the same dimensions!");
  502.             Matrix r = new Matrix(m1.rows, m1.cols);
  503.             for (int i = 0; i < r.rows; i++)
  504.                 for (int j = 0; j < r.cols; j++)
  505.                     r[i, j] = m1[i, j] + m2[i, j];
  506.             return r;
  507.         }
  508.  
  509.         //   O P E R A T O R S
  510.  
  511.         public static Matrix operator -(Matrix m)
  512.         { return Matrix.Multiply(-1, m); }
  513.  
  514.         public static Matrix operator +(Matrix m1, Matrix m2)
  515.         { return Matrix.Add(m1, m2); }
  516.  
  517.         public static Matrix operator -(Matrix m1, Matrix m2)
  518.         { return Matrix.Add(m1, -m2); }
  519.  
  520.         public static Matrix operator *(Matrix m1, Matrix m2)
  521.         { return Matrix.StrassenMultiply(m1, m2); }
  522.  
  523.         public static Matrix operator *(double n, Matrix m)
  524.         { return Matrix.Multiply(n, m); }
  525.     }
  526.  
  527.     //  The class for exceptions
  528.  
  529.     public class MException : Exception
  530.     {
  531.         public MException(string Message)
  532.             : base(Message)
  533.         { }
  534.     }
  535. }