src/Matrix/n_matrix.cc
// ______________________________________________________________________
// I              CLASS N_MATRIX SIMPLE BINARY FUNCTIONS
// ______________________________________________________________________
I              CLASS N_MATRIX SIMPLE BINARY FUNCTIONS
I CLASS N_MATRIX SIMPLE BINARY FUNCTIONS
complex n_matrix::trace(_matrix* mx)
  {
  if(!CheckDim(mx,1,1) || !CheckDim(mx,3,1))	
    NMxfatal(24);				
  complex z(0);				
  switch(mx->stored_type())
    {
    case n_matrix_type: 			
	{
        complex *n00 = data;			
        complex *m00 = ((n_matrix*)mx)->data;	
        complex *nz = n00 + rows_*cols_; 	
        complex *mz = m00 + rows_*cols_; 	
        complex *nij, *mji, *m0i;		
	for(nij=n00,m0i=m00; nij<nz; m0i++)	
	  for(mji=m0i; mji<mz; nij++,mji+=rows_)
	    z += (*nij) * (*mji); 		
	return z;
	}
	break;
    case h_matrix_type: 			
	{
        complex *n00 = data;			
        complex *h00 = ((h_matrix*)mx)->data;	
        complex *nij, *hij, *hji, *h0i;		
        int i=0, j;
        hij = h00;                              
        z = complex0;                           
	for(h0i=h00,nij=n00; i<rows_; i++,h0i++)
          {					
	  for(hji=h0i,j=0; j<=i;		
                j++,nij++,hij=hji,hji+=cols_-j)	
	    z += (*nij) * (*hji);		
          for(hij++; j<cols_; j++,nij++,hij++)
	    z += (*nij) * conj(*hij); 		
          }
	return z;
	}
	break;
    case d_matrix_type: 			
	{
        complex *n00 = data;			
        complex *d00 = ((d_matrix*)mx)->data;	
        complex *dend = d00 + rows_;		
        complex *nii, *dii;			
        int incd = cols_+1;			
	for(nii=n00,dii=d00; dii<dend;
                             nii+=incd, dii++)
          z += (*nii) * (*dii);			
	return z;
	}
	break;
    case i_matrix_type: 			
	{
        for(int i=0; i<rows_; i++)		
          z += data[i*cols_+i];			
        return z; 				
	}
	break;
    default:                 			
	{
        complex *n00 = data;			
        complex *nij;				
        int i,j;
	for(i=0, nij=n00; i<rows_; i++)
	  for(j=0; j<cols_; j++, nij++)
	    z += (*nij) * (*mx)(j,i); 		
	return z;
	}
    }
  }
complex n_matrix :: trace ( _matrix * mx ) { if ( ! CheckDim ( mx , 1 , 1 ) || ! CheckDim ( mx , 3 , 1 ) ) NMxfatal ( 24 ) ; complex z ( 0 ) ; switch ( mx -> stored_type ( ) ) { case n_matrix_type : { complex * n00 = data ; complex * m00 = ( ( n_matrix * ) mx ) -> data ; complex * nz = n00 + rows_ * cols_ ; complex * mz = m00 + rows_ * cols_ ; complex * nij , * mji , * m0i ; for ( nij = n00 , m0i = m00 ; nij < nz ; m0i ++ ) for ( mji = m0i ; mji < mz ; nij ++ , mji += rows_ ) z += ( * nij ) * ( * mji ) ; return z ; } break ; case h_matrix_type : { complex * n00 = data ; complex * h00 = ( ( h_matrix * ) mx ) -> data ; complex * nij , * hij , * hji , * h0i ; int i = 0 , j ; hij = h00 ; z = complex0 ; for ( h0i = h00 , nij = n00 ; i < rows_ ; i ++ , h0i ++ ) { for ( hji = h0i , j = 0 ; j <= i ; j ++ , nij ++ , hij = hji , hji += cols_ - j ) z += ( * nij ) * ( * hji ) ; for ( hij ++ ; j < cols_ ; j ++ , nij ++ , hij ++ ) z += ( * nij ) * conj ( * hij ) ; } return z ; } break ; case d_matrix_type : { complex * n00 = data ; complex * d00 = ( ( d_matrix * ) mx ) -> data ; complex * dend = d00 + rows_ ; complex * nii , * dii ; int incd = cols_ + 1 ; for ( nii = n00 , dii = d00 ; dii < dend ; nii += incd , dii ++ ) z += ( * nii ) * ( * dii ) ; return z ; } break ; case i_matrix_type : { for ( int i = 0 ; i < rows_ ; i ++ ) z += data [ i * cols_ + i ] ; return z ; } break ; default : { complex * n00 = data ; complex * nij ; int i , j ; for ( i = 0 , nij = n00 ; i < rows_ ; i ++ ) for ( j = 0 ; j < cols_ ; j ++ , nij ++ ) z += ( * nij ) * ( * mx ) ( j , i ) ; return z ; } } }
______________________________________________________________________
I              CLASS N_MATRIX SIMPLE BINARY FUNCTIONS
______________________________________________________________________
______________________________________________________________________ I CLASS N_MATRIX SIMPLE BINARY FUNCTIONS ______________________________________________________________________
// Input            nmx : An n_matrix (this) //                   mx : A second matrix // Output            z  : A complex value which is the trace of //                        the product of two input arrays //                         z = Tr{nmx*mx} = Tr{mx*nmx} // Note                 : This is faster than taking the matrix //                        product and then using unary trace! // Note                 : Uses internal structure of other matrix types //       ---                 --- --- //   z = \   <i|nmx*mx|i>  = \   \   <i|nmx|j><j|mx|i> //       /                   /   / //       ---                 --- --- //        i                   i   j // Insure result is square // If not, cannot do trace // Initialize the result at zero // Tr{n_matrix * n_matrix}     // Start of nmx: n00 = <0|nmx|0> // Start of mx: m00 = <0|mx|0> // End of nmx: nz = <rows_|nmx|cols_+1> // End of mx: mz = <cols_|mx|rows_+1> // These are the matrix elements // Effective loop over i // Effective loop over j // Tr{nmx*mx} += <i|nmx|j><j|mx|i> // Tr{n_matrix * h_matrix}     // Start of nmx: n00 = <0|nmx|0> // Start of hmx: h00 = <0|hmx|0> // These are the matrix elements // Keeps the compiler happy // Keeps the compiler happy // Loop rows of nmx, cols of hmx // Split loop of nmx cows, hmx rows // to be over upper triangle of hmx // (i.e. hmx when j<=i, hmx* j>i) // Tr{nmx*hmx} += <i|nmx|j><j|hmx|i>, j<i // Tr{nmx*hmx} += <i|nmx|j><i|hmx*|j>, j>=i // Tr{n_matrix * d_matrix}     // Start of nmx: n00 = <0|nmx|0> // Start of dmx: d00 = <0|dmx|0> // End of dmx: dend = <rows_|dmx|rows+1> // These are the matrix elements // <i|nii|i>+incd -> <i+1|nii|i+1> // Tr{nmx*dmx} += <i|nmx|i><i|dmx|i> // Tr{n_matrix * i_matrix}     // Sum over the diagonal elements // Tr{nmx*imx} += <i|nmx|i><i|imx|i> //             += <i|nmx|i> // Tr{nmx * unknown matrix} // Start of nmx: n00 = <0|nmx|0> // This is matrix element <i|nmx|j> // Tr{nmx*mx} += <i|nmx|j><j|mx|i>