src/Matrix/n_matrix.cc
// ____________________________________________________________________________
// F              CLASS N_MATRIX BINARY ARITHMETIC FUNCTIONS
// ____________________________________________________________________________
F              CLASS N_MATRIX BINARY ARITHMETIC FUNCTIONS
F CLASS N_MATRIX BINARY ARITHMETIC FUNCTIONS
_matrix* n_matrix::add(_matrix* mx)
  {
  if((rows_!=mx->rows()) || (cols_!=mx->cols()))
    {
    NMxerror(51,1);				
    NMxfatal(20);				
    return mx;
    }
  else
    switch(mx->stored_type())
      {
      case n_matrix_type:				
	{ 
	n_matrix* sum = new n_matrix(rows_,cols_);	
	complex *mij = &((n_matrix*)mx)->data[size-1];	
	complex *nij = &data[size-1];		    	
	complex *sij = &(sum->data[size-1]);	        
	for(; nij>=data; sij--,nij--,mij--)		
	  (*sij) = (*mij)+(*nij);			
	return sum;
	}
	break;
      case h_matrix_type:				
	{
  	n_matrix* sum = new n_matrix(rows_,cols_);	
        complex *hij = ((h_matrix*)mx)->data;		
        complex *nij = data;				
        complex *mij = ((n_matrix*)sum)->data;		
        complex *nji, *mji;				
        int i,j; 
        for(i=0; i<rows_; i++, nij+=i, mij+=i)		
          {
    	  (*mij) = (*nij) + (*hij);			
          mji = mij + cols_;				
          nji = nij + cols_;				
          mij++;					
          nij++;					
          hij++;					
  	  for(j=i+1; j<cols_; j++,hij++,mij++,nij++)	
            {
    	    (*mij) = (*nij) + (*hij);			
  	    (*mji) = (*nji) + conj(*hij);		
            mji += cols_;				
            nji += cols_;				
            }
          } 
	return sum;
	}
	break;
      case d_matrix_type:				
	{
	n_matrix* sum = new n_matrix(*this);		
	for(int i=0; i<rows_; i++)			
	  sum->data[i*rows_+i] +=			
                             ((d_matrix*)mx)->data[i];	
        sum->unitary = false;
	return sum;
	}
	break;
      case i_matrix_type:				
	{
	n_matrix* sum = new n_matrix(*this);		
	for(int i=0; i<rows_; i++)			
	  sum->data[i*rows_+i] += complex1;		
        sum->unitary = false;
	return sum;					
	}
	break;
      default:						
	{ 
	n_matrix* sum = new n_matrix(rows_,cols_);
	int pos = 0;
	for(int i=0; i<rows_; i++)
	  for(int j=0; j<cols_; j++,pos++)
	    sum->data[pos] = data[pos]+(*mx)(i,j);
	return sum;
	}
      }
    }
_matrix * n_matrix :: add ( _matrix * mx ) { if ( ( rows_ != mx -> rows ( ) ) || ( cols_ != mx -> cols ( ) ) ) { NMxerror ( 51 , 1 ) ; NMxfatal ( 20 ) ; return mx ; } else switch ( mx -> stored_type ( ) ) { case n_matrix_type : { n_matrix * sum = new n_matrix ( rows_ , cols_ ) ; complex * mij = & ( ( n_matrix * ) mx ) -> data [ size - 1 ] ; complex * nij = & data [ size - 1 ] ; complex * sij = & ( sum -> data [ size - 1 ] ) ; for ( ; nij >= data ; sij -- , nij -- , mij -- ) ( * sij ) = ( * mij ) + ( * nij ) ; return sum ; } break ; case h_matrix_type : { n_matrix * sum = new n_matrix ( rows_ , cols_ ) ; complex * hij = ( ( h_matrix * ) mx ) -> data ; complex * nij = data ; complex * mij = ( ( n_matrix * ) sum ) -> data ; complex * nji , * mji ; int i , j ; for ( i = 0 ; i < rows_ ; i ++ , nij += i , mij += i ) { ( * mij ) = ( * nij ) + ( * hij ) ; mji = mij + cols_ ; nji = nij + cols_ ; mij ++ ; nij ++ ; hij ++ ; for ( j = i + 1 ; j < cols_ ; j ++ , hij ++ , mij ++ , nij ++ ) { ( * mij ) = ( * nij ) + ( * hij ) ; ( * mji ) = ( * nji ) + conj ( * hij ) ; mji += cols_ ; nji += cols_ ; } } return sum ; } break ; case d_matrix_type : { n_matrix * sum = new n_matrix ( * this ) ; for ( int i = 0 ; i < rows_ ; i ++ ) sum -> data [ i * rows_ + i ] += ( ( d_matrix * ) mx ) -> data [ i ] ; sum -> unitary = false ; return sum ; } break ; case i_matrix_type : { n_matrix * sum = new n_matrix ( * this ) ; for ( int i = 0 ; i < rows_ ; i ++ ) sum -> data [ i * rows_ + i ] += complex1 ; sum -> unitary = false ; return sum ; } break ; default : { n_matrix * sum = new n_matrix ( rows_ , cols_ ) ; int pos = 0 ; for ( int i = 0 ; i < rows_ ; i ++ ) for ( int j = 0 ; j < cols_ ; j ++ , pos ++ ) sum -> data [ pos ] = data [ pos ] + ( * mx ) ( i , j ) ; return sum ; } } }
____________________________________________________________________________
F              CLASS N_MATRIX BINARY ARITHMETIC FUNCTIONS
____________________________________________________________________________
____________________________________________________________________________ F CLASS N_MATRIX BINARY ARITHMETIC FUNCTIONS ____________________________________________________________________________
// Input                nmx  : Input normal matrix (this) //			 mx  : Second matrix // Output               mx1  : New matrix which is the addition of //			       the two input arrays; mx1 = nmx + mx // Note			     : Uses internal structure of other matrix types // Insure sizes compatible  //    Mismatched dimensions //    Cannot perform addition // Add two n_matrices // Construct new empty matrix sum // Element of mx:  mij = <i|mx|j>  -> <rows_|mx|cols_> // Element of nmx: nij = <i|nmx|j> -> <rows_|nmx|cols_> // Element of sum: sij = <i|sum|j> -> <rows_|sum|cols_> // Perform the addition element-wise // <i|sum|j> = <i|nmx|j> + <i|mx|j> // Add h_matrix to n_matrix // Construct new empty matrix sum // Element of hmx: hij = <i|hmx|j> -> <0|hmx|0> // Element of nmx: nij = <i|nmx|j> -> <0|nmx|0> // Element of sum: mij = <i|sum|j> -> <0|sum|0> // Elements nji=<j|nmx|i>, mji=<j|sum|i> // Loop over the rows // <i|sum|j> = <i|nmx|j> + <i|hmx|j>, i=j // <j|sum|i> -> <i+1|sum|i> // <j|nmx|i> -> <i+1|nmx|i> // <i|sum|j> -> <i|sum|i+1> // <i|nmx|j> -> <i|nmx|i+1> // <i|hmx|j> -> <i|hmx|i+1> // Loop over the columns // <i|sum|j> = <i|nmx|j> + <i|hmx|j>, i<j // <j|sum|i> = <j|nmx|i> + <i|hmx*|j>, i<j // <j|sum|i> -> <j+1|sum|i> // <j|nmx|i> -> <j+1|nmx|i> // Add d_matrix to n_matrix // New n_matrix equal to this // Deal with only the diagonals // <i|sum|j> = <i|nmx|j>+del   <i|dmx|j> //                          i,j // Add i_matrix to n_matrix // New n_matrix equal to this // Add 1 to its diagonal // <i|sum|j> = <i|nmx|j>+del   <i|imx|j> //                          i,j // Add generic matrix to n_matrix