matrix.cpp

摄影测量课程实习程序：空间后方交会小程序

//Matrix.cpp

#include<iostream.h>
#include"Matrix.h"

Matrix::Matrix()
{
    r=0;c=0;
}

Matrix::Matrix(int h,int l)
{
    float * m[h][l];
    r=h;
    c=l;
    M=m;
}
Matrix::Matrix(int h,int l,float * m)
{
    r=h;
    c=l;
    M=m;
}
//构造对角阵 
Matrix::Matrix(int hl,float n)
{
    float m[hl][hl];
    int i,j;
    //float n;
    //n=float(*N);
    for(i=0;i<hl;i++)
    {
        for(j=0;j<hl;j++)
        {
            if(j==i) m[i][j]=n;
            else m[i][j]=0.0;
        }
    }
}

Matrix::~Matrix()
{
}

void Matrix::exchange(Matrix * X,int h1,int h2)
{
     float tM[*X.c];
     int i;
     for(i=0;i<X.c;i++)
     {
         tM[i]=*(X.M+h1)+i;
         *(X.M+h1)+i=*(X.M+h2)+i;
         *(X.M+h2)+i=tM[i];
     }
}

Matrix Matrix::Trans()
{
       int i,j,t;
       for(i=0;i<r-1;i++)
           for(j=i+1;j<c;j++)
           {
               t=M[i][j];
               M[i][j]=M[j][i];
               M[j][i]=t;
           }
       t=r;
       r=c;
       c=t;
}

Matrix Matrix::Multi(Matrix A,Matrix B)
{
       int h,l;
       if(A.c!=B.r)
       {
           cout<<"违反矩阵相乘法则!";
           exit(1); 
       }
       Matrix C(A.r,B.c);
       for(i=0;i<A.r;i++)
       {
           for(j=0;j<B.c;j++)
           {
               C.M[i][j]=0;
               for(l=0,h=0;l<A.c;l++,h++)
                 C.M[i][j]+=A.M[i][l]*B.M[h][j];
           }
       }
       return C;
}
Matrix Matrix::operator*(Matrix B)
{
       return Multi(*this,B);
}

Matrix Matrix::Invers()
{
       int l,i,flag;
       if(c!=r)
       {
           cout<<"此矩阵不可逆!";
           exit(1); 
       }
       Matrix E(c,1);
       for(l=0;l<c;l++)
       {
           i=l;flag=0;
           while(M[i][l]==0)
           {
               i++;
               if(i==r) {i-=r;  flag=1;}
           }
           if(i==r&&flag==1)
           {
               cout<<"此矩阵不可求逆!";
               exit(1); 
           }
           if(i!=l) exchange(this,i,l);
           //整行除一数 
           for(i=0;i<c;i++)
           {    M[l][i]/=M[l][l]; E.M[l][i]/=M[l][l];    }
           for(i=0;i<r;i++)
           {
               if(i==l||M[i][l]==0) continue;
               for(int j=0;j<c;j++)
               {
                   M[i][j]/=M[i][l]; E.M[i][j]/=M[i][l];//化为1 
                   M[i][j]-=M[l][j]; E.M[i][j]-=E.M[l][j];
               }
           }
       }
       return E;
}