song.yz@foxmail.com wechat: math-box

计算数学

附录:矩阵及线性代数类设计方法



大部分高级语言并不是直接针对数值算法设计。比如C语言中,实现矩阵运算通常通过数组, 而二维动态数组实际通过指针实现。用户需要非常小心的分配与释放内存,否则内存泄漏造成程序异常。 这里以C++和C#为例,设计矩阵与线性代数类,通过类的设计,可以在析构函数中自动释放内存。 还可以通过索引器,使得C++或C#类似于MATLAB一样,实现矩阵运算。在java中不支持运算符重载, 因此java矩阵类,可以把各功能定义成函数。在Python等语言中,已经有第三方数值代数库支持。

C++矩阵向量类与线性代数

关于C++的语法书,请参考Stoustrup编写的《the C++ programming language》等专著或教材。

头文件如下


/*------------------------------------------------------
   Versions and Changes :
      v1.0----2015-7-18
          矩阵向量类与数值代数函数库
          c++  matrix and vector library     
      
          主要包括
          1. 简单的矩阵向量加减乘法运算,坐标旋转矩阵
          2. 常用的矩阵分解 LDL, QR(MGS正交化,householder变换)
          3. 最小二乘问题求解
          4. 一般线性方程求解(正交化方法)
          5. 矩阵求逆
          6. 拉格朗日插值
         
          Notes1: 
          1. 矩阵类的模板设计问题
             矩阵向量类可以设计为模板,这样运算时候可以是的对各种数据类型具有通用性。
             但这里仅把该类设计为double型,因为在数值运算中基本不可能使用单精度。
             另外矩阵分解更不可能针对整型。实际运算中需要用整型矩阵运算的情况较少,
             如果出现可以使用C++标准库vector构造,方法见note2。
		  2. 用c++ STL标准库中vector构造矩阵方法 
  		     构造方法(int ,double ,float),以double为例
			 (1) vector > A(M,vector(N))
			 (2) vector > A(M)
			     for(int i=0;i(N)
			 (3) vector > A(N)
			     for(int i=0;i
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/ 

#ifndef cmatlib_h
#define cmatlib_h
#include 
#include 
#include   //format control 
using namespace std; 
class VEC
/*------------------------------------------------------
   Versions and Changes :
      v1.0----2015-7-18
          vector class  
          下标索引算符为()
--------------------------------------------------------
   Public Paras:
           int size1 ---------  rows  of matrix
           int size2 ---------  cols  of matri
   Funcitons :       
          VEC (M)  ---- creat a vector with lenght of M
          operator() --- access the component    
--------------------------------------------------------
   Author      : Song Yezhi  
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/ 
{
   private:
      double *x ;
      int M1;
   public:
      // Constructors
      VEC (int M)
      //creat a vector
      { 
        M1=M;
        x= new double [M] ;
        for(int i=0;i< M;i++)
          x[i]=0.0;
      }      
      // Destructor
      ~ VEC ()
      { delete [] x; }
      
      double& operator() (int i)
      // Component access 
      { return x[i]; }           
      double  operator () (int i) const { return x[i]; };      
      int size() const { return M1; };      
      void  output(int width,int precision);     
	  //output to the screen 
	  void  setv(double value);
	  //set the elements of the vector to a double value 
	    
     //+,-,*,/ for the same index of two vectors
     //which means that the size of the vectors must be the same 
     friend VEC operator + (const VEC& V1, const VEC& V2);        //V1+V2   
     friend VEC operator - (const VEC& V1, const VEC& V2);        //V1-V2
     friend VEC operator * (const VEC& V1, const VEC& V2);        //V1.*V2, for the matlab notation  
     friend VEC operator / (const VEC& V1, const VEC& V2);        //V1./V2, for the matlab notation  
     //V1./V2
     //+,-,*,/ for the vector and a scala
     // the operation works on every elements of the vector
     friend VEC operator + (const VEC& V1, double a);            //V1+a
     friend VEC operator - (const VEC& V1, double a);            //V1-a
     friend VEC operator * (const VEC& V1, double a);            //V1*a
     friend VEC operator / (const VEC& V1, double a);            //V1/a                   
} ;
 
class MAT
/*------------------------------------------------------
   Versions and Changes :
      v1.0----2015-7-18
         matrix class
          在数组基础上构造,下标算符重载()
--------------------------------------------------------
   Public Paras:
           int size1 ---------  rows  of matrix
           int size2 ---------  cols  of matrix
   Functions   :          
           MAT(M,N)  ---- creat a vector with the size of (M,N)
           operator() ---  caccess the component 
--------------------------------------------------------
   Author      : Song Yezhi  
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/ 
{
   private:
     double ** A;
     int M1,N1;
   public:
     // Constructors
     MAT(int M,int N)
     {
        M1=M;
        N1=N;  
        
        int i,j;
        
        A= new double *[M];
        for(i=0;i< M;i++)
          A[i]=new double [N] ;
     
            
        // set the initial value to zero 
        for (i=0;i< M;i++)
         for(j=0;j< N;j++)
            A[i][j]=0.0;  
     }
     // Destructor
     ~MAT()
     {
         int i;
         for(i=0;i< M1;i++)
           delete[] A[i];
         delete[] A;
     }
    
    //Component access
     double   operator() (int i,int j) const  { return A[i][j]; };             
     double & operator() (int i,int j)        { return A[i][j]; };           
     
     int size1() const { return M1; };
     int size2() const { return N1; };
        
     void output(int width,int precision);
	 //output to the screen 
	 void setv(double value);
	 //set the elements of the matrix to a double value 

     //+,-,*,/ for the same index of two matrixs
     //which means that the size of the matrixs must be the same 
     friend MAT operator + (const MAT& A1, const MAT& A2);   //A1+A2
     friend MAT operator - (const MAT& A1, const MAT& A2);   //A1-A2
     friend MAT operator * (const MAT& A1, const MAT& A2);   //A1.*A2
     friend MAT operator / (const MAT& A1, const MAT& A2);   //A1./A2
     
     //+,-,*,/ for the matrix and a scala
     // the operation works on every elements of the matrix
     friend MAT operator + (const MAT& A1, double a);      //A1+a
     friend MAT operator - (const MAT& A1, double a);      //A1-a
     friend MAT operator * (const MAT& A1, double a);      //A1*a
     friend MAT operator / (const MAT& A1, double a);      //A1/a
              
};

//--------------基本向量运算
void veccopy(VEC &b,const VEC &a);
//按值复制一个向量
double vecdot(const VEC &v1,const VEC & v2);
//向量内积
double norm(const VEC &v);
//向量2范数
//-----------------基本矩阵运算
void matcopy(MAT &B,const MAT &A);
//矩阵复制
void transpose(MAT &AT,const MAT &A);
//矩阵转置
void matmul(VEC &b,const MAT &A,const VEC &x);
// 矩阵乘以向量 
void matmul(MAT &C,const MAT &A,const MAT &B);
//矩阵乘以矩阵
//------------------矩阵分解与求逆
void LDL(const MAT &A,MAT &L,VEC &D);
// LDL 分解
void MGS(const MAT &A, MAT &Q, MAT &R);
//修正的Gram-Schmidt正交化方法
void householder(const MAT &A,MAT &Q,MAT &R);
//hoseholder正交变换
void invlowtri(const MAT &R,MAT &S);
//inverse of lower triangula matrix
void invuptri(const MAT &U,MAT &R);
//inverse of upper triangula matrix
void invmat(const MAT &A,MAT &invA);
//对称正定矩阵逆矩阵
void inv(const MAT &A,MAT &iA);
//一般矩阵的逆矩阵 采用MGS分解方法
//--------------------线性代数方程
void LS_LDL(VEC &x,const MAT &A,const VEC &b);
//基于不开平方的cholesky分解计算最小二乘问题, A为对称正定矩阵
//void LS_hous(VEC &x,const MAT &A,const VEC &b);
//基于Householder变换求解最新二乘问题或适定问题
void LS_MGS( VEC &x,const MAT &A,const VEC &b);
//通过修正的Gram-Schimdt正交化求解最小二乘问题或适定问题
void LS_hous( VEC &x,const MAT &A,const VEC &b);
//least square solution or general linear equation 

void uptri(const MAT &A,const VEC &b,VEC &x);
//上三角矩阵方程计算
void lowtri(const MAT &A, const VEC &b, VEC &x);
//下三角矩阵方程计算

void downtri(const MAT &A, const VEC &b, VEC &x);

//------------------插值算法
void lagrange(const VEC &x,const VEC &y,const VEC &xx,VEC &yy);
//lagrange interp
void rot_mat(MAT &Rmat,double angle,char axID);
//rotation matrix 
double robustweight(double v,double sigma);
//used for robust estimation 

#endif
 

C++代码如下:

/*------------------------------------------------------
   Versions and Changes :
      v1.0----2015-7-18
       c++ matrix library
       the comment is in  cmatlib.h
        
      v2.0---2015-11-05
        the second editon of cmatlib

                        
--------------------------------------------------------
   Notes:          
        1. 头文件为 matlib.h
        2. 调用该函数库需要加入头文件 
           #include "matlib.h"
         
        3. 编译方式为  $(cc) main.cpp matlib.cpp  

--------------------------------------------------------
   Author      : Song Yezhi  
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/ 

//#include 
//#include 
//#include   //
#include "cmatlib.h" 
//
//using namespace std;

void VEC::output(int width,int precision)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-09-11 
      output the vector to the screen
--------------------------------------------------------
   Input Parameters   : 
        width   ---  the width of the element value
        precision--  presicion of the element value
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
	int i; 	
	for(i=0;i< M1;i++)
	 	cout<
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
	int i;
	for(i=0;i< M1;i++)
		x[i]=value;
}


void MAT::output(int width,int precision) 
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-09-11 
      output the matrix to the screen
      in many case this function is used for debug
--------------------------------------------------------
   Input Parameters   :
        width   ---  the width of the element value
        precision--  presicion of the element value
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
   int i,j;	
   for(i=0;i< M1;i++)
   {
   	 for(j=0;j< N1;j++) cout<
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
	int i,j;
	for(i=0;i< M1;i++)
	{
		for(j=0;j< N1;j++)
			A[i][j]=value;
	}
}

VEC operator + (const VEC& V1, const VEC& V2)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      addition of the the two vectors
      V1+V2
--------------------------------------------------------
   Input Parameters   :
        V1
        V2
   Output Parameters  :
        V3=
   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
	int N=V1.size();
	int i;
	VEC VT(N);
	for(i=0;i< N;i++)
		VT(i)=V1(i)+V2(i);	
	return VT;
}

VEC operator - (const VEC& V1, const VEC& V2)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
    vector subtraction  V1-V2
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
	int N=V1.size();
	int i;
	VEC VT(N);
	for(i=0;i< N;i++)
		VT(i)=V1(i)-V2(i);
	
	return VT;
}

VEC operator * (const VEC& V1, const VEC& V2)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      V1.*V2
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
	int N=V1.size();
	int i;
	VEC VT(N);
	for(i=0;i< N;i++)
		VT(i)=V1(i)*V2(i);
	
	return VT;
}

VEC operator / (const VEC& V1, const VEC& V2)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      V1./V2
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
	int N=V1.size();
	int i;
	VEC VT(N);
	for(i=0;i< N;i++)
		VT(i)=V1(i)/V2(i);
	
	return VT;
}


VEC operator + (const VEC& V1,double a)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      V1+a
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
	int N=V1.size();
	int i;
	VEC VT(N);
	for(i=0;i< N;i++)
		VT(i)=V1(i)+a;
	return VT;
}

VEC operator - (const VEC& V1,double a)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      V1-a
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
	int N=V1.size();
	int i;
	VEC VT(N);
	for(i=0;i< N;i++)
		VT(i)=V1(i)-a;
	return VT;
}

VEC operator * (const VEC& V1,double a)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      V1*a
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
	int N=V1.size();
	int i;
	VEC VT(N);
	for(i=0;i< N;i++)
		VT(i)=V1(i)*a;
	return VT;
}

VEC operator / (const VEC& V1,double a)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      V1/a
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
	int N=V1.size();
	int i;
	VEC VT(N);
	for(i=0;i< N;i++)
		VT(i)=V1(i)/a;
	return VT;
}
void veccopy(VEC &b,const VEC &a)
/*------------------------------------------------------
   Versions and Changes :
      v1.0----2015-8-5
         按值复制一个向量    
--------------------------------------------------------
   Input Parameters   :
        a
   Output Parameters  :
        b
--------------------------------------------------------
   Author      : Song Yezhi  
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/ 
{
	int N=a.size(); 
	int M=b.size();	
	if(N!=M)
     {cout<<"the dimention of two vectors must be the same"<
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)            2015
------------------------------------------------------- */ 
{ 
    int N=v1.size();     
    double sum=0; 
    for(int i=0;i< N;i++)
       sum=sum+v1(i)*v2(i);        
    return sum;        
}

double norm(const VEC &v)
/*------------------------------------------------------
   Versions and Changes :
      v1.0----2015-7-30
          norm-2 of a vector  
--------------------------------------------------------
   Input Parameters   :
        v1---  vector
        v2---  vector 
   Output Parameters  :
        a --- double
 
--------------------------------------------------------
   Author      : Song Yezhi  
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)            2015
------------------------------------------------------- */ 
{ 
    int N=v.size();     
    double sum=0; 
    for(int i=0;i< N;i++)
       sum=sum+v(i)*v(i);        
    sum=sqrt(sum);
    return sum;        
}

void matcopy(MAT &B,const MAT &A)
/* ------------------------------------------------------
   Versions and Changes :
      v1.0----2015-8-5
          按值复制矩阵
          MAT matcopy(const MAT &A)
      v1.1--- 改为void类型函数
          通过地址传递变了,因为MAT的()算符重载问题
--------------------------------------------------------
   Author      : Song Yezhi  
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
------------------------------------------------------- */
{
    int M=A.size1();
	int N=A.size2();    
    int i,j;
    for(i=0;i< M;i++)
     for(j=0;j< N;j++)
      B(i,j)=A(i,j);  
}

void transpose(MAT &AT,const MAT &A)
/* ------------------------------------------------------
   Versions and Changes :
      v1.0----2015-7-16
          矩阵转置
--------------------------------------------------------
   Author      : Song Yezhi  
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
------------------------------------------------------- */
 {
     int M, N;
     M = A.size1();
     N = A.size2();     
     int P=AT.size1();
     int Q=AT.size2();     
     if(N!=P||M!=Q)   
     {cout<<"the dimention of two vectors must be the same"<
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
   int row=A1.size1();
   int col=A1.size2();   
   MAT A3(row,col);      
   for(int i=0;i< row;i++)
    for(int j=0;j< col;j++)
      A3(i,j)=A1(i,j)+A2(i,j);  
   return A3;
	
}

MAT operator -(const MAT& A1,const MAT& A2)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      A1-A2
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
   int row=A1.size1();
   int col=A1.size2();   
   MAT A3(row,col);      
   for(int i=0;i< row;i++)
    for(int j=0;j< col;j++)
      A3(i,j)=A1(i,j)-A2(i,j);  
   return A3;	
}

MAT operator *(const MAT& A1,const MAT& A2)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      A1.*A2
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
   int row=A1.size1();
   int col=A1.size2();   
   MAT A3(row,col);      
   for(int i=0;i< row;i++)
    for(int j=0;j< col;j++)
      A3(i,j)=A1(i,j)*A2(i,j);  
   return A3;	
}

MAT operator /(const MAT& A1,const MAT& A2)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      A1./A2
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
   int row=A1.size1();
   int col=A1.size2();
   
   MAT A3(row,col);
      
   for(int i=0;i< row;i++)
    for(int j=0;j< col;j++)
      A3(i,j)=A1(i,j)/A2(i,j);  
   return A3;	
}

MAT operator +(const MAT& A1,double a)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      A1+a
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
   int row=A1.size1();
   int col=A1.size2();
   
   MAT A3(row,col);
      
   for(int i=0;i< row;i++)
    for(int j=0;j< col;j++)
      A3(i,j)=A1(i,j)+a;  
   return A3;	
}

MAT operator -(const MAT& A1,double a)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
      A1-a
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
   int row=A1.size1();
   int col=A1.size2();
   
   MAT A3(row,col);
      
   for(int i=0;i< row;i++)
    for(int j=0;j< col;j++)
      A3(i,j)=A1(i,j)-a;  
   return A3;	
}

MAT operator *(const MAT& A1,double a)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
     A1*a  
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
   int row=A1.size1();
   int col=A1.size2();
   
   MAT A3(row,col);
      
   for(int i=0;i< row;i++)
    for(int j=0;j< col;j++)
      A3(i,j)=A1(i,j)*a;  
   return A3;	
}

MAT operator /(const MAT& A1,double a)
/*------------------------------------------------------
   Versions and Changes :
    v1.0---- 2015-11-05 
     A1/a 
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

   Notes :         
--------------------------------------------------------
   Author      : Song Yezhi           
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/
{
   int row=A1.size1();
   int col=A1.size2();
   
   MAT A3(row,col);
      
   for(int i=0;i< row;i++)
    for(int j=0;j< col;j++)
      A3(i,j)=A1(i,j)/a;  
   return A3;	
}

void matmul(VEC &b,const MAT &A,const VEC &x)
/*-----------------------------------------------------
   Versions and Changes :
      v1.0----2015-7-17
          b=Ax     
--------------------------------------------------------
   Input Parameters   :
        A(M,N)---  matrix
        x(N)---  vector
   Output Parameters  :
        b(M) --- vector 
--------------------------------------------------------
   Author      : Song Yezhi  
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)             2015
-------------------------------------------------------*/ 
{
   int M,N;
   M=A.size1();
   N=A.size2(); 
   double sum;
   int i,j;      
   for (i=0;i< M;i++)
   {
       sum=0.0;
       for(j=0;j< N;j++)
        sum=sum+A(i,j)*x(j);      
       b(i)=sum; 
   }     
}

void matmul(MAT &C,const MAT &A,const MAT &B)
/* ------------------------------------------------------
   Versions and Changes :
      v1.0----2015-7-17
          C=AB     
--------------------------------------------------------
   Input Parameters   :
        A(M,N)---  matrix
        B(N,P)---  matrix
   Output Parameters  :
        C(M,P) --- matrix
 
--------------------------------------------------------
   Author      : Song Yezhi      
   Copyrigt(C) : Shanghai Astronomical Observatory, CAS
                (All rights reserved)            2015
-------------------------------------------------------*/
{
	int M=A.size1();
	int N=A.size2();
	int P=C.size2();
	int i,j,k;			
	double tmp1,tmp2;
	double sum;
	for(i=0;i< M;i++)
	 for(j=0;j< P;j++)
	 {
	 	sum=0.0;
	 	for(k=0;k< N;k++) sum=sum+A(i,k)*B(k,j);	 	  	  
	 	C(i,j)=sum;	 	
	 }
	return;	
}