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

第一次搬进自己的blog，打算放点东东到里面。在工程软件开发过程中，会碰到很多有关矩阵的运算。这是一年前的源代码

  1using System;
  2using System.IO;
  3using System.Diagnostics;
  4
  5
  6namespace 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            }
306StreamWriter 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                }
342sw.WriteLine("tmp=\r\n"+tmp);
343sw.WriteLine("ret=\r\n"+ret);
344            }//end for
345
346
347sw.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
360sw.WriteLine("tmp=\r\n"+tmp);
361sw.WriteLine("ret=\r\n"+ret);
362sw.WriteLine("***=\r\n"+ this*ret);
363sw.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}

上面的矩阵类，用起来很方便。

例子：

//Square n*n Matrix 
Matrix m1 = new Matrix(n);
//m*n Matrix
Matrix m2 = new Matrix(m,n):
//Matrix Add Operattion
Mtrix result = m_1 + m_2;
//Display Matrix Method 1
for(int i=0;i<m.Row;i++)
  {
    for(int j=0;j<m.Col;j++)
      Console.WriteLine(m[i,j]);
    Console.Writeline();
   }
//Display Matrix Method 2
Console.WriteLine(m.ToString()); 

说明：类实现了矩阵的所有基本操作，包括单双目的＋，－，×、/(Inversion)，初等变换（三种数学上的初等变换），还有很多，希望读者自己体会。另有标准C＋＋源代码，实现了上面的所有的功能。