矩阵基本操作的实现(C# 源代码）

 1 using  System;
  2 using  System.IO;
  3 using  System.Diagnostics;
  4
  5
  6 namespace  Adjust
  7 {
  8      /**/ ///   <summary>
  9      ///  Matrix 的摘要说明。
 10      ///  实现矩阵的基本运算
 11      ///   </summary>
 12      public   class  Matrix
 13      {
 14     
 15          // 构造方阵
 16          public   Matrix( int  row)
 17          {
 18             m_data  =   new   double [row,row];
 19
 20         }
 21          public  Matrix( int  row, int  col)
 22          {
 23             m_data  =   new   double [row,col];
 24         }
 25          // 复制构造函数
 26          public  Matrix(Matrix m)
 27          {
 28              int  row  =  m.Row;
 29              int  col  =  m.Col;
 30             m_data  =   new   double [row,col];
 31
 32              for ( int  i = 0 ;i < row;i ++ )
 33                  for ( int  j = 0 ;j < col;j ++ )
 34                     m_data[i,j]  =  m[i,j];
 35
 36         }
 37
 38          /**/ /*
 39         //分配方阵的大小
 40         //对于已含有内存的矩阵，将清空数据
 41         public void SetSize(int row)
 42         {
 43             m_data = new double[row,row];
 44         }
 45
 46         
 47         //分配矩阵的大小
 48         //对于已含有内存的矩阵，将清空数据
 49         public void SetSize(int row,int col)
 50         {
 51             m_data = new double[row,col];
 52         }
 53          */
 54
 55          // unit matrix:设为单位阵
 56          public   void  SetUnit()
 57          {
 58              for ( int  i = 0 ;i < m_data.GetLength( 0 );i ++ )
 59                  for ( int  j = 0 ;j < m_data.GetLength( 1 );j ++ )
 60                     m_data[i,j]  =  ((i == j) ? 1 : 0 );
 61         }
 62
 63          // 设置元素值
 64          public   void  SetValue( double  d)
 65          {
 66              for ( int  i = 0 ;i < m_data.GetLength( 0 );i ++ )
 67                  for ( int  j = 0 ;j < m_data.GetLength( 1 );j ++ )
 68                     m_data[i,j]  =  d;
 69         }
 70
 71          //  Value extraction：返中行数
 72          public   int  Row
 73          {
 74              get
 75              {
 76
 77                  return  m_data.GetLength( 0 );
 78             }
 79         }
 80
 81          // 返回列数
 82          public   int  Col
 83          {
 84              get
 85              {
 86                  return  m_data.GetLength( 1 );
 87             }
 88         }
 89
 90          // 重载索引
 91          // 存取数据成员
 92          public   double   this [ int  row, int  col]
 93          {
 94              get
 95              {
 96                  return  m_data[row,col];
 97             }
 98              set
 99              {
100                 m_data[row,col]  =  value;
101             }
102         }
103
104          // primary change
105          // 　初等变换　对调两行：ri<-->rj
106          public  Matrix  Exchange( int  i, int  j)
107          {
108              double  temp;
109
110              for ( int  k = 0 ;k < Col;k ++ )
111              {
112                 temp  =  m_data[i,k];
113                 m_data[i,k]  =  m_data[j,k];
114                 m_data[j,k]  =  temp;
115             }
116              return   this ;
117         }
118
119
120          // 初等变换　第index 行乘以mul
121         Matrix Multiple( int  index, double  mul)    
122          {
123              for ( int  j = 0 ;j < Col;j ++ )
124              {
125                 m_data[index,j]  *=  mul;
126             }
127              return   this ;
128         }
129                 
130
131          // 初等变换 第src行乘以mul加到第index行
132         Matrix MultipleAdd( int  index, int  src, double  mul)
133          {
134              for ( int  j = 0 ;j < Col;j ++ )
135              {
136                 m_data[index,j]  +=  m_data[src,j] * mul;
137             }
138
139              return   this ;
140         }
141
142          // transpose 转置
143          public  Matrix Transpose()
144          {
145             Matrix ret  =   new  Matrix(Col,Row);
146
147              for ( int  i = 0 ;i < Row;i ++ )
148                  for ( int  j = 0 ;j < Col;j ++ )
149                  {
150                     ret[j,i]  =  m_data[i,j];
151                 }
152              return  ret;
153         }
154         
155          // binary addition 矩阵加
156          public   static  Matrix  operator +  (Matrix lhs,Matrix rhs)
157          {
158              if (lhs.Row  !=  rhs.Row)     // 异常
159              {
160                 System.Exception e  =   new  Exception( " 相加的两个矩阵的行数不等 " );
161                  throw  e;
162             }
163              if (lhs.Col  !=  rhs.Col)      // 异常
164              {
165                 System.Exception e  =   new  Exception( " 相加的两个矩阵的列数不等 " );
166                  throw  e;
167             }
168
169              int  row  =  lhs.Row;
170              int  col  =  lhs.Col;
171             Matrix ret = new  Matrix(row,col);
172
173              for ( int  i = 0 ;i < row;i ++ )
174                  for ( int  j = 0 ;j < col;j ++ )
175                  {
176                      double  d  =  lhs[i,j]  +  rhs[i,j];
177                     ret[i,j]  =  d;
178                 }
179              return  ret;
180
181         }
182
183          // binary subtraction 矩阵减
184          public   static  Matrix  operator -  (Matrix lhs,Matrix rhs)
185          {
186              if (lhs.Row  !=  rhs.Row)     // 异常
187              {
188                 System.Exception e  =   new  Exception( " 相减的两个矩阵的行数不等 " );
189                  throw  e;
190             }
191              if (lhs.Col  !=  rhs.Col)      // 异常
192              {
193                 System.Exception e  =   new  Exception( " 相减的两个矩阵的列数不等 " );
194                  throw  e;
195             }
196
197              int  row  =  lhs.Row;
198              int  col  =  lhs.Col;
199             Matrix ret = new  Matrix(row,col);
200
201              for ( int  i = 0 ;i < row;i ++ )
202                  for ( int  j = 0 ;j < col;j ++ )
203                  {
204                      double  d  =  lhs[i,j]  -  rhs[i,j];
205                     ret[i,j]  =  d;
206                 }
207              return  ret;
208         }
209
210
211          // binary multiple 矩阵乘
212          public   static  Matrix  operator *  (Matrix lhs,Matrix rhs)
213          {
214              if (lhs.Col  !=  rhs.Row)     // 异常
215              {
216                 System.Exception e  =   new  Exception( " 相乘的两个矩阵的行列数不匹配 " );
217                  throw  e;
218             }
219
220             Matrix ret  =   new  Matrix (lhs.Row,rhs.Col);
221              double  temp;
222              for ( int  i = 0 ;i < lhs.Row;i ++ )
223              {
224                  for ( int  j = 0 ;j < rhs.Col;j ++ )
225                  {
226                     temp  =   0 ;
227                      for ( int  k = 0 ;k < lhs.Col;k ++ )
228                      {
229                         temp  +=  lhs[i,k]  *  rhs[k,j];
230                     }
231                     ret[i,j]  =  temp;
232                 }
233             }
234
235              return  ret;
236         }
237
238
239          // binary division 矩阵除
240          public   static  Matrix  operator /  (Matrix lhs,Matrix rhs)
241          {
242              return  lhs  *  rhs.Inverse();
243         }
244
245          // unary addition单目加
246          public   static  Matrix  operator +  (Matrix m)
247          {
248             Matrix ret  =   new  Matrix(m);
249              return  ret;
250         }
251
252          // unary subtraction 单目减
253          public   static  Matrix  operator -  (Matrix m)
254          {
255             Matrix ret  =   new  Matrix(m);
256              for ( int  i = 0 ;i < ret.Row;i ++ )
257                  for ( int  j =   0 ;j < ret.Col;j ++ )
258                  {
259                     ret[i,j]  =   - ret[i,j];
260                 }
261
262              return  ret;
263         }
264
265          // number multiple 数乘
266          public   static  Matrix  operator *  ( double  d,Matrix m)
267          {
268             Matrix ret  =   new  Matrix(m);
269              for ( int  i = 0 ;i < ret.Row;i ++ )
270                  for ( int  j = 0 ;j < ret.Col;j ++ )
271                     ret[i,j]  *=  d;
272
273              return  ret;
274         }
275
276          // number division 数除
277          public   static  Matrix  operator /  ( double  d,Matrix m)
278          {
279              return  d * m.Inverse();
280         }
281
282          // 功能：返回列主元素的行号
283          // 参数：row为开始查找的行号
284          // 说明：在行号[row,Col)范围内查找第row列中绝对值最大的元素，返回所在行号
285          int  Pivot( int  row)
286          {
287              int  index = row;
288
289              for ( int  i = row + 1 ;i < Row;i ++ )
290              {
291                  if (m_data[i,row]  >  m_data[index,row])
292                     index = i;
293             }
294
295              return  index;
296         }
297
298          // inversion 逆阵：使用矩阵的初等变换，列主元素消去法
299          public  Matrix Inverse()
300          {
301              if (Row  !=  Col)     // 异常,非方阵
302              {
303                 System.Exception e  =   new  Exception( " 求逆的矩阵不是方阵 " );
304                  throw  e;
305             }
306 StreamWriter sw  =   new  StreamWriter( " ..//annex//close_matrix.txt " );
307             Matrix tmp  =   new  Matrix( this );
308             Matrix ret  = new  Matrix(Row);     // 单位阵
309             ret.SetUnit();
310
311              int  maxIndex;
312              double  dMul;
313
314              for ( int  i = 0 ;i < Row;i ++ )
315              {
316                 maxIndex  =  tmp.Pivot(i);
317     
318                  if (tmp.m_data[maxIndex,i] == 0 )
319                  {
320                     System.Exception e  =   new  Exception( " 求逆的矩阵的行列式的值等于0, " );
321                      throw  e;
322                 }
323
324                  if (maxIndex  !=  i)     // 下三角阵中此列的最大值不在当前行，交换
325                  {
326                     tmp.Exchange(i,maxIndex);
327                     ret.Exchange(i,maxIndex);
328
329                 }
330
331                 ret.Multiple(i, 1 / tmp[i,i]);
332
333                 tmp.Multiple(i, 1 / tmp[i,i]);
334
335                  for ( int  j = i + 1 ;j < Row;j ++ )
336                  {
337                     dMul  =   - tmp[j,i] / tmp[i,i];
338                     tmp.MultipleAdd(j,i,dMul);
339                     ret.MultipleAdd(j,i,dMul);
340     
341                 }
342 sw.WriteLine( " tmp=/r/n " + tmp);
343 sw.WriteLine( " ret=/r/n " + ret);
344             } // end for
345
346
347 sw.WriteLine( " **=/r/n " +   this * ret);
348
349              for ( int  i = Row - 1 ;i > 0 ;i -- )
350              {
351                  for ( int  j = i - 1 ;j >= 0 ;j -- )
352                  {
353                     dMul  =   - tmp[j,i] / tmp[i,i];
354                     tmp.MultipleAdd(j,i,dMul);
355                     ret.MultipleAdd(j,i,dMul);
356                 }
357             } // end for
358
359
360 sw.WriteLine( " tmp=/r/n " + tmp);
361 sw.WriteLine( " ret=/r/n " + ret);
362 sw.WriteLine( " ***=/r/n " +   this * ret);
363 sw.Close();
364         
365              return  ret;
366
367         } // end Inverse
368
369          #region
370 /**/ /*
371         //inversion 逆阵：使用矩阵的初等变换，列主元素消去法
372         public Matrix Inverse()
373         {
374             if(Row != Col)    //异常,非方阵
375             {
376                 System.Exception e = new Exception("求逆的矩阵不是方阵");
377                 throw e;
378             }
379             ///
380             StreamWriter sw = new StreamWriter("..//annex//matrix_mul.txt");
381             
382             ///    
383             Matrix tmp = new Matrix(this);
384             Matrix ret =new Matrix(Row);    //单位阵
385             ret.SetUnit();
386
387             int maxIndex;
388             double dMul;
389
390             for(int i=0;i<Row;i++)
391             {
392
393                 maxIndex = tmp.Pivot(i);
394     
395                 if(tmp.m_data[maxIndex,i]==0)
396                 {
397                     System.Exception e = new Exception("求逆的矩阵的行列式的值等于0,");
398                     throw e;
399                 }
400
401                 if(maxIndex != i)    //下三角阵中此列的最大值不在当前行，交换
402                 {
403                     tmp.Exchange(i,maxIndex);
404                     ret.Exchange(i,maxIndex);
405
406                 }
407
408                 ret.Multiple(i,1/tmp[i,i]);
409
410                 /
411                 //sw.WriteLine("nul /t"+tmp[i,i]+"/t"+ret[i,i]);
412                 
413                 tmp.Multiple(i,1/tmp[i,i]);
414                 //sw.WriteLine("mmm /t"+tmp[i,i]+"/t"+ret[i,i]);
415                 sw.WriteLine("111111111 tmp=/r/n"+tmp);
416                 for(int j=i+1;j<Row;j++)
417                 {
418                     dMul = -tmp[j,i];
419                     tmp.MultipleAdd(j,i,dMul);
420                     ret.MultipleAdd(j,i,dMul);
421     
422                 }
423                 sw.WriteLine("222222222222=/r/n"+tmp);
424
425             }//end for
426
427
428
429             for(int i=Row-1;i>0;i--)
430             {
431                 for(int j=i-1;j>=0;j--)
432                 {
433                     dMul = -tmp[j,i];
434                     tmp.MultipleAdd(j,i,dMul);
435                     ret.MultipleAdd(j,i,dMul);
436                 }
437             }//end for
438         
439             //
440
441
442             sw.WriteLine("tmp = /r/n" + tmp.ToString());
443
444             sw.Close();
445             ///
446             ///
447             return ret;
448
449         }//end Inverse
450
451 */
452
453          #endregion
454
455          // determine if the matrix is square:方阵
456          public   bool  IsSquare()
457          {
458              return  Row == Col;
459         }
460
461          // determine if the matrix is symmetric对称阵
462          public   bool  IsSymmetric()
463          {
464
465              if (Row  !=  Col) 
466                  return   false ;
467                              
468              for ( int  i = 0 ;i < Row;i ++ )
469                  for ( int  j = i + 1 ;j < Col;j ++ )
470                      if ( m_data[i,j]  !=  m_data[j,i])
471                          return   false ;
472
473              return   true ;
474         }
475
476          // 一阶矩阵->实数
477          public   double  ToDouble()
478          {
479             Trace.Assert(Row == 1   &&  Col == 1 );
480
481              return  m_data[ 0 , 0 ];
482         }
483
484          // conert to string
485          public   override   string  ToString()
486          {
487         
488              string  s = "" ;
489              for ( int  i = 0 ;i < Row;i ++ )
490              {
491                  for ( int  j = 0 ;j < Col;j ++ )
492                     s  +=   string .Format( " {0}  " ,m_data[i,j]);
493
494                 s  +=   " /r/n " ;
495             }
496              return  s;
497
498         }
499
500
501          // 私有数据成员
502          private   double [,] m_data;
503         
504     } // end class Matrix
505 }