矩阵运算.txt

关于矩阵计算的一些代码

290
291    for(int k=0;k<m_col;k++)
292    {
293        swap(m_data[i][k],m_data[j][k]);
294    }
295
296    return *this;
297}
298
299
300Matrix&  Matrix::Multiple(int index,double mul)    //初等变换　第index 行乘以mul
301{
302    if(index <0 || index >= m_row)
303    {
304        ostringstream oss;
305        oss<<"In void Matrix::Multiple(int index,double k) index "<<index<<" out of bounds.Please check it.";
306        throw oss.str();
307    }
308
309    transform(m_data[index],m_data[index]+m_col,m_data[index],bind2nd(multiplies<double>(),mul));
310    return *this;
311}
312
313
314Matrix& Matrix::MultipleAdd(int index,int src,double mul)    //初等变换 第src行乘以mul加到第index行
315{
316    if(index <0 || index >= m_row)
317    {
318        ostringstream oss;
319        oss<<"In void MultipleAdd(int index,int src,double k) index "<<index<<" out of bounds.Please check it.";
320        throw oss.str();
321    }
322
323    if(src < 0 || src >= m_row)
324    {
325        ostringstream oss;
326        oss<<"In void MultipleAdd(int index,int src,double k) src "<<src<<" out of bounds.Please check it.";
327        throw oss.str();
328    }
329        
330    for(int j=0;j<m_col;j++)
331    {
332        m_data[index][j] += m_data[src][j] * mul;
333    }
334
335    return *this;
336
337}
338
339        
340//inverse 逆阵：使用矩阵的初等变换，列主元素消去法
341Matrix Matrix::Inverse() const
342{
343    if(m_row != m_col)    //非方阵
344    {
345        ostringstream oss;
346        oss<<"In Matrix Matrix::Invert() const. m_row "<<m_row<<" != m_col "<<m_col<<" Please check it";
347        throw oss.str();
348    }
349
350    Matrix tmp(*this);
351    Matrix ret(m_row);    //单位阵
352    ret.SetUnit();
353
354    int maxIndex;
355    double dMul;
356
357    for(int i=0;i<m_row;i++)
358    {
359
360        maxIndex = tmp.Pivot(i);
361    
362        if(tmp.m_data[maxIndex][i]==0)
363        {
364            ostringstream oss;
365            oss<<"In Matrix Matrix::Invert() const 行列式的值等于0.Please check it";
366            throw oss.str();
367        }
368
369        if(maxIndex != i)    //下三角阵中此列的最大值不在当前行，交换
370        {
371            tmp.Exchange(i,maxIndex);
372            ret.Exchange(i,maxIndex);
373
374        }
375
376        ret.Multiple(i,1/tmp.m_data[i][i]);
377        tmp.Multiple(i,1/tmp.m_data[i][i]);
378
379
380        for(int j=i+1;j<m_row;j++)
381        {
382            dMul = -tmp.m_data[j][i];
383            tmp.MultipleAdd(j,i,dMul);
384            ret.MultipleAdd(j,i,dMul);
385    
386        }
387
388    }//end for
389
390    for(int i=m_row-1;i>0;i--)
391    {
392        for(int j=i-1;j>=0;j--)
393        {
394            dMul = -tmp.m_data[j][i];
395            tmp.MultipleAdd(j,i,dMul);
396            ret.MultipleAdd(j,i,dMul);
397        }
398    }//end for
399        
400    return ret;
401}
402
403    
404
405// assignment operator　賦值运算符
406Matrix& Matrix::operator= (const Matrix& m)
407{
408    if(m_data != 0)
409    {
410        for(int i=0;i<m_row;i++)
411        {
412            delete[] m_data[i];
413        }
414        delete[] m_data;
415    }
416
417    m_row = m.m_row;
418    m_col = m.m_col;
419
420    m_data = new double*[m_row];
421    for(int i=0;i< m_row;i++)
422    {
423        m_data[i] = new double[m_col];
424    }
425
426    for(int i=0;i<m_row;i++)
427    {
428        copy(m.m_data[i],m.m_data[i]+m_col,m_data[i]);
429    }
430    return *this;
431}
432
433
434//binary addition 矩阵加
435Matrix Matrix::operator+ (const Matrix& m) const
436{
437    if(m_row != m.m_row)
438    {
439        ostringstream oss;
440        oss<<"In Matrix::operator+(const Matrix& m) this->m_row "<<m_row<<" != m.m_row "<<m.m_row<<" Please check it.";
441        throw oss.str();
442    }
443    if(m_col != m.m_col)
444    {
445        ostringstream oss;
446        oss<<" Matrix::operator+ (const Matrix& m) this->m_col "<<m_col<<" != m.m_col "<<m.m_col<<" Please check it.";
447        throw oss.str();
448    }
449
450    Matrix ret(m_row,m_col);
451    for(int i=0;i<m_row;i++)
452    {
453        transform(m_data[i],m_data[i]+m_col,m.m_data[i],ret.m_data[i],plus<double>());
454    }
455    return ret;
456}
457
458//unary addition 求正
459Matrix Matrix::operator+() const
460{
461    return *this;
462}
463
464//binary subtraction 矩阵减
465Matrix Matrix::operator- (const Matrix& m) const
466{
467    if(m_row != m.m_row)
468    {
469        ostringstream oss;
470        oss<<"In Matrix::operator-(const Matrix& m) this->m_row "<<m_row<<" != m.m_row "<<m.m_row<<" Please check it.";
471        throw oss.str();
472    }
473    if(m_col != m.m_col)
474    {
475        ostringstream oss;
476        oss<<"In Matrix::operator- (const Matrix& m) this->m_col "<<m_col<<" != m.m_col "<<m.m_col<<" Please check it.";
477        throw oss.str();
478    }
479
480    Matrix ret(m_row,m_col);
481    for(int i=0;i<m_row;i++)
482    {
483        transform(m_data[i],m_data[i]+m_col,m.m_data[i],ret.m_data[i],minus<double>());
484    }
485    return ret;
486}
487
488//unary substraction 求负
489Matrix Matrix::operator-() const
490{
491    Matrix ret(*this);
492
493    for(int i=0;i<m_row;i++)
494        for(int j=0;j<m_col;j++)
495            ret.m_data[i][j] *= -1;
496
497    return ret;
498
499}
500
501
502//binary multiple 矩阵乘
503Matrix Matrix::operator*(const Matrix& m) const
504{
505    if(m_col != m.m_row)
506    {
507        ostringstream oss;
508        oss<<"In Matrix::operator*(const Matrix& m) this->m_col "<<m_col<<" != m.m_row "<<m.m_row<<" Please check it.";
509        throw oss.str();
510    }
511
512    Matrix ret(m_row,m.m_col);
513    Matrix tmp=m.Transpose();
514
515    for(int i=0;i<m_row;i++)
516    {
517        for(int j=0;j<m.m_col;j++)
518        {
519            ret.m_data[i][j]=inner_product(m_data[i],m_data[i]+m_col,tmp.m_data[j],0.0);
520        }
521    }
522
523    return ret;;
524}
525
526//friend scalar multiple　数乘
527Matrix operator* (double d,const Matrix& m)
528{
529    Matrix ret(m);
530
531    for(int i=0;i<ret.m_row;i++)
532        for(int j=0;j<ret.m_col;j++)
533            ret.m_data[i][j] *= d;
534
535    return ret;
536
537}
538
539//binary matrix division equivalent to multiple inverse矩阵除，等价于乘以逆阵
540Matrix Matrix::operator/(const Matrix& m) const
541{
542    return *this * m.Inverse();
543}
544
545//binary scalar division equivalent to multiple reciprocal数除，等价于乘以此数的倒数
546Matrix Matrix::operator/(double d) const
547{
548    return 1/d * (*this);
549}
550
551//friend division
552Matrix operator/(double d,const Matrix& m)
553{
554    return d * m.Inverse();
555}
556
557
558// subscript operator to get individual elements　下标运算符
559double Matrix::operator()(int i,int j) const
560{
561    if(i<0 || i>= m_row)
562    {
563        ostringstream oss;
564        oss<<"In double Matrix::operator()(int i,int j) const.i "<<i<<" out of bounds.Please check it";
565        throw oss.str();
566    }
567    if(j<0 || j>= m_col)
568    {
569        ostringstream oss;
570        oss<<"In double Matrix::operator()(int i,int j) const.j "<<j<<" out of bounds.Please check it";
571        throw oss.str();
572    }
573
574    return m_data[i][j];
575}
576
577// subscript operator to set individual elements　下标运算符
578double& Matrix::operator ()(int i,int j)
579{
580    if(i<0 || i>= m_row)
581    {
582        ostringstream oss;
583        oss<<"In double Matrix::operator()(int i,int j) const.i "<<i<<" out of bounds.Please check it";
584        throw oss.str();
585    }
586    if(j<0 || j>= m_col)
587    {
588        ostringstream oss;
589        oss<<"In double Matrix::operator()(int i,int j) const.j "<<j<<" out of bounds.Please check it";
590        throw oss.str();
591    }
592
593    return m_data[i][j];
594}
595
596//to string　化为标准串
597string Matrix::ToString() const
598{
599    ostringstream oss;
600    for(int i=0;i<m_row;i++)
601    {
602        copy(m_data[i],m_data[i]+m_col,ostream_iterator<double>(oss," "));
603        oss<<"\n";
604    }
605
606    return oss.str();
607}
608
609// outputt stream function　输出
610
611ostream& operator<<(ostream& os,const Matrix& m)
612{
613    for(int i=0;i<m.m_row;i++)
614    {
615        copy(m.m_data[i],m.m_data[i]+m.m_col,ostream_iterator<double>(os,"  "));
616        os<<'\n';
617    }
618    return os;
619}
620
621// input  stream function　输入
622istream& operator>>(istream& is,Matrix& m)
623{
624    for(int i=0;i<m.m_row;i++)
625    {
626        for(int j=0;j<m.m_col;j++)
627        {
628            is>>m.m_data[i][j];
629        }
630    }
631    return is;
632}
633
634
635// reallocation method
636
637// public method for resizing matrix
638
639// logical equal-to operator
640
641// logical no-equal-to operator
642
643// combined addition and assignment operator
644
645// combined subtraction and assignment operator
646
647// combined scalar multiplication and assignment operator
648
649
650// combined matrix multiplication and assignment operator
651
652// combined scalar division and assignment operator
653
654// combined power and assignment operator
655
656// unary negation operator
657
658// binary addition operator
659
660// binary subtraction operator
661
662// binary scalar multiplication operator
663
664// binary scalar multiplication operator
665
666// binary matrix multiplication operator
667
668// binary scalar division operator
669
670
671// binary scalar division operator
672
673// binary matrix division operator
674
675// binary power operator
676
677// unary transpose operator
678
679// unary Inverse operator
680
681// Inverse function
682
683// solve simultaneous equation
684
685// set zero to all elements of this matrix
686
687// set this matrix to unity
688
689// private partial pivoting method
690
691// calculate the determinant of a matrix
692
693// calculate the norm of a matrix
694
695// calculate the condition number of a matrix
696
697// calculate the cofactor of a matrix for a given element
698
699// calculate adjoin of a matrix
700
701// Determine if the matrix is singular
702
703// Determine if the matrix is diagonal
704
705// Determine if the matrix is scalar
706
707// Determine if the matrix is a unit matrix
708
709// Determine if this is a null matrix
710
711// Determine if the matrix is symmetric
712
713// Determine if the matrix is skew-symmetric
714// Determine if the matrix is upper triangular
715
716// Determine if the matrix is lower triangular