class DMatrix5x2{
  public:
    DMatrix5x2(){
      for (unsigned int j=0;j<3;j++){
	mA[0].v[j]=_mm_setzero_pd();
	mA[1].v[j]=_mm_setzero_pd();
      }
    }
    DMatrix5x2(__m128d aa,__m128d ab,__m128d ba,__m128d bb,__m128d ca,__m128d cb){
      mA[0].v[0]=aa;
      mA[0].v[1]=ba;
      mA[0].v[2]=ca;
      mA[1].v[0]=ab;
      mA[1].v[1]=bb; 
      mA[1].v[2]=cb;
    }
    ~DMatrix5x2(){};

    __m128d GetV(int pair,int col) const{
      return mA[col].v[pair];
    }
    void SetV(int pair,int col,__m128d v){
      mA[col].v[pair]=v;
    }
      
    double &operator() (int row,int col){
      return mA[col].d[row];
    } 
    double operator() (int row,int col) const{
      return mA[col].d[row];
    }
    
    DMatrix5x1 operator*(const DMatrix2x1 &m2){
      __m128d a=_mm_set1_pd(m2(0));
      __m128d b=_mm_set1_pd(m2(1));
      return DMatrix5x1(_mm_add_pd(_mm_mul_pd(GetV(0,0),a),
				   _mm_mul_pd(GetV(0,1),b)),
			_mm_add_pd(_mm_mul_pd(GetV(1,0),a),
				   _mm_mul_pd(GetV(1,1),b)),
			_mm_add_pd(_mm_mul_pd(GetV(2,0),a),
				   _mm_mul_pd(GetV(2,1),b)));
    }

    DMatrix5x2 operator*(const DMatrix2x2 &m2){
      __m128d a11=_mm_set1_pd(m2(0,0)); // row,col
      __m128d a12=_mm_set1_pd(m2(0,1));
      __m128d a21=_mm_set1_pd(m2(1,0));
      __m128d a22=_mm_set1_pd(m2(1,1));
      
      return DMatrix5x2(_mm_add_pd(_mm_mul_pd(GetV(0,0),a11),
				   _mm_mul_pd(GetV(0,1),a21)),
			_mm_add_pd(_mm_mul_pd(GetV(0,0),a12),
				   _mm_mul_pd(GetV(0,1),a22)),
			_mm_add_pd(_mm_mul_pd(GetV(1,0),a11),
				   _mm_mul_pd(GetV(1,1),a21)),
			_mm_add_pd(_mm_mul_pd(GetV(1,0),a12),
				   _mm_mul_pd(GetV(1,1),a22)),	
			_mm_add_pd(_mm_mul_pd(GetV(2,0),a11),
				   _mm_mul_pd(GetV(2,1),a21)),
			_mm_add_pd(_mm_mul_pd(GetV(2,0),a12),
				   _mm_mul_pd(GetV(2,1),a22)));

    }




    void Print(){
      cout << "DMatrix5x2:" <<endl;
      cout << "     |      0    |      1    |" <<endl;
      cout << "----------------------------------" <<endl;
      for (unsigned int i=0;i<5;i++){
	cout <<"   "<<i<<" |"<< setw(11)<<setprecision(4) << mA[0].d[i] 
	     << setw(11)<<setprecision(4) << mA[1].d[i]<< endl;
      }      
    }
    
  private:
    union dvec{
      __m128d v[3];
      double d[6];
    };
    union dvec mA[2];
    
  };