VIS2 SS2013 CVD DVR
 All Classes Namespaces Functions Enumerations Properties
cMatrixLib.cs
1 //------------------------------------------------------------------------
2 //
3 // Author : Anas Abidi
4 // Date : 18 Dec 2004
5 // Version : 2.0
6 // Description : Matrix Operations Library
7 //
8 //------------------------------------------------------------------------
9 using System;
10 
11 namespace MatrixLibrary
12 {
13  #region "Exception in the Library"
14  class MatrixLibraryExceptions: ApplicationException
15  {public MatrixLibraryExceptions (string message) : base(message){}}
16 
17  // The Exceptions in this Class
18  class MatrixNullException: ApplicationException
19  {
20  public MatrixNullException():
21  base("To do this operation, matrix can not be null"){}
22  }
23  class MatrixDimensionException: ApplicationException
24  {
25  public MatrixDimensionException():
26  base("Dimension of the two matrices not suitable for this operation !"){}
27  }
28  class MatrixNotSquare: ApplicationException
29  {
30  public MatrixNotSquare():
31  base("To do this operation, matrix must be a square matrix !"){}
32  }
33  class MatrixDeterminentZero: ApplicationException
34  {
35  public MatrixDeterminentZero():
36  base("Determinent of matrix equals zero, inverse can't be found !"){}
37  }
38  class VectorDimensionException: ApplicationException
39  {
40  public VectorDimensionException():
41  base("Dimension of matrix must be [3 , 1] to do this operation !"){}
42  }
43  class MatrixSingularException: ApplicationException
44  {
45  public MatrixSingularException():
46  base("Matrix is singular this operation cannot continue !"){}
47  }
48  #endregion
49 
59  public class Matrix
60  {
61  private double[,] in_Mat;
62 
63  #region "Class Constructor and Indexer"
64  public Matrix(int noRows, int noCols)
71  {
72  this.in_Mat = new double[noRows, noCols];
73  }
74 
80  public Matrix(double [,] Mat)
81  {
82  this.in_Mat = (double[,])Mat.Clone();
83  }
84 
88  public double this[int Row, int Col]
89  {
90  get{return this.in_Mat[Row, Col];}
91  set{this.in_Mat[Row,Col] = value;}
92  }
93  #endregion
94 
95  #region "Public Properties"
96  public int NoRows
102  {
103  get{return this.in_Mat.GetUpperBound(0)+1;}
104  set{this.in_Mat = new double[value, this.in_Mat.GetUpperBound(0)];}
105  }
106 
112  public int NoCols
113  {
114  get{return this.in_Mat.GetUpperBound(1)+1;}
115  set{this.in_Mat = new double[this.in_Mat.GetUpperBound(0),value];}
116  }
117 
121  public double[,] toArray
122  {
123  get{return this.in_Mat;}
124  }
125  #endregion
126 
127  #region "Find Rows and Columns in a Matrix"
128  private static void Find_R_C(double[] Mat, out int Row)
129  {
130  Row = Mat.GetUpperBound(0);
131  }
132 
133  private static void Find_R_C(double[,] Mat, out int Row, out int Col)
134  {
135  Row = Mat.GetUpperBound(0);
136  Col = Mat.GetUpperBound(1);
137  }
138  #endregion
139 
140  #region "Change 1D array ([n]) to/from 2D array [n,1]"
141 
151  public static double[,] OneD_2_TwoD(double[] Mat)
152  {
153  int Rows;
154  //Find The dimensions !!
155  try{Find_R_C(Mat,out Rows);}
156  catch{throw new MatrixNullException();}
157 
158  double[,] Sol = new double[Rows+1,1];
159 
160  for (int i=0; i<=Rows;i++)
161  {
162  Sol[i,0] = Mat[i];
163  }
164 
165  return Sol;
166  }
167 
177  public static double[] TwoD_2_OneD(double[,] Mat)
178  {
179  int Rows;
180  int Cols;
181  //Find The dimensions !!
182  try{Find_R_C(Mat,out Rows, out Cols);}
183  catch{throw new MatrixNullException();}
184 
185  if (Cols!=0) throw new MatrixDimensionException();
186 
187  double[] Sol = new double[Rows+1];
188 
189  for (int i=0; i<=Rows;i++)
190  {
191  Sol[i] = Mat[i,0];
192  }
193  return Sol;
194  }
195  #endregion
196 
197  #region "Identity Matrix"
198  public static double[,] Identity(int n)
204  {
205  double[,] temp = new double[n,n];
206  for (int i=0; i<n;i++) temp[i,i] = 1;
207  return temp;
208  }
209  #endregion
210 
211  #region "Add Matrices"
212  public static double[,] Add(double[,] Mat1, double[,] Mat2)
221  {
222  double[,] sol;
223  int i, j;
224  int Rows1, Cols1;
225  int Rows2, Cols2;
226 
227  //Find The dimensions !!
228  try
229  {
230  Find_R_C(Mat1,out Rows1, out Cols1);
231  Find_R_C(Mat2,out Rows2, out Cols2);
232  }
233  catch
234  {
235  throw new MatrixNullException();
236  }
237 
238  if (Rows1 != Rows2 || Cols1 != Cols2) throw new MatrixDimensionException();
239 
240  sol = new double[Rows1+1,Cols1+1];
241 
242  for (i = 0;i <= Rows1; i++)
243  for (j = 0; j<= Cols1; j++)
244  {
245  sol[i, j] = Mat1[i, j] + Mat2[i, j];
246  }
247 
248  return sol;
249  }
250 
259  public static Matrix Add(Matrix Mat1, Matrix Mat2)
260  {return new Matrix(Add(Mat1.in_Mat,Mat2.in_Mat));}
261 
270  public static Matrix operator+(Matrix Mat1, Matrix Mat2)
271  {return new Matrix(Add(Mat1.in_Mat,Mat2.in_Mat));}
272  #endregion
273 
274  #region "Subtract Matrices"
275 
284  public static double[,] Subtract(double[,] Mat1, double[,] Mat2)
285  {
286  double[,] sol;
287  int i, j;
288  int Rows1, Cols1;
289  int Rows2, Cols2;
290 
291  //Find The dimensions !!
292  try
293  {
294  Find_R_C(Mat1,out Rows1, out Cols1);
295  Find_R_C(Mat2,out Rows2, out Cols2);
296  }
297  catch
298  {
299  throw new MatrixNullException();
300  }
301 
302  if (Rows1 != Rows2 || Cols1 != Cols2) throw new MatrixDimensionException();
303 
304  sol = new double[Rows1+1,Cols1+1];
305 
306  for (i = 0;i <= Rows1; i++)
307  for (j = 0; j<= Cols1; j++)
308  {
309  sol[i, j] = Mat1[i, j] - Mat2[i, j];
310  }
311 
312  return sol;
313  }
314 
323  public static Matrix Subtract(Matrix Mat1, Matrix Mat2)
324  {return new Matrix(Subtract(Mat1.in_Mat,Mat2.in_Mat));}
325 
334  public static Matrix operator-(Matrix Mat1, Matrix Mat2)
335  {return new Matrix(Subtract(Mat1.in_Mat,Mat2.in_Mat));}
336  #endregion
337 
338  #region "Multiply Matrices"
339 
348  public static double[,] Multiply(double[,] Mat1, double[,] Mat2)
349  {
350  double[,] sol;
351  int Rows1, Cols1;
352  int Rows2, Cols2;
353  double MulAdd = 0;
354 
355  try
356  {
357  Find_R_C(Mat1, out Rows1, out Cols1);
358  Find_R_C(Mat2, out Rows2, out Cols2);
359  }
360  catch
361  {
362  throw new MatrixNullException();
363  }
364  if (Cols1 != Rows2) throw new MatrixDimensionException();
365 
366  sol = new double[Rows1+1, Cols2+1];
367 
368  for (int i=0; i<=Rows1; i++)
369  for (int j = 0; j<=Cols2; j++)
370  {
371  for (int l = 0; l<=Cols1; l++)
372  {
373  MulAdd = MulAdd + Mat1[i, l] * Mat2[l, j];
374  }
375  sol[i, j] = MulAdd;
376  MulAdd = 0;
377  }
378  return sol;
379  }
380 
395  public static Matrix Multiply(Matrix Mat1, Matrix Mat2)
396  {
397  if ((Mat1.NoRows==3) && (Mat2.NoRows==3) &&
398  (Mat1.NoCols==1) && (Mat1.NoCols==1))
399  {return new Matrix(CrossProduct(Mat1.in_Mat,Mat2.in_Mat));}
400  else
401  {return new Matrix(Multiply(Mat1.in_Mat,Mat2.in_Mat));}
402  }
403 
418  public static Matrix operator*(Matrix Mat1, Matrix Mat2)
419  {
420  if ((Mat1.NoRows==3) && (Mat2.NoRows==3) &&
421  (Mat1.NoCols==1) && (Mat1.NoCols==1))
422  {
423  return new Matrix(CrossProduct(Mat1.in_Mat,Mat2.in_Mat));
424  }
425  else
426  {
427  return new Matrix(Multiply(Mat1.in_Mat,Mat2.in_Mat));
428  }
429  }
430  #endregion
431 
432  #region "Determinant of a Matrix"
433  public static double Det(double[,] Mat)
442  {
443  int S, k, k1, i, j;
444  double[,] DArray;
445  double save, ArrayK, tmpDet;
446  int Rows, Cols;
447 
448  try
449  {
450  DArray = (double[,])Mat.Clone();
451  Find_R_C(Mat , out Rows, out Cols);
452  }
453  catch{throw new MatrixNullException();}
454 
455  if (Rows != Cols) throw new MatrixNotSquare();
456 
457  S = Rows;
458  tmpDet = 1;
459 
460  for (k = 0; k <= S; k++)
461  {
462  if (DArray[k, k] == 0)
463  {
464  j = k;
465  while ((j < S) && (DArray[k, j] == 0)) j = j + 1;
466  if (DArray[k, j] == 0) return 0;
467  else
468  {
469  for (i = k; i <= S; i++)
470  {
471  save = DArray[i, j];
472  DArray[i, j] = DArray[i, k];
473  DArray[i, k] = save;
474  }
475  }
476  tmpDet = -tmpDet;
477  }
478  ArrayK = DArray[k, k];
479  tmpDet = tmpDet * ArrayK;
480  if (k < S)
481  {
482  k1 = k + 1;
483  for (i = k1; i <= S; i++)
484  {
485  for (j = k1; j <= S; j++)
486  DArray[i, j] = DArray[i, j] - DArray[i, k] * (DArray[k, j] / ArrayK);
487  }
488  }
489  }
490  return tmpDet;
491  }
492 
501  public static double Det(Matrix Mat)
502  {return Det(Mat.in_Mat);}
503  #endregion
504 
505  #region "Inverse of a Matrix"
506  public static double[,] Inverse(double[,] Mat)
516  {
517  double[,] AI, Mat1;
518  double AIN, AF;
519  int Rows, Cols;
520  int LL, LLM, L1, L2, LC, LCA, LCB;
521 
522  try
523  {
524  Find_R_C(Mat, out Rows, out Cols);
525  Mat1 = (double[,])Mat.Clone();
526  }
527  catch{throw new MatrixNullException();}
528 
529  if (Rows != Cols) throw new MatrixNotSquare();
530  if (Det(Mat) == 0) throw new MatrixDeterminentZero();
531 
532  LL = Rows;
533  LLM = Cols;
534  AI = new double[LL+1, LL+1];
535 
536  for (L2 = 0; L2 <= LL; L2++)
537  {
538  for (L1 = 0; L1 <= LL; L1++) AI[L1, L2] = 0;
539  AI[L2, L2] = 1;
540  }
541 
542  for (LC = 0; LC <= LL; LC++)
543  {
544  if (Math.Abs(Mat1[LC, LC]) < 0.0000000001)
545  {
546  for (LCA = LC + 1; LCA <= LL; LCA++)
547  {
548  if (LCA == LC) continue;
549  if (Math.Abs(Mat1[LC, LCA]) > 0.0000000001)
550  {
551  for (LCB = 0; LCB <= LL; LCB++)
552  {
553  Mat1[LCB, LC] = Mat1[LCB, LC] + Mat1[LCB, LCA];
554  AI[LCB, LC] = AI[LCB, LC] + AI[LCB, LCA];
555  }
556  break;
557  }
558  }
559  }
560  AIN = 1 / Mat1[LC, LC];
561  for (LCA = 0; LCA <= LL; LCA++)
562  {
563  Mat1[LCA, LC] = AIN * Mat1[LCA, LC];
564  AI[LCA, LC] = AIN * AI[LCA, LC];
565  }
566 
567  for (LCA = 0; LCA <= LL; LCA++)
568  {
569  if (LCA == LC) continue;
570  AF = Mat1[LC, LCA];
571  for (LCB = 0; LCB <= LL; LCB++)
572  {
573  Mat1[LCB, LCA] = Mat1[LCB, LCA] - AF * Mat1[LCB, LC];
574  AI[LCB, LCA] = AI[LCB, LCA] - AF * AI[LCB, LC];
575  }
576  }
577  }
578  return AI;
579  }
580 
590  public static Matrix Inverse(Matrix Mat)
591  {return new Matrix(Inverse(Mat.in_Mat));}
592  #endregion
593 
594  #region "Transpose of a Matrix"
595  public static double[,] Transpose(double[,] Mat)
602  {
603  double[,] Tr_Mat;
604  int i, j, Rows, Cols;
605 
606  try {Find_R_C(Mat, out Rows, out Cols);}
607  catch{throw new MatrixNullException();}
608 
609  Tr_Mat = new double[Cols+1, Rows+1];
610 
611  for (i = 0; i<=Rows;i++)
612  for (j = 0; j<= Cols;j++) Tr_Mat[j, i] = Mat[i, j];
613 
614  return Tr_Mat;
615  }
616 
623  public static Matrix Transpose(Matrix Mat)
624  {return new Matrix(Transpose(Mat.in_Mat));}
625  #endregion
626 
627  #region "Singula Value Decomposition of a Matrix"
628  public static void SVD(double[,] Mat_, out double[,] S_, out double[,] U_, out double[,] V_)
643  {
644  int Rows, Cols;
645  int m, MP, n, NP;
646  double[] w;
647  double[,] A, v;
648 
649  try {Find_R_C(Mat_, out Rows, out Cols);}
650  catch{throw new MatrixNullException();}
651 
652  m = Rows+1;
653  n = Cols+1;
654 
655  if (m < n)
656  {
657  m=n;
658  MP = NP = n;
659  }
660  else if (m>n)
661  {
662  n=m;
663  NP = MP = m;
664  }
665  else
666  {
667  MP = m;
668  NP = n;
669  }
670 
671  A = new double[m+1,n+1];
672 
673  for (int row = 1; row<=Rows+1;row++)
674  for(int col = 1;col<=Cols+1;col++)
675  {A[row,col] = Mat_[row-1,col-1];}
676 
677  const int NMAX = 100;
678  v = new double[NP+1, NP+1];
679  w = new double[NP+1];
680 
681  int k, l, nm;
682  int flag, i, its, j, jj;
683 
684  double[,] U_temp, S_temp, V_temp;
685  double anorm, c, f, g, h, s, scale, x, y, z;
686  double [] rv1 = new double[NMAX];
687 
688  l=0;
689  nm = 0;
690  g = 0.0;
691  scale = 0.0;
692  anorm = 0.0;
693 
694  for (i=1;i<=n;i++)
695  {
696  l=i+1;
697  rv1[i]=scale*g;
698  g=s=scale=0.0;
699  if (i <= m)
700  {
701  for (k=i;k<=m;k++) scale += Math.Abs(A[k,i]);
702  if (scale!=0)
703  {
704  for (k=i;k<=m;k++)
705  {
706  A[k,i] /= scale;
707  s += A[k,i]*A[k,i];
708  }
709  f=A[i,i];
710  g = -Sign(Math.Sqrt(s),f);
711  h=f*g-s;
712  A[i,i]=f-g;
713  if (i != n)
714  {
715  for (j=l;j<=n;j++)
716  {
717  for (s=0,k=i;k<=m;k++) s += A[k,i]*A[k,j];
718  f=s/h;
719  for (k=i;k<=m;k++) A[k,j] += f*A[k,i];
720  }
721  }
722  for (k=i;k<=m;k++) A[k,i] *= scale;
723  }
724  }
725  w[i]=scale*g;
726  g=s=scale=0.0;
727  if (i <= m && i != n)
728  {
729  for (k=l;k<=n;k++) scale += Math.Abs(A[i,k]);
730  if (scale!=0)
731  {
732  for (k=l;k<=n;k++)
733  {
734  A[i,k] /= scale;
735  s += A[i,k]*A[i,k];
736  }
737  f=A[i,l];
738  g = -Sign(Math.Sqrt(s),f);
739  h=f*g-s;
740  A[i,l]=f-g;
741  for (k=l;k<=n;k++) rv1[k]=A[i,k]/h;
742  if (i != m)
743  {
744  for (j=l;j<=m;j++)
745  {
746  for (s=0.0,k=l;k<=n;k++) s += A[j,k]*A[i,k];
747  for (k=l;k<=n;k++) A[j,k] += s*rv1[k];
748  }
749  }
750  for (k=l;k<=n;k++) A[i,k] *= scale;
751  }
752  }
753  anorm=Math.Max(anorm,(Math.Abs(w[i])+Math.Abs(rv1[i])));
754  }
755  for (i=n;i>=1;i--)
756  {
757  if (i < n)
758  {
759  if (g!=0)
760  {
761  for (j=l;j<=n;j++)
762  v[j,i]=(A[i,j]/A[i,l])/g;
763  for (j=l;j<=n;j++)
764  {
765  for (s=0.0,k=l;k<=n;k++) s += A[i,k]*v[k,j];
766  for (k=l;k<=n;k++) v[k,j] += s*v[k,i];
767  }
768  }
769  for (j=l;j<=n;j++) v[i,j]=v[j,i]=0.0;
770  }
771  v[i,i]=1.0;
772  g=rv1[i];
773  l=i;
774  }
775  for (i=n;i>=1;i--)
776  {
777  l=i+1;
778  g=w[i];
779  if (i < n)
780  for (j=l;j<=n;j++) A[i,j]=0.0;
781  if (g!=0)
782  {
783  g=1.0/g;
784  if (i != n)
785  {
786  for (j=l;j<=n;j++)
787  {
788  for (s=0.0,k=l;k<=m;k++) s += A[k,i]*A[k,j];
789  f=(s/A[i,i])*g;
790  for (k=i;k<=m;k++) A[k,j] += f*A[k,i];
791  }
792  }
793  for (j=i;j<=m;j++) A[j,i] *= g;
794  }
795  else
796  {
797  for (j=i;j<=m;j++) A[j,i]=0.0;
798  }
799  ++A[i,i];
800  }
801  for (k=n;k>=1;k--)
802  {
803  for (its=1;its<=30;its++)
804  {
805  flag=1;
806  for (l=k;l>=1;l--)
807  {
808  nm=l-1;
809  if (Math.Abs(rv1[l])+anorm == anorm)
810  {
811  flag=0;
812  break;
813  }
814  if (Math.Abs(w[nm])+anorm == anorm) break;
815  }
816  if (flag!=0)
817  {
818  c=0.0;
819  s=1.0;
820  for (i=l;i<=k;i++)
821  {
822  f=s*rv1[i];
823  if (Math.Abs(f)+anorm != anorm)
824  {
825  g=w[i];
826  h=PYTHAG(f,g);
827  w[i]=h;
828  h=1.0/h;
829  c=g*h;
830  s=(-f*h);
831  for (j=1;j<=m;j++)
832  {
833  y=A[j,nm];
834  z=A[j,i];
835  A[j,nm]=y*c+z*s;
836  A[j,i]=z*c-y*s;
837  }
838  }
839  }
840  }
841  z=w[k];
842  if (l == k)
843  {
844  if (z < 0.0)
845  {
846  w[k] = -z;
847  for (j=1;j<=n;j++) v[j,k]=(-v[j,k]);
848  }
849  break;
850  }
851  if (its == 30) Console.WriteLine("No convergence in 30 SVDCMP iterations");
852  x=w[l];
853  nm=k-1;
854  y=w[nm];
855  g=rv1[nm];
856  h=rv1[k];
857  f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
858  g=PYTHAG(f,1.0);
859  f=((x-z)*(x+z)+h*((y/(f+Sign(g,f)))-h))/x;
860  c=s=1.0;
861  for (j=l;j<=nm;j++)
862  {
863  i=j+1;
864  g=rv1[i];
865  y=w[i];
866  h=s*g;
867  g=c*g;
868  z=PYTHAG(f,h);
869  rv1[j]=z;
870  c=f/z;
871  s=h/z;
872  f=x*c+g*s;
873  g=g*c-x*s;
874  h=y*s;
875  y=y*c;
876  for (jj=1;jj<=n;jj++)
877  {
878  x=v[jj,j];
879  z=v[jj,i];
880  v[jj,j]=x*c+z*s;
881  v[jj,i]=z*c-x*s;
882  }
883  z=PYTHAG(f,h);
884  w[j]=z;
885  if (z!=0)
886  {
887  z=1.0/z;
888  c=f*z;
889  s=h*z;
890  }
891  f=(c*g)+(s*y);
892  x=(c*y)-(s*g);
893  for (jj=1;jj<=m;jj++)
894  {
895  y=A[jj,j];
896  z=A[jj,i];
897  A[jj,j]=y*c+z*s;
898  A[jj,i]=z*c-y*s;
899  }
900  }
901  rv1[l]=0.0;
902  rv1[k]=f;
903  w[k]=x;
904  }
905  }
906 
907  S_temp = new double[NP,NP];
908  V_temp = new double[NP, NP];
909  U_temp = new double[MP, NP];
910 
911  for (i = 1; i<= NP; i++) S_temp[i - 1, i - 1] = w[i];
912 
913  S_ = S_temp;
914 
915  for (i = 1; i<=NP;i++)
916  for (j = 1; j<=NP;j++) V_temp[i - 1, j - 1] = v[i, j];
917 
918  V_ = V_temp;
919 
920  for (i = 1; i<= MP;i++)
921  for (j = 1; j<=NP;j++) U_temp[i - 1, j - 1] = A[i, j];
922 
923  U_ = U_temp;
924  }
925 
926  private static double SQR(double a)
927  {
928  return a*a;
929  }
930 
931  private static double Sign (double a, double b)
932  {
933  if (b >= 0.0) {return Math.Abs(a);}
934  else {return -Math.Abs(a);}
935  }
936 
937  private static double PYTHAG(double a, double b)
938  {
939  double absa,absb;
940 
941  absa=Math.Abs(a);
942  absb=Math.Abs(b);
943  if (absa > absb) return absa*Math.Sqrt(1.0+SQR(absb/absa));
944  else return (absb == 0.0 ? 0.0 : absb*Math.Sqrt(1.0+SQR(absa/absb)));
945  }
946 
961  public static void SVD(Matrix Mat, out Matrix S, out Matrix U, out Matrix V)
962  {
963  double [,] s, u, v;
964  SVD(Mat.in_Mat, out s, out u, out v);
965  S = new Matrix(s);
966  U = new Matrix(u);
967  V = new Matrix(v);
968  }
969  #endregion
970 
971  #region "LU Decomposition of a matrix"
972  public static void LU(double[,]Mat , out double[,]L, out double[,] U, out double[,] P)
988  {
989  double[,] A;
990  int i,j,k, Rows, Cols;
991 
992  try
993  {
994  A = (double[,])Mat.Clone();
995  Find_R_C(Mat , out Rows, out Cols);
996  }
997  catch{throw new MatrixNullException();}
998 
999  if (Rows != Cols) throw new MatrixNotSquare();
1000 
1001  int IMAX = 0, N = Rows;
1002  double AAMAX, Sum, Dum, TINY = 1E-20;
1003 
1004  int[] INDX = new int[N+1];
1005  double[] VV = new double[N*10];
1006  double D = 1.0;
1007 
1008  for (i = 0; i<= N;i++)
1009  {
1010  AAMAX = 0.0;
1011  for (j = 0; j<= N;j++)
1012  if (Math.Abs(A[i, j]) > AAMAX) AAMAX = Math.Abs(A[i, j]);
1013  if (AAMAX == 0.0) throw new MatrixSingularException();
1014  VV[i] = 1.0 / AAMAX;
1015  }
1016  for (j = 0; j<= N; j++)
1017  {
1018  if (j > 0)
1019  {
1020  for (i = 0; i<= (j - 1); i++)
1021  {
1022  Sum = A[i, j];
1023  if (i > 0)
1024  {
1025  for (k = 0; k<= (i - 1); k++)
1026  Sum = Sum - A[i, k] * A[k, j];
1027  A[i, j] = Sum;
1028  }
1029  }
1030  }
1031  AAMAX = 0.0;
1032  for (i = j; i<= N; i++)
1033  {
1034  Sum = A[i, j];
1035  if (j > 0)
1036  {
1037  for (k = 0; k<= (j - 1); k++)
1038  Sum = Sum - A[i, k] * A[k, j];
1039  A[i, j] = Sum;
1040  }
1041  Dum = VV[i] * Math.Abs(Sum);
1042  if (Dum >= AAMAX)
1043  {
1044  IMAX = i;
1045  AAMAX = Dum;
1046  }
1047  }
1048  if (j != IMAX)
1049  {
1050  for (k = 0; k<= N; k++)
1051  {
1052  Dum = A[IMAX, k];
1053  A[IMAX, k] = A[j, k];
1054  A[j, k] = Dum;
1055  }
1056  D = -D;
1057  VV[IMAX] = VV[j];
1058  }
1059  INDX[j] = IMAX;
1060  if (j != N)
1061  {
1062  if (A[j, j] == 0.0) A[j, j] = TINY;
1063  Dum = 1.0 / A[j, j];
1064  for (i = j + 1; i<=N; i++)
1065  A[i, j] = A[i, j] * Dum;
1066 
1067  }
1068  }
1069 
1070  if (A[N, N] == 0.0) A[N, N] = TINY;
1071 
1072  int count=0;
1073  double[,] l = new double[N+1,N+1];
1074  double[,] u = new double[N+1,N+1];
1075 
1076  for (i = 0; i<=N; i++)
1077  {
1078  for (j = 0; j<=count; j++)
1079  {
1080  if (i!=0) l[i,j] = A[i,j];
1081  if (i==j) l[i,j] = 1.0;
1082  u[N-i,N-j] = A[N-i,N-j];
1083  }
1084  count++;
1085  }
1086 
1087  L = l;
1088  U = u;
1089 
1090  P = Identity(N+1);
1091  for (i=0;i<=N;i++)
1092  {
1093  SwapRows(P,i,INDX[i]);
1094  }
1095  }
1096 
1097  private static void SwapRows(double[,] Mat, int Row, int toRow)
1098  {
1099  int N = Mat.GetUpperBound(0);
1100  double[,] dumArray = new double[1,N+1];
1101  for (int i = 0; i <= N; i++)
1102  {
1103  dumArray[0, i] = Mat[Row, i];
1104  Mat[Row, i] = Mat[toRow, i];
1105  Mat[toRow, i] = dumArray[0, i];
1106  }
1107  }
1108 
1124  public static void LU(Matrix Mat , out Matrix L, out Matrix U, out Matrix P)
1125  {
1126  double [,] l, u, p;
1127  LU(Mat.in_Mat, out l, out u, out p);
1128  L = new Matrix(l);
1129  U = new Matrix(u);
1130  P = new Matrix(p);
1131  }
1132  #endregion
1133 
1134  #region "Solve system of linear equations"
1135  public static double[,] SolveLinear(double[,]MatA , double[,] MatB)
1150  {
1151  double[,] A;
1152  double[,] B;
1153  double SUM;
1154  int i,ii,j,k,ll, Rows, Cols;
1155 
1156  #region "LU Decompose"
1157  try
1158  {
1159  A = (double[,])MatA.Clone();
1160  B = (double[,])MatB.Clone();
1161  Find_R_C(A , out Rows, out Cols);
1162  }
1163  catch{throw new MatrixNullException();}
1164 
1165  if (Rows != Cols) throw new MatrixNotSquare();
1166  if ((B.GetUpperBound(0)!=Rows) || (B.GetUpperBound(1)!=0))
1167  throw new MatrixDimensionException();
1168 
1169  int IMAX = 0, N = Rows;
1170  double AAMAX, Sum, Dum, TINY = 1E-20;
1171 
1172  int[] INDX = new int[N+1];
1173  double[] VV = new double[N*10];
1174  double D = 1.0;
1175 
1176  for (i = 0; i<= N;i++)
1177  {
1178  AAMAX = 0.0;
1179  for (j = 0; j<= N;j++)
1180  if (Math.Abs(A[i, j]) > AAMAX) AAMAX = Math.Abs(A[i, j]);
1181  if (AAMAX == 0.0) throw new MatrixSingularException();
1182  VV[i] = 1.0 / AAMAX;
1183  }
1184  for (j = 0; j<= N; j++)
1185  {
1186  if (j > 0)
1187  {
1188  for (i = 0; i<= (j - 1); i++)
1189  {
1190  Sum = A[i, j];
1191  if (i > 0)
1192  {
1193  for (k = 0; k<= (i - 1); k++)
1194  Sum = Sum - A[i, k] * A[k, j];
1195  A[i, j] = Sum;
1196  }
1197  }
1198  }
1199  AAMAX = 0.0;
1200  for (i = j; i<= N; i++)
1201  {
1202  Sum = A[i, j];
1203  if (j > 0)
1204  {
1205  for (k = 0; k<= (j - 1); k++)
1206  Sum = Sum - A[i, k] * A[k, j];
1207  A[i, j] = Sum;
1208  }
1209  Dum = VV[i] * Math.Abs(Sum);
1210  if (Dum >= AAMAX)
1211  {
1212  IMAX = i;
1213  AAMAX = Dum;
1214  }
1215  }
1216  if (j != IMAX)
1217  {
1218  for (k = 0; k<= N; k++)
1219  {
1220  Dum = A[IMAX, k];
1221  A[IMAX, k] = A[j, k];
1222  A[j, k] = Dum;
1223  }
1224  D = -D;
1225  VV[IMAX] = VV[j];
1226  }
1227  INDX[j] = IMAX;
1228  if (j != N)
1229  {
1230  if (A[j, j] == 0.0) A[j, j] = TINY;
1231  Dum = 1.0 / A[j, j];
1232  for (i = j + 1; i<=N; i++)
1233  A[i, j] = A[i, j] * Dum;
1234 
1235  }
1236  }
1237  if (A[N, N] == 0.0) A[N, N] = TINY;
1238  #endregion
1239 
1240  ii=-1;
1241  for ( i=0; i<=N; i++)
1242  {
1243  ll=INDX[i];
1244  SUM=B[ll,0];
1245  B[ll,0]=B[i,0];
1246  if (ii!=-1)
1247  {
1248  for (j=ii; j<=i-1; j++) SUM=SUM-A[i,j]*B[j,0];
1249  }
1250  else if (SUM!=0) ii=i;
1251  B[i,0]=SUM;
1252  }
1253  for (i=N;i>=0;i--)
1254  {
1255  SUM=B[i,0];
1256  if (i < N)
1257  {
1258  for (j=i+1; j<=N; j++) SUM=SUM-A[i,j]*B[j,0];
1259  }
1260  B[i,0]=SUM/A[i,i];
1261  }
1262  return B;
1263  }
1264 
1279  public static Matrix SolveLinear(Matrix MatA , Matrix MatB)
1280  {return new Matrix(Matrix.SolveLinear(MatA.in_Mat, MatB.in_Mat));}
1281  #endregion
1282 
1283  #region "Rank of a matrix"
1284  public static int Rank(double[,] Mat)
1291  {
1292  int r = 0;
1293  double[,] S, U, V;
1294  try
1295  {
1296  int Rows,Cols;
1297  Find_R_C(Mat , out Rows, out Cols);
1298  }
1299  catch{throw new MatrixNullException();}
1300  double EPS = 2.2204E-16;
1301  SVD(Mat, out S, out U, out V);
1302 
1303  for (int i = 0; i<= S.GetUpperBound(0);i++)
1304  {if (Math.Abs(S[i, i]) > EPS) r++;}
1305 
1306  return r;
1307  }
1308 
1315  public static int Rank(Matrix Mat)
1316  {return Rank(Mat.in_Mat);}
1317  #endregion
1318 
1319  #region "Pseudoinverse of a matrix"
1320  public static double[,] PINV(double[,] Mat)
1329  {
1330  double[,] Matrix;
1331  int i, m, n, j, l;
1332  double[,] S, Part_I, Part_II;
1333  double EPS, MulAdd, Tol, Largest_Item=0;
1334 
1335  double[,] svdU, svdS, svdV;
1336 
1337  try
1338  {
1339  Matrix = (double[,])Mat.Clone();
1340  Find_R_C(Mat , out m, out n);
1341  }
1342  catch{throw new MatrixNullException();}
1343 
1344  SVD(Mat, out svdS, out svdU, out svdV);
1345 
1346  EPS = 2.2204E-16;
1347  m ++;
1348  n ++;
1349 
1350  Part_I = new double[m,n];
1351  Part_II = new double[m,n];
1352  S = new Double[n,n];
1353 
1354  MulAdd = 0;
1355  for (i = 0 ; i<= svdS.GetUpperBound(0);i++)
1356  {
1357  if (i == 0) Largest_Item = svdS[0, 0];
1358  if (Largest_Item < svdS[i, i]) Largest_Item = svdS[i, i];
1359  }
1360 
1361  if (m > n) Tol = m * Largest_Item * EPS;
1362  else Tol = n * Largest_Item * EPS;
1363 
1364  for (i = 0; i < n; i++) S[i, i] = svdS[i,i];
1365 
1366  for (i = 0;i < m; i++)
1367  {
1368  for (j = 0; j < n; j++)
1369  {
1370  for (l = 0; l< n; l++)
1371  {
1372  if (S[l, j] > Tol) MulAdd += svdU[i, l] * (1 / S[l, j]);
1373  }
1374  Part_I[i, j] = MulAdd;
1375  MulAdd = 0;
1376  }
1377  }
1378 
1379  for (i = 0;i < m; i++)
1380  {
1381  for (j = 0; j < n; j++)
1382  {
1383  for (l = 0; l< n; l++)
1384  {
1385  MulAdd += Part_I[i, l] * svdV[j, l];
1386  }
1387  Part_II[i, j] = MulAdd;
1388  MulAdd = 0;
1389  }
1390  }
1391  return Transpose(Part_II);
1392  }
1393 
1402  public static Matrix PINV(Matrix Mat)
1403  {return new Matrix(PINV(Mat.in_Mat));}
1404  #endregion
1405 
1406  #region "Eigen Values and Vactors of Symmetric Matrix"
1407  public static void Eigen(double[,] Mat, out double[,] d,out double[,] v)
1422  {
1423 
1424  double [,] a;
1425  int Rows, Cols;
1426  try
1427  {
1428  Find_R_C(Mat, out Rows, out Cols);
1429  a = (double[,])Mat.Clone();
1430  }
1431  catch{throw new MatrixNullException();}
1432 
1433  if (Rows != Cols) throw new MatrixNotSquare();
1434 
1435  int j,iq,ip,i,n, nrot;
1436  double tresh,theta,tau,t,sm,s,h,g,c;
1437  double[] b,z;
1438 
1439  n = Rows;
1440  d = new double[n+1,1];
1441  v = new double[n+1,n+1];
1442 
1443  b=new double[n+1];
1444  z=new double[n+1];
1445  for (ip=0;ip <= n;ip++)
1446  {
1447  for (iq=0;iq <= n;iq++) v[ip,iq]=0.0;
1448  v[ip,ip]=1.0;
1449  }
1450  for (ip=0;ip <= n;ip++)
1451  {
1452  b[ip]=d[ip,0]=a[ip,ip];
1453  z[ip]=0.0;
1454  }
1455 
1456  nrot=0;
1457  for (i=0;i <= 50;i++)
1458  {
1459  sm=0.0;
1460  for (ip=0;ip <= n-1;ip++)
1461  {
1462  for (iq=ip+1;iq <= n;iq++)
1463  sm += Math.Abs(a[ip,iq]);
1464  }
1465  if (sm == 0.0)
1466  {
1467  return;
1468  }
1469  if (i < 4)
1470  tresh=0.2*sm/(n*n);
1471  else
1472  tresh=0.0;
1473  for (ip=0;ip <= n-1;ip++)
1474  {
1475  for (iq=ip+1;iq <= n;iq++)
1476  {
1477  g=100.0*Math.Abs(a[ip,iq]);
1478  if (i > 4 && (double)(Math.Abs(d[ip,0])+g) == (double)Math.Abs(d[ip,0])
1479  && (double)(Math.Abs(d[iq,0])+g) == (double)Math.Abs(d[iq,0]))
1480  a[ip,iq]=0.0;
1481  else if (Math.Abs(a[ip,iq]) > tresh)
1482  {
1483  h=d[iq,0]-d[ip,0];
1484  if ((double)(Math.Abs(h)+g) == (double)Math.Abs(h))
1485  t=(a[ip,iq])/h;
1486  else
1487  {
1488  theta=0.5*h/(a[ip,iq]);
1489  t=1.0/(Math.Abs(theta)+Math.Sqrt(1.0+theta*theta));
1490  if (theta < 0.0) t = -t;
1491  }
1492  c=1.0/Math.Sqrt(1+t*t);
1493  s=t*c;
1494  tau=s/(1.0+c);
1495  h=t*a[ip,iq];
1496  z[ip] -= h;
1497  z[iq] += h;
1498  d[ip,0] -= h;
1499  d[iq,0] += h;
1500  a[ip,iq]=0.0;
1501  for (j=0;j <= ip-1;j++)
1502  {
1503  ROT(g,h,s,tau,a,j,ip,j,iq);
1504  }
1505  for (j=ip+1;j <= iq-1;j++)
1506  {
1507  ROT(g,h,s,tau,a,ip,j,j,iq);
1508  }
1509  for (j=iq+1;j <= n;j++)
1510  {
1511  ROT(g,h,s,tau,a,ip,j,iq,j);
1512  }
1513  for (j=0;j <= n;j++)
1514  {
1515  ROT(g,h,s,tau,v,j,ip,j,iq);
1516  }
1517  ++(nrot);
1518  }
1519  }
1520  }
1521  for (ip=0;ip <= n;ip++)
1522  {
1523  b[ip] += z[ip];
1524  d[ip,0]=b[ip];
1525  z[ip]=0.0;
1526  }
1527  }
1528  Console.WriteLine("Too many iterations in routine jacobi");
1529  }
1530 
1531  private static void ROT(double g ,double h ,double s ,double tau,
1532  double[,] a ,int i ,int j ,int k ,int l)
1533  {
1534  g=a[i,j]; h=a[k,l];
1535  a[i,j]=g-s*(h+g*tau);
1536  a[k,l]=h+s*(g-h*tau);
1537  }
1538 
1553  public static void Eigen(Matrix Mat, out Matrix d,out Matrix v)
1554  {
1555  double [,] D, V;
1556  Eigen(Mat.in_Mat, out D, out V);
1557  d = new Matrix(D);
1558  v = new Matrix(V);
1559  }
1560  #endregion
1561 
1562  #region "Multiply a matrix or a vector with a scalar quantity"
1563  public static double[,] ScalarMultiply(double Value, double[,] Mat)
1572  {
1573  int i, j, Rows, Cols;
1574  double[,] sol;
1575 
1576  try {Find_R_C(Mat, out Rows, out Cols);}
1577  catch{throw new MatrixNullException();}
1578  sol = new double[Rows+1, Cols+1];
1579 
1580  for (i = 0; i<=Rows;i++)
1581  for (j = 0; j<=Cols;j++)
1582  sol[i, j] = Mat[i, j] * Value;
1583 
1584  return sol;
1585  }
1586 
1595  public static Matrix ScalarMultiply(double Value, Matrix Mat)
1596  {return new Matrix(ScalarMultiply(Value,Mat.in_Mat));}
1597 
1609  public static Matrix operator*(Matrix Mat, double Value)
1610  {return new Matrix(ScalarMultiply(Value,Mat.in_Mat));}
1611 
1623  public static Matrix operator*(double Value,Matrix Mat)
1624  {return new Matrix(ScalarMultiply(Value,Mat.in_Mat));}
1625  #endregion
1626 
1627  #region "Divide a matrix or a vector with a scalar quantity"
1628  public static double[,] ScalarDivide(double Value, double[,] Mat)
1637  {
1638  int i, j, Rows, Cols;
1639  double[,] sol;
1640 
1641  try {Find_R_C(Mat, out Rows, out Cols);}
1642  catch{throw new MatrixNullException();}
1643 
1644  sol = new double[Rows+1, Cols+1];
1645 
1646  for (i = 0; i<=Rows;i++)
1647  for (j = 0; j<=Cols;j++)
1648  sol[i, j] = Mat[i, j] / Value;
1649 
1650  return sol;
1651  }
1652 
1661  public static Matrix ScalarDivide(double Value, Matrix Mat)
1662  {return new Matrix(ScalarDivide(Value,Mat.in_Mat));}
1663 
1674  public static Matrix operator/(Matrix Mat,double Value)
1675  {return new Matrix(ScalarDivide(Value,Mat.in_Mat));}
1676  #endregion
1677 
1678  #region "Vectors Cross Product"
1679  public static double[] CrossProduct(double[] V1, double[] V2)
1688  {
1689  double i, j, k;
1690  double[] sol = new double[2];
1691  int Rows1;
1692  int Rows2;
1693 
1694  try
1695  {
1696  Find_R_C(V1, out Rows1);
1697  Find_R_C(V2, out Rows2);
1698  }
1699  catch{throw new MatrixNullException();}
1700 
1701  if (Rows1 != 2) throw new VectorDimensionException();
1702 
1703  if (Rows2 != 2) throw new VectorDimensionException();
1704 
1705  i = V1[1] * V2[2] - V1[2] * V2[1];
1706  j = V1[2] * V2[0] - V1[0] * V2[2];
1707  k = V1[0] * V2[1] - V1[1] * V2[0];
1708 
1709  sol[0] = i ; sol[1] = j ; sol[2] = k;
1710 
1711  return sol;
1712  }
1713 
1722  public static double[,] CrossProduct(double[,] V1, double[,] V2)
1723  {
1724  double i, j, k;
1725  double[,] sol = new double[3,1];
1726  int Rows1, Cols1;
1727  int Rows2, Cols2;
1728 
1729  try
1730  {
1731  Find_R_C(V1, out Rows1, out Cols1);
1732  Find_R_C(V2, out Rows2, out Cols2);
1733  }
1734  catch{throw new MatrixNullException();}
1735 
1736  if (Rows1 != 2 || Cols1 != 0) throw new VectorDimensionException();
1737 
1738  if (Rows2 != 2 || Cols2 != 0) throw new VectorDimensionException();
1739 
1740  i = V1[1, 0] * V2[2, 0] - V1[2, 0] * V2[1, 0];
1741  j = V1[2, 0] * V2[0, 0] - V1[0, 0] * V2[2, 0];
1742  k = V1[0, 0] * V2[1, 0] - V1[1, 0] * V2[0, 0];
1743 
1744  sol[0, 0] = i ; sol[1, 0] = j ; sol[2, 0] = k;
1745 
1746  return sol;
1747  }
1748 
1757  public static Matrix CrossProduct(Matrix V1, Matrix V2)
1758  {return (new Matrix((CrossProduct(V1.in_Mat,V2.in_Mat))));}
1759  #endregion
1760 
1761  #region "Vectors Dot Product"
1762  public static double DotProduct(double[] V1, double[] V2)
1771  {
1772  int Rows1;
1773  int Rows2;
1774 
1775  try
1776  {
1777  Find_R_C(V1, out Rows1);
1778  Find_R_C(V2, out Rows2);
1779  }
1780  catch{throw new MatrixNullException();}
1781 
1782  if (Rows1 != 2) throw new VectorDimensionException();
1783 
1784  if (Rows2 != 2) throw new VectorDimensionException();
1785 
1786  return (V1[0] * V2[0] + V1[1] * V2[1] + V1[2] * V2[2]);
1787  }
1788 
1797  public static double DotProduct(double[,] V1, double[,] V2)
1798  {
1799  int Rows1, Cols1;
1800  int Rows2, Cols2;
1801 
1802  try
1803  {
1804  Find_R_C(V1, out Rows1, out Cols1);
1805  Find_R_C(V2, out Rows2, out Cols2);
1806  }
1807  catch{throw new MatrixNullException();}
1808 
1809  if (Rows1 != 2 || Cols1 != 0) throw new VectorDimensionException();
1810 
1811  if (Rows2 != 2 || Cols2 != 0) throw new VectorDimensionException();
1812 
1813  return (V1[0,0] * V2[0,0] + V1[1,0] * V2[1,0] + V1[2,0] * V2[2,0]);
1814  }
1815 
1824  public static double DotProduct(Matrix V1, Matrix V2)
1825  {return (DotProduct(V1.in_Mat, V2.in_Mat));}
1826  #endregion
1827 
1828  #region "Magnitude of a Vector"
1829  public static double VectorMagnitude(double[] V)
1836  {
1837  int Rows;
1838 
1839  try {Find_R_C(V, out Rows);}
1840  catch{throw new MatrixNullException();}
1841 
1842  if (Rows != 2) throw new VectorDimensionException();
1843 
1844  return Math.Sqrt(V[0] * V[0] + V[1] * V[1] + V[2] * V[2]);
1845  }
1846 
1853  public static double VectorMagnitude(double[,] V)
1854  {
1855  int Rows, Cols;
1856 
1857  try {Find_R_C(V, out Rows, out Cols);}
1858  catch{throw new MatrixNullException();}
1859 
1860  if (Rows != 2 || Cols != 0) throw new VectorDimensionException();
1861 
1862  return Math.Sqrt(V[0, 0] * V[0, 0] + V[1, 0] * V[1, 0] + V[2, 0] * V[2, 0]);
1863  }
1864 
1871  public static double VectorMagnitude(Matrix V)
1872  {return (VectorMagnitude(V.in_Mat));}
1873  #endregion
1874 
1875  #region "Two Matrices are equal or not"
1876  public static bool IsEqual (double[,] Mat1, double[,] Mat2)
1884  {
1885  double eps = 1E-14;
1886  int Rows1, Cols1;
1887  int Rows2, Cols2;
1888 
1889  //Find The dimensions !!
1890  try
1891  {
1892  Find_R_C(Mat1,out Rows1, out Cols1);
1893  Find_R_C(Mat2,out Rows2, out Cols2);
1894  }
1895  catch
1896  {
1897  throw new MatrixNullException();
1898  }
1899  if (Rows1 != Rows2 || Cols1 != Cols2) throw new MatrixDimensionException();
1900 
1901  for (int i=0; i<=Rows1; i++)
1902  {
1903  for (int j=0; j<=Cols1; j++)
1904  {
1905  if (Math.Abs(Mat1[i,j]-Mat2[i,j])>eps) return false;
1906  }
1907  }
1908  return true;
1909  }
1910 
1918  public static bool IsEqual (Matrix Mat1, Matrix Mat2)
1919  {return IsEqual(Mat1.in_Mat, Mat2.in_Mat);}
1920 
1928  public static bool operator ==(Matrix Mat1, Matrix Mat2)
1929  {return IsEqual(Mat1.in_Mat, Mat2.in_Mat);}
1930 
1938  public static bool operator !=(Matrix Mat1, Matrix Mat2)
1939  {return (! IsEqual(Mat1.in_Mat, Mat2.in_Mat));}
1940 
1947  public override bool Equals(Object obj)
1948  {
1949  try {return (bool) (this == (Matrix) obj);}
1950  catch {return false;}
1951  }
1952  #endregion
1953 
1954  #region "Print Matrix"
1955  public static string PrintMat(double[,] Mat)
1963  {
1964  int N_Rows, N_Columns, k, i, j, m;
1965  string StrElem;
1966  int StrLen; int[] Greatest;
1967  string LarString, OptiString, sol;
1968 
1969  //Find The dimensions !!
1970  try {Find_R_C(Mat, out N_Rows, out N_Columns);}
1971  catch
1972  {
1973  throw new MatrixNullException();
1974  }
1975 
1976  sol = "";
1977  LarString = "";
1978  OptiString = "";
1979  Greatest = new int[N_Columns+1];
1980 
1981  for (i = 0;i<= N_Rows;i++)
1982  {
1983  for (j = 0;j<= N_Columns;j++)
1984  {
1985  if (i == 0)
1986  {
1987  Greatest[j] = 0;
1988  for (m = 0;m<= N_Rows;m++)
1989  {
1990  StrElem =Mat[m, j].ToString("0.0000");
1991  StrLen = StrElem.Length;
1992  if (Greatest[j] < StrLen)
1993  {
1994  Greatest[j] = StrLen;
1995  LarString = StrElem;
1996  }
1997  }
1998  if (LarString.StartsWith("-")) Greatest[j] = Greatest[j];
1999  }
2000  StrElem = Mat[i, j].ToString("0.0000");
2001  if (StrElem.StartsWith( "-"))
2002  {
2003  StrLen = StrElem.Length;
2004  if (Greatest[j] >= StrLen)
2005  {
2006  for (k = 1;k<=(Greatest[j] - StrLen);k++) OptiString += " ";
2007  OptiString += " ";
2008  }
2009  }
2010  else
2011  {
2012  StrLen = StrElem.Length;
2013  if (Greatest[j] > StrLen)
2014  {
2015  for (k = 1;k<= (Greatest[j] - StrLen);k++)
2016  {
2017  OptiString += " ";
2018  }
2019  }
2020  }
2021  OptiString += " " + (Mat[i, j].ToString("0.0000"));
2022  }
2023 
2024  if (i != N_Rows)
2025  {
2026  sol += OptiString + "\n";
2027  OptiString = "";
2028  }
2029  sol += OptiString;
2030  OptiString = "";
2031  }
2032  return sol;
2033  }
2034 
2042  public static string PrintMat(Matrix Mat)
2043  {return (PrintMat(Mat.in_Mat));}
2044 
2051  public override string ToString()
2052  {return (PrintMat(this.in_Mat));}
2053 
2054  #endregion
2055  }
2056 }
Generated on Wed Jun 12 2013 23:20:54 for VIS2 SS2013 CVD DVR by   doxygen 1.8.4